1 Introduction

We will use two sets of Bone Marrow Mononuclear Cells (BMMC):

Fastq files were retrieved from publicly available archive (SRA and HCA).

Sequencing quality was assessed and visualised using fastQC and MultiQC.

Reads were aligned against GRCh38 and features counted using cellranger (v3.1.0).

We will now check the quality of the data further:

We will then:

2 Load packages

library(DropletUtils)
library(scater)
library(ensembldb)
library(AnnotationHub)
library(BiocParallel)
library(tidyverse)
library(patchwork)
library(ggvenn)

3 Reading CellRanger output into R

3.1 Sample meta data

We will load both the Caron and HCA data sets. We have already prepared a sample meta data table that relates the sample/run ID to the sample group.

samplesheet <- read_tsv("Data/sample_sheet.tsv")

NOTE: There are two PBMMC_1 samples. These are two libraries from the same sample material. Later on, we will combine these to form a single sample.

3.1.1 Parallelisation

scRNAseq data sets tend to be very large and processing them requires a large amount of computing power and can take time. Many of the commands we will use have the option to be run in parallels across multiple processors. By default they will only use a single processor, but parallelisation will greatly speed up the analysis.

We need to first set up some parallel parameters using the package BiocParallel.

bp.params <- MulticoreParam(workers = 7)

Here were are selecting to use forked processes with MulticoreParam and instructing the function to use 7 cores (our machines have 8, this leaves 1 to run the desktop etc.) Note that on Windows MulticoreParam does not work and it is necessary to use SnowParam - please refer to the BiocParallel vignettes for further information.

3.1.2 Loading a single sample

We will load the data for the SRR9264343. To load the data from the CellRanger outputs, we need to use the function read10xCounts from the DropletUtils package. We pass the function the location of the directory containing the counts matrix, cell barcodes and features (genes).

We could load the raw data, which includes counts for all cell barcodes detected in the sample, and use the emptyDrops function in DropletUtils to call cells, however, CellRanger has already called cells and so we are going to work with the filtered matrix, which only contains droplets called as cells by CellRanger.

sample.path <- "CellRanger_Outputs/SRR9264343/outs/filtered_feature_bc_matrix/"
sce.sing <- read10xCounts(sample.path, col.names=TRUE, BPPARAM = bp.params)
sce.sing
## class: SingleCellExperiment 
## dim: 36601 3094 
## metadata(1): Samples
## assays(1): counts
## rownames(36601): ENSG00000243485 ENSG00000237613 ... ENSG00000278817
##   ENSG00000277196
## rowData names(3): ID Symbol Type
## colnames(3094): AAACCTGAGACTTTCG-1 AAACCTGGTCTTCAAG-1 ...
##   TTTGTCACAGGCTCAC-1 TTTGTCAGTTCGGCAC-1
## colData names(2): Sample Barcode
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):

3.1.3 The SingleCellExperiment object

The data has been loaded as a SingleCellExperiment object. The details of the structure of the object are described here. In summary, it stores various data types in a single object. Currently it will contain:

  • the count matrix
  • feature (gene) metadata
  • cell (droplet) metadata

Later we will also add the outcomes of downstream analysis such as dimensionality reduction.

3.1.4 The counts matrix

Compared to bulk RNA-seq, Single-cell RNA-seq data is sparse, i.e. there many missing values or zeroes. This is particularly true with droplet-based methods such as 10X, mostly because:

  • any given cell does not express each gene
  • the library preparation does not capture all transcript the cell does express
  • the sequencing depth per cell is far lower and so fewer of the expressed genes are detected

We can access the counts matrix with counts. Given the large number of droplets in a sample, count matrices can be large.

dim(counts(sce.sing))
## [1] 36601  3094

They are however very sparse, that is, most of the entries are 0’s. To save memory the counts can be stored in a ‘sparse matrix’ that only stores non-zero values, for example a dgCMatrix object.

counts(sce.sing)[1:10, 1:10]
## 10 x 10 sparse Matrix of class "dgCMatrix"
##                                    
## ENSG00000243485 . . . . . . . . . .
## ENSG00000237613 . . . . . . . . . .
## ENSG00000186092 . . . . . . . . . .
## ENSG00000238009 . . . . . . . . . .
## ENSG00000239945 . . . . . . . . . .
## ENSG00000239906 . . . . . . . . . .
## ENSG00000241860 . . . . . . . . . .
## ENSG00000241599 . . . . . . . . . .
## ENSG00000286448 . . . . . . . . . .
## ENSG00000236601 . . . . . . . . . .

3.1.5 Features

Details about the “features” (in this case genes) can by accessed using the rowData function. Currently it contains the ensembl gene ID and the gene symbol, which have been derived from the 10x reference used by CellRanger. It also contains a “Type” column, which tells us what sort of data we are looking at; in this case it is “Expression” for gene expression. If we wish to, we can add further annotation to the features by adding extra columns to this data frame.

rowData(sce.sing)
## DataFrame with 36601 rows and 3 columns
##                              ID      Symbol            Type
##                     <character> <character>     <character>
## ENSG00000243485 ENSG00000243485 MIR1302-2HG Gene Expression
## ENSG00000237613 ENSG00000237613     FAM138A Gene Expression
## ENSG00000186092 ENSG00000186092       OR4F5 Gene Expression
## ENSG00000238009 ENSG00000238009  AL627309.1 Gene Expression
## ENSG00000239945 ENSG00000239945  AL627309.3 Gene Expression
## ...                         ...         ...             ...
## ENSG00000277836 ENSG00000277836  AC141272.1 Gene Expression
## ENSG00000278633 ENSG00000278633  AC023491.2 Gene Expression
## ENSG00000276017 ENSG00000276017  AC007325.1 Gene Expression
## ENSG00000278817 ENSG00000278817  AC007325.4 Gene Expression
## ENSG00000277196 ENSG00000277196  AC007325.2 Gene Expression

3.1.6 Droplet annotation

Details about the droplets can be accessed using colData. Currently it contains the sample names and droplet Barcodes. As with the feature data, we can add additional information about each droplet, e.g. counts of genes or the percentage of mitochondrial genes, to this data frame. The rows of this table correspond to the data in the columns of the count matrix; the row names of this table will match the column names of the counts matrix - currently these are the droplet barcodes.

colData(sce.sing)
## DataFrame with 3094 rows and 2 columns
##                                    Sample            Barcode
##                               <character>        <character>
## AAACCTGAGACTTTCG-1 CellRanger_Outputs/S.. AAACCTGAGACTTTCG-1
## AAACCTGGTCTTCAAG-1 CellRanger_Outputs/S.. AAACCTGGTCTTCAAG-1
## AAACCTGGTGCAACTT-1 CellRanger_Outputs/S.. AAACCTGGTGCAACTT-1
## AAACCTGGTGTTGAGG-1 CellRanger_Outputs/S.. AAACCTGGTGTTGAGG-1
## AAACCTGTCCCAAGTA-1 CellRanger_Outputs/S.. AAACCTGTCCCAAGTA-1
## ...                                   ...                ...
## TTTGGTTTCTTTAGGG-1 CellRanger_Outputs/S.. TTTGGTTTCTTTAGGG-1
## TTTGTCAAGAAACGAG-1 CellRanger_Outputs/S.. TTTGTCAAGAAACGAG-1
## TTTGTCAAGGACGAAA-1 CellRanger_Outputs/S.. TTTGTCAAGGACGAAA-1
## TTTGTCACAGGCTCAC-1 CellRanger_Outputs/S.. TTTGTCACAGGCTCAC-1
## TTTGTCAGTTCGGCAC-1 CellRanger_Outputs/S.. TTTGTCAGTTCGGCAC-1
colnames(counts(sce.sing))[1:6]
## [1] "AAACCTGAGACTTTCG-1" "AAACCTGGTCTTCAAG-1" "AAACCTGGTGCAACTT-1"
## [4] "AAACCTGGTGTTGAGG-1" "AAACCTGTCCCAAGTA-1" "AAACCTGTCGAATGCT-1"

4 Properties of scRNA-seq data

4.1 Number of genes detected per cell

The number and identity of genes detected in a cell vary greatly across cells: the total number of genes detected across all cells is far larger than the number of genes detected in each cell.

For the current set of samples the total number of genes detected across cells was 19938 out of 36601 gene in the reference, but if we look at the number of genes detected in each cell, we can see that this ranges from 35 to 4131, with a median of 1771.

genesPerCell <- colSums(counts(sce.sing) > 0)
summary(genesPerCell)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      35    1354    1771    1760    2181    4131
plot(density(genesPerCell), main="", xlab="Genes per cell")

4.2 Total UMI for a gene versus the number of times detected

If we compare the number of UMI’s assigned to an individual gene to the number of cells in which that gene is detected, we can see that highly expressed genes tend to be detected in a higher proportion of cells than lowly expressed genes.

tmpCounts <- counts(sce.sing)[,1:1000]

plot(rowSums(tmpCounts),
     rowMeans(tmpCounts > 0),
     log = "x",
     xlab="total number of UMIs",
     ylab="proportion of cells expressing the gene"
)

rm(tmpCounts)

4.3 Distribution of counts for a gene across cells

We could also look at the distribution of counts for individual genes across all cells. The plot below shows this distribution for the top 20 genes detected.

rel_expression <- t( t(counts(sce.sing)) / colSums(counts(sce.sing))) * 100
rownames(rel_expression) <- rowData(sce.sing)$Symbol
most_expressed <- sort(rowSums( rel_expression ),T)[20:1]
plot_data <- as.matrix(t(rel_expression[names(most_expressed),]))

boxplot(plot_data, cex=0.1, las=1, xlab="% total count per cell", horizontal=TRUE)

5 Quality control

The cell calling performed by CellRanger does not always retain only droplets containing cells. Poor-quality cells, or rather droplets, may be caused by cell damage during dissociation or failed library preparation. They usually have low UMI counts, few genes detected and/or high mitochondrial content. The presence of these droplets in the data set may affect normalisation, assessment of cell population heterogeneity, clustering and trajectory:

In order to remove or reduce the impact of poor-quality droplets on our downstream analysis we will attempt to filter them out using some QC metrics. The three principle means of doing this are to apply thresholds for inclusion on three characteristics:

The scater function addPerCellQC() will compute various per droplet QC metrics and will add this information as new columns in the droplet annotation (colData) of the single cell object.

5.1 Load multiple samples

We can load multiple samples at the same time using the read10xCounts command. This will create a single object containing the data for multiple samples. We can then QC and filter the samples in conjunction. As we will see later, this is not always optimal when samples have been processed in multiple batches.

As an example we will one sample from each sample group. Again we will start with the filtered counts matrix, which only contains cells called by CellRanger. We pass the read10xCounts a named vector containing the paths to the filtered counts matrices that we wish to load; the names of the vector will be used as the sample names in the Single Cell Experiment object.

samples_list <- samplesheet %>% 
    group_by(SampleGroup) %>%  
    slice(1) %>%  
    pull(SampleId)
list_of_files <- str_c("CellRanger_Outputs/", 
                       samples_list, 
                       "/outs/filtered_feature_bc_matrix")
names(list_of_files) <- samples_list
list_of_files
##                                                       MantonBM1 
##  "CellRanger_Outputs/MantonBM1/outs/filtered_feature_bc_matrix" 
##                                                      SRR9264343 
## "CellRanger_Outputs/SRR9264343/outs/filtered_feature_bc_matrix" 
##                                                      SRR9264347 
## "CellRanger_Outputs/SRR9264347/outs/filtered_feature_bc_matrix" 
##                                                      SRR9264351 
## "CellRanger_Outputs/SRR9264351/outs/filtered_feature_bc_matrix" 
##                                                      SRR9264349 
## "CellRanger_Outputs/SRR9264349/outs/filtered_feature_bc_matrix"
sce <- read10xCounts(list_of_files, col.names=TRUE, BPPARAM = bp.params)
sce
## class: SingleCellExperiment 
## dim: 36601 36447 
## metadata(1): Samples
## assays(1): counts
## rownames(36601): ENSG00000243485 ENSG00000237613 ... ENSG00000278817
##   ENSG00000277196
## rowData names(3): ID Symbol Type
## colnames(36447): 1_AAACCTGAGAAACCTA-1 1_AAACCTGAGACAAGCC-1 ...
##   5_TTTGTCAGTCATTAGC-1 5_TTTGTCAGTTGTGGCC-1
## colData names(2): Sample Barcode
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):

5.2 Modify the droplet annotation

Currently, the droplet annotation in colData slot of the sce object has two columns: “Sample” and “Barcode”. The “Sample” is the name of the sample as we provided it to read10xCounts, the “Barcode” is the barcode for the droplet (cell).

colData(sce)
## DataFrame with 36447 rows and 2 columns
##                           Sample            Barcode
##                      <character>        <character>
## 1_AAACCTGAGAAACCTA-1   MantonBM1 AAACCTGAGAAACCTA-1
## 1_AAACCTGAGACAAGCC-1   MantonBM1 AAACCTGAGACAAGCC-1
## 1_AAACCTGAGAGACTTA-1   MantonBM1 AAACCTGAGAGACTTA-1
## 1_AAACCTGAGCAACGGT-1   MantonBM1 AAACCTGAGCAACGGT-1
## 1_AAACCTGAGCAGCGTA-1   MantonBM1 AAACCTGAGCAGCGTA-1
## ...                          ...                ...
## 5_TTTGGTTTCACGACTA-1  SRR9264349 TTTGGTTTCACGACTA-1
## 5_TTTGTCACAACTGCGC-1  SRR9264349 TTTGTCACAACTGCGC-1
## 5_TTTGTCACATGTTGAC-1  SRR9264349 TTTGTCACATGTTGAC-1
## 5_TTTGTCAGTCATTAGC-1  SRR9264349 TTTGTCAGTCATTAGC-1
## 5_TTTGTCAGTTGTGGCC-1  SRR9264349 TTTGTCAGTTGTGGCC-1

The “Barcode” column contains the cell/droplet barcode and comprises the actual sequence and a ‘group ID’, e.g. AAACCTGAGAAACCAT-1. The ‘group ID’ helps distinguish cells from different samples that have identical barcode sequences, however, as each sample was processed separately with CellRanger, the group ID is set to 1 in all data sets. Some of our downstream tools will use this column to identify different droplets, so we will need modify these to be unique.

In order to distinguish droplets that originate from different samples but have the same barcode, read10XCounts has added the “index” of the sample in the list_of_files object (1-10) to the beginning of the cell barcode in the row names of the droplet annotation table and the column names of the count matrix. We will use this number to modify the “Barcode” column.

We will also switch the “Sample” column to be the sample name and add information from the sample sheet to the droplet annotation.

colData(sce) <- colData(sce) %>% 
    as.data.frame() %>%
    rownames_to_column("RowName") %>% 
    mutate(SampleNum = str_extract(RowName, "^[0-9]+")) %>%
    mutate(Barcode = str_replace(Barcode, "1$", SampleNum)) %>%
    left_join(samplesheet, by=c(Sample="SampleId")) %>%
    rename(SampleId=Sample) %>% 
    rename(Sample=SampleName) %>%    
#     mutate(Sample = case_when(
#         SampleId == "SRR9264351" ~ str_c(Sample, "a"),
#         SampleId == "SRR9264352" ~ str_c(Sample, "b"),
#         TRUE ~ Sample)) %>% 
    column_to_rownames("RowName") %>% 
    select(Sample, Barcode, SampleId, SampleGroup, DatasetName) %>%
    DataFrame()
colData(sce)
## DataFrame with 36447 rows and 5 columns
##                           Sample            Barcode    SampleId SampleGroup
##                      <character>        <character> <character> <character>
## 1_AAACCTGAGAAACCTA-1     ABMMC_1 AAACCTGAGAAACCTA-1   MantonBM1       ABMMC
## 1_AAACCTGAGACAAGCC-1     ABMMC_1 AAACCTGAGACAAGCC-1   MantonBM1       ABMMC
## 1_AAACCTGAGAGACTTA-1     ABMMC_1 AAACCTGAGAGACTTA-1   MantonBM1       ABMMC
## 1_AAACCTGAGCAACGGT-1     ABMMC_1 AAACCTGAGCAACGGT-1   MantonBM1       ABMMC
## 1_AAACCTGAGCAGCGTA-1     ABMMC_1 AAACCTGAGCAGCGTA-1   MantonBM1       ABMMC
## ...                          ...                ...         ...         ...
## 5_TTTGGTTTCACGACTA-1     PRE-T_1 TTTGGTTTCACGACTA-5  SRR9264349       PRE-T
## 5_TTTGTCACAACTGCGC-1     PRE-T_1 TTTGTCACAACTGCGC-5  SRR9264349       PRE-T
## 5_TTTGTCACATGTTGAC-1     PRE-T_1 TTTGTCACATGTTGAC-5  SRR9264349       PRE-T
## 5_TTTGTCAGTCATTAGC-1     PRE-T_1 TTTGTCAGTCATTAGC-5  SRR9264349       PRE-T
## 5_TTTGTCAGTTGTGGCC-1     PRE-T_1 TTTGTCAGTTGTGGCC-5  SRR9264349       PRE-T
##                      DatasetName
##                      <character>
## 1_AAACCTGAGAAACCTA-1         HCA
## 1_AAACCTGAGACAAGCC-1         HCA
## 1_AAACCTGAGAGACTTA-1         HCA
## 1_AAACCTGAGCAACGGT-1         HCA
## 1_AAACCTGAGCAGCGTA-1         HCA
## ...                          ...
## 5_TTTGGTTTCACGACTA-1       Caron
## 5_TTTGTCACAACTGCGC-1       Caron
## 5_TTTGTCACATGTTGAC-1       Caron
## 5_TTTGTCAGTCATTAGC-1       Caron
## 5_TTTGTCAGTTGTGGCC-1       Caron

5.3 Undetected genes

Although the count matrix has 36601 genes, many of these will not have been detected in any droplet.

detected_genes <- rowSums(counts(sce)) > 0
table(detected_genes)
## detected_genes
## FALSE  TRUE 
##  8984 27617

About a quarter of the genes have not been detected in any droplet. We can remove these before proceeding in order to reduce the size of the single cell experiment object.

sce <- sce[detected_genes,]

5.4 Annotate genes

In order to assess the percentage of mitochondrial UMIs, we will need to be able to identify mitochondrial genes. The simplest way to do this is to annotate the genes with their chromosome of origin.

There are many ways we could annotate our genes in R. We will use AnnotationHub. AnnotationHub has access to a large number of annotation databases. Our genes are currently annotated with Ensembl IDs, so we will use Ensembl human database. We will also specify that we want the database corresponding to Ensembl release 98 as this the release from which the CellRanger gene annotation was derived.

ah <- AnnotationHub()
ens.mm.98 <- query(ah, c("Homo sapiens", "EnsDb", 98))[[1]] 

genes <- rowData(sce)$ID
gene_annot <- AnnotationDbi::select(ens.mm.98, 
                                    keys = genes,
                                    keytype = "GENEID",
                                    columns = c("GENEID", "SEQNAME")) %>%
    set_names(c("ID", "Chromosome"))
rowData(sce) <- merge(rowData(sce), gene_annot, by = "ID", sort=FALSE)
rownames(rowData(sce)) <- rowData(sce)$ID

rowData(sce)
## DataFrame with 27617 rows and 4 columns
##                              ID      Symbol            Type  Chromosome
##                     <character> <character>     <character> <character>
## ENSG00000243485 ENSG00000243485 MIR1302-2HG Gene Expression           1
## ENSG00000238009 ENSG00000238009  AL627309.1 Gene Expression           1
## ENSG00000239945 ENSG00000239945  AL627309.3 Gene Expression           1
## ENSG00000241860 ENSG00000241860  AL627309.5 Gene Expression           1
## ENSG00000286448 ENSG00000286448  AP006222.2 Gene Expression           1
## ...                         ...         ...             ...         ...
## ENSG00000275869 ENSG00000275869  AC136612.1 Gene Expression  KI270728.1
## ENSG00000277761 ENSG00000277761  AC136616.2 Gene Expression  KI270728.1
## ENSG00000277836 ENSG00000277836  AC141272.1 Gene Expression  KI270728.1
## ENSG00000278633 ENSG00000278633  AC023491.2 Gene Expression  KI270731.1
## ENSG00000278817 ENSG00000278817  AC007325.4 Gene Expression  KI270734.1

5.5 Add per cell QC metrics

We can now add the per cell QC metrics to the droplet annotation using the function addPerCellQC. In order to get the metrics for the subset of mitochondrial genes, we need to pass the function a vector indicating which genes are mitochondrial.

is.mito <- which(rowData(sce)$Chromosome=="MT")

sce <- addPerCellQC(sce, subsets=list(Mito=is.mito), BPPARAM = bp.params)

The function has added six columns to the droplet annotation:

  • sum: total UMI count
  • detected: number of features (genes) detected
  • subsets_Mito_sum: number of UMIs mapped to mitochondrial transcripts
  • subsets_Mito_detected: number of mitochondrial genes detected
  • subsets_Mito_percent: percentage of UMIs mapped to mitochondrial transcripts
  • total: also the total UMI count

We will use sum, detected, and subsets_Mito_percent to further filter the cells.

colData(sce)
## DataFrame with 36447 rows and 11 columns
##                           Sample            Barcode    SampleId SampleGroup
##                      <character>        <character> <character> <character>
## 1_AAACCTGAGAAACCTA-1     ABMMC_1 AAACCTGAGAAACCTA-1   MantonBM1       ABMMC
## 1_AAACCTGAGACAAGCC-1     ABMMC_1 AAACCTGAGACAAGCC-1   MantonBM1       ABMMC
## 1_AAACCTGAGAGACTTA-1     ABMMC_1 AAACCTGAGAGACTTA-1   MantonBM1       ABMMC
## 1_AAACCTGAGCAACGGT-1     ABMMC_1 AAACCTGAGCAACGGT-1   MantonBM1       ABMMC
## 1_AAACCTGAGCAGCGTA-1     ABMMC_1 AAACCTGAGCAGCGTA-1   MantonBM1       ABMMC
## ...                          ...                ...         ...         ...
## 5_TTTGGTTTCACGACTA-1     PRE-T_1 TTTGGTTTCACGACTA-5  SRR9264349       PRE-T
## 5_TTTGTCACAACTGCGC-1     PRE-T_1 TTTGTCACAACTGCGC-5  SRR9264349       PRE-T
## 5_TTTGTCACATGTTGAC-1     PRE-T_1 TTTGTCACATGTTGAC-5  SRR9264349       PRE-T
## 5_TTTGTCAGTCATTAGC-1     PRE-T_1 TTTGTCAGTCATTAGC-5  SRR9264349       PRE-T
## 5_TTTGTCAGTTGTGGCC-1     PRE-T_1 TTTGTCAGTTGTGGCC-5  SRR9264349       PRE-T
##                      DatasetName       sum  detected subsets_Mito_sum
##                      <character> <numeric> <integer>        <numeric>
## 1_AAACCTGAGAAACCTA-1         HCA      4183       969              138
## 1_AAACCTGAGACAAGCC-1         HCA      2851       763               68
## 1_AAACCTGAGAGACTTA-1         HCA     27790      4681             1243
## 1_AAACCTGAGCAACGGT-1         HCA     20543      4121              829
## 1_AAACCTGAGCAGCGTA-1         HCA      2521       921              126
## ...                          ...       ...       ...              ...
## 5_TTTGGTTTCACGACTA-1       Caron     17960      4143              638
## 5_TTTGTCACAACTGCGC-1       Caron      1510       837              283
## 5_TTTGTCACATGTTGAC-1       Caron     13076      3041              385
## 5_TTTGTCAGTCATTAGC-1       Caron      5784      1961               90
## 5_TTTGTCAGTTGTGGCC-1       Caron      6675      2112              158
##                      subsets_Mito_detected subsets_Mito_percent     total
##                                  <integer>            <numeric> <numeric>
## 1_AAACCTGAGAAACCTA-1                    11              3.29907      4183
## 1_AAACCTGAGACAAGCC-1                    10              2.38513      2851
## 1_AAACCTGAGAGACTTA-1                    13              4.47283     27790
## 1_AAACCTGAGCAACGGT-1                    12              4.03544     20543
## 1_AAACCTGAGCAGCGTA-1                    11              4.99802      2521
## ...                                    ...                  ...       ...
## 5_TTTGGTTTCACGACTA-1                    13              3.55234     17960
## 5_TTTGTCACAACTGCGC-1                    12             18.74172      1510
## 5_TTTGTCACATGTTGAC-1                    13              2.94433     13076
## 5_TTTGTCAGTCATTAGC-1                    11              1.55602      5784
## 5_TTTGTCAGTTGTGGCC-1                    11              2.36704      6675

5.6 QC metric distribution

Before moving on to do the actual cell filtering, it is always a good idea to explore the distribution of the metrics across the droplets.

We can use the scater function plotColData to generate plots that provide a look at these distributions on a per sample basis.

plotColData(sce, x="Sample", y="sum",other_fields="SampleGroup") + 
    facet_wrap(~SampleGroup, nrow=1, scales = "free_x") + 
    scale_y_log10() + 
    ggtitle("Total count")

plotColData(sce, x="Sample", y="detected", other_fields="SampleGroup") + 
    facet_wrap(~SampleGroup, nrow=1, scales = "free_x") + 
    scale_y_log10() + 
    ggtitle("Detected features")

plotColData(sce, x="Sample", y="subsets_Mito_percent", other_fields="SampleGroup") + 
    facet_wrap(~SampleGroup, nrow=1, scales = "free_x") +
    ggtitle("Mito percent")

A scatter plot shows the extent to which library size and numbers of genes detected are correlated.

colData(sce) %>% 
    as.data.frame() %>% 
    arrange(subsets_Mito_percent) %>% 
    ggplot(aes(x = sum, y = detected)) +
      geom_point(aes(colour = subsets_Mito_percent > 10)) + 
      facet_wrap(vars(SampleGroup))

5.7 Identification of low-quality cells with adaptive thresholds

One could use hard threshold for the library size, number of genes detected and mitochondrial content based on the distributions seen above. These would need vary across runs and the decision making process is somewhat arbitrary. It may therefore be preferable to rely on outlier detection to identify cells that markedly differ from most cells.

We saw above that the distribution of the QC metrics is close to Normal. Hence, we can detect outliers using the median and the median absolute deviation (MAD) from the median (not the mean and the standard deviation which both are sensitive to outliers).

For a given metric, an outlier value is one that lies over some number of MADs away from the median. A cell will be excluded if it is an outlier in the part of the range to avoid, for example low gene counts, or high mitochondrial content. For a normal distribution, a threshold defined with a distance of 3 MADs from the median retains about 99% of values.

The scater function isOutlier can be used to detect outlier cells based on any metric in the colData table. It returns a boolean vector that identifies outliers. By default it will mark any cell that is 3 MADS in either direction from the median as an outlier.

5.7.1 Library size

With library size we wish to identify outliers that have very low library sizes, this indicates that the droplets either contain poor quality cells, perhaps damaged or dying, or do not contain a cell at all.

The library size distribution tends to have a long tail to the right (small numbers of cells with very high UMI counts). We therefore log transform the library size in order to the make the distribution closer to normal. This also improves the resolution of the smaller library sizes and ensures that we do not end up with negative threshold.

low_lib_size <- isOutlier(sce$sum, log=TRUE, type="lower")
table(low_lib_size)
## low_lib_size
## FALSE  TRUE 
## 35544   903

This has excluded 903 cells. We can view the threshold values to check that they seem reasonable.

attr(low_lib_size, "thresholds")
##    lower   higher 
## 781.7415      Inf

We can view the effect of the filtering using plotColData.

colData(sce)$low_lib_size <- low_lib_size
plotColData(sce, 
            x="Sample", 
            y="sum",
            other_fields="SampleGroup", 
            colour_by = "low_lib_size") + 
    facet_wrap(~SampleGroup, nrow=1, scales = "free_x") + 
    scale_y_log10() + 
    labs(y = "Total count", title = "Total count") +
    guides(colour=guide_legend(title="Discarded"))

5.7.2 Number of genes

As with the library size, we will log tranform the number of genes detected prior to filtering using the median absolute deviation.

low_n_features <- isOutlier(sce$detected, log=TRUE, type="lower")
table(low_n_features)
## low_n_features
## FALSE  TRUE 
## 36233   214

This has excluded out 214 cells. The threshold value was:

attr(low_n_features, "thresholds")[1]
##    lower 
## 286.9177

We can view the effect of the filtering using plotColData.

colData(sce)$low_n_features <- low_n_features
plotColData(sce, 
            x="Sample", 
            y="detected",
            other_fields="SampleGroup", 
            colour_by = "low_n_features") + 
    facet_wrap(~SampleGroup, nrow=1, scales = "free_x") + 
    scale_y_log10() + 
    labs(y = "Genes detected", title = "Genes detected") +
    guides(colour=guide_legend(title="Discarded"))

5.7.3 Mitochondrial content

For the mitochondrial content the exclusion zone is in the higher part of the distribution. For this reason we do not need to worry about log transforming the data as want to remove the long right hand tail anyway.

high_Mito_percent <- isOutlier(sce$subsets_Mito_percent, type="higher")
table(high_Mito_percent)
## high_Mito_percent
## FALSE  TRUE 
## 32863  3584

This has removed 3584 cells in total. The upper threshold value:

attr(high_Mito_percent, "thresholds")[2]
##   higher 
## 9.308917

We can view the effect of the filtering using plotColData.

colData(sce)$high_Mito_percent <- high_Mito_percent
plotColData(sce,  
            x="Sample",
            y="subsets_Mito_percent",
            other_fields="SampleGroup",
            colour_by = "high_Mito_percent") + 
    facet_wrap(~SampleGroup, nrow=1, scales = "free_x") + 
    labs(y = "Percentage mitochondrial UMIs",
         title = "Mitochondrial UMIs") +
    guides(colour=guide_legend(title="Discarded"))

5.7.4 Summary of discarded cells

Having applied each of the three thresholds separately, we can now combine them to see how many droplets in total we will be excluding.

data.frame(`Library Size` = sum(low_lib_size),
           `Genes detected` = sum(low_n_features),
           `Mitochondrial UMIs` = sum(high_Mito_percent),
           Total = sum(low_lib_size | low_n_features | high_Mito_percent))

5.7.5 All three filter steps at once

The three steps above may be run in one go using the quickPerCellQC function. This creates a DataFrame with 4 columns containing TRUE/FALSE - one for each filter metric and one called “discard” that combined the three logicals.

cell_qc_results <- quickPerCellQC(colData(sce),
              percent_subsets=c("subsets_Mito_percent"))
colSums(as.data.frame(cell_qc_results))
##              low_lib_size            low_n_features high_subsets_Mito_percent 
##                       903                       214                      3584 
##                   discard 
##                      3808

5.7.6 Assumptions

Data quality depends on the tissue analysed, some being difficult to dissociate, e.g. brain, so that one level of QC stringency will not fit all data sets.

Filtering based on QC metrics as done here assumes that these QC metrics are not correlated with biology. This may not necessarily be true in highly heterogenous data sets where some cell types represented by good-quality cells may have low RNA content or high mitochondrial content.

5.7.7 Considering experimental factors when filtering

The HCA and Caron data sets analysed here may have been obtained in experiments with different conditions, such as cell preparation or sequencing depth. Such differences between these two batches would affect the adaptive thesholds discussed above - that is, the distributions of the metrics may be different in each batch and so perhaps we should really apply the adaptive thresholding for each batch. The quickPerCellQC has a “batch” argument that allows us to specify with samples belong to which batches. The batches are then filtered independently.

batch.cell_qc_results <- quickPerCellQC(colData(sce),
                                percent_subsets=c("subsets_Mito_percent"),
                                batch=sce$DatasetName)
colSums(as.data.frame(batch.cell_qc_results))
##              low_lib_size            low_n_features high_subsets_Mito_percent 
##                         0                       143                      3001 
##                   discard 
##                      3103

The table below shows how the thresholds for each metric differ between the batch-wise analysis and the analysis using all samples.

all.thresholds <- tibble(`Batch`="All",
       `Library Size`=attr(cell_qc_results$low_lib_size, "thresholds")[1],
       `Genes detected`=attr(cell_qc_results$low_n_features, "thresholds")[1],
       `Mitochondrial UMIs`=attr(cell_qc_results$high_subsets_Mito_percent, "thresholds")[2])


tibble(`Batch`=names(attr(batch.cell_qc_results$low_lib_size, "thresholds")[1,]),
       `Library Size`=attr(batch.cell_qc_results$low_lib_size, "thresholds")[1,],
       `Genes detected`=attr(batch.cell_qc_results$low_n_features, "thresholds")[1,],
       `Mitochondrial UMIs`=attr(batch.cell_qc_results$high_subsets_Mito_percent, "thresholds")[2,]) %>% 
    bind_rows(all.thresholds) %>% 
    mutate(across(where(is.numeric), round, digits=2)) %>% 
    datatable(rownames = FALSE, options = list(dom="t"))

Let’s replace the columns in the droplet annotation with these new filters.

sce$low_lib_size <- batch.cell_qc_results$low_lib_size
sce$low_n_features <- batch.cell_qc_results$low_n_features
sce$high_Mito_percent <- batch.cell_qc_results$high_subsets_Mito_percent
sce$discard <- batch.cell_qc_results$discard

We can visualise how the new filters look using violin plots.

plotColData(sce, 
            x="Sample", 
            y="sum",
            other_fields="SampleGroup", 
            colour_by = "low_lib_size") + 
    facet_wrap(vars(SampleGroup), nrow=1, scales = "free_x") + 
    scale_y_log10() + 
    labs(y = "Total count", title = "Total count") +
    guides(colour=guide_legend(title="Discarded"))

plotColData(sce, 
            x="Sample", 
            y="detected",
            other_fields="SampleGroup", 
            colour_by = "low_n_features") + 
    facet_wrap(vars(SampleGroup), nrow=1, scales = "free_x") + 
    scale_y_log10() + 
    labs(y = "Genes detected", title = "Genes detected") +
    guides(colour=guide_legend(title="Discarded"))

plotColData(sce, 
        x="Sample", 
        y="subsets_Mito_percent",
        other_fields="SampleGroup", 
        colour_by = "high_Mito_percent") + 
    facet_wrap(vars(SampleGroup), nrow=1, scales = "free_x") + 
    labs(y = "Percentage mitochondrial UMIs",
         title = "Mitochondrial UMIs") +
    guides(colour=guide_legend(title="Discarded"))

There are some distinct differences, most noticeable is that there is now no filtering using library size. The venn diagrams below show how the number of discarded droplets in HCA and Caron have changed for each filter in comparison to when the MAD filtering was applied across all samples.

libDat <- tibble(`All together`=cell_qc_results$low_lib_size, 
                 `By batch`=batch.cell_qc_results$low_lib_size,
                 Batch=sce$Sample)
    
ph1 <- libDat %>% 
    dplyr::filter(Batch=="HCA") %>% 
    ggvenn(show_percentage = FALSE) +
        labs(title="Library Size - HCA")
pc1 <- libDat %>% 
    dplyr::filter(Batch=="Caron") %>% 
    ggvenn(show_percentage = FALSE) +
        labs(title="Library Size - Caron")

nGenDat <- tibble(`All together`=cell_qc_results$low_n_features, 
                  `By batch`=batch.cell_qc_results$low_n_features,
                 Batch=sce$DatasetName)
ph2 <- nGenDat %>% 
    dplyr::filter(Batch=="HCA") %>% 
        ggvenn(show_percentage = FALSE) +
            labs(title="Genes detected - HCA")
pc2 <- nGenDat %>% 
    dplyr::filter(Batch=="Caron") %>% 
           ggvenn(show_percentage = FALSE) +
            labs(title="Genes detected - Caron")


mitDat <- tibble(`All together`=cell_qc_results$high_subsets_Mito_percent, 
       `By batch`=batch.cell_qc_results$high_subsets_Mito_percent,
                 Batch=sce$DatasetName)
ph3 <- mitDat %>% 
    dplyr::filter(Batch=="HCA") %>% 
        ggvenn(show_percentage = FALSE) +
            labs(title="Mitochondrial UMIs - HCA")
pc3 <- mitDat %>% 
    dplyr::filter(Batch=="Caron") %>% 
           ggvenn(show_percentage = FALSE) +
            labs(title="Mitochondrial UMIs - Caron")

(pc1 + pc2 + pc3) / (ph1 + ph2 + ph3)

The most striking difference is in the filtering of the Caron data by library size. As we can see from the violin plots the ABMMC sample from the HCA has a radically different library size distribution to the Caron samples, with all cells having > 1000 UMIs. When we applied the adaptive filters across all samples, this sample caused the MADs to be distorted and resulted in a threshold that was inappropriately high for the Caron samples.

5.8 Filtering out poor quality cells

Now that we have identified poor quality cells we can filter them out before proceeding to do any further analysis.

sce.filtered <- sce[, !sce$discard]

Note: An important thing to note is that, now that we have filtered this object, some of the QC metrics that we calculated across all genes (for colData) and across all cells (for rowData) are no longer correct for the filtered data set. We need to remove them, and if necessary recalculate them.

colData(sce.filtered) <- colData(sce.filtered)[,1:3]
sce.filtered <- addPerCellQC(sce.filtered, BPPARAM = bp.params)

6 QC and filtering by combining the metrics.

In the previous approach we used the three metrics in isolation to filter droplets. Another approach is to combine the three (or more) metrics in a single filtering step by looking for outliers in the multi-dimensional space defined by the metrics.

As with the adaptive thresholds above, this method should not be applied across batches or samples with differing distributions in the metrics or it will exclude many good quality cells. To demonstrate these methods, we’ll just extract one sample from our SingleCellExperiment object.

sce.BM1 <- sce[ , sce$Sample == "PBMMC_1"]

6.1 Using “outlyingness”

Essentially we need to reduce our 3 metrics to a single metric, we can then use isOutlier to select outliers based on this metric. One way to do this is to use the function adjOutlyingness from the robustbase package. This function computes the “outlyingness” for each droplet.

Here we will use the same three metrics as before: library size, the number of genes detected and the mitochondrial content. Remember that for “sum” (total UMIs) and “detected” (number of genes detected), we want to use the log10 value.

library(robustbase)
stats <- cbind(log10(sce.BM1$sum),
               log10(sce.BM1$detected),
               sce.BM1$subsets_Mito_percent)

outlying <- adjOutlyingness(stats, only.outlyingness = TRUE)
multi.outlier <- isOutlier(outlying, type = "higher")
summary(multi.outlier)
##    Mode   FALSE    TRUE 
## logical     898      82

6.2 Using PCA

Another approach is to perform a principal component analysis (PCA) on the table of metrics, apply adjOutlyingness to the metrics table and use this to detect outliers. The scater function runColDataPCA can be used to perform the PCA and detect outliers. We’ll need to add a couple of columns to the colData for the log10 metrics first.

sce.BM1$log10sum <- log10(sce.BM1$sum)
sce.BM1$log10detected <- log10(sce.BM1$detected)
sce.BM1 <- runColDataPCA(sce.BM1, 
                     variables=list("log10sum", "log10detected", "subsets_Mito_percent"),
                     outliers=TRUE,
                     BPPARAM = bp.params)

This has added the results of the principal component analysis into a new slot in the SingleCellExperiment object specifically for holding the results of dimension reduction transformations such as PCA, t-SNE and UMAP. The results can be accessed using the reducedDim function.

head(reducedDim(sce.BM1))
##                             PC1         PC2
## 4_AAACCTGCACTTCGAA-1 -0.4059768  0.54959341
## 4_AAACCTGCAGACGCAA-1 -0.3242857  0.52845630
## 4_AAACCTGTCATCACCC-1 -0.3979909  0.09660094
## 4_AAAGATGCATAAAGGT-1  0.2074401 -0.07273329
## 4_AAAGATGTCTCGAGTA-1 -2.3969937  1.34212551
## 4_AAAGCAAAGCAGCGTA-1 -0.1425594  0.27437045

It has also added a column “outlier” to the colData, which specifies the droplets that have been identified as outliers.

summary(sce.BM1$outlier)
##    Mode   FALSE    TRUE 
## logical     917      63

6.3 A note on multi-dimensional filtering

These types of approach can provide more power for detecting outliers as they are looking at patterns across multiple metrics, however, it can be difficult to interpret the reason why any particular droplet has been excluded.

7 Mitochondrial content versus library size

A useful diagnostic plot for assessing the impact of the filtering is to do a scatter plot of the mitochondrial content against the library size. We can overlay our final filter metric using the point colour.

plotColData(sce, 
            x="sum", 
            y="subsets_Mito_percent", 
            other_fields="Sample",
            colour_by="discard") +
    facet_wrap(~Sample, ncol=5, scale="free_x")

8 QC and Filtering based on sparsity

The approach above identified poor-quality using thresholds on the number of genes detected and mitochondrial content. We will here specifically look at the sparsity of the data, both at the gene and cell levels.

8.1 Sparsity plots

We will compute:

  • the cell sparsity: for each cell, the proportion of genes that are not detected
  • the gene sparsity: for each gene, the proportion of cells in which it is not detected

To help calculate the gene sparsity we can generate QC metrics for genes with addPerFeatureQC. This adds two columns to the gene annotation (rowData):

  • mean - the mean UMI count for the gene across all cells
  • detected - the percentage of cells in which the gene was detected
sce <- addPerFeatureQC(sce, BPPARAM = bp.params)
rowData(sce)
## DataFrame with 27617 rows and 6 columns
##                              ID      Symbol            Type  Chromosome
##                     <character> <character>     <character> <character>
## ENSG00000243485 ENSG00000243485 MIR1302-2HG Gene Expression           1
## ENSG00000238009 ENSG00000238009  AL627309.1 Gene Expression           1
## ENSG00000239945 ENSG00000239945  AL627309.3 Gene Expression           1
## ENSG00000241860 ENSG00000241860  AL627309.5 Gene Expression           1
## ENSG00000286448 ENSG00000286448  AP006222.2 Gene Expression           1
## ...                         ...         ...             ...         ...
## ENSG00000275869 ENSG00000275869  AC136612.1 Gene Expression  KI270728.1
## ENSG00000277761 ENSG00000277761  AC136616.2 Gene Expression  KI270728.1
## ENSG00000277836 ENSG00000277836  AC141272.1 Gene Expression  KI270728.1
## ENSG00000278633 ENSG00000278633  AC023491.2 Gene Expression  KI270731.1
## ENSG00000278817 ENSG00000278817  AC007325.4 Gene Expression  KI270734.1
##                        mean   detected
##                   <numeric>  <numeric>
## ENSG00000243485 2.74371e-05 0.00274371
## ENSG00000238009 7.68239e-04 0.07682388
## ENSG00000239945 2.74371e-05 0.00274371
## ENSG00000241860 2.79858e-03 0.27711471
## ENSG00000286448 1.92060e-04 0.01920597
## ...                     ...        ...
## ENSG00000275869 2.74371e-05 0.00274371
## ENSG00000277761 2.74371e-05 0.00274371
## ENSG00000277836 8.23113e-05 0.00823113
## ENSG00000278633 5.48742e-05 0.00548742
## ENSG00000278817 3.31989e-03 0.32924521

Now we can calculate sparsity using the “detected” columns in the colData and the rowData.

colData(sce)$cell_sparsity <- 1 - (colData(sce)$detected / nrow(sce))
rowData(sce)$gene_sparsity <- (100 - rowData(sce)$detected) / 100

We now plot the distribution of these two metrics.

The cell sparsity plot shows that most cells have between 85% and 99% 0’s, which is typical.

hist(sce$cell_sparsity, breaks=50, col="grey80", xlab="Cell sparsity", main="")

The gene sparsity plot shows that a large number of genes are almost never detected, which is also regularly observed.

hist(rowData(sce)$gene_sparsity, breaks=50, col="grey80", xlab="Gene sparsity", main="")

8.2 Filter by sparsity

We could remove cells with sparsity higher than 0.99, and/or mitochondrial content higher than 10%.

Genes detected in a few cells only are unlikely to be informative and would hinder normalisation. We will remove genes that are expressed in fewer than 20 cells.

sparse.cells <- sce$cell_sparsity > 0.99
mito.cells <- sce$subsets_Mito_percent > 10

min.cells <- 1 - (20 / ncol(sce))
sparse.genes <- rowData(sce)$gene_sparsity > min.cells

Number of genes removed:

table(sparse.genes)
## sparse.genes
## FALSE  TRUE 
## 19317  8300

Number of cells removed:

table(sparse.cells, mito.cells)
##             mito.cells
## sparse.cells FALSE  TRUE
##        FALSE 33108  3151
##        TRUE     70   118

8.3 Session information

sessionInfo()
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] robustbase_0.93-8           ggvenn_0.1.9               
##  [3] patchwork_1.1.1             forcats_0.5.1              
##  [5] stringr_1.4.0               dplyr_1.0.7                
##  [7] purrr_0.3.4                 readr_2.0.1                
##  [9] tidyr_1.1.3                 tibble_3.1.4               
## [11] tidyverse_1.3.1             BiocParallel_1.26.2        
## [13] AnnotationHub_3.0.1         BiocFileCache_2.0.0        
## [15] dbplyr_2.1.1                ensembldb_2.16.4           
## [17] AnnotationFilter_1.16.0     GenomicFeatures_1.44.2     
## [19] AnnotationDbi_1.54.1        scater_1.20.1              
## [21] ggplot2_3.3.5               scuttle_1.2.1              
## [23] DropletUtils_1.12.2         SingleCellExperiment_1.14.1
## [25] SummarizedExperiment_1.22.0 Biobase_2.52.0             
## [27] GenomicRanges_1.44.0        GenomeInfoDb_1.28.2        
## [29] IRanges_2.26.0              S4Vectors_0.30.0           
## [31] BiocGenerics_0.38.0         MatrixGenerics_1.4.3       
## [33] matrixStats_0.60.1          knitr_1.33                 
## [35] DT_0.18                     rmarkdown_2.10             
## [37] nvimcom_0.9-82             
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1                  backports_1.2.1              
##   [3] lazyeval_0.2.2                crosstalk_1.1.1              
##   [5] digest_0.6.27                 htmltools_0.5.2              
##   [7] viridis_0.6.1                 fansi_0.5.0                  
##   [9] magrittr_2.0.1                memoise_2.0.0                
##  [11] ScaledMatrix_1.0.0            tzdb_0.1.2                   
##  [13] limma_3.48.3                  Biostrings_2.60.2            
##  [15] modelr_0.1.8                  vroom_1.5.4                  
##  [17] R.utils_2.10.1                prettyunits_1.1.1            
##  [19] colorspace_2.0-2              rvest_1.0.1                  
##  [21] blob_1.2.2                    rappdirs_0.3.3               
##  [23] haven_2.4.3                   xfun_0.25                    
##  [25] crayon_1.4.1                  RCurl_1.98-1.4               
##  [27] jsonlite_1.7.2                glue_1.4.2                   
##  [29] gtable_0.3.0                  zlibbioc_1.38.0              
##  [31] XVector_0.32.0                DelayedArray_0.18.0          
##  [33] BiocSingular_1.8.1            Rhdf5lib_1.14.2              
##  [35] DEoptimR_1.0-9                HDF5Array_1.20.0             
##  [37] scales_1.1.1                  DBI_1.1.1                    
##  [39] edgeR_3.34.0                  Rcpp_1.0.7                   
##  [41] viridisLite_0.4.0             xtable_1.8-4                 
##  [43] progress_1.2.2                dqrng_0.3.0                  
##  [45] bit_4.0.4                     rsvd_1.0.5                   
##  [47] htmlwidgets_1.5.3             httr_1.4.2                   
##  [49] ellipsis_0.3.2                farver_2.1.0                 
##  [51] pkgconfig_2.0.3               XML_3.99-0.7                 
##  [53] R.methodsS3_1.8.1             sass_0.4.0                   
##  [55] locfit_1.5-9.4                utf8_1.2.2                   
##  [57] labeling_0.4.2                tidyselect_1.1.1             
##  [59] rlang_0.4.11                  later_1.3.0                  
##  [61] cellranger_1.1.0              munsell_0.5.0                
##  [63] BiocVersion_3.13.1            tools_4.1.1                  
##  [65] cachem_1.0.6                  cli_3.0.1                    
##  [67] generics_0.1.0                RSQLite_2.2.8                
##  [69] broom_0.7.9                   evaluate_0.14                
##  [71] fastmap_1.1.0                 yaml_2.2.1                   
##  [73] fs_1.5.0                      bit64_4.0.5                  
##  [75] KEGGREST_1.32.0               sparseMatrixStats_1.4.2      
##  [77] mime_0.11                     R.oo_1.24.0                  
##  [79] xml2_1.3.2                    biomaRt_2.48.3               
##  [81] rstudioapi_0.13               compiler_4.1.1               
##  [83] beeswarm_0.4.0                filelock_1.0.2               
##  [85] curl_4.3.2                    png_0.1-7                    
##  [87] interactiveDisplayBase_1.30.0 reprex_2.0.1                 
##  [89] bslib_0.2.5.1                 stringi_1.7.4                
##  [91] highr_0.9                     lattice_0.20-45              
##  [93] ProtGenerics_1.24.0           Matrix_1.3-4                 
##  [95] vctrs_0.3.8                   pillar_1.6.2                 
##  [97] lifecycle_1.0.0               rhdf5filters_1.4.0           
##  [99] BiocManager_1.30.16           jquerylib_0.1.4              
## [101] BiocNeighbors_1.10.0          cowplot_1.1.1                
## [103] bitops_1.0-7                  irlba_2.3.3                  
## [105] httpuv_1.6.2                  rtracklayer_1.52.1           
## [107] R6_2.5.1                      BiocIO_1.2.0                 
## [109] promises_1.2.0.1              gridExtra_2.3                
## [111] vipor_0.4.5                   codetools_0.2-18             
## [113] assertthat_0.2.1              rhdf5_2.36.0                 
## [115] rjson_0.2.20                  withr_2.4.2                  
## [117] GenomicAlignments_1.28.0      Rsamtools_2.8.0              
## [119] GenomeInfoDbData_1.2.6        hms_1.1.0                    
## [121] beachmat_2.8.1                DelayedMatrixStats_1.14.3    
## [123] Cairo_1.5-12.2                lubridate_1.7.10             
## [125] shiny_1.6.0                   ggbeeswarm_0.6.0             
## [127] restfulr_0.0.13