Single Cells with edgeR

This vignette demonstrates usage of Glimma on a single cell dataset. The data here comes from brain cells from the Zeisel A et al. (2015) study on mouse brain single cells. We will use the MDS plot to perform unsupervised clustering of the cells. A pseudo-bulk single cell aggregation approach with edgeR will be used to test for differential expression, and two styles of MA plots will be used to investigate the results.

This is a simplified workflow not intended to represent best practice, but to produce reasonable looking plots in the minimal amount of code. Please refer to a resource such Orchestrating Single-Cell Analysis with Bioconductor (OSCA) for appropriate workflows for analysis.

We start by loading in the data using the scRNAseq package.

library(scRNAseq)
library(scater)
library(scran)
library(Glimma)
library(edgeR)

sce <- ZeiselBrainData(ensembl=TRUE)
#> snapshotDate(): 2022-04-19
#> see ?scRNAseq and browseVignettes('scRNAseq') for documentation
#> loading from cache
#> see ?scRNAseq and browseVignettes('scRNAseq') for documentation
#> loading from cache
#> see ?scRNAseq and browseVignettes('scRNAseq') for documentation
#> loading from cache
#> snapshotDate(): 2022-04-19
#> see ?scRNAseq and browseVignettes('scRNAseq') for documentation
#> loading from cache
#> snapshotDate(): 2022-04-21
#> loading from cache
#> require("ensembldb")
#> Warning: Unable to map 1565 of 20006 requested IDs.

Once the data is loaded in we follow the OSCA procedure for identifying highly variable genes for creating a multi-dimensional scaling (MDS) plot. We use the functions provided by scran to identify the most highly variable genes rather than the algorithm within glimmaMDS, as scran is tailored towards single cells.

sce <- logNormCounts(sce)

var_mod <- modelGeneVar(sce)
hvg_genes <- getTopHVGs(var_mod, n=500)
hvg_sce <- sce[hvg_genes, ]

hvg_sce <- logNormCounts(hvg_sce)

Choosing to colour the MDS plot using level1class reveals separation between cell types.

glimmaMDS(
    exprs(hvg_sce),
    groups = colData(hvg_sce)
)

To demonstrate the MA plot we will perform a differential expression analysis using the pseudo-bulk approach. This involves creating pseudo-bulk samples by aggregating single cells as an analogue of biological replicates. Here the pseudo-bulk samples will be generated from combinations of level1class and level2class, the cells belonging to unique combinations of the two factors will be aggregated into samples.

colData(sce)$pb_group <-
    paste0(colData(sce)$level1class,
           "_",
           colData(sce)$level2class)

sce_counts <- counts(sce)
pb_counts <- t(rowsum(t(sce_counts), colData(sce)$pb_group))

pb_samples <- colnames(pb_counts)
pb_samples <- gsub("astrocytes_ependymal", "astrocytes-ependymal", pb_samples)
pb_split <- do.call(rbind, strsplit(pb_samples, "_"))
pb_sample_anno <- data.frame(
    sample = pb_samples,
    cell_type = pb_split[, 1],
    sample_group = pb_split[, 2]
)

With the pseudo-bulk annotations and counts we can construct a DGEList object.

pb_dge <- DGEList(
    counts = pb_counts,
    samples = pb_sample_anno,
    group = pb_sample_anno$cell_type
)

pb_dge <- calcNormFactors(pb_dge)

With this we perform differential expression analysis between “pyramidal SS” and “pyramidal CA1” samples using edgeR’s generalised linear models.

design <- model.matrix(~0 + cell_type, data = pb_dge$samples)
colnames(design) <- make.names(gsub("cell_type", "", colnames(design)))

pb_dge <- estimateDisp(pb_dge, design)

contr <- makeContrasts("pyramidal.SS - pyramidal.CA1", levels = design)

pb_fit <- glmFit(pb_dge, design)
pb_lrt <- glmLRT(pb_fit, contrast = contr)

The results of this analysis can be visualised using glimmaMA() as it would be for bulk RNA-seq.

glimmaMA(pb_lrt, dge = pb_dge)

An alternative view of the data can be constructed using the single cells in the expression plot rather than the pseudo-bulk samples. Since the MA plot is related to the expressions by only the genes in the rows, another expression matrix containing the same genes can be substituted in as below.

We construct a new DGE list from the raw single cell counts, then filter it down to just the cells used in our comparison and further down-sampled to 100 cells. This is done because Glimma does not handle a large number of cells well, the limit being a few hundred for most computers. Sampling still provides an approximate representation of the data without computation strain.

The code is not evaluated here to keep the vignette compact.

sc_dge <- DGEList(
    counts = sce_counts,
    group = colData(sce)$level1class
)

sc_dge <- sc_dge[, colData(sce)$level1class %in% c("pyramidal CA1", "pyramidal SS")]

glimmaMA(
    pb_lrt,
    dge = sc_dge[, sample(1:ncol(sc_dge), 100)]
)

Session Info

sessionInfo()
#> R version 4.2.0 RC (2022-04-19 r82224)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 20.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.15-bioc/R/lib/libRblas.so
#> LAPACK: /home/biocbuild/bbs-3.15-bioc/R/lib/libRlapack.so
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_GB              LC_COLLATE=C              
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] ensembldb_2.20.0            AnnotationFilter_1.20.0    
#>  [3] GenomicFeatures_1.48.0      AnnotationDbi_1.58.0       
#>  [5] AnnotationHub_3.4.0         BiocFileCache_2.4.0        
#>  [7] dbplyr_2.1.1                scran_1.24.0               
#>  [9] scater_1.24.0               ggplot2_3.3.5              
#> [11] scuttle_1.6.0               scRNAseq_2.9.2             
#> [13] SingleCellExperiment_1.18.0 DESeq2_1.36.0              
#> [15] SummarizedExperiment_1.26.0 Biobase_2.56.0             
#> [17] MatrixGenerics_1.8.0        matrixStats_0.62.0         
#> [19] GenomicRanges_1.48.0        GenomeInfoDb_1.32.0        
#> [21] IRanges_2.30.0              S4Vectors_0.34.0           
#> [23] BiocGenerics_0.42.0         edgeR_3.38.0               
#> [25] limma_3.52.0                Glimma_2.6.0               
#> 
#> loaded via a namespace (and not attached):
#>   [1] igraph_1.3.1                  lazyeval_0.2.2               
#>   [3] splines_4.2.0                 BiocParallel_1.30.0          
#>   [5] digest_0.6.29                 htmltools_0.5.2              
#>   [7] viridis_0.6.2                 fansi_1.0.3                  
#>   [9] magrittr_2.0.3                memoise_2.0.1                
#>  [11] ScaledMatrix_1.4.0            cluster_2.1.3                
#>  [13] Biostrings_2.64.0             annotate_1.74.0              
#>  [15] prettyunits_1.1.1             colorspace_2.0-3             
#>  [17] blob_1.2.3                    rappdirs_0.3.3               
#>  [19] ggrepel_0.9.1                 xfun_0.30                    
#>  [21] dplyr_1.0.8                   crayon_1.5.1                 
#>  [23] RCurl_1.98-1.6                jsonlite_1.8.0               
#>  [25] genefilter_1.78.0             survival_3.3-1               
#>  [27] glue_1.6.2                    gtable_0.3.0                 
#>  [29] zlibbioc_1.42.0               XVector_0.36.0               
#>  [31] DelayedArray_0.22.0           BiocSingular_1.12.0          
#>  [33] scales_1.2.0                  DBI_1.1.2                    
#>  [35] Rcpp_1.0.8.3                  viridisLite_0.4.0            
#>  [37] xtable_1.8-4                  progress_1.2.2               
#>  [39] dqrng_0.3.0                   bit_4.0.4                    
#>  [41] rsvd_1.0.5                    metapod_1.4.0                
#>  [43] htmlwidgets_1.5.4             httr_1.4.2                   
#>  [45] RColorBrewer_1.1-3            ellipsis_0.3.2               
#>  [47] pkgconfig_2.0.3               XML_3.99-0.9                 
#>  [49] sass_0.4.1                    locfit_1.5-9.5               
#>  [51] utf8_1.2.2                    tidyselect_1.1.2             
#>  [53] rlang_1.0.2                   later_1.3.0                  
#>  [55] munsell_0.5.0                 BiocVersion_3.15.2           
#>  [57] tools_4.2.0                   cachem_1.0.6                 
#>  [59] cli_3.3.0                     generics_0.1.2               
#>  [61] RSQLite_2.2.12                ExperimentHub_2.4.0          
#>  [63] evaluate_0.15                 stringr_1.4.0                
#>  [65] fastmap_1.1.0                 yaml_2.3.5                   
#>  [67] knitr_1.38                    bit64_4.0.5                  
#>  [69] purrr_0.3.4                   KEGGREST_1.36.0              
#>  [71] sparseMatrixStats_1.8.0       mime_0.12                    
#>  [73] xml2_1.3.3                    biomaRt_2.52.0               
#>  [75] compiler_4.2.0                beeswarm_0.4.0               
#>  [77] filelock_1.0.2                curl_4.3.2                   
#>  [79] png_0.1-7                     interactiveDisplayBase_1.34.0
#>  [81] statmod_1.4.36                tibble_3.1.6                 
#>  [83] geneplotter_1.74.0            bslib_0.3.1                  
#>  [85] stringi_1.7.6                 bluster_1.6.0                
#>  [87] lattice_0.20-45               ProtGenerics_1.28.0          
#>  [89] Matrix_1.4-1                  vctrs_0.4.1                  
#>  [91] pillar_1.7.0                  lifecycle_1.0.1              
#>  [93] BiocManager_1.30.17           jquerylib_0.1.4              
#>  [95] BiocNeighbors_1.14.0          bitops_1.0-7                 
#>  [97] irlba_2.3.5                   httpuv_1.6.5                 
#>  [99] rtracklayer_1.56.0            R6_2.5.1                     
#> [101] BiocIO_1.6.0                  promises_1.2.0.1             
#> [103] gridExtra_2.3                 vipor_0.4.5                  
#> [105] assertthat_0.2.1              rjson_0.2.21                 
#> [107] withr_2.5.0                   GenomicAlignments_1.32.0     
#> [109] Rsamtools_2.12.0              GenomeInfoDbData_1.2.8       
#> [111] parallel_4.2.0                hms_1.1.1                    
#> [113] grid_4.2.0                    beachmat_2.12.0              
#> [115] rmarkdown_2.14                DelayedMatrixStats_1.18.0    
#> [117] shiny_1.7.1                   ggbeeswarm_0.6.0             
#> [119] restfulr_0.0.13