dataesr/scientific-paragraphs-categorization huggingface dataset:
ZINB-WaVE is more accurate than state-of-the-art methods We next compared ZINB-WaVE to PCA and ZIFA, in terms of their ability to recover the true underlying low-dimensional signal and clustering structure. For datasets simulated from a ZINB model, our estimation procedure with correctly specified K led to almost perfect correlation between distances in the true and estimated low-dimensional space (Fig. 7a, b ). The correlation remained high even for misspecified K, in most cases higher than for PCA and ZIFA. ZINB-WaVE performed consistently well across a range of simulation scenarios, including different numbers of cells n, different zero fractions, and varying cluster tightness (Supplementary Fig. S28 ). We observed a consistent ranking, although noisier, when the methods were compared in terms of silhouette widths (Fig.
7c, d ). As with the real datasets, the choice of normalization influenced the simulation results. Critically, there was not an overall best normalization method; rather, the performance of the normalization methods depended on the dimensionality reduction method and on intrinsic characteristics of the data, such as the fraction of zero counts and the number and tightness of the clusters (Fig. 7  and Supplementary Fig. S28 ). For instance, TC normalization worked best for PCA on data simulated from the V1 dataset (Fig. 7a, c ), while FQ and TMM normalization worked best for PCA and ZIFA, respectively, on data simulated from the S1/CA1 dataset (Fig. 7b, d ). It is important to note that the previous results were obtained for data simulated from the ZINB-WaVE model underlying our estimation procedure.
It is hence not surprising that ZINB-WaVE outperformed its competitors. To provide a fairer comparison, we also assessed the methods on data simulated from the model For each gene and cell, bias and MSE were averaged over B = 10 datasets simulated from our ZINB-WaVE model, based on the S1/CA1 dataset and with n = 1, 000 cells, J = 1, 000 genes, scaling of one for the ratio of within to between-cluster sums of squares (b 2 = 1), K = 2 unknown factors, zero fraction of about 80%, X = 1 n , cell-level intercept (V = 1 J ), and genewise dispersion. The lower and upper hinges of the boxplots correspond to the first and third quartiles (the 25th and 75th percentiles), respectively. The upper/lower whisker extends from the hinge to the largest/smallest value no further than 1. 5×IQR from the hinge (where IQR is the inter-quartile range or difference between the third and first quartiles).
However, there is no guarantee that the low-dimensional signal extracted by ZINB-WaVE is biologically relevant: If unwanted technical variation affects the data and is not accounted for in the model (or in prior normalization), the low-rank matrix W inferred by ZINB-WaVE will capture such confounding effects. It is therefore important to explore the correlation between the latent factors estimated by our procedure and known measures of quality control that can be computed for scRNA-seq libraries, using, for instance, the Bioconductor R package scater  [43]  (see Fig. 2f ). If one observes high correlation between one or more latent factors and some QC measures, it may be beneficial to include these QC measures as covariates in the model. Several authors have recognized that high-dimensional genomic data are affected by a variety of unwanted technical effects (e.
Such a low-dimensional representation can be used in downstream analyses, such as clustering or pseudotime ordering of the cells  [15]. Visualization of high-dimensional datasets is an equally important area of research and many algorithms are available, among which t-SNE  [25]  has become the most popular for scRNA-seq data. Recently, Wang et al. [48]  have proposed a novel visualization algorithm that can account for zero inflation and showed improvement over t-SNE. As t-SNE takes as input a matrix of cell pairwise distances, which may be noisy in high dimensions, a typical pipeline involves computing such distances in PCA space, selecting, for example, the first 50 PCs. An alternative approach is to derive such distances from the low-dimensional space defined by the factors inferred by ZINB-WaVE.
This strategy was used effectively in Figure  3c  to visualize the PBMC dataset. In this article, we have focused on an unsupervised setting, where the goal is to extract a low-dimensional signal from noisy zero-inflated data. However, our proposed ZINB model is more general and can be used, in principle, for supervised differential expression analysis, where the parameters of interest are regression coefficients β corresponding to known sample-level covariates in the matrix X (e. g. , cell type, treatment/control status). Differentially expressed genes may be identified via likelihood ratio tests or Wald tests, with standard errors of estimators of β obtained from the Hessian matrix of the likelihood function. Additionally, posterior dropout probabilities can be readily derived from the model and used as weights to unlock standard bulk RNA-seq methods  [49] , such as edgeR  [41].
Real datasets V1 dataset. Tasic et al. [3]  characterized more than 1,600 cells from the primary visual cortex (V1) in adult male mice, using a set of established Cre lines. Single cells were isolated by FACS into 96-well plates and RNA was reverse transcribed and amplified using the SMARTer kit. Sequencing was performed using the Illumina HiSeq platform, yielding 100 bp-long reads. We selected a subset of three Cre lines, Ntsr1-Cre, Rbp4-Cre, and Scnn1a-Tg3-Cre, that label Layer 4, Layer 5, and Layer 6 excitatory neurons, respectively. This subset consists of 379 cells, grouped by the authors into 17 clusters; we excluded the cells that did not pass the authors' quality control filters and that were classified by the authors as "intermediate" cells between two clusters, retaining a total of 285 cells.
Gene expression was quantified by gene-level read counts. Raw gene-level read counts and QC metrics (see below) are available as part of the scRNAseq Bioconductor R package (https://bioconductor. org/packages/scRNAseq). We applied the dimensionality reduction methods to the 1,000 most variable genes. S1/CA1 dataset. Zeisel et al. [4]  characterized 3,005 cells from the primary somatosensory cortex (S1) and the hippocampal CA1 region, using the Fluidigm C1 microfluidics cell capture platform followed by Illumina sequencing. Gene expression was quantified by unique molecular identifier (UMI) counts. In addition to gene expression measures, we have access to metadata that can be used to assess the methods: batch, sex, number of mRNA molecules. Raw UMI counts and metadata were downloaded from http://linnarssonlab.
org/ cortex/. mESC dataset. Kolodziejczyk et al. [37]  sequenced the transcriptome of 704 mouse embryonic stem cells (mESCs), across three culture conditions (serum, 2i, and a2i), using the Fluidigm C1 microfluidics cell capture platform followed by Illumina sequencing. We selected only the cells from the second and third batch, after excluding the samples that did not pass the authors' QC filtering. This allowed us to have cells from each culture condition in each batch and resulted in a total of 169 serum cells, 141 2i cells, and 159 a2i cells. In addition to gene expression measures, we have access to batch and plate information that can be included as covariates in our model. Raw gene-level read counts were downloaded from http: //www. ebi. ac. uk/teichmann-srv/espresso/. Batch and plate information was extracted from the sample names, as done in Lun & Marioni  [42].
We applied the dimensionality reduction methods to the 1,000 most variable genes. Glioblastoma dataset. Patel et al. [6]  collected 672 cells from five dissociated human glioblastomas. Transcriptional profiles were generated using the SMART-Seq protocol. We analyzed only the cells that passed the authors' QC filtering. The raw data were downloaded from the NCBI GEO database (accession GSE57872). Reads were aligned using TopHat with the following parameters: -rg-library Illumina -rgplatform Illumina -keep-fasta-order -G -N 3 -read-edit-dist 3 -no-coverage-search -x 1 -M -p 12. Counts were obtained using htseq-count with the following parameters (http://www-huber. embl. de/HTSeq/doc/ count. html): -a 10 -q -s no -m union. We applied the dimensionality reduction methods to the 1,000 most variable genes.
OE dataset. Fletcher et al. [29]  characterized 849 FACS-purified cells from the mouse olfactory epithelium (OE), using the Fluidigm C1 microfluidics cell capture platform followed by Illumina sequencing. Gene-level read counts were downloaded from GEO (GSE95601; file GSE95601_oeHBCdiff_Cufflinks_eSet_counts_table. txt. gz). As done in Perraudeau et al. [31] , we filtered the cells that exhibited poor sample quality using SCONE  [36]  (v. 1. 1. 2). A total of 747 cells passed this filtering procedure. To compare with the original results, we also re-analyze the final repertoire of 13 stable clusters found in Fletcher et al. [29] , consisting of 616 cells, downloaded from https://github. com/rufletch/p63-HBC-diff. See Fletcher et al. [29]  for details on the original analysis and Perraudeau et al.
[31]  for details on the ZINB-WaVE based workflow. 10x Genomics 68k PBMCs dataset. Zheng et al. [32]  characterized 68,579 peripheral blood mononuclear cells (PBMCs) from a healthy donor. The gene-level UMI counts were downloaded from https: //www. 10xgenomics. com/single-cell/ using the cellrangerRkit R package (Version 1. 1. 0). We applied the dimensionality reduction methods to the 1,000 most variable genes.
Simulating from the ZINB-WaVE model Setting the simulation parameters. In order to simulate realistic data, we fitted our ZINB-WaVE model to two real datasets (V1 and S1/CA1) and used the resulting parameter estimates as the truth to be estimated in the simulation. Genes that did not have at least 5 reads in at least 5 cells were filtered out and J = 1, 000 genes were then sampled at random for each dataset. The ZINB-WaVE model was fit to the count matrix Y with the number of unknown cell-level covariates set to K = 2, genewise dispersion (ζ = ln θ = -ln φ), X µ , X π , V µ , and V π as columns of ones (i. e. , intercept only), and no offset matrices, to get estimates for W , α µ , α π , β µ , β π , γ µ , γ π , and ζ. The parameters which were varied in the simulations are as follows. • The number of cells: n = 100, 1, 000, 10, 000.
Then, 2 the within and between-cluster sums of squares, respectively. The level of difficulty of the clustering problem can be controlled by the ratio of within to between-cluster sums of squares. However, we want to keep the overall mean W and overall variance (i. e. , T SS) constant, so that the simulated values of W stay in the same range as the estimated W from the real dataset; this prevents us from simulating an unrealistic count matrix Y. T SS = C c=1 nc i=1 (W ic -W ) 2 = C c=1 nc i=1 (W ic -Wc ) 2 + C c=1 n c ( Wc -W ) 2 = W SS + BSS, with W SS = C c=1 nc i=1 (W ic -Wc ) 2 and BSS = C c=1 n c ( Wc -W ) Let us scale the between-cluster sum of squares by a 2 and the within-cluster sum of squares by a 2 b 2 , i. e. , replace ( Wc -W ) by a( Wc -W ) and (W ic -Wc ) by ab(W ic -Wc ), with a ≥ 0 and b ≥ 0 such that T SS and W are constant.
The total sum of squares T SS remains constant, i. e. , T SS = a 2 b 2 W SS + a 2 BSS, provided a 2 = T SS b 2 W SS + BSS. Requiring the overall mean W to remain constant implies that W * c = (1 -a) W + a Wc , where the * superscript refers to the transformed W. Thus, W * ic = W * c + ab(W ic -Wc ) = (1 -a) W + a Wc + ab(W ic -Wc ) = (1 -a) W + a(1 -b) Wc + abW ic. ( 13 ) The above transformation results in a scaling of the ratio of within to between-cluster sums of squares by b 2 , while keeping the overall mean and variance constant. In our simulations, we fixed C = 3 clusters and considered three values for b 2 , where the same value of b 2 is applied to each of the K columns of W : b 2 = 1, corresponding to the case where the within and between-cluster sums of squares are the same as the ones of the fitted W from the real dataset; b 2 = 5, corresponding to a larger ratio of within to between-cluster sums of squares and hence a harder clustering problem; b 2 = 10, corresponding to a very large ratio of within to between-cluster sums of squares and hence almost no clustering.
Overall, 2 (real datasets) × 3 (n) × 3 (zf rac) × 3 (clustering) = 54 scenarios were considered in the simulation. Simulating the datasets. For each of the 54 scenarios, we simulated B = 10 datasets, resulting in a total of 54 × 10 = 540 datasets. Using the fitted W , α µ , α π , β µ , β π , γ µ , γ π , and ζ from one of the two real datasets, the datasets were simulated according to the following steps. 1. Simulate W with desired clustering strength. First fit a K-variate Gaussian mixture distribution to W inferred from one of the real datasets using the R function Mclust from the mclust package and specifying the number of clusters C. Then, for each of B = 10 datasets, simulate W cluster by cluster from K-variate Gaussian distributions using the mvrnorm function from the MASS package, with the cluster means, covariance matrices, and frequencies output by Mclust.
Transform W as in Equation (  13 ) to get the desired ratio of within to between-cluster sums of squares. 2. Simulate γ µ and γ π to get the desired zero fraction. We only considered cell-level intercept n-vectors γ µ and γ π , i. e. , L = 1 and a matrix V of gene-level covariates consisting of a single column of ones. As the fitted γ µ and γ π from the original datasets are correlated n-vectors, fit a bivariate Gaussian distribution to γ µ and γ π using the function Mclust from the mclust package with C = 1 cluster. Then, for each of B = 10 datasets, simulate γ µ and γ π from a bivariate Gaussian distribution using the mvrnorm function from the MASS package, with the mean and covariance matrix output by Mclust. To increase/decrease the zero fraction, increase/decrease each element of the mean vector for γ π inferred from Mclust (shifts of {0, 2, 5} for V1 dataset and {-1.
Simulating from the Lun & Marioni [42] model To simulate datasets from a different model than our ZINB-WaVE model, we simulated counts using the procedure described in Lun & Marioni  [42]  (details in Supplementary Materials of original publication and code available from the Github repository https://github. com/MarioniLab/PlateEffects2016). Although the Lun & Marioni  [42]  model is also based on a zero-inflated negative binomial distribution, the distribution is parameterized differently and also fit differently, gene by gene. In particular, the negative binomial mean is parameterized as a product of the expression level of interest and two nuisance technical effects, a gene-level effect assumed to have a log-normal distribution and a cell-level effect (cf. library size) whose distribution is empirically derived.
The zero inflation probability is assumed to be constant across cells for each gene and is estimated independently of the negative binomial mean. The ZINB distribution is fit gene by gene using the zeroinfl function from the R package pscl. We used the raw gene-level read counts for the mESC dataset as input to the script reference/submitter. sh, to create a simulation function constructed by fitting a ZINB distribution to these counts. The script simulations/submitter. sh was then run to simulate counts based on the estimated parameters (negative binomial mean, zero inflation probability, and genewise dispersion parameter). We simulated C = 3 clusters with equal number of cells per cluster. The parameters which were varied in the simulations are as follows. • The number of cells: n = 100, 1, 000, 10, 000.
• The proportion of zero counts, zf rac = i,j 1(Y ij = 0)/nJ, via the zero inflation probability: zf rac ≈ 0. 4, 0. 6, 0. 8. For zf rac = 0. 4, we did not modify the code in Lun & Marioni  [42]. However, to simulate datasets with greater zero fractions, namely zf rac = 0. 6 and zf rac = 0. 8, we added respectively 0. 3 and 0. 6 to the zero inflation probability (p i , in their notation). For each of the 3 (n) × 3 (zf rac) = 9 scenarios, we simulated B (B = 10) datasets, resulting in 9 × 10 = 90 simulated datasets in total.
ZINB-WaVE We applied the ZINB-WaVE procedure using the function zinbFit from our R package zinbwave, with the following parameter choices. • Number of unknown cell-level covariates K: K = 1, 2, • Gene-level covariate matrix V : not included or set to a column of ones 1 J. • Cell-level covariate matrix X: set to a column of ones 1 n. For the mESC dataset, we also considered including batch covariates in X. • Dispersion parameter ζ: the same for all genes (common dispersion) or specific to each gene (genewise dispersion).
Zero-inflated factor analysis We used the zero-inflated factor analysis (ZIFA) method  [26] , as implemented in the ZIFA python package (Version 0. 1) available at https://github. com/epierson9/ZIFA, with the block algorithm (function block_ZIFA. fitModel, with default parameters). The output of ZIFA is an n × K matrix corresponding to a projection of the counts onto a latent low-dimensional space of dimension K.
Normalization methods As normalization is essential, especially for zero-inflated distributions, PCA and ZIFA were applied to both raw and normalized counts. The following normalization methods were used. • Total-count normalization (TC). Counts are divided by the total number of reads sequenced for each sample and multiplied by the mean total count across all the samples. This method is related to the popular transcripts per million (TPM)  [52]  and fragments per kilobase million (FPKM)  [53]  methods. • Full-quantile normalization (FQ)  [54]. The quantiles of the distributions of the gene-level read counts are matched across samples. We used the function betweenLaneNormalization from the Bioconductor R package EDASeq. • Trimmed mean of M values (TMM)  [55]. The TMM global-scaling factor is computed as the weighted mean of log-ratios between each sample and a reference sample.
Clustering methods Clustering of the OE dataset. We used the resampling-based sequential ensemble clustering (RSEC) framework implemented in the RSEC function from the Bioconductor R package clusterExperiment  [56]. Briefly, RSEC implements a consensus clustering algorithm which generates and aggregates a collection of clusterings, based on resampling cells and using a sequential tight clustering algorithm  [57]. See Fletcher et al. [29]  for details on the parameters used for RSEC in the original analysis and Perraudeau et al. [31]  for details on the parameters used for RSEC in the ZINB-WaVE workflow. Clustering of the 10x Genomics 68k PBMCs dataset. We used two different clustering methods to assess the ability of ZINB-WaVE to extract biologically meaningful signal from the data. First, we used a a clustering procedure similar to the one implemented in the R package Seurat [34] (Version 2.
0. 1). In particular, we used ZINB-WaVE (K = 10) instead of PCA for dimensionality reduction. The clustering is based on shared nearest neighbor modularity  [33]. We used the following parameters: k. param = 10, k. scale=10, resolution = 0. 6. All the other parameters were left at their default values. We also clustered the data using a sequential k-means clustering approach  [57] , implemented in the clusterSingle function of the clusterExperiment package  [56]. We used the following parameters: sequential = TRUE, subsample = FALSE, k0 = 15, beta = 0. 95, clusterFunction = "kmeans".
Models for count data We compared the goodness-of-fit of our ZINB-WaVE model to that of two other models for count data: a standard negative binomial model, that does not account for zero inflation, and the MAST hurdle model, that is specifically designed for scRNA-seq data  [23]. Goodness-of-fit was assessed on the V1 dataset using mean-difference plots (MD-plots) of estimated vs. observed mean count and zero probability, as well as plots of the estimated dispersion parameter against the observed zero frequency.
ZINB-WaVE model For our ZINB-WaVE model, the overall mean and zero probability are E[Y ij ] = (1 -π ij )µ ij , P (Y ij = 0) = π ij + (1 -π ij )(1 + φ j µ ij ) 1 φ j. ZINB-WaVE was fit to the V1 dataset using the function zinbFit from our R package zinbwave, with the following parameter choices: K = 0 unknown cell-level covariates, gene-level intercept (X = 1 n ), cell-level intercept (V = 1 J ), and genewise dispersion.
Negative binomial model The negative binomial (NB) distribution is a special case of the ZINB distribution that does not account for zero inflation, i. e. , for which π = 0. Thus, E[Y ij ] = µ ij , P (Y ij = 0) = (1 + φ j µ ij ) 1 φ j. A NB distribution was fit gene by gene to the V1 dataset, after full-quantile normalization, using the Bioconductor R package edgeR  [41]  (Version 3. 16. 5) with only an intercept (i. e. , default value for the design argument as a single column of ones) and genewise dispersion.
MAST is implemented in the Bioconductor R package MAST  [58]  (which allows different covariate vectors X i for the continuous and discrete components). We use the function zlm to fit MAST with an intercept and a covariate for the cellular detection rate (as recommended in the MAST vignette for the MAIT data analysis) for both the discrete and continuous parts. Note that in contrast to the NB and ZINB models, the MAST hurdle model is for log2(TPM+1) instead of counts and has no dispersion parameter. To be able to compare the fits of the three models, the MAST goodness-of-fit plots display estimated vs. observed mean log2(TPM+1) and estimated variance σ 2 j vs. observed zero frequency (Supplementary Fig. S21 ). Although not a direct comparison, we think this allows a fair assessment of the goodness-of-fit of the different models.
Evaluation criteria Clustering: Silhouette width Given a set of labels for the cells (e. g. , biological condition, batch), silhouette widths provide a measure of goodness of the clustering of the cells with respect to these labels. Silhouette widths may be averaged within clusters or across all observations to assess clustering strength at the level of clusters or overall, respectively. The silhouette width s i of a sample i is defined as follows: s i = b i -a i max{a i , b i } , where a i = d(i, C cl(i) ), b i = min l =cl(i) d(i, C l ), C cl(i) is the cluster to which i belongs, and d(i, C l ) is the average distance between sample i and the samples in cluster C l. Average silhouette widths were used to compare ZINB-WaVE, PCA, and ZIFA on both real and simulated datasets. For simulated data, the cluster labels correspond to the true simulated W and, for each scenario, silhouette widths were computed and averaged over B datasets.
Clustering: Precision and recall Two different clusterings (i. e. , partitions) may be compared quantitatively using the precision and recall coefficients. These measures involve assessing whether pairs of cells cluster together in each of the two clusterings. The four different possibilities are enumerated in the following table. Then, using Clustering 1 as a reference, the precision coefficient is defined as the proportion of cells clustered together in Clustering 2 which are also (correctly) clustered together in Clustering 1 P recision = Y Y Y Y + N Y. Similarly, the recall coefficient is defined as the proportion of cells (correctly) clustered together in Clustering 1 which are also clustered together in Clustering 2 Recall = Y Y Y Y + Y N. Precision and recall were used to compare ZINB-WaVE, PCA, and ZIFA on simulated datasets, where the reference Clustering 1 corresponds to the true simulated labels and Clustering 2 was performed using the R function kmeans in the reduced-dimensional space: the first two columns of W for ZINB-WaVE, the first two principal components for PCA, and the first two latent variables for ZIFA.