Created by: Ahmed Mahfouz
In this practical, we will walk through a pipeline to analyze single cell RNA-sequencing (scRNA-seq) data. Starting from a count matrix, we will cover the following steps of the analysis: 1. Quality control 2. Normalization 3. Feature selection
For this tutorial we will use 3 different PBMC datasets from the 10x Genomics website (https://support.10xgenomics.com/single-cell-gene-expression/datasets).
- 1k PBMCs using 10x v2 chemistry
- 1k PBMCs using 10x v3 chemistry
- 1k PBMCs using 10x v3 chemistry in combination with cell surface proteins, but disregarding the protein data and only looking at gene expression.
The datasets are available in this repository.
Load required packages:
suppressMessages(require(Seurat))
suppressMessages(require(scater))
suppressMessages(require(scran))
suppressMessages(require(Matrix))Here, we use the function Read10X_h5 to read in the expression matrices in R.
v3.1k <- Read10X_h5("pbmc_1k_v3_filtered_feature_bc_matrix.h5", use.names = T)
v2.1k <- Read10X_h5("pbmc_1k_v2_filtered_feature_bc_matrix.h5", use.names = T)
p3.1k <- Read10X_h5("pbmc_1k_protein_v3_filtered_feature_bc_matrix.h5", use.names = T)## Genome matrix has multiple modalities, returning a list of matrices for this genome
# select only gene expression data from the CITE-seq data.
p3.1k <- p3.1k$`Gene Expression`First, create Seurat objects for each of the datasets, and then merge into one large seurat object.
sdata.v2.1k <- CreateSeuratObject(v2.1k, project = "v2.1k")
sdata.v3.1k <- CreateSeuratObject(v3.1k, project = "v3.1k")
sdata.p3.1k <- CreateSeuratObject(p3.1k, project = "p3.1k")
# merge into one single seurat object. Add cell ids just in case you have overlapping barcodes between the datasets.
alldata <- merge(sdata.v2.1k, c(sdata.v3.1k,sdata.p3.1k), add.cell.ids=c("v2.1k","v3.1k","p3.1k"))
# also add in a metadata column that indicates v2 vs v3 chemistry
chemistry <- rep("v3",ncol(alldata))
chemistry[Idents(alldata) == "v2.1k"] <- "v2"
alldata <- AddMetaData(alldata, chemistry, col.name = "Chemistry")
alldata## An object of class Seurat
## 33538 features across 2931 samples within 1 assay
## Active assay: RNA (33538 features, 0 variable features)
Check number of cells from each sample, is stored in the orig.ident slot of metadata and is autmatically set as active ident.
table(Idents(alldata))##
## p3.1k v2.1k v3.1k
## 713 996 1222
Seurat automatically calculates some QC-stats, like number of UMIs and features per cell. Stored in columns nCount_RNA & nFeature_RNA of the metadata.
head(alldata@meta.data)## orig.ident nCount_RNA nFeature_RNA Chemistry
## v2.1k_AAACCTGAGCGCTCCA-1 v2.1k 6631 2029 v2
## v2.1k_AAACCTGGTGATAAAC-1 v2.1k 2196 881 v2
## v2.1k_AAACGGGGTTTGTGTG-1 v2.1k 2700 791 v2
## v2.1k_AAAGATGAGTACTTGC-1 v2.1k 3551 1183 v2
## v2.1k_AAAGCAAGTCTCTTAT-1 v2.1k 3080 1333 v2
## v2.1k_AAAGCAATCCACGAAT-1 v2.1k 5769 1556 v2
We will manually calculate the proportion of mitochondrial reads and add to the metadata table.
percent.mito <- PercentageFeatureSet(alldata, pattern = "^MT-")
alldata <- AddMetaData(alldata, percent.mito, col.name = "percent.mito")In the same manner we will calculate the proportion gene expression that comes from ribosomal proteins.
percent.ribo <- PercentageFeatureSet(alldata, pattern = "^RP[SL]")
alldata <- AddMetaData(alldata, percent.ribo, col.name = "percent.ribo")Now have another look at the metadata table
head(alldata@meta.data)## orig.ident nCount_RNA nFeature_RNA Chemistry
## v2.1k_AAACCTGAGCGCTCCA-1 v2.1k 6631 2029 v2
## v2.1k_AAACCTGGTGATAAAC-1 v2.1k 2196 881 v2
## v2.1k_AAACGGGGTTTGTGTG-1 v2.1k 2700 791 v2
## v2.1k_AAAGATGAGTACTTGC-1 v2.1k 3551 1183 v2
## v2.1k_AAAGCAAGTCTCTTAT-1 v2.1k 3080 1333 v2
## v2.1k_AAAGCAATCCACGAAT-1 v2.1k 5769 1556 v2
## percent.mito percent.ribo
## v2.1k_AAACCTGAGCGCTCCA-1 5.172674 25.84829
## v2.1k_AAACCTGGTGATAAAC-1 4.143898 20.81056
## v2.1k_AAACGGGGTTTGTGTG-1 3.296296 51.55556
## v2.1k_AAAGATGAGTACTTGC-1 5.885666 29.25936
## v2.1k_AAAGCAAGTCTCTTAT-1 2.987013 17.53247
## v2.1k_AAAGCAATCCACGAAT-1 2.010747 45.69249
Now we can plot some of the QC-features as violin plots
VlnPlot(alldata, features = "nFeature_RNA", pt.size = 0.1) + NoLegend()
VlnPlot(alldata, features = "nCount_RNA", pt.size = 0.1) + NoLegend()
VlnPlot(alldata, features = "percent.mito", pt.size = 0.1) + NoLegend()
VlnPlot(alldata, features = "percent.ribo", pt.size = 0.1) + NoLegend()
As you can see, the v2 chemistry gives lower gene detection, but higher detection of ribosomal proteins. As the ribosomal proteins are highly expressed they will make up a larger proportion of the transcriptional landscape when fewer of the lowly expressed genes are detected.
We can also plot the different QC-measures as scatter plots.
FeatureScatter(alldata, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
FeatureScatter(alldata, feature1 = "nFeature_RNA", feature2 = "percent.mito")
FeatureScatter(alldata, feature1="percent.ribo", feature2="nFeature_RNA")
We can also subset the data to only plot one sample.
FeatureScatter(alldata, feature1 = "nCount_RNA", feature2 = "nFeature_RNA",
cells = WhichCells(alldata, expression = orig.ident == "v3.1k") )
We have quite a lot of cells with high proportion of mitochondrial reads. It could be wise to remove those cells, if we have enough cells left after filtering. Another option would be to either remove all mitochondrial reads from the dataset and hope that the remaining genes still have enough biological signal. A third option would be to just regress out the percent.mito variable during scaling.
In this case we have as much as 99.7% mitochondrial reads in some of the cells, so it is quite unlikely that there is much celltype signature left in those.
Looking at the plots, make resonable decisions on where to draw the cutoff. In this case, the bulk of the cells are below 25% mitochondrial reads and that will be used as a cutoff.
#select cells with percent.mito < 25
idx <- which(alldata$percent.mito < 25)
selected <- WhichCells(alldata, cells = idx)
length(selected)## [1] 2703
# and subset the object to only keep those cells
data.filt <- subset(alldata, cells = selected)
# plot violins for new data
VlnPlot(data.filt, features = "percent.mito")
As you can see, there is still quite a lot of variation in percent mito, so it will have to be dealt with in the data analysis step.
Extremely high number of detected genes could indicate doublets. However, depending on the celltype composition in your sample, you may have cells with higher number of genes (and also higher counts) from one celltype.
In these datasets, there is also a clear difference between the v2 vs v3 10x chemistry with regards to gene detection, so it may not be fair to apply the same cutoffs to all of them.
Also, in the protein assay data there is a lot of cells with few detected genes giving a bimodal distribution. This type of distribution is not seen in the other 2 datasets. Considering that they are all pbmc datasets it makes sense to regard this distribution as low quality libraries.
Filter the cells with high gene detection (putative doublets) with cutoffs 4100 for v3 chemistry and 2000 for v2.
#start with cells with many genes detected.
high.det.v3 <- WhichCells(data.filt, expression = nFeature_RNA > 4100)
high.det.v2 <- WhichCells(data.filt, expression = nFeature_RNA > 2000 & orig.ident == "v2.1k")
# remove these cells
data.filt <- subset(data.filt, cells=setdiff(WhichCells(data.filt),c(high.det.v2,high.det.v3)))
# check number of cells
ncol(data.filt)## [1] 2631
Filter the cells with low gene detection (low quality libraries) with less than 1000 genes for v2 and < 500 for v2.
#start with cells with many genes detected.
low.det.v3 <- WhichCells(data.filt, expression = nFeature_RNA < 1000 & orig.ident != "v2.1k")
low.det.v2 <- WhichCells(data.filt, expression = nFeature_RNA < 500 & orig.ident == "v2.1k")
# remove these cells
data.filt <- subset(data.filt, cells=setdiff(WhichCells(data.filt),c(low.det.v2,low.det.v3)))
# check number of cells
ncol(data.filt)## [1] 2531
Lets plot the same qc-stats another time.
VlnPlot(data.filt, features = "nFeature_RNA", pt.size = 0.1) + NoLegend()
VlnPlot(data.filt, features = "nCount_RNA", pt.size = 0.1) + NoLegend()
VlnPlot(data.filt, features = "percent.mito", pt.size = 0.1) + NoLegend()
VlnPlot(data.filt, features = "percent.ribo", pt.size = 0.1) + NoLegend()
# and check the number of cells per sample before and after filtering
table(Idents(alldata))##
## p3.1k v2.1k v3.1k
## 713 996 1222
table(Idents(data.filt))##
## p3.1k v2.1k v3.1k
## 526 933 1072
Seurat has a function for calculating cell cycle scores based on a list of know S-phase and G2/M-phase genes.
data.filt <- CellCycleScoring(
object = data.filt,
g2m.features = cc.genes$g2m.genes,
s.features = cc.genes$s.genes
)## Warning: The following features are not present in the object: MLF1IP, not
## searching for symbol synonyms
## Warning: The following features are not present in the object: FAM64A, HN1, not
## searching for symbol synonyms
VlnPlot(data.filt, features = c("S.Score","G2M.Score"))
In this case it looks like we only have a few cycling cells in the datasets.
options(stringsAsFactors = FALSE)
set.seed(32546)
To speed things up, we will continue working with the v3.1k dataset only. We will convert the Seurat object to a SCE object to work with the scater package. You can read more about SCE objects here.
Note: to create an SCE object directly from the count matrices, have a look at their tutorial at: https://bioconductor.org/packages/release/bioc/vignettes/scater/inst/doc/overview.html.
pbmc.sce <- SingleCellExperiment(assays = list(counts = as.matrix(v3.1k)))
pbmc.sce <- pbmc.sce[rowSums(counts(pbmc.sce) > 0) > 2,]
isSpike(pbmc.sce, "MT") <- grepl("^MT-", rownames(pbmc.sce))## Warning: 'isSpike<-' is deprecated.
## Use 'isSpike<-' instead.
## See help("Deprecated")
## Warning: 'spikeNames' is deprecated.
## See help("Deprecated")
## Warning: 'isSpike' is deprecated.
## See help("Deprecated")
pbmc.sce <- calculateQCMetrics(pbmc.sce)## Warning: 'calculateQCMetrics' is deprecated.
## Use 'perCellQCMetrics' or 'perFeatureQCMetrics' instead.
print(colnames(colData(pbmc.sce)))## [1] "is_cell_control"
## [2] "total_features_by_counts"
## [3] "log10_total_features_by_counts"
## [4] "total_counts"
## [5] "log10_total_counts"
## [6] "pct_counts_in_top_50_features"
## [7] "pct_counts_in_top_100_features"
## [8] "pct_counts_in_top_200_features"
## [9] "pct_counts_in_top_500_features"
## [10] "total_features_by_counts_endogenous"
## [11] "log10_total_features_by_counts_endogenous"
## [12] "total_counts_endogenous"
## [13] "log10_total_counts_endogenous"
## [14] "pct_counts_endogenous"
## [15] "pct_counts_in_top_50_features_endogenous"
## [16] "pct_counts_in_top_100_features_endogenous"
## [17] "pct_counts_in_top_200_features_endogenous"
## [18] "pct_counts_in_top_500_features_endogenous"
## [19] "total_features_by_counts_feature_control"
## [20] "log10_total_features_by_counts_feature_control"
## [21] "total_counts_feature_control"
## [22] "log10_total_counts_feature_control"
## [23] "pct_counts_feature_control"
## [24] "pct_counts_in_top_50_features_feature_control"
## [25] "pct_counts_in_top_100_features_feature_control"
## [26] "pct_counts_in_top_200_features_feature_control"
## [27] "pct_counts_in_top_500_features_feature_control"
## [28] "total_features_by_counts_MT"
## [29] "log10_total_features_by_counts_MT"
## [30] "total_counts_MT"
## [31] "log10_total_counts_MT"
## [32] "pct_counts_MT"
## [33] "pct_counts_in_top_50_features_MT"
## [34] "pct_counts_in_top_100_features_MT"
## [35] "pct_counts_in_top_200_features_MT"
## [36] "pct_counts_in_top_500_features_MT"
Filter out poor quality cells to avoid negative size factors. These steps are very similar to what we have already done on the combined Seurat object but now we perform them on one dataset only using the Scater package.
pbmc.sce <- pbmc.sce[, pbmc.sce$pct_counts_MT < 20]
pbmc.sce <- pbmc.sce[, (pbmc.sce$total_features_by_counts > 1000 & pbmc.sce$total_features_by_counts < 4100)]Create a new assay with unnormalized counts for comparison to post-normalization.
assay(pbmc.sce, "logcounts_raw") <- log2(counts(pbmc.sce) + 1)
plotRLE(pbmc.sce[,1:50], exprs_values = "logcounts_raw", style = "full")
Run PCA and save the result in a new object, as we will overwrite the PCA slot later.
raw.sce <- runPCA(pbmc.sce, exprs_values = "logcounts_raw")
scater::plotPCA(raw.sce, colour_by = "total_counts")
Plot the expression of the B cell marker MS4A1.
plotReducedDim(raw.sce, dimred = "PCA", by_exprs_values = "logcounts_raw",
colour_by = "MS4A1")
In the default normalization method in Seurat, counts for each cell are divided by the total counts for that cell and multiplied by the scale factor 10,000. This is then log transformed.
Here we use the filtered data from the counts slot of the SCE object to create a Seurat object. After normalization, we convert the result back into a SingleCellExperiment object for comparing plots.
pbmc.seu <- CreateSeuratObject(counts(pbmc.sce), project = "PBMC")
pbmc.seu <- NormalizeData(pbmc.seu)
pbmc.seu.sce <- as.SingleCellExperiment(pbmc.seu)
pbmc.seu.sce <- calculateQCMetrics(pbmc.seu.sce)## Warning: 'calculateQCMetrics' is deprecated.
## Use 'perCellQCMetrics' or 'perFeatureQCMetrics' instead.
Perform PCA and examine the normalization results with plotRLE and plotReducedDim. This time, use “logcounts” as the expression values to plot (or omit the parameter, as “logcounts” is the default value). Check some marker genes, for example GNLY (NK cells) or LYZ (monocytes).
plotRLE(pbmc.seu.sce[,1:50], style = "full")
pbmc.seu.sce <- runPCA(pbmc.seu.sce)
scater::plotPCA(pbmc.seu.sce, colour_by = "total_counts")
plotReducedDim(pbmc.seu.sce, dimred = "PCA", colour_by = "MS4A1")
The normalization procedure in scran is based on the deconvolution method by Lun et al (2016). Counts from many cells are pooled to avoid the drop-out problem. Pool-based size factors are then “deconvolved” into cell-based factors for cell-specific normalization. Clustering cells prior to normalization is not always necessary but it improves normalization accuracy by reducing the number of DE genes between cells in the same cluster.
qclust <- scran::quickCluster(pbmc.sce)## Warning: handling of spike-ins via 'isSpike()' is deprecated.
## Store spike-ins in 'altExps' instead.
pbmc.sce <- scran::computeSumFactors(pbmc.sce, clusters = qclust)## Warning: handling of spike-ins via 'isSpike()' is deprecated.
## Store spike-ins in 'altExps' instead.
summary(sizeFactors(pbmc.sce))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.3095 0.6559 0.8292 1.0000 1.1831 2.7118
pbmc.sce <- normalize(pbmc.sce)## Warning: 'normalizeSCE' is deprecated.
## Use 'logNormCounts' instead.
## See help("Deprecated")
## Warning: 'centreSizeFactors' is deprecated.
## See help("Deprecated")
## Warning in .get_all_sf_sets(object): spike-in set 'MT' should have its own size
## factors
Examine the results and compare to the log-normalized result. Are they different?
plotRLE(pbmc.sce[,1:50], exprs_values = "logcounts", exprs_logged = FALSE,
style = "full")
pbmc.sce <- runPCA(pbmc.sce)
scater::plotPCA(pbmc.sce, colour_by = "total_counts")
plotReducedDim(pbmc.sce, dimred = "PCA", colour_by = "MS4A1")
## 3. Feature selection
In the scran method for finding HVGs, a trend is first fitted to the technical variances. In the absence of spike-ins, this is done using the whole data, assuming that the majority of genes are not variably expressed. Then, the biological component of the variance for each endogenous gene is computed by subtracting the fitted value of the trend from the total variance. HVGs are then identified as those genes with the largest biological components. This avoids prioritizing genes that are highly variable due to technical factors such as sampling noise during RNA capture and library preparation. see the scran vignette for details.
dec <- modelGeneVar(pbmc.sce)
dec <- dec[!is.na(dec$FDR),]
top.hvgs <- order(dec$bio, decreasing = TRUE)
head(dec[top.hvgs,])## DataFrame with 6 rows and 6 columns
## mean total tech bio
## <numeric> <numeric> <numeric> <numeric>
## S100A9 2.19635761188417 9.9921611783313 0.812459482093853 9.17970169623745
## S100A8 1.96333931256695 8.9185913683661 0.825600780035344 8.09299058833076
## LYZ 2.144445380119 8.871163376052 0.815450514367692 8.05571286168431
## HLA-DRA 2.25495356643387 5.46022312130995 0.808059791165969 4.65216333014398
## CD74 2.83875237521701 4.48816089915459 0.772655780129843 3.71550511902475
## IGKC 1.01035988303756 4.4247505696434 0.727408602628664 3.69734196701473
## p.value FDR
## <numeric> <numeric>
## S100A9 0 0
## S100A8 1.48282482742041e-248 7.53077302352581e-245
## LYZ 2.11182408062784e-252 1.60878758462229e-248
## HLA-DRA 3.03427863892219e-87 1.15575673356546e-83
## CD74 1.60410027553219e-61 4.07334529966808e-58
## IGKC 1.7442340055978e-68 5.31502986185763e-65
dec$HVG <- (dec$FDR<0.00001)
hvg_genes <- rownames(dec[dec$FDR < 0.00001, ])
# plot highly variable genes
plot(dec$mean, dec$total, pch=16, cex=0.6, xlab="Mean log-expression",
ylab="Variance of log-expression")
o <- order(dec$mean)
lines(dec$mean[o], dec$tech[o], col="dodgerblue", lwd=2)
points(dec$mean[dec$HVG], dec$total[dec$HVG], col="red", pch=16)
## save the decomposed variance table and hvg_genes into metadata for safekeeping
metadata(pbmc.sce)$hvg_genes <- hvg_genes
metadata(pbmc.sce)$dec_var <- decWe choose genes that have a biological component that is significantly greater than zero, using a false discovery rate (FDR) of 5%.
plotExpression(pbmc.sce, features = rownames(dec[top.hvgs[1:10],]))
The default method in Seurat 3 is variance-stabilizing transformation. A trend is fitted to to predict the variance of each gene as a function of its mean. For each gene, the variance of standardized values is computed across all cells and used to rank the features. By default, 2000 top genes are returned.
pbmc.seu <- FindVariableFeatures(pbmc.seu, selection.method = "vst")
top10 <- head(VariableFeatures(pbmc.seu), 10)
vplot <- VariableFeaturePlot(pbmc.seu)
LabelPoints(plot = vplot, points = top10, repel = TRUE, xnudge = 0, ynudge = 0)## Warning: Transformation introduced infinite values in continuous x-axis
How many of the variable genes detected with scran are included in VariableFeatures in Seurat?
table(hvg_genes %in% VariableFeatures(pbmc.seu))##
## FALSE TRUE
## 5 70
sessionInfo()## R version 3.6.1 (2019-07-05)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Dutch_Netherlands.1252 LC_CTYPE=Dutch_Netherlands.1252
## [3] LC_MONETARY=Dutch_Netherlands.1252 LC_NUMERIC=C
## [5] LC_TIME=Dutch_Netherlands.1252
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] Matrix_1.2-18 scran_1.14.6
## [3] scater_1.14.6 ggplot2_3.3.2
## [5] SingleCellExperiment_1.8.0 SummarizedExperiment_1.16.1
## [7] DelayedArray_0.12.3 BiocParallel_1.20.1
## [9] matrixStats_0.56.0 Biobase_2.46.0
## [11] GenomicRanges_1.38.0 GenomeInfoDb_1.22.1
## [13] IRanges_2.20.2 S4Vectors_0.24.4
## [15] BiocGenerics_0.32.0 Seurat_3.2.0
##
## loaded via a namespace (and not attached):
## [1] plyr_1.8.6 igraph_1.2.5 lazyeval_0.2.2
## [4] splines_3.6.1 listenv_0.8.0 digest_0.6.25
## [7] htmltools_0.5.0 viridis_0.5.1 magrittr_1.5
## [10] tensor_1.5 cluster_2.1.0 ROCR_1.0-11
## [13] limma_3.42.2 globals_0.12.5 colorspace_1.4-1
## [16] rappdirs_0.3.1 ggrepel_0.8.2 xfun_0.17
## [19] dplyr_1.0.2 crayon_1.3.4 RCurl_1.95-4.12
## [22] jsonlite_1.7.0 spatstat_1.64-1 spatstat.data_1.4-3
## [25] survival_3.2-3 zoo_1.8-8 ape_5.4-1
## [28] glue_1.4.1 polyclip_1.10-0 gtable_0.3.0
## [31] zlibbioc_1.32.0 XVector_0.26.0 leiden_0.3.3
## [34] BiocSingular_1.2.2 future.apply_1.6.0 abind_1.4-5
## [37] scales_1.1.1 edgeR_3.28.1 miniUI_0.1.1.1
## [40] Rcpp_1.0.5 viridisLite_0.3.0 xtable_1.8-4
## [43] reticulate_1.16 dqrng_0.2.1 bit_1.1-14
## [46] rsvd_1.0.3 htmlwidgets_1.5.1 httr_1.4.2
## [49] RColorBrewer_1.1-2 ellipsis_0.3.1 ica_1.0-2
## [52] pkgconfig_2.0.3 farver_2.0.3 uwot_0.1.8
## [55] deldir_0.1-28 locfit_1.5-9.4 tidyselect_1.1.0
## [58] labeling_0.3 rlang_0.4.7 reshape2_1.4.4
## [61] later_1.1.0.1 munsell_0.5.0 tools_3.6.1
## [64] generics_0.0.2 ggridges_0.5.2 evaluate_0.14
## [67] stringr_1.4.0 fastmap_1.0.1 yaml_2.2.1
## [70] goftest_1.2-2 knitr_1.29 bit64_0.9-7
## [73] fitdistrplus_1.1-1 purrr_0.3.4 RANN_2.6.1
## [76] pbapply_1.4-3 future_1.18.0 nlme_3.1-148
## [79] mime_0.9 hdf5r_1.3.3 compiler_3.6.1
## [82] beeswarm_0.2.3 plotly_4.9.2.1 png_0.1-7
## [85] spatstat.utils_1.17-0 tibble_3.0.3 statmod_1.4.34
## [88] stringi_1.4.6 lattice_0.20-41 vctrs_0.3.2
## [91] pillar_1.4.6 lifecycle_0.2.0 lmtest_0.9-37
## [94] RcppAnnoy_0.0.16 BiocNeighbors_1.4.2 data.table_1.13.0
## [97] cowplot_1.1.0 bitops_1.0-6 irlba_2.3.3
## [100] httpuv_1.5.4 patchwork_1.0.1 R6_2.4.1
## [103] promises_1.1.1 KernSmooth_2.23-17 gridExtra_2.3
## [106] vipor_0.4.5 codetools_0.2-16 MASS_7.3-52
## [109] withr_2.2.0 sctransform_0.2.1 GenomeInfoDbData_1.2.2
## [112] mgcv_1.8-32 grid_3.6.1 rpart_4.1-15
## [115] tidyr_1.1.1 rmarkdown_2.3 DelayedMatrixStats_1.8.0
## [118] Rtsne_0.15 shiny_1.5.0 ggbeeswarm_0.6.0