dataesr/scientific-paragraphs-categorization huggingface dataset:
Identification of Differentially Expressed Genes Between Self-Renewing and Differentiating Progenitors In order to identify a pool of genes potentially relevant in the differentiation process, we analyzed the transcriptome of self-renewing and differentiating primary chicken erythrocytic progenitor cells (T2EC) using RNA-Seq. We sequenced two independent libraries from selfrenewing T2EC and two independent libraries from T2EC induced to differentiate for 48 h. For each condition, we first verified that read counts between replicates were reproducible (S3A and S3B Fig). We then identified 424 significantly differentially expressed genes (p-value < 0. 05, S3C Fig). Gene ontology analysis using the DAVID database  [60]  revealed a clear over-representation of genes involved in sterol biosynthesis in this list (not shown).
Based on the result of Kmeans clustering, we selected around thirteen genes per group to represent each cluster equally. This left us with 92 genes for further analysis (S1 Table ). We then used STRING database to search for known connections among these genes. The result confirmed the existence of a strongly connected subnetwork associated with sterol synthesis (S5B Fig). Moreover, this analysis also revealed the presence of another highly connected subnetwork mostly composed of genes involved in signaling cascades and two transcription factors (BATF and RUNX2). Those two main networks are linked by the gene HSP90AA1 which encodes the molecular chaperone HSP90alpha. Its activity is not only involved in stress response but also in many different molecular and biological processes because of its important interactome.
Cell-to-Cell Heterogeneity Blurred Cell Differentiation Process We measured the expression level of the selected 92 genes by single-cell RT-qPCR using 96 cells isolated from the most informative time-points of the differentiation sequence. Based upon preliminary experiments, we decided to analyze cells from six time-points during differentiation. After data cleaning (see Materials and Methods), we obtained the expression level of 90 genes in 55,  73, 72, 70, 68 , and 51 single cells from 0, 8, 24, 33, 48, and 72 h of differentiation, respectively. One should note that the variability we observed at the single-cell level originates from two types of sources: biological sources and experimental sources. We therefore tested the technical reproducibility of different RT-qPCR steps liable to generate such experimental noise (see Materials and Methods).
9% of the variability in the dataset. Hence, unlike in the population-averaged data, the differentiation process did not represent the main source of variability at the single-cell level. The application of HCA further confirmed that the classification became more complex for single-cell data (Fig  2B ). Contrary to bulk analysis, individual cells from the same time-point were not necessarily more similar to each other than to cells from neighboring time-points. Consequently, the clustering of individual cells into groups became complicated. The picture of cell differentiation process that emerged from the single-cell analysis thus far was more complex than the one obtained from the population level analysis. This difference between single-cell and population-level analysis arises from the unraveling of cell-to-cell heterogeneity in the single-cell data, which could have been hidden by the averaging effect of the population (see below).
PCA is a linear method for dimensionality reduction of single-cell data. In view of non-linear relationships of cell states in state space, recently nonlinear techniques like t-SNE  [55]  or diffusion maps  [63]  have been applied in single-cell data analysis. t-SNE is a variation of Stochastic Neighbor Embedding deemed capable of capturing more local structures than classical PCA, while also revealing global structure such as the presence of clusters at several scales. Diffusion maps use a non-linear distance metric (referred to as diffusion distance), which is deemed conceptually relevant in view of noisy diffusion-like dynamics during differentiation  [63]. We therefore applied these algorithms on our datasets, as well as another non-linear version of PCA, called Kernel PCA  [64] , not previously applied to single-cell gene expression data (Fig 2C  to 2E ).
Single-Cell Data Embed Population Information and Reveal New Discriminating Genes Involved in the Differentiation Process In order to assess to what extent the differentiation process was still visible in the single-cell data, we performed PCA on datasets from the two extreme time-points, 0 and 72 h (Fig  3A ). The result showed a clear separation of both time-points with only a few cells intermingled. We also performed HCA on datasets from the same time-points (Fig  3B ). Again, the segregation of the cells was still not perfect, but cells were not as mixed as before. Here, there exist two clusters of self-renewing and differentiating cells. When compared to the analysis of the entire time series, the separation between cells from the two extreme time-points looked clearer. Therefore, the analysis of single-cell data confirmed that part of the information present in the single-cell data is linked to the differentiation process.
The idea that shared information was present in single-cell and population-based data was reinforced by the analysis of the correlation matrices within and between the two datasets (S6 Fig). It was apparent that (1) the global intensity of the correlations was higher with population-based data and (2) there existed a co-structure between the two datasets. At the population level, we showed that the set of genes selected was relevant to analyze the differentiation process (Fig  1 ). The cross-correlation analysis strengthened this view and demonstrated that when looking at the single-cell scale, the information held by these genes was not totally erased by cell-to-cell variability. We then looked at the genes that contributed the most to the PCA outcome (Fig  3C ). Among the genes that discriminate the most self-renewing cells, one could highlight LDHA (Lactate deshydrogenase A), CRIP2, and Sca2.
The introduction of asynchrony induced a flatter time profile of the entropy (Fig  9F ). This finding did not, however, prove that our cells are synchronously differentiating, but only demonstrated the effect of asynchrony: in the background of bursty gene transcriptional process, asynchrony will tend to smoothen (and not augment) the entropy of the system. Therefore the observed surge in entropy can not be attributed to the asynchrony of the process. The rise-and-fall of entropy in our data is in line was examined in a different setting, namely a reprogramming process  [58]. The authors stated, "The initial transcriptional response is relatively homogeneous," offering the opportunity to examine the entropy time profile in such a homogeneous process. Our analysis of this dataset produced a similar behavior for entropy which significantly increased initially, before returning to lower values (S8 Fig).
However, when assessed at the single-cell level, our analyses revealed a much higher cellular heterogeneity. Despite this heterogeneity, the selected genes were still effective in separating the two most extreme time-points in T2EC differentiation, confirming that information associated with the differentiation process is embodied in the gene expression data at the single-cell level. From the dataset that we generated at the single-cell level, two main results could be obtained: (i) regarding the biology of the erythroid differentiation, we identified previously unidentified genes as being important components of the self-renewal and differentiation of erythroid progenitors, and (ii) on a larger perspective, our results fully supported a dynamical view where differentiation can be seen as a critical phase transition driven by stochasticity.
5 min at 65˚C. RNA were reversetranscribed in presence of SuperScript III Reverse Transcriptase enzyme, from the Super-Script VILO cDNA Synthesis kit (Invitrogen), and T4 gene 32 protein (New England Biolabs) to improve reverse transcription efficiency. The reaction thermal cycling conditions were 5 min at 25˚C, 30 min at 50˚C, 25 min at 55˚C, 5 min at 60˚C and 10 min at 70˚C. Specific target amplification of cDNA. Primers were designed using the Ensembl database (http://www. ensembl. org/Gallus_Gallus/Info/Index/) and Primer3Plus software (http:// www. bioinformatics. nl/primer3plus/). For information about the primers sequences used, please contact the authors. The cDNA pre-amplification was performed using the TaqMan PreAmpMaster (Applied Biosystems) mixed with all primer pairs of the genes of interest (Sigma-Aldrich), diluted at 500 M.
Pre-amplified cDNA were mixed with Sso Fast EvaGreen Supermix With Low ROX (Bio-Rad) and DNA binding dye sample loading reagent (Fluidigm). Primer pairs of the genes of interest were diluted at 5 μM with the Assay Loading Reagent (Fluidigm) and low EDTA buffer. First, the 96. 96 DynamicArray IFC chip (Fluidigm) was primed. Then, prepared cDNA and primer pairs were loaded in the inlets of this device. To avoid chip-linked variability, when analyzing single-cell data we were careful to represent every time-point in each of the four microfluidic-based chip analyzed. The prime step and transfer of cDNA samples and primers from the inlets into the chip were performed using the IFC Controller HX (BioMark HD system). The chip was analyzed using the BioMark HD reader according to the GE 96 × 96 PCR + Melt v2.
pcl program, thanks to the data collection software. Then, raw data were analyzed with the Fluidigm Real-Time PCR Analysis software. Positive exogenous controls (RNA spikes) were used to validate the RT-qPCR experiment as recommended by Fluidigm Company. We also used the RNA spikes to normalize the data (see below). To determine qPCR efficiency of every primer pairs used, serial dilution scales of bulk-cell cDNA were performed. PCR efficiencies were calculated as follows: E = 10 -1/slope. Primer pairs presenting PCR efficiency less than 80% or more than 120% were removed from subsequent analyses. RT-qPCR: low-level data analysis. First, a manual examination was performed regarding data quality. RTqPCR data were exported from the BioMark HD data collection software. On every microfluidic-based chip, each gene was controlled in a qualitative manner in order to keep only reliable and good quality data.
After this quality control, Cq values of data tagged as "1" were replaced with UD (for "undefined") in the raw data file, since they would not be taken into account in later analysis. Then the new table underwent an automatic formatting consisting in a second multiple-criteria cleaning process using an in-house R script. Cq values were converted into (approximately) absolute numbers of molecules according to the following steps. First, we selected cells with at least one valid spike measurement (i. e. , whose Cq is different from UD and 999). Then, we normalized the raw value c Cq i;j for cell i and RNA j according to the cell mean spike value Cq i (or the only available spike if one is invalid), with the global mean spike value Cq 0 as reference. That is, the normalized value Cq i,j for cell i and RNA j is defined by Cq i;j ¼ c Cq i;j À Cq i À Cq 0 À Á : After removing cells with abnormally important amount of genes with low expression (high Cq i,j values, suggesting the absence of a cell in the well), the numbers of mRNA molecules were estimated, considering the following: a maximum Cq equal to 30 as the measurement of 1 molecule in the well after 22 cycles of pre-amplification, a dilution factor corresponding to 1 cell extract diluted in 96 wells, and a sampling of 1/45 for PCR measurement.
Statistical Analysis Software. Most of the statistical analyses were performed using R  [46]. The k-means clustering was performed using the stats R library. PCAs were performed using the ade4 package  [47]. All PCAs were centered (mean subtraction) and normalized (dividing by the standard deviation). All PCAs were displayed according to PC1 and PC2, which are the first and second axis of the PCA respectively. Hierarchical cluster analysis was performed applying the R hclust function, using the complete linkage method on Euclidean distances. Dendrograms were built and plotted using the dendextend R library. Correlation analysis was performed using rcorr from the Hmisc R library. The p-value was corrected for multiple testing using the Bonferroni method  [48]. Networks were computed using Cytoscape  [49].
Cross-correlation analysis was performed using the matcor function from the CCA R library. Normality of the distributions was tested using the shapiro. test function. The variances were compared using the F test with the var. test function. Wilcoxon test was performed using the wilcox. test function. t-SNE and diffusion maps were computed using the Matlab Toolbox for Dimensionality Reduction (http://lvdmaaten. github. io/drtoolbox/). The t-SNE analysis was performed on a normalized version of the data, using zscore function. Kernel PCA was computed using the Matlab kPCA script  [50]  applying polynomial with fractional power 0. 1. All linear analysis methods (PCA, HCA and correlation analysis) were performed after applying the transformation m 7 ! ln(m + 1) to the data, which gives access to the more linear Cq structure.
For each time-point combination, a Mann-Whitney U test was performed using the ranksum Matlab command to detect significant variation. Wave membership was based on time variations. By definition a gene belongs to the wave at time T if there is at least one variation detected between time T and a previous time-point and if the gene does not belong to a previous wave. Only genes identified in a wave that displayed a significant variation in more than 90% of boot-straped samples were kept in this wave. Estimation of entropy. We estimated the Shannon entropy of each gene j at each timepoint t as follows: we computed basic histograms of the genes with N = N c /2 bins, where N c is the number of cells, which provided the probabilities p t j;k of each class k. Finally, the entropies were defined by E t j ¼ À X N k¼1 p t j;k log 2 ðp t j;k Þ: When all cells express the same amount of a given gene, this gene's entropy will be null.
ncbi. nlm. nih. gov/geo/download/?acc=GSE67310. The original (Log2 [FKPM]) data were transformed into FKPM data for analysis using the BPglm algorithm  [59]. Running the algorithm with an FDR value of less than 0. 00005 and using the Bonferroni correction method for multiple testing led us to a list of 776 differentially expressed genes, on which entropy was computed. Statistical significance was computed using the Wilcoxon non parametric test.
Acknowledgments We would like to thank the following colleagues for helpful advice and discussions throughout this work:  Vincent Lacroix  ( LBBE/UCBL ),  Marie Se  ´mon (IGFL/ENSL), Gae ¨l Yvert (LBMC/ENSL),  Didier Auboeuf  ( LBMC/ENSL ),  Isabelle Durand  (CLB),  David Cox  (CLB),  Nicole Dalla-Venezia  (CLB),  Jordan C. Moore  (Fluidigm),  Frederic Moret (INMG/UCBL ),  Julien Falk  ( INMG/UCBL ), and all the members of the Stochagène and Iceberg projects. We thank the Genoscope, and especially  Carole Dossat , for their invaluable help in sequencing of the RNAseq libraries. We would like to thank  Gerard Benoit (LBMC/ENSL) ,  Marieke von Lindern  (Sanquin, Amsterdam), and  Sui Huang  ( ICSB, Seattle ) for critical reading of the manuscript. We thank  Geneviève Fourel  for pointing our attention to the [58] paper.
This work was supported by funding from the  Institut Rhonalpin des Systèmes Complexes (IXXI) ,  La Ligue contre le Cancer  ( comitede Haute-Savoie ) and by two grants from the  French agency ANR  ( Stochagene ;  ANR 2011 BSV6 014 01  and  ICEBERG ;  ANR-IABI-3096 ). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
NET v3. 5. 10 (http:// www. getpaint. net, accessed 12 August 2015). Using the ImageJ software package (National Institutes of Health, Bethesda, MD, USA; http://imagej. nih. gov/ij/, accessed 12 August 2015), a set of morphological shape indices was obtained from physical measurements of each otolith image, as follows: (1) circularity, calculated as (4p Â (area/perimeter 2 )); (2) aspect ratio, the ratio of the major and minor axes of the ellipse that binds the outline; and (3) roundness, calculated as (4 Â area/(p Â major axis 2 )). Using the TPSdig utility (http://life. bio. sunysb. edu/morph/ software. html, accessed 12 August 2015 ), images were converted to binary and otolith outlines were traced using edge detection and saved as a series of x,y coordinates. Elliptical Fourier harmonics were extracted from smoothed otolith outlines (using 200 smoothing iterations) using the momocs package in R (http://cran.
Eleven coefficients showed significant deviation from normality (based on visual inspection of the probability distribution) that could not be corrected by transformation, leaving a total of 37 shape descriptors in the subsequent analysis. Principal component analysis (PCA) was conducted to reduce the dataset to a manageable number of orthogonal descriptor variables that summarised the variability in otolith shape. The scree plot (eigenvalues against principal component number) and cumulative variance values were used to identify the principal components that explained most of the variability in otolith shape. Correlations between these principal components and FL were tested statistically using Pearson's correlation coefficient. A series of general linear models (GLMs) compared principal component scores between locations and years (fixed factors), with FL included as a covariate.
Principal components that were significantly correlated with length and showed no spatial heterogeneity in the size-shape relationship (i. e. length Â location interaction P. 0. 05) were standardised using the common within-group slope, according to the following equation: Y c ¼ Y À b Â L where Y c is the corrected variable, Y is the original variable, b is the common within-group slope of the shape-size relationship and L is the measurement of fish size (FL, cm). In order to visually represent how otolith shape changes with length, the otolith outlines were grouped into 60-cm length classes  (,60, 61-120, 121-180, 181-240 and 241-300 cm)  and average outlines for each length class were reconstructed from the first 12 elliptical Fourier harmonics using the mean Shapes function in the momocs package in R.
Characterisation of Atlantic bluefin tuna for the east and west Atlantic using otolith shape From the full dataset, a subset was selected for which stock origin could be assumed with a reasonable degree of certainty. Samples of adult bluefin from the western Atlantic spawning grounds in the Gulf of Mexico, which would constitute the ideal baseline for this stock, were not available. A previous study that measured stable isotopes in the otoliths of Atlantic bluefin tuna collected from the Gulf of St Lawrence and adjacent areas over three decades showed that over 99% the fish in this area originate from nurseries in the western Atlantic  (Schloesser et al. 2010). Therefore, the samples collected from Canadian fishery in September 2013 were assumed to be representative of the west Atlantic stock.
3) and showed no spatial heterogeneity in the size-shape relationship (i. e. length Â location interaction P. 0. 05) were standardised using the common within-group slope as described above. A series of GLMs was used to identify which shape descriptors showed significant variation between the east Atlantic and west Atlantic baseline samples. These shape descriptors were included in a second PCA to reduce the dimensionality of the dataset. The principal components that captured the majority of the variation in otolith shape were used in a stepwise discriminant function analysis (DFA) to distinguish between fish from the western and eastern baseline samples.
Estimation of stock composition in mixed samples of unknown stock origin A mixed sample was selected from the dataset that included Atlantic bluefin. 200 cm FL from the east Atlantic (Portugal, Morocco and the Straits of Gibraltar) and the central Atlantic (Table  1 ). All images used were of right otoliths. The principal components that were selected in the DFA to distinguish between the two baseline groups were used in a Bayesian stock mixture analysis to estimate the proportion of fish in each sample that originated from the east and west Atlantic stocks. The analysis was conducted using the package mixFish in R as described by  Smith and Campana (2010). When using the Bayesian approach, the observations from the mixed samples classified to the base populations (the western and eastern stocks in this case) are used to update the parameter estimation of the base population (unconditional estimation).
Analysis of spatial, temporal and size-based sources of shape variability The PCA of the full dataset was based on 34 harmonics and three shape indices. The first 11 principal components explained 70% of the total variance in the otolith shape descriptors. Principal component (PC) 1, PC2, PC3 and PC4 were significantly correlated with length and were therefore standardised for the length effect using the common within-group slope. For the fish ,160 cm FL, there were small (R 2 values ranged from 1. 7 to 6. 5%) but significant (P , 0. 05) differences between sampling locations in the east Atlantic and Mediterranean for PC2, PC4, PC6, PC8, PC10 and PC11, whereas PC3 and PC5 showed significant differences between locations and sampling years. For the principal components that showed no interannual variation, further examination of the mean and 95% confidence limits (Fig.
Characterisation of Atlantic bluefin tuna for the east and west Atlantic using otolith shape In all, 25 elliptical Fourier coefficients and three shape indices showed significant variation between the east and west Atlantic (GLM P , 0. 05) and were not significantly correlated with length (in some cases after standardisation). The first 12 principal components from a PCA of these shape descriptors explained 84% of the variability in the dataset. The loading plot for the first two principal components (Fig. 5 ) highlights the shape descriptors that make major contributions to the variability in the dataset (e. g. roundness, aspect ratio, Elliptical Fourier coefficients D1, A11, C11, D9, A12 and B10). Six principal components were retained in the DFA by stepwise selection (PC1, PC2, PC5, PC6, PC8, PC11) producing one canonical function that distinguished between otoliths from east Atlantic and west Atlantic fish (Wilk's l ¼ 0.
The simultaneous collection of otolith shape and chemistry data is logistically very feasible because otoliths can be photographed without interfering with subsequent chemical analysis. It is recommended that this approach be taken in future sampling programs for Atlantic bluefin tuna. The estimates of stock composition in the mixed samples that were obtained using otolith shape descriptors are within the ranges reported by  Rooker et al. (2014)  for Atlantic bluefin collected from the same locations (albeit in different years and from a wider size range of fish). According to our estimates, the sample from the central north Atlantic collected in 2012 comprised 94 AE 7% (mean AE s. d. ) eastern origin fish. Although this is higher than the 63. 9 AE 9. 6% reported by  Rooker et al. (2014)  for the same location in 2010, it is consistent with their estimates of 90.
Acknowledgements This work was performed under the provision of the  ICCAT Atlantic Wide Research Program for Bluefin Tuna (GBYP) , funded by the  European Community  (Grant  SI2/542789 ),  Canada ,  Croatia, Japan, Norway, Turkey, US (NMFS NA11NMF4720107 ),  Chinese Taipei  and the  ICCAT Secretariat. The contents of the paper do not necessarily reflect the point of view of ICCAT or of the other funders. This work was also supported by grant  CTM2011-27505  from the  Spanish Ministry of Economy and Competitiveness (MINECO).