14.1 Basic usage. The volcano plot that is being produced after this analysis is wierd and seems not to be correct. For each method, we compared the permutation P-values to the P-values directly computed by each method, which we define as the method P-values. Generally, the NPV values were more similar across methods. However, a better approach is to avoid using p-values as quantitative / rankable results in plots; they're not meant to be used in that way. . In that case, the number of modes in the expression distribution in the CF group (bimodal) and the non-CF group (unimodal) would be different, but the pseudobulk method may not detect a difference, because it is only able to detect differences in mean expression. Cons: Standard normalization, scaling, clustering and dimension reduction were performed using the R package Seurat version 3.1.1 (Butler et al., 2018; Satija et al., 2015; Stuart et al., 2019). We evaluated the performance of our tested approaches for human multi-subject DS analysis in health and disease. Another interactive feature provided by Seurat is being able to manually select cells for further investigation. Cons: Specifically, we considered a setting in which there were two groups of subjects to compare, containing four and three subjects, respectively with 21 731 genes. Volcano plots represent a useful way to visualise the results of differential expression analyses. Here, we introduce a mathematical framework for modeling different sources of biological variation introduced in scRNA-seq data, and we provide a mathematical justification for the use of pseudobulk methods for DS analysis. Supplementary Figure S12b shows the top 50 genes for each method, defined as the genes with the 50 smallest adjusted P-values. As in Section 3.5, in the bulk RNA-seq, genes with adjusted P-values less than 0.05 and at least a 2-fold difference in gene expression between healthy and IPF are considered true positives and all others are considered true negatives. Published by Oxford University Press. True positives were identified as those genes in the bulk RNA-seq analysis with FDR<0.05 and |log2(CD66+/CD66)|>1. The volcano plots for the three scRNA-seq methods have similar shapes, but the wilcox and mixed methods have inflated adjusted P-values relative to subject (Fig. (a) Volcano plots and (b) heatmaps of top 50 genes for 7 different DS analysis methods. These approaches will likely yield better type I and type II error rate control, but as we saw for the mixed method in our simulation, the computation times can be substantially longer and the computational burden of these methods scale with the number of cells, whereas the pseudobulk method scales with the number of subjects. ## [7] crosstalk_1.2.0 listenv_0.9.0 scattermore_0.8 Visualize single cell expression distributions in each cluster, # Violin plot - Visualize single cell expression distributions in each cluster, # Feature plot - visualize feature expression in low-dimensional space, # Dot plots - the size of the dot corresponds to the percentage of cells expressing the, # feature in each cluster. Multiple methods and bioinformatic tools exist for initial scRNA-seq data processing, including normalization, dimensionality reduction, visualization, cell type identification, lineage relationships and differential gene expression (DGE) analysis (Chen et al., 2019; Hwang et al., 2018; Luecken and Theis, 2019; Vieth et al., 2019; Zaragosi et al., 2020). Help! 6b). FindMarkers from Seurat returns p values as 0 for highly significant genes. Figure 4a shows volcano plots summarizing the DS results for the seven methods. All seven methods identify two distinct groups of genes: those with higher average expression in large airways and those with higher average expression in small airways. Performance measures for DS analysis of simulated data. We performed DS analysis using the same seven methods as Section 3.1. Consider a purified cell type (PCT) study design, in which many cells from a cell type of interest could be isolated and profiled using bulk RNA-seq. Four of the methods were applications of the FindMarkers function in the R package Seurat (Butler et al., 2018; Satija et al., 2015; Stuart et al., 2019) with different options for the type of test performed: for the method wilcox, cell counts were normalized, log-transformed and a Wilcoxon rank sum test was performed for each gene; for the method NB, cell counts were modeled using a negative binomial generalized linear model; for the method MAST, cell counts were modeled using a hurdle model based on the MAST software (Finak et al., 2015) and for the method DESeq2, cell counts were modeled using the DESeq2 software (Love et al., 2014). Supplementary Figure S14(cd) show that generally the shapes of the volcano plots are more similar between the subject and mixed methods than the wilcox method. If zjc1,zjc2,,zjcL are L cell-level covariates, then a log-linear regression model could take the form logijc=lzjclijl. Create volcano plot. The volcano plot for the subject method shows three genes with adjusted P-value <0.05 (log10(FDR) > 1.3), whereas the other six methods detected a much larger number of genes. Was this translation helpful? In this comparison, many genes were detected by all seven methods. In addition to returning a vector of cell names, CellSelector() can also take the selected cells and assign a new identity to them, returning a Seurat object with the identity classes already set. Although, in this work, we only consider the simple model presented above, the model could be extended to allow for systematic variation between cells by imposing a regression model in stage ii. For each subject, the number of cells and numbers of UMIs per cell were matched to the pig data. For each of these two cell types, the expression profiles are compared to all other cells as in traditional marker detection analysis. In addition to the inference reports and the associated Volcano plot views that allow users to visualize the distribution of fold change of all genes from say, one cluster to another, or one cluster to all cells, users can also visualize the normalized read . With this data you can now make a volcano plot. ADD REPLY link 18 months ago by Kevin Blighe 84k 0. I understand a little bit more now. S14e), we find that the subject and wilcox methods produce ranked gene lists with higher frequencies of marker genes than the mixed method, with subject having a slightly higher detection of known markers than wilcox. Our study highlights user-friendly approaches for analysis of scRNA-seq data from multiple biological replicates. Importantly, although these results specifically target differences in small airway secretory cells and are not directly comparable with other transcriptome studies, previous bulk RNA-seq (Bartlett et al., 2016) and microarray (Stoltz et al., 2010) studies have suggested few gene expression differences in airway epithelial tissues between CF and non-CF pigs; true differential gene expression between genotypes at birth is therefore likely to be small, as detected by the subject method. EnhancedVolcano (Blighe, Rana, and Lewis 2018) will attempt to fit as many labels in the plot window as possible, thus avoiding 'clogging' up the . In terms of identifying the true positives, wilcox and mixed had better performance (TPR = 0.62 and 0.56, respectively) than subject (TPR = 0.34). The volcano plots for subject and mixed show a stronger association between effect size (absolute log2-transformed fold change) and statistical significance (negative log10-transformed adjusted P-value). Seurat utilizes Rs plotly graphing library to create interactive plots. ## loaded via a namespace (and not attached): ## [1] systemfonts_1.0.4 plyr_1.8.8 igraph_1.4.1, ## [4] lazyeval_0.2.2 sp_1.6-0 splines_4.2.0, ## [7] crosstalk_1.2.0 listenv_0.9.0 scattermore_0.8, ## [10] digest_0.6.31 htmltools_0.5.5 fansi_1.0.4, ## [13] magrittr_2.0.3 memoise_2.0.1 tensor_1.5, ## [16] cluster_2.1.3 ROCR_1.0-11 limma_3.54.1, ## [19] globals_0.16.2 matrixStats_0.63.0 pkgdown_2.0.7, ## [22] spatstat.sparse_3.0-1 colorspace_2.1-0 rappdirs_0.3.3, ## [25] ggrepel_0.9.3 textshaping_0.3.6 xfun_0.38, ## [28] dplyr_1.1.1 crayon_1.5.2 jsonlite_1.8.4, ## [31] progressr_0.13.0 spatstat.data_3.0-1 survival_3.3-1, ## [34] zoo_1.8-11 glue_1.6.2 polyclip_1.10-4, ## [37] gtable_0.3.3 leiden_0.4.3 future.apply_1.10.0, ## [40] abind_1.4-5 scales_1.2.1 spatstat.random_3.1-4, ## [43] miniUI_0.1.1.1 Rcpp_1.0.10 viridisLite_0.4.1, ## [46] xtable_1.8-4 reticulate_1.28 ggmin_0.0.0.9000, ## [49] htmlwidgets_1.6.2 httr_1.4.5 RColorBrewer_1.1-3, ## [52] ellipsis_0.3.2 ica_1.0-3 farver_2.1.1, ## [55] pkgconfig_2.0.3 sass_0.4.5 uwot_0.1.14, ## [58] deldir_1.0-6 utf8_1.2.3 tidyselect_1.2.0, ## [61] labeling_0.4.2 rlang_1.1.0 reshape2_1.4.4, ## [64] later_1.3.0 munsell_0.5.0 tools_4.2.0, ## [67] cachem_1.0.7 cli_3.6.1 generics_0.1.3, ## [70] ggridges_0.5.4 evaluate_0.20 stringr_1.5.0, ## [73] fastmap_1.1.1 yaml_2.3.7 ragg_1.2.5, ## [76] goftest_1.2-3 knitr_1.42 fs_1.6.1, ## [79] fitdistrplus_1.1-8 purrr_1.0.1 RANN_2.6.1, ## [82] pbapply_1.7-0 future_1.32.0 nlme_3.1-157, ## [85] mime_0.12 formatR_1.14 compiler_4.2.0, ## [88] plotly_4.10.1 png_0.1-8 spatstat.utils_3.0-2, ## [91] tibble_3.2.1 bslib_0.4.2 stringi_1.7.12, ## [94] highr_0.10 desc_1.4.2 lattice_0.20-45, ## [97] Matrix_1.5-3 vctrs_0.6.1 pillar_1.9.0, ## [100] lifecycle_1.0.3 spatstat.geom_3.1-0 lmtest_0.9-40, ## [103] jquerylib_0.1.4 RcppAnnoy_0.0.20 data.table_1.14.8, ## [106] cowplot_1.1.1 irlba_2.3.5.1 httpuv_1.6.9, ## [109] R6_2.5.1 promises_1.2.0.1 KernSmooth_2.23-20, ## [112] gridExtra_2.3 parallelly_1.35.0 codetools_0.2-18, ## [115] MASS_7.3-56 rprojroot_2.0.3 withr_2.5.0, ## [118] sctransform_0.3.5 parallel_4.2.0 grid_4.2.0, ## [121] tidyr_1.3.0 rmarkdown_2.21 Rtsne_0.16, ## [124] spatstat.explore_3.1-0 shiny_1.7.4, Fast integration using reciprocal PCA (RPCA), Integrating scRNA-seq and scATAC-seq data, Demultiplexing with hashtag oligos (HTOs), Interoperability between single-cell object formats. For higher numbers of differentially expressed genes (pDE > 0.01), the subject method had lower NPV values when = 0.5 and similar or higher NPV values when > 0.5. The vertical axes give the performance measures, and the horizontal axes label each method. ## Running under: Ubuntu 20.04.5 LTS When only 1% of genes were differentially expressed (pDE = 0.01), all methods had NPV values near 1. Increasing sequencing depth can reduce technical variation and achieve more precise expression estimates, and collecting samples from more subjects can increase power to detect differentially expressed genes. As scRNA-seq costs have decreased, collecting data from more than one biological replicate has become more feasible, but careful modeling of different layers of biological variation remains challenging for many users. We propose an extension of the negative binomial model to scRNA-seq data by introducing an additional stage in the model hierarchy. The negative binomial distribution has a convenient interpretation as a hierarchical model, which is particularly useful for sequencing studies. In summary, here we (i) suggested a modeling framework for scRNA-seq data from multiple biological sources, (ii) showed how failing to account for biological variation could inflate the FDR of DS analysis and (iii) provided a formal justification for the validity of pseudobulking to allow DS analysis to be performed on scRNA-seq data using software designed for DS analysis of bulk RNA-seq data (Crowell et al., 2020; Lun et al., 2016; McCarthy et al., 2017). The main idea of the theorem is that if gene counts are summed across cells and the number of cells grows large for each subject, the influence of cell-level variation on the summed counts is negligible. However, a better approach is to avoid using p-values as quantitative / rankable results in plots; they're not meant to be used in that way. It is helpful to inspect the proposed model under a simplifying assumption. Red and blue dots represent genes with a log 2 FC (fold . Among the other five methods, when the number of differentially expressed genes was small (pDE = 0.01), the mixed method had the highest PPV values, whereas for higher numbers of differentially expressed genes (pDE > 0.01), the DESeq2 method had the highest PPV values. In stage iii, technical variation in counts is generated from a Poisson distribution. Four of the methods were applications of the FindMarkers function in the R package Seurat (Butler et al., 2018; . Step 3: Create a basic volcano plot. Define the aggregated countsKij=cKijc, and let sj=csjc. A software package, aggregateBioVar, is freely available on Bioconductor (https://www.bioconductor.org/packages/release/bioc/html/aggregateBioVar.html) to accommodate compatibility with upstream and downstream methods in scRNA-seq data analysis pipelines. True positives were identified as those genes in the bulk RNA-seq analysis with FDR<0.05 and |log2(IPF/healthy)|>1. ## [64] later_1.3.0 munsell_0.5.0 tools_4.2.0 Marker detection methods allow quantification of variation between cells and exploration of expression heterogeneity within tissues. . So, If I change the assay to "RNA", how we can trust that the DEGs are not due . In addition to simulated data, we analysed an animal model dataset containing large and small airway epithelia from CF and non-CF pigs (Rogers et al., 2008). For the AM cells (Fig. make sure label exists on your cells in the metadata corresponding to treatment (before- and after-), You will be returned a gene list of pvalues + logFc + other statistics. These analyses provide guidance on strengths and weaknesses of different methods in practice. ## [100] lifecycle_1.0.3 spatstat.geom_3.1-0 lmtest_0.9-40 If a gene was not differentially expressed, the value of i2 was set to 0. Results for alternative performance measures, including receiver operating characteristic (ROC) curves, TPRs and false positive rates (FPRs) can be found in Supplementary Figures S7 and S8. Third, the proposed model also ignores many aspects of the gene expression distribution in favor of simplicity. Alternatively, batch correction methods have been proposed to remove inter-individual differences prior to DS analysis, however, this increases type I error rates and disturbs the rank-order of results as explained in Zimmerman et al. d Volcano plots showing DE between T cells from random groups of unstimulated controls drawn . To generate such a plot, one can use SCpubr::do_VolcanoPlot (), which needs as input the Seurat object and the result of running Seurat::FindMarkers () choosing two groups. Future work with mixed models for scRNA-seq data should focus on maintaining scalable and computationally efficient implementation in software. 1 Answer. ## [55] pkgconfig_2.0.3 sass_0.4.5 uwot_0.1.14 ## [10] digest_0.6.31 htmltools_0.5.5 fansi_1.0.4 We have found this particularly useful for small clusters that do not always separate using unbiased clustering, but which look tantalizingly distinct. The null and alternative hypotheses for the i-th gene are H0i:i2=0 and H0i:i20, respectively. The difference between these formulas is in the mean calculation. ## [46] xtable_1.8-4 reticulate_1.28 ggmin_0.0.0.9000 However, the plot does not look well volcanic. (Lahnemann et al., 2020). The scRNA-seq data for the analysis of human lung tissue were obtained from GEO accession GSE122960, and the bulk RNA-seq of purified AT2 and AM fractions were shared by the authors immediately upon request. To use, simply make a ggplot2-based scatter plot (such as DimPlot() or FeaturePlot()) and pass the resulting plot to HoverLocator(). PR curves for DS analysis methods. ## [79] fitdistrplus_1.1-8 purrr_1.0.1 RANN_2.6.1 Carver College of Medicine, University of Iowa, Seq-Well: a sample-efficient, portable picowell platform for massively parallel single-cell RNA sequencing, Newborn cystic fibrosis pigs have a blunted early response to an inflammatory stimulus, Controlling the false discovery rate: a practical and powerful approach to multiple testing, The dynamics of gene expression in vertebrate embryogenesis at single-cell resolution, Integrating single-cell transcriptomic data across different conditions, technologies, and species, Comprehensive single-cell transcriptional profiling of a multicellular organism, Single-cell reconstruction of human basal cell diversity in normal and idiopathic pulmonary fibrosis lungs, Single-cell RNA-seq technologies and related computational data analysis, Muscat detects subpopulation-specific state transitions from multi-sample multi-condition single-cell transcriptomics data, Discrete distributional differential expression (D3E)a tool for gene expression analysis of single-cell RNA-seq data, MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data, PanglaoDB: a web server for exploration of mouse and human single-cell RNA sequencing data, Highly multiplexed single-cell RNA-seq by DNA oligonucleotide tagging of cellular proteins, Data Analysis Using Regression and Multilevel/Hierarchical Models, Seq-Well: portable, low-cost RNA sequencing of single cells at high throughput, SINCERA: a pipeline for single-cell RNA-seq profiling analysis, baySeq: empirical Bayesian methods for identifying differential expression in sequence count data, Single-cell RNA sequencing technologies and bioinformatics pipelines, Multiplexed droplet single-cell RNA-sequencing using natural genetic variation, Bayesian approach to single-cell differential expression analysis, Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells, A statistical approach for identifying differential distributions in single-cell RNA-seq experiments, Eleven grand challenges in single-cell data science, EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments, Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2, Current best practices in single-cell RNA-seq analysis: a tutorial, A step-by-step workflow for low-level analysis of single-cell RNA-seq data with bioconductor, Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets, Scater: pre-processing, quality control, normalization and visualization of single-cell RNA-seq data in R, DEsingle for detecting three types of differential expression in single-cell RNA-seq data, Comparative analysis of sequencing technologies for single-cell transcriptomics, Single-cell mRNA quantification and differential analysis with Census, Reversed graph embedding resolves complex single-cell trajectories, Single-cell transcriptomic analysis of human lung provides insights into the pathobiology of pulmonary fibrosis, edgeR: a Bioconductor package for differential expression analysis of digital gene expression data, Disruption of the CFTR gene produces a model of cystic fibrosis in newborn pigs, Single-cell profiling of the developing mouse brain and spinal cord with split-pool barcoding, Spatial reconstruction of single-cell gene expression data, Single-cell transcriptomes of the human skin reveal age-related loss of fibroblast priming, Cystic fibrosis pigs develop lung disease and exhibit defective bacterial eradication at birth, Comprehensive integration of single-cell data, The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells, RNA sequencing data: Hitchhikers guide to expression analysis, A systematic evaluation of single cell RNA-seq analysis pipelines, Sequencing thousands of single-cell genomes with combinatorial indexing, Comparative analysis of differential gene expression analysis tools for single-cell RNA sequencing data, SigEMD: A powerful method for differential gene expression analysis in single-cell RNA sequencing data, Using single-cell RNA sequencing to unravel cell lineage relationships in the respiratory tract, Comparative analysis of droplet-based ultra-high-throughput single-cell RNA-seq systems, Comparative analysis of single-cell RNA sequencing methods, A practical solution to pseudoreplication bias in single-cell studies. (d) ROC and PR curves for subject, wilcox and mixed methods using bulk RNA-seq as a gold standard. The authors thank Michael J. Welsh, Joseph Zabner, Kai Wang and Keyan Zarei for careful reading of the manuscript and helpful feedback that improved the clarity and content in the final draft. We will also label the top 10 most significant genes with their . The analyses presented here have illustrated how different results could be obtained when data were analysed using different units of analysis. Yes, you can use the second one for volcano plots, but it might help to understand what it's implying. As an example, were going to select the same set of cells as before, and set their identity class to selected. In your last function call, you are trying to group based on a continuous variable pct.1 whereas group_by expects a categorical variable. ## [94] highr_0.10 desc_1.4.2 lattice_0.20-45 Furthermore, guidelines for library complexity in bulk RNA-seq studies apply to data with heterogeneity between cell types, so these recommendations should be sufficient for both PCT and scRNA-seq studies, in which data have been stratified by cell type. The subject and mixed methods are composed of genes that have high inter-group (CF versus non-CF) and low intra-group (between subject) variability, whereas the wilcox, NB, MAST, DESeq2 and Monocle methods tend to be sensitive to a highly variable gene expression pattern from the third CF pig. Next, we applied our approach for marker detection and DS analysis to published human datasets. ## locale: Next, we used subject, wilcox and mixed to test for differences in expression between healthy and IPF subjects within the AT2 and AM cell populations. Generally, tests for marker detection, such as the wilcox method, are sufficient if type I error rate control is less of a concern than type II error rate and in circumstances where type I error rate is most important, methods like subject and mixed can be used. Further, they used flow cytometry to isolate alveolar type II (AT2) cell and alveolar macrophage (AM) fractions from the lung samples and profiled these PCTs using bulk RNA-seq. ## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3 data("pbmc_small") # Find markers for cluster 2 markers <- FindMarkers(object = pbmc_small, ident.1 = 2) head(x = markers) # Take all cells in cluster 2, and find markers that separate cells in the 'g1' group (metadata # variable 'group') markers <- FindMarkers(pbmc_small, ident.1 = "g1", group.by = 'groups', subset.ident = "2") head(x = markers) # Pass 'clustertree' or an object of class . (b) AT2 cells and AM express SFTPC and MARCO, respectively. Returns a volcano plot from the output of the FindMarkers function from the Seurat package, which is a ggplot object that can be modified or plotted. Supplementary Table S2 contains performance measures derived from the ROC and PR curves. I have successfully installed ggplot, normalized my datasets, merged the datasets, etc., but what I do not understand is how to transfer the sequencing data to the ggplot function. ## 13714 features across 2638 samples within 1 assay, ## Active assay: RNA (13714 features, 2000 variable features), ## 2 dimensional reductions calculated: pca, umap, # Ridge plots - from ggridges. (2019) used scRNA-seq to profile cells from the lungs of healthy subjects and those with pulmonary fibrosis disease subtypes, including hypersensitivity pneumonitis, systemic sclerosis-associated and myositis-associated interstitial lung diseases and IPF (Reyfman et al., 2019). This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, https://doi.org/10.1093/bioinformatics/btab337, https://www.bioconductor.org/packages/release/bioc/html/aggregateBioVar.html, https://creativecommons.org/licenses/by/4.0/, Receive exclusive offers and updates from Oxford Academic, Academic Pulmonary Sleep Medicine Physician Opportunity in Scenic Central Pennsylvania, MEDICAL MICROBIOLOGY AND CLINICAL LABORATORY MEDICINE PHYSICIAN, CLINICAL CHEMISTRY LABORATORY MEDICINE PHYSICIAN. In Supplementary Figure S14(ef), we quantify the ability of each method to correctly identify markers of T cells and macrophages from a database of known cell type markers (Franzen et al., 2019). Theorem 1 implies that when the number of cells per subject is large, the aggregated counts follow a distribution with the same mean and variance structure as the negative binomial model used in many software packages for DS analysis of bulk RNA-seq data. # Particularly useful when plotting multiple markers, # Visualize co-expression of two features simultaneously, # Split visualization to view expression by groups (replaces FeatureHeatmap), # Violin plots can also be split on some variable. (Crowell et al., 2020) provides a thorough comparison of a variety of DGE methods for scRNA-seq with biological replicates including: (i) marker detection methods, (ii) pseudobulk methods, where gene counts are aggregated between cells from different biological samples and (iii) mixed models, where models for gene expression are adjusted for sample-specific or batch effects. As we observed in Figure 2, the subject method had a larger area under the curve than the other six methods in all simulation settings, with larger differences for higher signal-to-noise ratios. The FindAllMarkers () function has three important arguments which provide thresholds for determining whether a gene is a marker: logfc.threshold: minimum log2 fold change for average expression of gene in cluster relative to the average expression in all other clusters combined. In scRNA-seq studies, where cells are collected from multiple subjects (e.g. Here, we propose a statistical model for scRNA-seq gene counts, describe a simple method for estimating model parameters and show that failing to account for additional biological variation in scRNA-seq studies can inflate false discovery rates (FDRs) of statistical tests. ## [124] spatstat.explore_3.1-0 shiny_1.7.4. Here, we present the DS results comparing CF and non-CF pigs only in secretory cells from the small airways. The wilcox, MAST and Monocle methods had intermediate performance in these nine settings. When samples correspond to different experimental subjects, the first stage characterizes biological variation in gene expression between subjects.
General Scott Miller Wife,
Coyote Lake Ending Explained,
Cuando Un Hombre Te Pregunta Si Llegaste Bien,
Articles F