dataesr/scientific-paragraphs-categorization huggingface dataset:
18  of the original paper  20  ). Taken together, STREAM recovers the correct trajectories for the wild-type cells and, using the mapping feature, also predicts and effectively visualizes the consequences of the genetic perturbation as validated in the original study. Delineating hierarchies of the zebrafish hematopoiesis. To test the robustness and scalability of STREAM, we next explored data derived from different platforms and organisms. To this end, we used two recently published zebrafish datasets obtained with single-cell qPCR  32  and inDrop  33  assays. These data provided the first comprehensive model of the zebrafish hematopoiesis system without biases introduced by FACS sorting subpopulations. The first dataset from Moore et al. 32  provided a first model of the zebrafish hematopoiesis system using a carefully designed panel of 96 genes.
In this study, 166 cells were profiled from the wild-type (WT) whole-kidney marrow (WKM). STREAM analysis uncovered four cell lineages trajectories (Fig. 3a left ) and based on the automatic gene detection module, uncovered marker genes for each trajectory (Fig. 3a  right), which includes T cell marker gene TCR-alpha, B-cell marker gene CD79a, myeloid marker gene nccrp-1, and erythroid marker gene band3. Based on this analysis, we hypothesize that the inferred four branches correspond to T cell, B cell, myeloid, and erythroid lineages (Fig. 3a  middle). To test this hypothesis, we used the STREAM mapping feature to map fluorescent-labeled and FACS-sorted cells from WKM: 20 erythroid cells from peripheral blood cells (per RBC), 24 erythroid cells Tg(gata1:dsRed), 48 myeloid cells Tg(mpx:GFP), 49 marrow-derived B cells Tg (rag2:dsRed), 83 mature T cells Tg (lck:GFP)cells, 85 HSPCs Tg(CD41:GFP)low.
lck + cells were mapped to the T cell branch, mpx + cells were mapped to the myeloid cell branch, rag2 + cells were mapped to the B cell branch, both gata1+ and peripheral RBCs were mapped to the erythroid branch, while the majority of HSPCs were mapped to the proposed starting state as expected (Fig. 3b ). This result provides additional support for our hypothesis that the four branches corresponded to well defined lineages, highlighting the utility of the STREAM mapping feature. To test the scalability and robustness of STREAM on a larger and more challenging scRNA-seq dataset, we next analyzed 9628 unlabeled cells from the zebrafish whole-kidney marrow generated by Tang et al. 33  using the inDrop protocol  2. The original study, based on dimensionality reduction and clustering, uncovered and annotated 10 different and imbalanced subpopulations (some of which were validated by the authors using sorting of fluorescent transgenic cell sub-populations) (Fig.
3c ). STREAM correctly recapitulated the hierarchy of the different lineages and unbiasedly recovered four main hematopoietic cellular trajectories: starting from HSCs, through blood progenitor cells, cells differentiate into erythroid, macrophage, neutrophil, and lymphoid lineages (Fig. 3d ). Importantly, we rediscovered well-known marker genes: hbaa1 for the erythroid branch, grna for the macrophage branch, mpx for the neutrophil branch, and igl3c3 for the lymphoid branch (Fig. 3e ). However, we noticed that B and T cells were not separated and were assigned to the same lineage branch. Therefore, we derived an improved seeding strategy that is well suited to learn complex trajectories in high dimensions and that well recapitulates the known lineage for this dataset as presented in Supplementary Note 2 and Supplementary Figs.
4 5 6. This new strategy is generalizable to other datasets and described in detail in the method section. In summary, these analyses highlight some important points of our approach: (1) STREAM is able to identify more refined trajectories increasing the number of dimensions, (2) we can recover trajectories using unsorted populations, (3) the trajectory inference is robust to subpopulation imbalance, (4) our gene expression analysis is a powerful tool to discover marker genes, and (5) our method is scalable to currently available large-scale single-cell assays. Comparison with other methods. Several methods have been proposed for pseudotime inference or trajectory reconstructions. In fact, more than 50 methods have been proposed for this task, making a systematic comparison unfeasible for the scope of this manuscript.
For each method, the analyses were performed with standard parameters when possible (following the guidelines provided in the documentation) otherwise the parameters were obtained by contacting the respective authors (Supplementary Note 3). To evaluate the ability of each method in recovering the correct topology we used a proposed gold-standard synthetic dataset by Rizvi et al. 34  with known topology and pseudotime: two bifurcation events and three different time points (Fig. 4a ). First, we started by quantitatively evaluating the number of correctly detected branching events and the pseudotime accuracy (Online methods). For topology correctness, STREAM and other five methods, including scTDA, SLICER, Monocle2, Mpath, GPfates, PHATE, successfully identified two bifurcation events. Second, for pseudotime accuracy, we calculated for each method the correlation between true pseudotime and inferred pseudotime as proposed before  34.
Four metrics were used to evaluate this correlation, including rank-based Pearson correlation, distancebased Pearson correlation, Spearman's rank correlation, and Kendall's tau coefficient. We use four different correlation metrics since some methods (scTDA, TSCAN, and Mpath) only return a simple ordering, i. e. , the ranks of cells, and do not provide the actual pseudotime defined as the distance of each cell from the origin in the proposed embedding (for Spearman's rank correlation and Kendall's tau coefficient, ranks-based and distance-based correlations are the same). STREAM has the best performance for two out of four metrics (and importantly when using distance-based pseudotime) and second-best performance for the other two rank-based metrics (following scTDA in which this synthetic dataset was proposed) (Fig.
4b ). Finally, we assessed the qualitative output of each method using their proposed visualization. STREAM is the only tool that provides a densitylevel visualization to study the composition of different cell types in different branches. (Fig. 4c , Supplementary Note 3). To compare the different methods on real datasets, we first used one of the most commonly used scRNA-seq datasets for this task, originally generated by Trapnell et al. 39. This dataset contains human skeletal muscle myoblasts (HSMM) cells differentiating along a linear trajectory. In this analysis, we were able to evaluate only methods capable to detect the correct bifurcation event. Regardless, the visual outputs of all the methods are presented for completeness (Supplementary Fig. 7 ). The original study proposed a single bifurcation, which leads to myoblast cells or separate potentially contaminating cells (Fig.
8 ). In Wishbone and SLICER, distance-based pseudotime shows a set of unexpected discrete segments. Neither Mpath nor TSCAN can generate distance-based pseudotime. In addition, Mpath does not recover a monotonically increasing trend (Supplementary Fig. 8 ). STREAM has also the highest average coefficient on ENO3 based on the four different metrics (Fig. 5c ). When combining all three marker genes, STREAM has the overall best performance (calculated as the average rank for the four proposed metrics) (Fig. 5d ). Finally, we analyzed a high-quality single-cell qPCR dataset containing ~270 blood cells sorted from six different populations: HSC, MPP, CMP, GMP, MEP, and common lymphoid progenitor cells (CLPs) profiled for ~170 key transcription factors important in mouse hematopoiesis  40. The output of each method is shown in Supplementary Fig.
These assays allow to study an important additional layer in gene regulation. In fact, accessible sites in the genome can be used to profile the activity of important cell type specific non-coding regulatory elements such as enhancers. These regions harbor transcription factor binding sites and can control the activity of relatively proximal genes via long-range interactions  7. The analysis of human scATAC-seq data is particularly challenging because sparsity is intrinsic to these assays. In fact, the signal is limited by the DNA copy number, which only consists of 0, 1, or 2 reads within a diploid genome. In addition, compared to previously published datasets in other model organisms (mouse, Drosophila melanogaster), the human genome is larger (respectively, 1. 2X and 27X) making this problem even more accentuated.
Although some initial efforts have been made to adapt existing trajectory inference methods originally developed for transcriptomic data to scATAC-seq data analysis  42, 43  , to our knowledge STREAM is the only documented end-to-end pipeline that provides to users the specific functions to analyze and visualize scATAC-seq data starting from raw count data and based on an unbiased approach to model important DNA sequence features associated with chromatin-accessibility. We present below how STREAM can be used to infer trajectories from single-cell epigenomic data and show its application to a recently published dataset, where a total of 3072 cells were profiled from the human bone marrow and isolated by FACS into nine different cellular populations, including HSC, MPP, CMP, CLP, LMPP, GMP, MEP, mono and plasmacytoid dendritic cells (pDCs)  44.
Discussion Large-scale single-cell measurements have opened unprecedented opportunities to study dynamic processes such as differentiation or response to stimuli. Trajectory inference methods are important tools to describe those processes based on snapshot measurements. In this study we presented STREAM, a trajectory inference tool capable of reliably reconstructing trajectories and inferring pseudotime from different data types and when multiple branching points are present. Our method does not require prior knowledge such as time points, start cell, or the number of branching events to reconstruct trajectories, and does not need extensive bioinformatics knowledge thanks to a user-friendly and interactive web interface. Additionally STREAM introduces four innovations compared to other existing methods: (1) a densitylevel trajectory visualization useful to study subpopulation composition and cell-fate genes along branching trajectories, (2) a documented end-to-end pipeline to reconstruct trajectories from chromatin-accessibility data, (3) an interactive website that can be used to not only compute trajectories but also host a database to readily visualize and explore precomputed trajectories for several published studies  [30] [31] [32] [33] 39, 44  , and (4) a trajectory mapping procedure to readily map new cells to precomputed structures without pooling data and re-computing trajectories.
STREAM has been extensively tested using several published datasets from different organisms (zebrafish, mouse, human) and single-cell technologies (qPCR, scRNA-seq, scATAC-seq). Our comparison of STREAM with 10 other methods shows that it is among the top scoring methods on both synthetic and real datasets (best average ranking across several metrics and most balanced method for branching detection). However, there are still several general limitations for the current trajectory inference methods, including STREAM. First, most of the methods are still based on linear or tree-structure models even though we know that developmental processes may involve more sophisticated topologies, e. g. , cyclic or disconnected graph. Therefore, further improvement is necessary to adapt STREAM to properly describing these more general structures.
MLLE addresses the regularization problem of standard LLE by introducing multiple weight vectors in each neighborhood rather than a single one, which makes it more stable and robust. MLLE minimizes the following embedding cost function: min X N i¼1 X s i l¼1 t i À X j2J i w ðlÞ ji t j 2 ; s:t:TT T ¼ I ð3Þ w ð1Þ i ; ; w s i i (s i ≤ k i ) are linearly independent weight vectors, which are obtained from a matrix of the difference between data point x i and its k i neighbors. In STREAM, the neighborhood size is chosen based on the number of cells and is set by default to 10% of the total number of cells. The number of MLLE components depends on the number of branches and on the complexity of the structure to learn. Typically, three components capture the main structure for most datasets (and this space can be readily visualized); however, increasing them may recover finer structures as discussed in main text under the section "STREAM trajectory inference in high dimensional spaces".
g. , in the set of all possible trees). Finally, the use of ElPiGraph results in an explicit definition of the principal tree embedded into the data space. The explicit tree structure can be studied independently on the data, or the data can be mapped onto the principal tree and studied in its intrinsic, geodesic, coordinates. A detailed description of ElPiGraph and related elastic principal graph approaches is available elsewhere  [27] [28] [29]. Concerning the choice of the parameter values, in STREAM these default values are used for the principal graph inference with ElPiGraph: R 0 = ∞, α = 0. 02, μ = 0. 1, λ = 0. 02. These values worked well across all the single-cell datasets tested. However, if necessary, these parameters can be easily modified by the user directly from the STREAM package.
The ElPiGraph. R R package, used as part of STREAM to fit principal trees to the data, is available at https://github. com/sysbio-curie/ ElPiGraph. R. Domain-specific optimizations to model single-cell data. Although ElPiGraph is a general approach to construct principal trees (and other topologies), the obtained structures may not optimally describe biologically relevant trajectories or accurately capture pseudotime information based on single-cell data. Therefore, in addition to the described seeding strategy, in STREAM several single-cell specific optimizations were introduced to the core algorithm of ElPiGraph (Supplementary Fig. 1 ): Control over-branching: A regularization parameter α with range (0,1] was introduced to explicitly control the complexity of the resulting graph structure.
Diverging gene detection. For each pair of branches B i and B j , and for the gene E, the gene expression values across cells from both branches are scaled to the range [0,1]. For gene expression E i fromB i and gene expression E j from B j , we first calculate their mean values. Then, we check the fold change between mean values to make sure it is above a specified threshold (the default log2 fold change value is >0. 25). Mann-Whitney U test is then used to test whether E i is greater than E j or E i is less than E j. Since the statisticU could be approximated by a normal distribution for large samples, and U depends on specific datasets, we standardize Uto Z-score to make it comparable between different datasets. For small samples where this test is underpowered (<20 cells per branch), we report only the fold change to qualitatively evaluate the difference between E i and E j.
STREAM analysis on scATAC-seq data. For the scATAC-seq analysis, a total of 3,072 cells were profiled using FACS to isolate 9 distinct populations from CD34+ human bone marrow, which encompassed progenitors for four welldefined lineages  44. Two thousand thirty-four high-quality cells passed quality control filtering and were used in the downstream analysis with STREAM. Specifically, cells were filtered so that 1000 unique nuclear fragments were observed for each cell and at least 60% of these reads aligned in open chromatin peaks. After filtering low quality cells, the mean intensity and GC content for each peak that was called for this dataset was computed using the addGCBias function for the hg19 genome using the BSgenome. Hsapiens. UCSC. hg19 package available through chromVAR  45. These two coordinates were used to infer an empirically-defined set of background peaks to compute accessibility deviations, which have been described elsewhere  7, 44.
As features we used an unbiased k-mer scoring, which is agnostic to any known transcription factor motif and thus generalizable to other systems. We used the matchKmers function in chromVAR with parameters k = 7 and genome = BSgenome. Hsapiens. UCSC. hg19, which returns a matrix of dimension number of peaks by number of k-mers where a 1 indicates that the peak contains the k-mer sequence. The output of this function was then included in the computeDeviations function to compute chromatinaccessibility z-scores for each of the k-mers in our dataset. This matrix of cells by k-mer accessibility z-scores serves as a data-driven dimensionality reduction of the chromatin-accessibility profiles of these cells. Based on the z-score matrix of k-mer DNA sequences, all the 7-mer features are standardized to have zero mean and unit variance.
Since the z-score matrix of k-mers has both positive and negative values, the variable gene selection method based on gene expression values is not directly applicable. As such, PCA is performed on the scaled matrix to convert z-score to principal components. According to the variance ratio elbow plot we selected the top 15 PCs but excluded the first component since it captured technical noise. Then the selected PCs are used as features for MLLE to reduce dimensionality. In the reduced MLLE space, the same strategy is used to reconstruct trajectories as previously described. Diverging and transition k-mers were selected with the same procedures used for gene selection. Finally, detected k-mers were mapped to known transcription factors using Tomtom 58 (http:// meme-suite. org/tools/tomtom) and a motif database previously assembled (github.
com/buenrostrolab/chromVARmotifs)  44. Comparison on simulated datasets. Given a set of n cells and assuming we know their developmental/sampling time and topological organization, i. e. , how they are organized in branches, we can easily evaluate a generic reconstruction method with the following two metrics: Difference between the number of inferred and true branches. Correlation between the true sampling time X and the inferred pseudotime Y. For the pseudotime we use either the proposed ranking or the actual distance from the starting point as provided by each method. We used three different measure of correlation: Pearson correlation r, Spearman correlation ρ, and Kendall's tau correlation τ, calculated as follow: r pseudotime ¼ P n i¼1 x i À x ð Þy i À y ð Þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P n i¼1 x i À x ð Þ 2 q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi P n i¼1 y i À y ð Þ 2 q ð9Þ r rank ¼ P n i¼1 x i À x ð Þ rg y i À rg y ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P n i¼1 x i À x ð Þ 2 q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi P n i¼1 rg y i À rg y 2 r ð10Þ ρ ¼ cov rg X ; rg Y ð Þ σ rgX σ rgY ð11Þ τ ¼ 1 n n À 1 ð Þ X i≠j sgn x i À x j sgn y i À y j ð12Þ Where rg X and rg Y are the ranks of cells, cov (rg X , rg Y ) is the covariance of rank variables, σ rgx and σ rgy are the standard deviations of rank variables.
Note that since both Spearman correlation ρ and Kendall's tau correlation τ are rank-based methods, the correlation between X and Y and the correlation between X and rg Y are the same, so we consider only the correlation between X and Y. Comparison on real datasets. To evaluate the quality of reconstruction in real datasets in which we do not have the real developmental time and topological information, we used the following two metrics: (1) Path-specific marker gene correlation analysis: In real datasets oftentimes, we do not have the sampling time along a branch. In this case, instead, it is helpful to evaluate how the inferred pseudotime recapitulates the progressive activation or repression of an important gene along that branch. The main idea here is that ordering cells based on a marker gene, which is important in defining a developmental trajectory, as a reasonable surrogate for the correct pseudotime ordering.
Acknowledgements This project has been made possible in part by grant number  2018-182734  to L. P. from the  Chan Zuckerberg Initiative DAF , an advised fund of  Silicon Valley Community Foundation. L. P. is also partially supported by a  National Human Genome Research Institute  ( NHGRI) Career Development Award  ( R00HG008399 ). G. Y. 's research was supported by a  Claudia Adams Barr Award  and a  Bridge Award. J. D. B. acknowledges support from the  Harvard Society of Fellows and Broad Institute Fellowship. J. D. B. also acknowledges the  Allen Distinguished Investigator Program , through The  Paul G. Allen Frontiers Group  for funding. A. Z. and L. A. were supported by  ITMO Cancer SysBio program (MOSAIC)  and  INCa   PLBIO program  (CALYS,  INCA_11692 ). D. M. L was supported by  R24OD016761  and  R01CA211734.
Code availability STREAM is available as a user-friendly open-source software and can be used interactively as a web-application at http://stream. pinellolab. org (Supplementary Fig. 11 , Supplementary Note 4), a bioconda package 'stream' for step-by-step analysis https:// bioconda. github. io/recipes/stream/README. html (Supplementary Note 5), or as a standalone command-line tool: https://github. com/pinellolab/STREAM (Supplementary Note 6). All the analyses presented in this manuscript can be reproduced using the bioconda package and the provided Jupyter notebooks in Supplementary Data 1 and 2.
Blood collection and laboratory analysis A standardized protocol for blood collection and processing was followed and fasting was not required (details published elsewhere  15  ). All plasma samples (citrate anticoagulant) were assayed at IARC, using the AbsoluteIDQ p180 Kit (Biocrates Life Sciences AG, Innsbruck, Austria) and following the procedure recommended by the vendor. Samples from matched case-control sets were assayed in the same analytical batch along with quality control samples from pooled plasma; laboratory personnel was blinded to sample category, that is, case, study control or quality control. A triple quadrupole mass spectrometer (Triple Quad 4500; AB Sciex, Framingham, MA) was used to quantify 148 metabolites. After exclusions of metabolites and men, 3,057 matched case-control sets had data for 119 metabolites (details in Supporting Information Methods and Table  S1 ).
For the purpose of comparison with previous studies, we also computed the associations of each individual metabolite and risk of prostate cancer, while accounting for multiple testing. This analysis was done for the full data set (n case = 3,057), and separately for the data not previously published  6  (n case = 2,018). All tests of statistical significance were two-sided, and 0. 05 was considered as the nominal level for statistical significance in analyses of metabolite patterns and prostate cancer risk. All analyses were conducted in Stata Statistical Software Package, version 15 (Stata Corporation, College Station, TX). The "tt" add-on for STATA is available from http://www. gorst. dk. See Supporting Information Methods for further details of statistical analyses.
Analysis of individual metabolites Results of the associations between individual metabolite concentrations and prostate cancer risk in the full dataset (n case = 3,057) and in the data not previously published  6  (n case = 2,018) are shown in Supporting Information Tables S10-S23. Results from the new dataset were in line with those of the full dataset, although slightly weaker, with two obvious exceptions; acylalkyl phosphatidylcholines were more often associated with overall prostate cancer in the new dataset than in the full dataset (Supporting Information Tables  S10  and S17 ), and the same was the case for diacyl phosphatidylcholines and risk of lowintermediate grade prostate cancer (Supporting Information Tables  S15  and S22 ).