1 Installation

To install the package, please use the following.

if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
BiocManager::install("miQC")

2 Introduction

This vignette provides a basic example of how to run miQC, which allows users to perform cell-wise filtering of single-cell RNA-seq data for quality control. Single-cell RNA-seq data is very sensitive to tissue quality and choice of experimental workflow; it’s critical to ensure compromised cells and failed cell libraries are removed. A high proportion of reads mapping to mitochondrial DNA is one sign of a damaged cell, so most analyses will remove cells with mtRNA over a certain threshold, but those thresholds can be arbitrary and/or detrimentally stringent, especially for archived tumor tissues. miQC jointly models both the proportion of reads mapping to mtDNA genes and the number of detected genes with mixture models in a probabilistic framework to identify the low-quality cells in a given dataset.

You’ll need the following packages installed to work this tutorial:

suppressPackageStartupMessages({
    library(SingleCellExperiment)
    library(scRNAseq)
    library(scater)
    library(flexmix)
    library(splines)
    library(biomaRt)
    library(miQC)
})

2.1 Example data

To demonstrate how to run miQC on a single-cell RNA-seq dataset, we’ll use data from mouse brain cells, originating from an experiment by Zeisel et al (Zeisel et al. 2015), and available through the Bioconductor package scRNAseq.

sce <- ZeiselBrainData()
sce
## class: SingleCellExperiment 
## dim: 20006 3005 
## metadata(0):
## assays(1): counts
## rownames(20006): Tspan12 Tshz1 ... mt-Rnr1 mt-Nd4l
## rowData names(1): featureType
## colnames(3005): 1772071015_C02 1772071017_G12 ... 1772066098_A12
##   1772058148_F03
## colData names(10): tissue group # ... level1class level2class
## reducedDimNames(0):
## mainExpName: endogenous
## altExpNames(2): ERCC repeat

2.2 Scater preprocessing

In order to calculate the percent of reads in a cell that map to mitochondrial genes, we first need to establish which genes are mitochondrial. For genes listed as HGNC symbols, this is as simple as searching for genes starting with mt-. For other IDs, we recommend using a biomaRt query to map to chromosomal location and identify all mitochondrial genes.

mt_genes <- grepl("^mt-",  rownames(sce))
feature_ctrls <- list(mito = rownames(sce)[mt_genes])

feature_ctrls
## $mito
##  [1] "mt-Tl1"  "mt-Tn"   "mt-Tr"   "mt-Tq"   "mt-Ty"   "mt-Tk"   "mt-Ta"  
##  [8] "mt-Tf"   "mt-Tp"   "mt-Tc"   "mt-Td"   "mt-Tl2"  "mt-Te"   "mt-Tv"  
## [15] "mt-Tg"   "mt-Tt"   "mt-Tw"   "mt-Tm"   "mt-Ti"   "mt-Nd3"  "mt-Nd6" 
## [22] "mt-Nd4"  "mt-Atp6" "mt-Nd2"  "mt-Nd5"  "mt-Nd1"  "mt-Co3"  "mt-Cytb"
## [29] "mt-Atp8" "mt-Co2"  "mt-Co1"  "mt-Rnr2" "mt-Rnr1" "mt-Nd4l"

miQC is designed to be run with the Bioconductor package scater, which has a built-in function addPerCellQC to calculate basic QC metrics like number of unique genes detected per cell and total number of reads. When we pass in our list of mitochondrial genes, it will also calculate percent mitochondrial reads.

sce <- addPerCellQC(sce, subsets = feature_ctrls)
head(colData(sce))
## DataFrame with 6 rows and 22 columns
##                     tissue   group # total mRNA mol      well       sex
##                <character> <numeric>      <numeric> <numeric> <numeric>
## 1772071015_C02    sscortex         1           1221         3         3
## 1772071017_G12    sscortex         1           1231        95         1
## 1772071017_A05    sscortex         1           1652        27         1
## 1772071014_B06    sscortex         1           1696        37         3
## 1772067065_H06    sscortex         1           1219        43         3
## 1772071017_E02    sscortex         1           1378         5         1
##                      age  diameter        cell_id  level1class level2class
##                <numeric> <numeric>    <character>  <character> <character>
## 1772071015_C02         2         1 1772071015_C02 interneurons       Int10
## 1772071017_G12         1       353 1772071017_G12 interneurons       Int10
## 1772071017_A05         1        13 1772071017_A05 interneurons        Int6
## 1772071014_B06         2        19 1772071014_B06 interneurons       Int10
## 1772067065_H06         6        12 1772067065_H06 interneurons        Int9
## 1772071017_E02         1        21 1772071017_E02 interneurons        Int9
##                      sum  detected subsets_mito_sum subsets_mito_detected
##                <numeric> <numeric>        <numeric>             <numeric>
## 1772071015_C02     22354      4871              774                    23
## 1772071017_G12     22869      4712             1121                    27
## 1772071017_A05     32594      6055              952                    27
## 1772071014_B06     33525      5852              611                    28
## 1772067065_H06     21694      4724              164                    23
## 1772071017_E02     25919      5427             1122                    19
##                subsets_mito_percent altexps_ERCC_sum altexps_ERCC_detected
##                           <numeric>        <numeric>             <numeric>
## 1772071015_C02             3.462468             6706                    43
## 1772071017_G12             4.901832             6435                    46
## 1772071017_A05             2.920783             6335                    47
## 1772071014_B06             1.822521             7032                    43
## 1772067065_H06             0.755969             5931                    39
## 1772071017_E02             4.328871             6671                    43
##                altexps_ERCC_percent altexps_repeat_sum altexps_repeat_detected
##                           <numeric>          <numeric>               <numeric>
## 1772071015_C02              18.0070               8181                     419
## 1772071017_G12              15.6349              11854                     480
## 1772071017_A05              11.1238              18021                     582
## 1772071014_B06              12.8999              13955                     512
## 1772067065_H06              17.1908               6876                     363
## 1772071017_E02              13.3543              17364                     618
##                altexps_repeat_percent     total
##                             <numeric> <numeric>
## 1772071015_C02                21.9677     37241
## 1772071017_G12                28.8012     41158
## 1772071017_A05                31.6435     56950
## 1772071014_B06                25.5999     54512
## 1772067065_H06                19.9299     34501
## 1772071017_E02                34.7600     49954

3 miQC

3.1 Basic usage

Using the QC metrics generated via scater, we can use the plotMetrics function to visually inspect the quality of our dataset.

plotMetrics(sce)

We can see that most cells have a fairly low proportion of mitochondrial reads, given that the graph is much denser at the bottom. We likely have many cells that are intact and biologically meaningful. There are also a few cells that have almost half of their reads mapping to mitochondrial genes, which are likely broken or otherwise compromised and we will want to exclude from our downstream analysis. However, it’s not clear what boundaries to draw to separate the two groups of cells. With that in mind, we’ll generate a linear mixture model using the mixtureModel function.

model <- mixtureModel(sce)

This function is a wrapper for flexmix, which fits a mixture model on our data and returns the parameters of the two lines that best fit the data, as well as the posterior probability of each cell being derived from each distribution.

We can look at the parameters and posterior values directly with the functions

parameters(model)
##                       Comp.1        Comp.2
## coef.(Intercept) 26.89445262  9.5036977900
## coef.detected    -0.00375561 -0.0008344042
## sigma             6.48969202  3.3055823842
head(posterior(model))
##            [,1]      [,2]
## [1,] 0.10549891 0.8945011
## [2,] 0.09961569 0.9003843
## [3,] 0.12855554 0.8714445
## [4,] 0.14700508 0.8529949
## [5,] 0.14392177 0.8560782
## [6,] 0.11502350 0.8849765

Or we can visualize the model results using the plotModel function:

plotModel(sce, model)

As expected, the cells at the very top of the graph are almost certainly compromised, most likely to have been derived from the distribution with fewer unique genes and higher baseline mitochondrial expression.

We can use these posterior probabilities to choose which cells to keep, and visualize the consequences of this filtering with the plotFiltering function.

plotFiltering(sce, model)

To actually perform the filtering and remove the indicated cells from our SingleCellExperiment object, we can run the filterCells parameter.

sce <- filterCells(sce, model)
## Removing 359 out of 3005 cells.
sce
## class: SingleCellExperiment 
## dim: 20006 2646 
## metadata(0):
## assays(1): counts
## rownames(20006): Tspan12 Tshz1 ... mt-Rnr1 mt-Nd4l
## rowData names(1): featureType
## colnames(2646): 1772071015_C02 1772071017_G12 ... 1772066097_D04
##   1772066098_A12
## colData names(23): tissue group # ... total prob_compromised
## reducedDimNames(0):
## mainExpName: endogenous
## altExpNames(2): ERCC repeat

3.2 Other model types

In most cases, a linear mixture model will be satisfactory as well as simpler, but miQC also supports some non-linear mixture models: currently polynomials and b-splines. A user should only need to change the model_type parameter when making the model, and all visualization and filtering functions will work the same as with a linear model.

sce <- ZeiselBrainData()
sce <- addPerCellQC(sce, subsets = feature_ctrls)

model2 <- mixtureModel(sce, model_type = "spline")
plotModel(sce, model2)

plotFiltering(sce, model2)

model3 <- mixtureModel(sce, model_type = "polynomial")
plotModel(sce, model3)

plotFiltering(sce, model3)

miQC is designed to combine information about mitochondrial percentage and library complexity (number of genes discovered) to make filtering decisions, but if an even simpler model is preferred, miQC can make a model based only on mitochondrial information. This can be done by changing the model_type parameter to “one_dimensional”, which runs a 1D gaussian mixture model.

model4 <- mixtureModel(sce, model_type = "one_dimensional")
plotModel(sce, model4)

plotFiltering(sce, model4)

3.3 Extras

3.3.1 Changing posterior cutoff

miQC defaults to removing any cell with 75% or greater posterior probability of being compromised, but if we want to be more or less stringent, we can alter the posterior_cutoff parameter, like so:

plotFiltering(sce, model4, posterior_cutoff = 0.9)

Note that when performing miQC multiple times on different samples for the same experiment, it’s recommended to select the same posterior_cutoff for all, to give consistency in addition to the flexibility of sample-specific models.

3.3.2 Preventing exclusion of low-mito cells

miQC includes two parameters to accomodate unusual and undesired behavior in the linear distributions. These issues are especially visible in some cancer datasets with a stringent posterior cutoff. We’ve included a bare-bones version of QC data for a high-grade serous ovarian tumor (full version of the data is available at GEO, accession GSM4816047).

data("hgsoc_metrics")
sce <- SingleCellExperiment(colData = metrics)
model <- mixtureModel(sce)
plotFiltering(sce, model, posterior_cutoff = 0.6, enforce_left_cutoff = FALSE,
                keep_all_below_boundary = FALSE)

The first issue is the group of cells at the bottom of the distribution getting marked for removal. These cells happen to be near the x-intercept of the compromised cell line, which increases their posterior probability of being compromised. But since they have decent library complexity and a low mitochondrial percentage, so it doesn’t make biological sense to exclude them. When the keep_all_below_boundary parameter is set to True, as is the default, any cells below the intact cell line are kept:

plotFiltering(sce, model, posterior_cutoff = 0.6, enforce_left_cutoff = FALSE,
                keep_all_below_boundary = TRUE)

3.3.3 Preventing U-shaped boundary

The second issue is the U-shape in the boundary between kept and discarded cells. When this occurs, there will be cells at the bottom of the trough that are discarded, but some cells with less library complexity (farther left) and higher percentage of mitochondrial reads (higher) – meaning they are worse in both of our QC metrics – will be kept. To avoid this happening, the parameter enforce_left_cutoff will identify the cell marked for removal with the lowest mitochondrial percentage, determine its library complexity, and discard all cells with both lower complexity and higher mitochondrial percentage:

plotFiltering(sce, model, posterior_cutoff = 0.6, enforce_left_cutoff = TRUE,
                keep_all_below_boundary = TRUE)

This will make a de facto mitochondrial percentage cutoff for all cells with low library complexity, but will be more permissive for cells with high library complexity and high mitochondrial percentage, which are more likely to be intact cells with a biological reason for high mitochondrial expression than their low library complexity counterparts.

3.4 When not to use miQC

The miQC model is based on the assumption that there are a non-trivial number of compromised cells in the dataset, which is not true in all datasets. We recommend running plotMetrics on a dataset before running miQC to see if the two-distribution model is appropriate. Look for the distinctive triangular shape where cells have a wide variety of mitochondrial percentages at lower gene counts and taper off to lower mitochondrial percentage at higher gene counts, as can be seen in the Zeisel data.

For example of a dataset where there’s not a significant number of compromised cells, so the two-distribution assumption is not met, look at another dataset from the scRNAseq package, mouse data from Lun et al (Lun et al. 2017). Note that this dataset uses ensembl IDs, so we will use biomaRt to determine which genes are mitochondrial and calculate mitochondrial percentage.

sce <- LunSpikeInData()
mouse_mart <- useEnsembl("ENSEMBL_MART_ENSEMBL", "mmusculus_gene_ensembl")
id_map <- getBM(values = rownames(sce),
                filters = "ensembl_gene_id",
                attributes = c("ensembl_gene_id", "chromosome_name"),
                mart = mouse_mart)
mt_genes <- subset(id_map,id_map$chromosome_name=="MT")$ensembl_gene_id
feature_ctrls <- list(mito = mt_genes)
sce <- addPerCellQC(sce, subsets = feature_ctrls)

plotMetrics(sce)

The mixtureModel function will throw a warning if only one distribution is found, in which case no cells would be filtered. In these cases, we recommend using other filtering methods, such as a cutoff on mitochondrial percentage or identifying outliers using median absolute deviation (MAD).

4 Session Information

## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.14-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.14-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_GB              LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] splines   stats4    stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ensembldb_2.18.0            AnnotationFilter_1.18.0     GenomicFeatures_1.46.0      AnnotationDbi_1.56.0       
##  [5] miQC_1.2.0                  biomaRt_2.50.0              flexmix_2.3-17              lattice_0.20-45            
##  [9] scater_1.22.0               ggplot2_3.3.5               scuttle_1.4.0               scRNAseq_2.7.2             
## [13] SingleCellExperiment_1.16.0 SummarizedExperiment_1.24.0 Biobase_2.54.0              GenomicRanges_1.46.0       
## [17] GenomeInfoDb_1.30.0         IRanges_2.28.0              S4Vectors_0.32.0            BiocGenerics_0.40.0        
## [21] MatrixGenerics_1.6.0        matrixStats_0.61.0          BiocStyle_2.22.0           
## 
## loaded via a namespace (and not attached):
##   [1] AnnotationHub_3.2.0           BiocFileCache_2.2.0           lazyeval_0.2.2               
##   [4] BiocParallel_1.28.0           digest_0.6.28                 htmltools_0.5.2              
##   [7] magick_2.7.3                  viridis_0.6.2                 fansi_0.5.0                  
##  [10] magrittr_2.0.1                memoise_2.0.0                 ScaledMatrix_1.2.0           
##  [13] Biostrings_2.62.0             prettyunits_1.1.1             colorspace_2.0-2             
##  [16] blob_1.2.2                    rappdirs_0.3.3                ggrepel_0.9.1                
##  [19] xfun_0.27                     dplyr_1.0.7                   crayon_1.4.1                 
##  [22] RCurl_1.98-1.5                jsonlite_1.7.2                glue_1.4.2                   
##  [25] gtable_0.3.0                  zlibbioc_1.40.0               XVector_0.34.0               
##  [28] DelayedArray_0.20.0           BiocSingular_1.10.0           scales_1.1.1                 
##  [31] DBI_1.1.1                     Rcpp_1.0.7                    viridisLite_0.4.0            
##  [34] xtable_1.8-4                  progress_1.2.2                bit_4.0.4                    
##  [37] rsvd_1.0.5                    httr_1.4.2                    ellipsis_0.3.2               
##  [40] modeltools_0.2-23             pkgconfig_2.0.3               XML_3.99-0.8                 
##  [43] farver_2.1.0                  nnet_7.3-16                   sass_0.4.0                   
##  [46] dbplyr_2.1.1                  utf8_1.2.2                    tidyselect_1.1.1             
##  [49] labeling_0.4.2                rlang_0.4.12                  later_1.3.0                  
##  [52] munsell_0.5.0                 BiocVersion_3.14.0            tools_4.1.1                  
##  [55] cachem_1.0.6                  generics_0.1.1                RSQLite_2.2.8                
##  [58] ExperimentHub_2.2.0           evaluate_0.14                 stringr_1.4.0                
##  [61] fastmap_1.1.0                 yaml_2.2.1                    knitr_1.36                   
##  [64] bit64_4.0.5                   purrr_0.3.4                   KEGGREST_1.34.0              
##  [67] sparseMatrixStats_1.6.0       mime_0.12                     xml2_1.3.2                   
##  [70] compiler_4.1.1                beeswarm_0.4.0                filelock_1.0.2               
##  [73] curl_4.3.2                    png_0.1-7                     interactiveDisplayBase_1.32.0
##  [76] tibble_3.1.5                  bslib_0.3.1                   stringi_1.7.5                
##  [79] highr_0.9                     ProtGenerics_1.26.0           Matrix_1.3-4                 
##  [82] vctrs_0.3.8                   pillar_1.6.4                  lifecycle_1.0.1              
##  [85] BiocManager_1.30.16           jquerylib_0.1.4               BiocNeighbors_1.12.0         
##  [88] bitops_1.0-7                  irlba_2.3.3                   httpuv_1.6.3                 
##  [91] rtracklayer_1.54.0            R6_2.5.1                      BiocIO_1.4.0                 
##  [94] bookdown_0.24                 promises_1.2.0.1              gridExtra_2.3                
##  [97] vipor_0.4.5                   assertthat_0.2.1              rjson_0.2.20                 
## [100] withr_2.4.2                   GenomicAlignments_1.30.0      Rsamtools_2.10.0             
## [103] GenomeInfoDbData_1.2.7        parallel_4.1.1                hms_1.1.1                    
## [106] grid_4.1.1                    beachmat_2.10.0               rmarkdown_2.11               
## [109] DelayedMatrixStats_1.16.0     shiny_1.7.1                   ggbeeswarm_0.6.0             
## [112] restfulr_0.0.13

References

Lun, Aaron T. L., Fernando J. Calero-Nieto, Liora Haim-Vilmovsky, Berthold Göttgens, and John C. Marioni. 2017. “Assessing the Reliability of Spike-in Normalization for Analyses of Single-Cell RNA Sequencing Data.” Genome Research 27 (11): 1795–1806. https://doi.org/10.1101/gr.222877.117.

Zeisel, Amit, Ana B. Muñoz-Manchado, Simone Codeluppi, Peter Lönnerberg, Gioele La Manno, Anna Juréus, Sueli Marques, et al. 2015. “Brain Structure. Cell Types in the Mouse Cortex and Hippocampus Revealed by Single-Cell RNA-Seq.” Science (New York, N.Y.) 347 (6226). https://doi.org/10.1126/science.aaa1934.