dataesr/scientific-paragraphs-categorization huggingface dataset:
Using our posterior probabilities as observation-level weights in EDGER recovers the original BCV plot and mean-variance trend (Fig. 2c ), illustrating the ability of our method to account for zero inflation. Hence, observation weights provide the key to unlocking standard bulk RNA-seq tools for zero-inflated data. The BCV plot for the Islam et al. [16]  scRNA-seq dataset (Fig. 2d )  [35]  shows similar striped patterns as for zeroinflated bulk RNA-seq data. Such patterns are observed in many single-cell datasets (Additional file 1: Figure  S2 ). ZINB-WaVE identifies many zeros to be excess for the Islam dataset. It also provides good classification power for excess zeros for data simulated from the Islam dataset (Fig. 2e ). Incorporating the ZINB-WaVE weights in an EDGER analysis removes the striped patterns and yields a BCV plot that is similar to that for bulk RNA-seq data (Fig.
High power and false positive control on simulated (sc)RNA-seq data We provide a scRNA-seq data simulation paradigm that retains gene-specific characteristics as well as global associations across all genes (see "Methods" for details). More specifically, we first estimate dataset-specific associations between zero abundance, sequencing depth, and average log counts per million (CPM), and then explicitly account for these associations in our simulation model (Additional file 1: Figures  S4, S5 ). The scRNA-seq simulation study is based on three datasets: the Islam et al. [16]  dataset  [35] , comparing 48 embryonic stem cells to 44 embryonic fibroblasts in mouse; a subset of the Trapnell et al. [36]  dataset, comparing differentiating human myoblasts at the 48 h (85 cells) and 72 h (64 cells) timepoints; and a 10x Genomics peripheral blood mononuclear cell (PBMC) dataset (see "Real datasets" in "Methods" for details).
The datasets differ in throughput, sequencing depth, and extent of zero inflation. For example, Additional file 1: Figure  S6  shows a higher proportion of excess zeros in the Islam dataset compared to the Trapnell dataset, an observation further supported because the Islam and Trapnell datasets contain ∼65% and ∼48% zeros, respectively. 10x Genomics datasets are known to contain even more zeros. The evaluated subset of the PBMC dataset contains ∼87% zeros. The simulated datasets successfully mimic the characteristics of the original datasets, as evaluated with the R package COUNTSIMQC  [37]  (Additional files 2, 3, and 4). This diverse range of datasets is, therefore, representative of scRNA-seq datasets that occur in practice and it is a suitable basis for method evaluation and comparison.
We evaluate the performance of the method in terms of sensitivity and false positive control using false discovery proportion vs. true positive rate (FDP-TPR) curves. Figure  3  (Additional file 1: Figure  S7 ) illustrates that many methods break down on the simulated Islam dataset due to a high degree of zero inflation. Surprisingly, even methods specifically developed to deal with excess zeros, like SCDE and METAGENOMESEQ, suffer from poor performance, with MAST being a notable exception. The DESEQ2 methods, however, are able to cope with the high degree of zero inflation. Note that, in general, it is a good strategy to disable the outlier imputation step in DESEQ2, since it deteriorates performance on scRNA-seq data (Additional file 1: Figure  S8 ). SEURAT, limma-voom, and SCDE have very low sensitivity.
The methods based on ZINB-WaVE weights dominate all competitors in terms of sensitivity and specificity, providing high power, good false discovery rate (FDR) control, and sensible p-value distributions (Additional file 1: Figure  S9 ). Note that the remaining methods also suffer from poor FDR control. Since zero inflation is fairly modest for the Trapnell dataset, most methods perform better than for the Islam simulation (Fig. 3 ). The ZINB-WaVE-based methods and DESEQ2 outperform the remaining methods in terms of sensitivity and provide good FDR control. EDGER is their closest competitor. The remaining methods provide much lower sensitivity and/or very liberal FDR control. Note how bespoke scRNA-seq methods seem to break down on datasets with a lower degree of zero inflation, often providing too liberal or too conservative p-value distributions, while ZINB-WaVE-based methods, in general, show a reasonable p-value distribution, with an enrichment of low p-values and approximately uniformly distributed larger p-values (Additional file 1: Figure  S10 ).
Typical 10x Genomics datasets contain a high number of cells with shallow sequencing depth, due to the extreme multiplexing of libraries. As a result, counts and hence, estimated NB means are lower, making zeros more plausible according to the NB distribution and excess zeros, thus, harder to identify. This is picked up by the simulation framework, where only ∼8% of the genes were simulated to have at least one excess zero in n = 1200 samples. Bulk RNA-seq methods can, hence, be expected to be among the top performers. Figure  4  shows FDP-TPR curves for the 10x Genomics simulation study, demonstrating the good performance of bulk RNA-seq methods EDGER and DESEQ2. ZINB-WaVE EDGER and ZINB-WaVE DESEQ2 are among the top performers, having comparable or slightly lower performance compared to their unweighted counterparts.
False positive rate control We compared our ZINB-WaVE-weight-based method to commonly used DE methods for mock comparisons based on two publicly available real scRNA-seq datasets. We assessed performance based on the per-comparison error rate (PCER), defined as the proportion of false positives (i. e. , type I errors) among all genes being considered for DE, where a gene is declared DE if its nominal unadjusted p-value is less than or equal to 0. 05. The first dataset, referred to as the Usoskin  [11]  dataset, is for 622 mouse neuronal cells from the dorsal root ganglion, classified into 11 categories. The authors acknowledge the existence of a batch effect related to the picking session for the cells. We find that the batch effect is not only associated with expression measures, but also influences the relationship between sequencing depth and zero abundance (Fig.
For the Usoskin dataset, we assessed false positive control by comparing the actual vs. the nominal PCER for mock null datasets where none of the genes are expected to be DE. Specifically, we generated 30 mock datasets where, for each dataset, two groups of 45 cells each were created by sampling 15 cells at random without replacement from each of the three picking sessions. Sampling cells within batch allows us to control for potential confounding by the batch variable. For each of the 30 mock datasets, we considered seven methods to identify genes that are DE between the two groups and declared a gene DE if its nominal unadjusted p-value was less than or equal to 0. 05. For these mock datasets, any gene declared DE between the two groups is a false positive. Thus, for each method, the nominal PCER of 0.
5d ). We also replicated the original analysis of Usoskin et al. [11] , by performing one-against-others tests of DE for each cell type (Additional file 1: Figure  S13 ). limma-voom found a high number of DE genes, confirming our results from the mock evaluations where it was too liberal. The ZINB-WaVE methods tended to find a high number of DE genes, which is promising combined with the good PCER control seen in the mock comparisons. While introducing ZINB-WaVE weights in DESEQ2 leads to a higher number of significant genes on average, the effect is less clear with EDGER and seems to depend on the contrast. Similar results were observed for a 10x Genomics PBMC dataset comprising 2700 single cells sequenced on an Illumina NextSeq 500 (Additional file 1: Figure  S14 ), with the distinction that we found a conservative p-value distribution for ZINB-WaVE-weighted DESEQ2.
Since no information was provided about potential batch effects, we did not consider batch covariates for this dataset. Additionally, we examined the PCER and p-value distributions on mock comparisons while varying the regularization parameter ( ) for the ZINB-WaVE estimation procedure. Not surprisingly, we observed that the PCER decreases with increasing , i. e. , as the parameters of the ZINB-WaVE model are subjected to more "shrinking" (Additional file 1: Figures  S15  and S16  for the Usoskin and 10x Genomics PBMC datasets, respectively).
Biologically meaningful clustering and DE results To analyze the 2700 cells from the 10x Genomics PBMC dataset (see "Methods"), we followed the tutorial available at http://satijalab. org/seurat/pbmc3k_tutorial. html and used the R package SEURAT  [39]. The major steps of the pipeline were quality control, data filtering, identification of high-variance genes, dimensionality reduction using the first ten components from principal component analysis (PCA), and graph-based clustering. The final step of the pipeline was to identify genes that are DE between clusters, to derive cell-type signatures. Two different parameterizations were used for the SEURAT clustering. With one parameterization, a single cluster was identified for CD4+ T cells, while with another, two CD4+ T-cell subclusters were identified, corresponding to CD4+ naive T cells and CD4+ memory T cells (gold and red To compute the weights for differential expression analysis, ZINB-WaVE was fitted with intercept and cell-type covariate in X, V = 1 J , K = 0 for W, common dispersion, and = 10 12.
6c  and Additional file 1: Figure  S17 ). As we do not have prior knowledge about the cells in the different clusters, we are unable to say whether the cluster is more representative of the CD4+ T-effector memory cell type or if our method missed the enrichment in the CD4+ T-central memory cell type. However, it is interesting that using ZINB-WaVE to account for zero inflation in the clustering allowed EDGER to find results that seem more biologically meaningful than without accounting for zero inflation. Finally, using a Benjamini and Hochberg  [42]  adjusted p-value cut-off of 0. 05, limma-voom declared 433 and 194 DE genes and weighted EDGER 371 and 151, for clustering based on, respectively, the first ten PCs and W from ZINB-WaVE. We additionally showed on mock comparisons for the same 10x Genomics PBMC dataset that limma-voom had a greater actual PCER than weighted EDGER (Additional file 1: Figure  S14 ), suggesting that some of the DE genes found by limma-voom are likely to be false positives.
However, for the Usoskin dataset, ZINGER seems to find a higher number of DE genes than ZINB-WaVE and its bulk RNA-seq counterparts (Additional file 1: Figure  S22 ), while also controlling the PCER in mock evaluations (Additional file 1: Figure  S23 ). Due to the computational burden of the ZINGER method, we were unable to apply it to large-scale datasets, such as those from the 10x Genomics platform, thus limiting our comparison.
Computational time The better performance of our ZINB-WaVE-weighted DE method comes at a computational cost, since we first fit ZINB-WaVE to the entire cells-by-genes matrix of read counts to compute the weights and then use a weighted version of DESEQ2 or EDGER for inferring DE. To give the reader an idea of how different methods scale in terms of computation time, we benchmarked three different datasets: the Islam dataset (92 cells), one of the mock null Usoskin datasets used in Fig. 5  (90 cells), and the CD4+ T-cell cluster of the 10x Genomics PBMC dataset (1151 cells). For each dataset, 10,000 genes were sampled at random and the two cell types were used as covariates. For the Usoskin dataset, batch was added as a covariate for all methods. For all datasets, the fastest method was limma-voom followed by EDGER (Additional file 1: Figure  S24 ).
Our results complement the findings of Jaakkola et al. [29]  and Soneson and Robinson  [30]  that bespoke scRNAseq tools do not systematically improve upon bulk RNAseq tools. Although MAST, METAGENOMESEQ, and SCDE were explicitly developed to handle excess zeros, they suffer from poor performance in a high zero-inflation setting, as demonstrated in the simulation study. The value of our method was demonstrated for scRNAseq protocols relying on both standard (Islam, Usoskin, and Trapnell datasets) and unique molecular identifier (UMI) (10x Genomics PBMC dataset) read counting. UMIs were recently proposed to reduce measurement variability across samples  [15]. In UMI-based protocols, transcripts are labeled with a small random UMI barcode prior to amplification. After amplification and sequencing, one enumerates the unique UMIs found for every transcript, which correspond to individual sequenced UMI-labeled transcripts.
There is some evidence in the literature that zero inflation is less of a problem for UMIbased than for full-length protocols and that UMI read counts could follow a NB distribution  [44, 45]. Hence, our method also provides good results for UMI-based data with limited zero inflation, demonstrating its broad applicability. In the simulation study, power to detect DE was generally lower for 10x Genomics UMI datasets (Fig. 4 ) than for full-length protocol datasets (Fig. 3 ). While the 10x Genomics platform has the advantage of an extremely high throughput, allowing many cells to be characterized, the resulting datasets often have the disadvantage of low library sizes, a logical consequence of UMI counting and of the trade-off between sequencing depth and number of cells to be sequenced in one sequencing run.
As a result, the sequencing depth of these datasets is much lower than that of bulk RNA-seq datasets, making it harder to identify excess zeros and assess DE, even in large sample size settings. Although the 10x Genomics platform may be well suited for hypothesis generation, e. g. , through cell-type discovery or lineage trajectory studies, full-length protocols may be more appropriate for discovering marker genes between inferred cell types or trajectories, an approach that has also been adopted in previous studies  [46]. We have used ZINB-WaVE in conjunction with either EDGER or DESEQ2. However, the ZINB-WaVE posterior probabilities could be used as weights to unlock other standard RNA-seq workflows in zero-inflation situations. Additional file 1: Figure  S7  shows that ZINB-WaVE weights combined with heteroscedastic weights in limma-voom also increase power in a scRNA-seq context, although this may be at the expense of type I error control.
The ZINB-WaVE method penalizes the L2 norm of the parameter estimates for regularization purposes. It requires a penalty parameter that is rescaled differently for gene-specific parameters, cell-specific parameters, and dispersion parameters  [23]. All analyses in this manuscript were performed with = 10 12 , to provide consistently comparable results. However, the optimal value of is dataset-specific and further research is needed to provide a data-driven approach for selecting an optimal. Indeed, based on our simulations, the value of the penalty parameter can have a profound influence on the results (Additional file 1: Figure  S25 ), but we found = 10 12 to have generally good performance. ZINB-WaVE has an option to infer latent variables W, which may correspond to either unmeasured confounding covariates or unmeasured covariates of interest.
The observational weights were computed with the number of unknown covariates K = 0, i. e. , no latent variables were inferred. To cluster the real datasets, we inferred an optimal choice of K using the AIC (Additional file 1: Figure  S18 ). However, further investigation is needed to confirm that the AIC is appropriate for selecting K. In principle, our proposed ZINB-WaVE model could also be used to identify DE genes both in terms of the NB mean and the zero-inflation probability, reflecting, respectively, a continuum in DE and a more binary (i. e. , presence or absence) DE pattern. In this context, the parameters of interest are regression coefficients β corresponding to known sample-level covariates in the matrix X used in either μ or π (Eq. 4). DE genes may be identified via likelihood ratio tests or Wald tests, with the standard errors of estimators of β obtained from the inverse of the Hessian matrix of the likelihood function.
DESEQ2 We extended the DESEQ2 package  [1]  to accommodate zero inflation by providing the option to use observationlevel weights in the parameter estimation steps. In this case, the ZINB-WaVE weights are supplied in the weight slot of an object of class DESeqDataSet. DESEQ2's default normalization procedure is based on geometric means of counts, which are zero for genes with at least one zero count. This greatly limits the number of genes that can be used for normalization in scRNAseq applications  [51]. We, therefore, use the normalization method suggested in the PHYLOSEQ package  [52] , which calculates the geometric mean for a gene using only its positive counts, so that genes with zero counts could still be used for normalization. The PHYLOSEQ normalization procedure can now be applied by setting the argument type equal to poscounts in the DESEQ2 function estimateSizeFactors.
limma-voom For the limma-voom approach  [3] , implemented in the voom function from the LIMMA package, heteroscedastic weights are estimated based on the mean-variance relationship of the log-transformed counts. We adapt limma-voom to zero-inflated situations by multiplying the heteroscedastic weights by the ZINB-WaVE weights and using the resulting weights in weighted linear regression. To account for the downweighting of zeros, the residual degrees of freedom of the linear model are adjusted, such as with EDGER, before the empirical Bayes variance shrinkage and are, therefore, also propagated to the moderated statistical tests. Both the standard and ZINB-WaVE-weighted versions of limma-voom were considered in the simulation study. The latter was not considered for the real datasets due to its poor performance in the simulation study.
DE method comparison We compared our weighted DE approach to state-of-theart bulk RNA-seq methods implemented in the packages EDGER (v3. 20. 1) [2, 50] , DESEQ2 (v1. 19. 8)  [1] , and LIMMA (v3. 34. 0)  [3]. We also considered dedicated scRNAseq tools from the packages SCDE (v2. 6. 0) [28] , MAST (v1. 4. 0)  [19] , and NODES (v0. 0. 0. 9010)  [55] , as well as METAGENOMESEQ (v1. 18. 0)  [48] , which was developed to account for zero inflation in metagenomics applications. A ZINB model is also implemented in SHRINKBAYES  [56] , but the method does not scale to the typical sample sizes encountered in scRNA-seq and has many tuning parameters, so we did not include it in our comparison. In DESEQ2, we disable the outlier imputation step and allow for shrinkage of fold-changes by default. In addition, for large 3 -end sequencing datasets like the Usoskin and 10x Genomics PBMC datasets, we set the minimum expected count estimated by DESEQ2 to 10 -6 , allowing the method to cope with large sample sizes and low counts.
scRNA-seq data simulation We extended the framework of Zhou et al. [57]  for scRNA-seq applications. In the GitHub repository linked to this manuscript (https://github. com/statOmics/ zinbwaveZinger), we provide user-friendly R code to simulate scRNA-seq read counts. The user can input a real scRNA-seq dataset to infer gene-level parameters for read count distributions. Library sizes for the simulated samples are by default resampled from the real dataset, but can also be user-specified. The simulation paradigm randomly resamples parameters estimated from the original dataset, where all parameters of a given gene are resampled jointly to retain gene-specific characteristics present in the original dataset. In scRNA-seq, dropouts and bursting lead to bias in parameter estimation if not properly accounted for.
We then compare, for every gene, the estimated probability of zero counts based on the model in Eq. 6 to the corresponding NB-based probability f NB 0; μij , θj with μij = λj N i , and randomly add excess zeros whenever the former probability is higher than the latter. The model in Eq. 6 is motivated by datasetspecific associations observed in real scRNA-seq datasets (Additional file 1: Figures  S4, S5 ). This framework acknowledges both gene-specific characteristics as well as broad dataset-specific associations across all genes and provides realistic scRNA-seq data for evaluating methods. We assessed the performance of various DE methods using data simulated based on the Islam et al. [16]  dataset, a subset of the Trapnell et al. [36]  dataset, and a 10x Genomics PBMC dataset. See "Real datasets" for information on these datasets.
Gene set enrichment analysis To identify cell types corresponding to the two CD4+ Tcell subclusters of the 10x Genomics PBMC dataset, we used GSEA with the function fgsea from the Bioconductor R package FGSEA (v1. 4. 0)  [40]  and gene sets for 64 immune and stroma cell types from the R package XCELL (v1. 1. 0) [41]. For each DE method, the input to fgsea is a list of genes ranked by a test statistic comparing expression in the two CD4+ T-cell subclusters. To facilitate comparison between DE methods, the test statistic used here is a transformation of the unadjusted p-values (p) with the sign of the log-fold-change (lfc): -1 (1p/2) sign(lfc), where (•) denotes the standard Gaussian cumulative distribution function. As suggested by FGSEA, all genes were used for the analysis. To assess the enrichment/depletion of one cluster compared to the other cluster, we used the normalized enrichment score.
Real datasets Usoskin dataset This dataset is for mouse neuronal cells from the dorsal root ganglion, sequenced on either an Illumina Genome Analyzer IIx or HiSeq 2000  [11]. The cells were robotically picked in three separate sessions and the 5 end of the transcripts sequenced. The expression measures were downloaded from supplementary data accompanying the original manuscript (http://linnarssonlab. org/drg/). After quality control and sample filtering (removal of non-single cells and non-neuronal cells), the authors considered 622 cells, which were classified into 11 neuronal cell-type categories. Only genes with more than 20 non-zero counts were retained, for a total of 12,132 genes. There is a batch effect related to the picking session for the cells. For the DE analysis, the picking session was, therefore, included as a batch covariate in all models.
To mimic a null dataset with no DE, we created two groups of 45 cells each, where, for each group, 15 cells were sampled at random, without replacement (over all cell types) from each picking session. For each of 30 such mock null datasets, we considered seven methods to identify genes that are DE between the two groups and declared a gene DE if its nominal unadjusted p-value was less than or equal to 0. 05. For these mock datasets, any gene declared DE between the two groups is a false positive. The nominal PCER of 0. 05 for each method is compared to its actual PCER, which is simply the proportion of genes declared DE.
10x Genomics PBMC dataset We analyzed a dataset of PBMCs that is freely available from 10x Genomics (https://support. 10xgenomics. com/single-cell-gene-expression/datasets/1. 1. 0/pbmc3k) [62]. We downloaded the data from https://s3-uswest-2. amazonaws. com/10x. files/samples/cell/pbmc3k/ pbmc3k_filtered_gene_bc_matrices. tar. gz, which correspond to 2700 single cells sequenced on an Illumina NextSeq 500 using UMIs. We clustered cells following the tutorial available at http://satijalab. org/seurat/pbmc3k_ tutorial. html and using the R package SEURAT (v2. 1. 0) [39]. The major steps of the pipeline are quality control, data filtering, identification of high-variance genes, dimensionality reduction using the first ten components from PCA, and graph-based clustering. To identify cluster markers, we used our ZINB-WaVE-weighted DE method instead of the method implemented in SEURAT.