dataesr/scientific-paragraphs-categorization huggingface dataset:
Since exploring all possible scenarios would be computationally too costly, IGoR restricts its exploration to the reasonably likely ones. Scenario exploration takes from 1 ms up to less than a second per sequence on a single CPU core, depending on the chain (see full distributions of runtimes in Fig. S1 ). Different recombination architectures and dependencies can be configured within IGoR by specifying dependencies between elementary events (gene choices, deletions, insertions, hypermutations) through an acyclic directed graph, or Bayesian network, as illustrated in Fig. 1b  for the case of TCR β chains (see Online Methods for the other used structures). IGoR functions according to three modes: learning, analysis, and generation (Fig. 1c ). In the learning mode, IGoR infers the recombination statistics of large datasets arXiv:1705.
(c) IGoR's pipeline includes three modes. In the learning mode, IGoR learns recombination statistics from data sequences. In the analysis mode, IGoR outputs detailed recombination scenario statistics for each sequence. In the generation mode, IGoR produces synthetic sequences with specified recombination statistics. of sequences using a Sparse Expectation-Maximization algorithm (see Online Methods). In the analysis mode, IGoR assigns recombination events to sequences in a probabilistic way, by outputing the most likely scenarios ranked by their probabilities, as well as the overall generation probability of the sequence. In the generation mode, IGoR outputs random sequences with specified statistics, e. g. learned from real datasets.
Inference of V(D)J recombination We used IGoR's learning mode to infer the accurate statistics of V(D)J recombination from four datasets comprised of unique sequences of non-productive rearrangements of three different chains, sequenced either at the levels of mRNA (TCRα chain or TRA, and TCRβ chain or TRB  [14] ) or DNA (TRB  [13] , BCR heavy chain or IGH from naive cells  [12] ), generalizing earlier methods  [15] [16] [17]. Restricting to nonproductive unique sequences allowed us to avoid biases introduced by functional selection. The Expectation-Maximization algorithm converged within a few iterations (see Fig. S2  for convergence of parameters, and Fig. S3  for the case of IGH). The same TRB insertion and deletion distributions were inferred regardless of the individual, laboratory of origin, or sequencing protocol, and of whether DNA  [13]  (light blue distributions in Fig.
2 ) or mRNA  [14]  (dark blue) was used. By contrast, V and J gene usage varied moderately but significantly across individuals, and even more across sequencing technologies, suggesting possible primer-dependent biases (Fig. S4 , see also Fig. S17  for IGH D-J gene usage). Insertions at the TRA V-J junction, and at the TRB V-D and D-J junctions have similar distributions (Fig. 2a ), as previously reported  [17]. IGH have significantly more insertions at the junctions than TCRs, consistent with previous observations  [16]. We then validated the learning algorithm on synthetic datasets. Sequences were generated in batches of 10 3 to 10 5 by IGoR with a variable error rate, using statistics inferred from 60bp DNA TRB data. IGoR's learning algorithm was then run on these raw sequences, and the resulting statistics compared to the known ground truth.
We found that the inference was highly accurate for datasets of 10 5 sequences and an error rate set to its typical experimental value, 10 -3 (Fig. 3a  and b ), and was not affected by overfitting. However, not all high-throughput sequencing datasets reach this depth, especially when restricted to unique non-productive sequences. In addition, hypermutation rates in BCRs, which IGoR treats in the same way as errors, can reach 1-10%. To assess how these limitations affect accuracy, we calculated the Kullback-Leibler divergence (a non-parametric measure of difference between probability distributions, see Online Methods) between the true distributions and the inferred ones, for varying sizes of datasets and error rates. For an er-   [12] , TRB at the VD and DJ junction from both DNA  [13]  and mRNA data  [14] , and TRA at the VJ insertion site from mRNA data  [14].
Comparison to other methods We compared our method to two representative stateof-the-art algorithms: MiXCR  [8] , an efficient assignment tool that finds the best matching germline genes, and Partis  [10] , a BCR-specific tool that uses maximum likelihood to find the most plausible scenario. 130 base-pair IGH sequences were synthetized in silico from a datainferred model using IGoR's generation mode. We then assigned recombination scenarios using MiXCR, Partis and IGoR, and compared them to the true scenarios with which sequences were generated. In IGoR's and Partis' case, the model parameters were learned from the generated dataset to mimick the analysis of real data. Fig. 4c  shows the performance of the three methods in assigning the correct scenario of recombination. IGoR performs about 2.
S8 ). However, the dependency between D an J usage in TRB is correctly captured by IGoR but not by the other methods (Fig. 4d ). TRB D and J genes are organised in two clusters, one containing D1 followed by genes of the J1 family, the other containing D2 followed by genes of the J2 family. Because of this organisation, D2 cannot be recombined with genes from the J1 family  [18]. MiXCR assigns 20% of impossible D2-J1 recombination events to sequences (note that Partis does not process TCRs). By constrast, IGoR correcly learns the rule by assigning zero frequency to these impossible D-J pairs. The same results are obtained directly on real data (see Fig. S9 ). Finally, IGoR accurately reconstructs the distribution of insertions, while the other methods systematically overestimate the probability of zero insertions (Fig.
(c) Frequency with which IGoR, MiXCR and Partis call the correct scenario of recombination as the most likely one ('scenario'), as well as each separate feature of the scenario ('V gene,' etc. ). 'Failed' corresponds to sequences for which the algorithm did not output an assignment. (c) Usage frequency of TRB D gene conditioned on the J gene, inferred by the IGoR and MiXCR (Partis does not handle TCR sequences). IGoR recovers the physiological exclusion between D2 and J1, while MiXCR does not. thetic data with hypermutations according to the model learned from the real dataset, and re-learned its parameters using IGoR, finding excellent agreement (Fig. S10 ). The PWM prediction for the position-dependent probability of hypermutations correlated well with that actually observed in the sequences (r = 0.
DISCUSSION By treating alignments of immune receptors to the germline probabilistically  [15] , IGoR corrects for systematic biases in the estimate of V(D)J recombination statistics, and predicts recombination scenarios more accurately than previous methods. Its detailed analysis of recombination scenarios further reveals that, even with a perfect estimator, the scenario is incorrectly called in more than 70% of sequences, suggesting caution when interpreting results from deterministic assignments. Although we demonstrated its functions on human TRB and IGH, IGoR's flexible structure makes it applicable to any variable lymphocyte receptor (TCR or immunoglobulin) and species for which genomic data is available. Unlike Hidden Markov Model based methods (e. g. [10, 17] ), it can include a wide array of possible dependencies between the recombination events.
It can also be adapted to handle unusual or incomplete rearrangements (D-J rearrangments, DD2/DD3 rearrangements in TCR δ chains, hybrid TRA/TRD recombinations, etc. ). IGoR can also help detect unusual rearrangement features by using its syntheticaly generated sequences as a control. For instance, rearrangements with tandem Ds have been reported  [12] , but distinguishing them from random insertions can be challenging. To test this, we counted sequences with two ≥10-nt D segments in the data, and compared it with predictions from IGoR's synthetic sequences generated with a single D segment (see Online Methods). We found 5 times more double-D assignments in IGH data than in the control, validating the findings of  [12]. In contrast, the same analysis performed on TRB showed no significant presence of tandem Ds.
Alignment to germline and scenario listing In the analysis and learning modes, each sequence is first aligned to all possible germline genes retrieved from genomic databases (e. g. IMGT), using the Smith-Waterman algorithm  [29]. Only germline genes with alignment scores higher than an adjustable threshold are considered for further analysis (see SI Text for details). Possible scenarios are then listed by picking germline genes with an above-threshold alignment score, and by choosing a number of base pairs to further delete from the ends of their aligned parts. The base pairs located between the germline segments trimmed in this manner are called insertions, and alignment mismatches to the germline are called errors or hypermutations. When the palindromic end of germline genes is not entirely deleted, the number of remaining palindromic base pairs are described as negative deletions.
Sequence analysis For each sequence in the dataset, the probability of possible scenarios are computed using the recombination probability of Eqs. 1 or 2, multiplied by the probability of errors or hypermutations P err : P scenario = P recomb × P err. Scenarios are then listed in order of decreasing probability. The sum of probabilities P recomb × P err of possible recombination and hypermutation events gives the probability of observation of that particular sequence read, P read. The probability that the pre-mutation sequence was generated by recombination, P gen , is defined as the sum of the probabilities P recomb of scenarios leading to that sequence. Since the pre-mutation sequence is not known with certainty, we calculated an approximate generation probability P gen as the geometric mean of P gen of all possible unmutated sequences consistent with the read, weighted by their posterior probabilities, P gen × P err /P read.
Learning algorithm The learning algorithm infers the parameters of Eqs. 1 or 2, as well as the error or hypermutation model parameters of Eq. 3, from a large datasets of unique sequences. It relies on the sequence analysis module, and follows an Expectation-Maximization procedure. Starting from an arbitrary (but reasonable) set of parameters, all sequences in the dataset are analysed as described above, producing a long list of scenarios associated with each sequence. We define the pseudo-loglikelihood as the weighted sum of the log-likelihoods of all scenarios of all sequences, where the weights are given by the conditional probabilities of scenarios given the sequence, P recomb /P read (Expectation step). This pseudolog-likelihood is then maximized with respect to the parameters of the log-likelihoods (Eqs.
Datasets We applied the learning algorithm on the following publicly available datasets: TCR alpha and beta chains RNA datasets from  [14]  are available on Sequence Read Archive (SRP078490); TCR beta chains 60 bp DNA datasets from  [15]  are available at http://physics. princeton. edu/ ~ccallan/TCRPaper/data/; naive and Memory BCR heavy chains DNA datasets from  [16, 30]  are available at http://physics. princeton. edu/ ~ccallan/BCRPaper/data/.
The dependence on θ reflects the fact that θ also includes the parameters of the error or hypermutation model. We write the joint probability of producing a given sequence S and observing a given read R as: P (R, S|θ) = P gen (S|θ)P err (R|S, θ). (8) Summing over all possible recombination products, the likelihood of a sequencing read is: P read (R|θ) = S P (R, S|θ) = E P recomb (E|θ)P err (R|E, θ), (9) and the total likelihood of the model given a dataset of reads (R 1 ,. , R N ) is given by: L total (θ) = N a=1 P read (R a |θ). (10) B Expectation-maximization 1 General scheme The recombination machinery is degenerate, as several scenarios of recombination and hypermutations can lead to the same sequence, and the recombination scenario E from which the sequencing read R comes from is in general unknown.
E Data and software 1 Genomic templates We used custom genomic templates derived from the IMGT database  [33]. TCR alpha V and J genomic templates were taken from the IMGT human database. For TCR beta V, D and J genes we used curated genomic templates from  [15]. BCR heavy chain V, D and J genes were taken from the customized genomic templates used in  [16]. For software comparison we used default genomic templates provided with Partis and MiXCR.
Comparison to other software We benchmarked our method against MiXCR 2. 0. 2  [8]  -a commonly used deterministic alignment method. We used the MiXCR sequence assignment to compute the frequency of gene usage, insertion length, deletions and obtain the distributions shown in Fig. S8. We also compared to Partis  [10]  -a recent HMM based model of recombination. Since Partis uses a Viterbi learning algorithm, we used the most likely assignments it outputs to compute the corresponding probability distribution shown in Fig. S8. Since Partis is designed to handle BCRs we assessed its performance on the BCR dataset only. We benchmarked IGoR's performance for evaluating possible recombination scenarios on real data sequences used to infer the models presented in the main text. We used 60bp TCR β sequences for benchmarking since the difficulty for finding the correct V and J for alignment is higher.
C. FIG. S4 : Gene usage in TRB mRNA vs DNA data. We plot the marginal gene usage averaged over conditional dependencies for V, D and J genes respectively inferred using IGoR from mRNA 100bp (red) and DNA 60bp (blue) technology datasets. We observe a higher inter-method than inter-individual variability. FIG. S6:  A probabilistic assignment approach is crucial for TCRs. Equivalent of main text Fig. 4b  for 30000 60bp TCRs. This figure shows the distribution of the number of scenarios that need to be enumerated (from most to least likely) to include the true scenario with 50% (blue), 75% (green), 90% (red), or 95% (cyan) confidence. The shorter read length compared to 130bp BCRs entail a higher uncertainty on the V gene identity, for which a higher number of scenarios must be considered. and Partis Viterbi learning we used the output assignments to compute the corresponding recombination statistics.
We plot a for base pairs that have been observed at least 2000 times on a 100 000 sequences dataset, for which we can compute a reliable mutation frequency, and the mutation frequency predicted by our model. The two top panels show good predictive power for the gene on which the model was learned. However the two bottom panels show a lesser ability to predict the correct mutation frequencies on the whole locus, hence suggesting that differences observed in inferred position weight matrices (Fig. S13 ) are of biological relevance.
Satellite Imagery A series of nine optical Sentinel-2 and nine SAR Sentinel-1 images were acquired from autumn to spring (August 2016 to May 2017) from the data hub of the European Space Agency and Centre National d'Etudes Spatiales  [41]. Sentinel-1 images were acquired in Single Look Complex (SLC) mode (delivered with VH and VV polarization states) with an incidence angle of 45-47 •. The associated range and azimuth spatial resolutions were respectively 2. 3 and 13. 9 m. Sentinel-2 level 2A images were acquired (i. e. , corrected for geometric and atmospheric effects) with a spatial resolution of 10-60 m and a spectral resolution of 13 bands. Characteristics of the optical and SAR images are summarized in Table  2.
Backscattering Coefficients Sentinel-1 images were first radiometrically calibrated (Figure  3 ) using SNAP (Sentinel Application Platform) v5. 0 software with the following equation  [42] : value(i) = |DN i | 2 A i 2 , ( 1 ) where DN is the digital number of each pixel (amplitude of backscattering signal) and A is the information necessary to convert SAR reflectivity into physical units provided in the Calibration Annotation Data Set in the image metadata. This equation transforms the DN of each pixel into a backscattering coefficient on a linear scale.
Backscattering Coefficients Sentinel-1 images were first radiometrically calibrated (Figure  3 ) using SNAP (Sentinel Application Platform) v5. 0 software with the following equation  [42] : ( ) = | |² ² , (1) where DN is the digital number of each pixel (amplitude of backscattering signal) and A is the information necessary to convert SAR reflectivity into physical units provided in the Calibration Annotation Data Set in the image metadata. This equation transforms the DN of each pixel into a backscattering coefficient on a linear scale. A refined Lee filter  [43, 44]  was then applied with a window of 7 × 7 pixels to reduce speckle noise using SNAP v5. 0 software. This window size was selected to decrease speckle noise while preserving a suitable spatial scale, which was necessary to ensure identification of winter land use.
Sentinel-1 images were then geocoded using Shuttle Radar Topography Mission 3s data to correct topographic deformations. The accuracy of geometric correction was less than 10 m per pixel. A backscattering ratio was calculated by dividing by. This ratio highlights differing scattering mechanisms of each target. The backscattering coefficients and and the backscattering ratio : were then converted into decibels (db) using the following equation: ( ) = 10 × ( ) , (2) Polarimetric Parameters A 2 × 2 covariance matrix ( ) was first extracted from the scattering matrix S of each SLC image using PolSARpro v5. 1. 1 software  [45]. The elements of the matrix, which are independent of the A refined Lee filter  [43, 44]  was then applied with a window of 7 × 7 pixels to reduce speckle noise using SNAP v5.
0 software. This window size was selected to decrease speckle noise while preserving a suitable spatial scale, which was necessary to ensure identification of winter land use. Sentinel-1 images were then geocoded using Shuttle Radar Topography Mission 3s data to correct topographic deformations. The accuracy of geometric correction was less than 10 m per pixel. A backscattering ratio was calculated by dividing σ 0 V H by σ 0 VV. This ratio highlights differing scattering mechanisms of each target. The backscattering coefficients σ 0 V H and σ 0 VV and the backscattering ratio σ 0 V H : σ 0.
VV were then converted into decibels (db) using the following equation: σ 0 (db) = 10 × log 10 (σ 0 ) , (2) Polarimetric Parameters A 2 × 2 covariance matrix (C 2 ) was first extracted from the scattering matrix S of each SLC image using PolSARpro v5. 1. 1 software  [45]. The elements of the matrix, which are independent of the polarimetric absolute phase  [46] , were directly geocoded using SNAP v5. 0 software. A refined Lee filter was then applied using a window of 7 × 7 pixels to reduce speckle noise. The second (C12) and fourth (C22) elements of the diagonal were compared to the backscattering coefficients. The SPAN, which corresponds to the total scattered power and the Shannon Entropy (SE), which corresponds to the sum of two contributions related to the intensity (SE i ) and the degree of polarization (SE p )  [46] , were then calculated from the matrix.
Calculation of Vegetation Indices and Biophysical Parameters Based on the literature, 2 vegetation indices, 1 water index and 3 biophysical parameters were computed using SNAP v5. 0 software. The two vegetation indices were derived from the near infrared and red bands: The Normalized Difference Vegetation Index (NDVI;  [49] ) and the Soil Adjusted Vegetation Index (SAVI;  [50] ) that have demonstrated their relevance to study land use  [51, 52]. The Normalized Difference Water Index (NDWI;  [53] ) was calculated from the near and middle infrared bands often used to estimate water in vegetation  [54]. Biophysical parameters (Leaf Area Index, fraction of photosynthetically active radiation and fractional vegetation cover) were calculated using the PROSAIL radiative transfer model  [55, 56]  implemented in SNAP v5.
Processing Sentinel-1 and -2 Time-Series The method for processing Sentinel-1 and -2 time-series we developed had three steps (Figure  4 ): (i) extracting features to yield one parameter value per field, (ii) reducing the optical and SAR parameter dataset using a data-mining algorithm and (iii) classifying vegetation cover using RF and SVM algorithms. In total, the number of parameters to be processed was 297: 162 for Sentinel-2 data (18 parameters X 9 dates) and 135 for Sentinel-1 (15 parameters X 9 dates).
Reduction of the Parameter Dataset A two-step approach was used to remove noise and inconsistent data from the parameter datasets used for the classification. First, a correlation matrix was computed and a threshold of 0. 95/-0. 95 was applied to find the most correlated parameters (268 out of 297 parameters). Then, the choice of parameters to be removed was made based on a current state-of-the-art to select the most relevant parameters for the study of inter land use  [27, 51, 52]. At this stage, the dataset consisted of 144 parameters (90 optical and 54 radar parameters). In a second step, an analysis of the parameter relevance was performed using the RF significance function to further reduce the number of parameters to be classified. A measure of variable importance was provided for each candidate predictor using this heuristic method based on the Gini Index  [36, 58, 59 ].
Land Use Classification The RF and SVM algorithms, which are supervised classification methods, were used to classify land use during winter 2016-2017. These algorithms were chosen for their consistently strong performance and the accuracy with which they classify LULC  [60, 61]. The RF algorithm is an ensemble algorithm that uses a set of classification and regression trees to make a prediction  [36]. The package randomForest developed by  [62]  and implemented in R (v. 3. 3. 2) was used to perform winter   3  for the definition of abbreviations.
Land Use Classification The RF and SVM algorithms, which are supervised classification methods, were used to classify land use during winter 2016-2017. These algorithms were chosen for their consistently strong performance and the accuracy with which they classify LULC  [60, 61]. The RF algorithm is an ensemble algorithm that uses a set of classification and regression trees to make a prediction  [36]. The package randomForest developed by  [62]  and implemented in R (v. 3. 3. 2) was used to perform winter land use classifications. Two Random Forest parameters, namely the number of trees (ntree), which was created by randomly selecting samples from the training dataset  [28]  and the number of variables used for tree nodes splitting (mtry) were tuned and randomly determined using the tune function implemented in the randomForest package.
The package e1071 (v 1. 7-0) developed by  [64]  and implemented in R (v. 3. 3. 2)  was used to perform SVM classifications. A set of four SVM parameters (gamma, cost, degree, nu) was randomly tuned using the tune function integrated in the e1071 package. Then, several tests were carried out to determine the optimal kernel for winter use classification. At the end of the selection process the polynomial kernel was selected. Results obtained with RF and SVM algorithm were compared to evaluate their suitability for classifying vegetation cover into land use classes. Classification performance was estimated using a cross-validation test. The classification was applied to a varying subset of 257 fields (500 pixels): two-thirds (171 fields or 332 pixels) were used for training, one-third (86 fields or 168 pixels) for validation  [19].