Carver College of Medicine, University of Iowa. As you can see, there are four major groups of genes: - Genes that surpass our p-value and logFC cutoffs (blue). ## [103] jquerylib_0.1.4 RcppAnnoy_0.0.20 data.table_1.14.8 To avoid confounding the results by disease, this analysis is confined to data from six healthy subjects in the dataset. ## [106] cowplot_1.1.1 irlba_2.3.5.1 httpuv_1.6.9 The method subject treated subjects as the units of analysis, and statistical tests were performed according to the procedure outlined in Sections 2.2 and 2.3. It enables quick visual identification of genes with large fold changes that are also statistically significant. 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. Session Info 1. Yes, you can use the second one for volcano plots, but it might help to understand what it's implying. PR curves for DS analysis methods. ## [13] SeuratData_0.2.2 SeuratObject_4.1.3 If mi is the sample mean of {Eij} over j, vi is the sample variance of {Eij} over j, mij is the sample mean of {Eijc} over c, and vij is the sample variance of {Eijc} over c, we fixed the subject-level and cell-level variance parameters to be i=vi/mi2 and ij2=vij/mij2, respectively. First, the adjusted P-values for each method are sorted from smallest to largest. Our study highlights user-friendly approaches for analysis of scRNA-seq data from multiple biological replicates. 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. #' @param plot.adj.pvalue logical specifying whether adjusted p-value should by plotted on the y-axis. ## [7] crosstalk_1.2.0 listenv_0.9.0 scattermore_0.8 Four of the methods were applications of the FindMarkers function in the R package Seurat (Butler et al., 2018; . We will call genes significant here if they have FDR < 0.01 and a log2 fold change of 0.58 (equivalent to a fold-change of 1.5). 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. The implemented methods are subject (red), wilcox (blue), NB (green), MAST (purple), DESeq2 (orange), monocle (gold) and mixed (brown). In practice, often only one cutoff value for the adjusted P-value will be chosen to detect genes. In a scRNA-seq experiment with multiple subjects, we assume that the observed data consist of gene counts for G genes drawn from multiple cells among n subjects. FindMarkers: Finds markers (differentially expressed genes) for identified clusters. ## [121] tidyr_1.3.0 rmarkdown_2.21 Rtsne_0.16 The subject method had the shortest average computation times, typically <1 min. The observed counts for the PCT study are analogous to the aggregated counts for one cell type in a scRNA-seq study. Because these assumptions are difficult to validate in practice, we suggest following the guidelines for library complexity in bulk RNA-seq studies. ## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C Supplementary Figure S12b shows the top 50 genes for each method, defined as the genes with the 50 smallest adjusted P-values. 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. This creates a data.frame with gene names as rows, and includes avg_log2FC, and adjusted p-values. . In recent years, the reagent and effort costs of scRNA-seq have decreased dramatically as novel techniques have been developed (Aicher et al., 2019; Briggs et al., 2018; Cao et al., 2017; Chen et al., 2019; Gehring et al., 2020; Gierahn et al., 2017; Klein et al., 2015; Macosko et al., 2015; Natarajan et al., 2019; Rosenberg et al., 2018; Vitak et al., 2017; Zhang et al., 2019; Ziegenhain et al., 2017), so that biological replication, meaning data collected from multiple independent biological units such as different research animals or human subjects, is becoming more feasible; biological replication allows generalization of results to the population from which the sample was drawn. A volcano plot is a type of scatterplot that shows statistical significance (P value) versus magnitude of change (fold change). Nine simulation settings were considered. 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 each panel, PR curves are plotted for each of seven DS analysis methods: subject (red), wilcox (blue), NB (green), MAST (purple), DESeq2 (orange), Monocle (gold) and mixed (brown). 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. See ?FindMarkers in the Seurat package for all options. ## [55] pkgconfig_2.0.3 sass_0.4.5 uwot_0.1.14 Improvements in type I and type II error rate control of the DS test could be considered by modeling cell-level gene expression adjusted for potential differences in gene expression between subjects, similar to the mixed method in Section 3. Plots a volcano plot from the output of the FindMarkers function from the Seurat package or the GEX_cluster_genes function alternatively. If the ident.2 parameter is omitted or set to NULL, FindMarkers () will test for differentially expressed features between the group specified by ident.1 and all other cells. Infinite p-values are set defined value of the highest -log(p) + 100. The recall, also known as the true positive rate (TPR), is the fraction of differentially expressed genes that are detected. If a gene was not differentially expressed, the value of i2 was set to 0. #' @return 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. This is the model used in DESeq2 (Love et al., 2014). With Seurat, all plotting functions return ggplot2-based plots by default, allowing one to easily capture and manipulate plots just like any other ggplot2-based plot. In order to contrast DS analysis with cells as units of analysis versus subjects as units of analysis, we analysed both simulated and experimental data. With this data you can now make a volcano plot; Repeat for all cell clusters/types of interest, depending on your research questions. 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. In extreme cases, where only a few cells have been collected for some subjects, interpretation of gene expression differences should be handled with caution. The expression level of gene i for group 1, i1, was matched to the pig data by setting ei1=jcKijc/i'jcKi'jc. We have developed the software package aggregateBioVar (available on Bioconductor) to facilitate broad adoption of pseudobulk-based DE testing; aggregateBioVar includes a detailed vignette, has low code complexity and minimal dependencies and is highly interoperable with existing RNA-seq analysis software using Bioconductor core data structures (Fig. These were the values used in the original paper for this dataset. 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. A more powerful statistical test that yields well-controlled FDR could be constructed by considering techniques that estimate all parameters of the hierarchical model. ## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C Supplementary data are available at Bioinformatics online. ## [9] panc8.SeuratData_3.0.2 ifnb.SeuratData_3.1.0 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. (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. For the T cells, (Supplementary Fig. 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. Give feedback. ## [115] MASS_7.3-56 rprojroot_2.0.3 withr_2.5.0 More conventional statistical techniques for hierarchical models, such as maximum likelihood or Bayesian maximum a posteriori estimation, could produce less noisy parameter estimates and hence, lead to a more powerful DS test (Gelman and Hill, 2007). In contrast, single-cell experiments contain an additional source of biological variation between cells. It furthers the University's objective of excellence in research, scholarship, and education by publishing worldwide, This PDF is available to Subscribers Only. ## [46] xtable_1.8-4 reticulate_1.28 ggmin_0.0.0.9000 For clarity of exposition, we adopt and extend notations similar to (Love et al., 2014). ## [52] ellipsis_0.3.2 ica_1.0-3 farver_2.1.1 The subject method had the highest PPV, and the NB method had the lowest PPV in all nine simulation settings. #' @param output_dir The relative directory that will be used to save results. ## [11] hcabm40k.SeuratData_3.0.0 bmcite.SeuratData_0.3.0 . The other two methods were Monocle, which utilized a negative binomial generalized additive model to test for differences in gene expression using the R package Monocle (Qiu et al., 2017a, b; Trapnell et al., 2014) and mixed, which modeled counts using a negative binomial generalized linear mixed model with a random effect to account for differences in gene expression between subjects and DS testing was performed using a Wald test. If subjects are composed of different proportions of types A and B, DS results could be due to different cell compositions rather than different mean expression levels. I change the test.use but did not work. Each panel shows results for 100 simulated datasets in one simulation setting. ## [94] highr_0.10 desc_1.4.2 lattice_0.20-45 (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). sessionInfo()## R version 4.2.0 (2022-04-22) In general, the method subject had lower area under the ROC curve and lower TPR but with lower FPR. As scRNA-seq studies grow in scope, due to technological advances making these studies both less labor-intensive and less expensive, biological replication will become the norm. ## [64] later_1.3.0 munsell_0.5.0 tools_4.2.0 Pseudobulking has been tested in real scRNA-seq studies (Kang et al., 2018) and benchmarked extensively via simulation (Crowell et al., 2020). Step-by-step guide to create your volcano plot. Aggregation technique accounting for subject-level variation in DS analysis. To measure heterogeneity in expression among different groups, we assume that mean expression for gene iin subject j is influenced by R subject-specific covariates xj1,,xjR. (a) t-SNE plot shows CD66+ (turquoise) and CD66- (salmon) basal cells from single-cell RNA-seq profiling of human trachea. ## [37] gtable_0.3.3 leiden_0.4.3 future.apply_1.10.0 Volcano plot in R with seurat and ggplot. "t" : Student's t-test. For each method, the computed P-values for all genes were adjusted to control the FDR using the BenjaminiHochberg procedure (Benjamini and Hochberg, 1995). ## [9] LC_ADDRESS=C LC_TELEPHONE=C ## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3 Analysis of AT2 cells and AMs from healthy and IPF lungs. FindMarkers from Seurat returns p values as 0 for highly significant genes. (b) CD66+ basal cells were identified via detection of CEACAM5 or CEACAM6. Supplementary Figure S14 shows the results of marker detection for T cells and macrophages. Along with new functions add interactive functionality to plots, Seurat provides new accessory functions for manipulating and combining plots. ## [5] ssHippo.SeuratData_3.1.4 pbmcsca.SeuratData_3.0.0 You signed in with another tab or window. Comparison of methods for detection of CD66+ and CD66- basal cell markers from human trachea. 6b). Then, we consider the top g genes for each method, which are the g genes with the smallest adjusted P-values, and find what percentage of these top genes are known markers. We have found this particularly useful for small clusters that do not always separate using unbiased clustering, but which look tantalizingly distinct. disease and intervention), (ii) variation between subjects, (iii) variation between cells within subjects and (iv) technical variation introduced by sampling RNA molecules, library preparation and sequencing. Andrew L Thurman, Jason A Ratcliff, Michael S Chimenti, Alejandro A Pezzulo, Differential gene expression analysis for multi-subject single-cell RNA-sequencing studies with aggregateBioVar, Bioinformatics, Volume 37, Issue 19, 1 October 2021, Pages 32433251, https://doi.org/10.1093/bioinformatics/btab337. In order to determine the reliability of the unadjusted P-values computed by each method, we compared them to the unadjusted P-values obtained from a permutation test. To use, simply make a ggplot2-based scatter plot (such as DimPlot() or FeaturePlot()) and pass the resulting plot to HoverLocator(). 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 . Figure 3(b and c) show the PPV and negative predictive value (NPV) for each method and simulation setting under an adjusted P-value cutoff of 0.05. 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). In a study in which a treatment has the effect of altering the composition of cells, subjects in the treatment and control groups may have different numbers of cells of each cell type. 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). ## [100] lifecycle_1.0.3 spatstat.geom_3.1-0 lmtest_0.9-40 . ## [7] pbmcMultiome.SeuratData_0.1.2 pbmc3k.SeuratData_3.1.4 In the second stage, the observed data for each gene, measured as a count, is assumed to follow a Poisson distribution with mean equal to the product of a size factor, such as sequencing depth, and gene expression generated in the first stage. ## In stage ii, we assume that we have not measured cell-level covariates, so that variation in expression between cells of the same type occurs only through the dispersion parameter ij2. 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. The general process for detecting genes then would be: Repeat for all cell clusters/types of interest, depending on your research questions. True positives were identified as those genes in the bulk RNA-seq analysis with FDR<0.05 and |log2(CD66+/CD66)|>1. ## I used ggplot to plot the graph, but my graph is blank at the center across Log2Fc=0. The second stage represents technical variation introduced by the processes of sampling from a population of RNAs, building a cDNA library and sequencing. We then compare multiple differential expression testing methods on scRNA-seq datasets from human samples and from animal models. 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. 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. (c) Volcano plots show results of three methods (subject, wilcox and mixed) used to identify CD66+ and CD66- basal cell marker genes. Developed by Paul Hoffman, Satija Lab and Collaborators. Supplementary Figure S10 shows concordance between adjusted P-values for each method. # Calculate feature-specific contrast levels based on quantiles of non-zero expression. Search for other works by this author on: Iowa Institute of Human Genetics, Roy J. and Lucille A. (b) AT2 cells and AM express SFTPC and MARCO, respectively. As a gold standard, results from bulk RNA-seq of isolated AT2 cells and AM comparing IPF and healthy lungs (bulk). The null and alternative hypotheses for the i-th gene are H0i:i2=0 and H0i:i20, respectively. 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). First, we identified the AT2 and AM cells via clustering (Fig. The vertical axis gives the precision (PPV) and the horizontal axis gives recall (TPR). The volcano plot for the subject method shows three genes with adjusted P-value <0.05 (-log 10 (FDR) > 1.3), whereas the other six methods detected a much larger number of genes. (a) t-SNE plot shows AT2 cells (red) and AM (green) from single-cell RNA-seq profiling of human lung from healthy subjects and subjects with IPF. ## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 1 Answer. In our simulation study, we also found that the pseudobulk method was conservative, but in some settings, mixed models had inflated FDR. The subject and mixed methods show the highest ratios of inter-group to intra-group variation in gene expression, whereas the other five methods have substantial intra-group variation. Infinite p-values are set defined value of the highest . The value of pDE describes the relative number of differentially expressed genes in a simulated dataset, and the value of controls the signal-to-noise ratio. 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. 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 = Seurat utilizes Rs plotly graphing library to create interactive plots. First, it is assumed that prerequisite steps in the bioinformatic pipeline produced cells that conform to the assumptions of the proposed model. The computations for each method were performed on the high-performance computing cluster at the University of Iowa.
Vetassess Outcome Letter,
Structure Of Greek Tragedy,
What Does Owa Stand For In Alabama,
Articles F