In this tutorial, we will run all tutorials with a set of 6 PBMC 10x
datasets from 3 covid-19 patients and 3 healthy controls, the samples
have been subsampled to 1500 cells per sample. They are part of the
github repo and if you have cloned the repo they should be available in
folder: labs/data/covid_data_GSE149689. Instructions on how
to download them can also be found in the Precourse material.
webpath <- "https://raw.githubusercontent.com/NBISweden/workshop-scRNAseq/new_dataset/labs/data/covid_data_GSE149689/sub/"
dir.create("./data/raw", recursive = T)## Warning in dir.create("./data/raw", recursive = T): './data/raw' already existsfile_list <- c("Normal_PBMC_13.h5", "Normal_PBMC_14.h5", "Normal_PBMC_5.h5", "nCoV_PBMC_15.h5",
"nCoV_PBMC_17.h5", "nCoV_PBMC_1.h5")
for (i in file_list) {
download.file(url = paste0(webpath, i), destfile = paste0("./data/raw/", i))
}With data in place, now we can start loading libraries we will use in this tutorial.
suppressMessages(require(scater))
suppressMessages(require(scran))
suppressMessages(require(cowplot))
suppressMessages(require(org.Hs.eg.db))
if (!require(DoubletFinder)) {
remotes::install_github("chris-mcginnis-ucsf/DoubletFinder", upgrade = F, dependencies = F)
}## Loading required package: DoubletFindersuppressMessages(require(DoubletFinder))We can first load the data individually by reading directly from HDF5 file format (.h5).
cov.15 <- Seurat::Read10X_h5(filename = "data/raw/nCoV_PBMC_15.h5", use.names = T)
cov.1 <- Seurat::Read10X_h5(filename = "data/raw/nCoV_PBMC_1.h5", use.names = T)
cov.17 <- Seurat::Read10X_h5(filename = "data/raw/nCoV_PBMC_17.h5", use.names = T)
ctrl.5 <- Seurat::Read10X_h5(filename = "data/raw/Normal_PBMC_5.h5", use.names = T)
ctrl.13 <- Seurat::Read10X_h5(filename = "data/raw/Normal_PBMC_13.h5", use.names = T)
ctrl.14 <- Seurat::Read10X_h5(filename = "data/raw/Normal_PBMC_14.h5", use.names = T)We can now load the expression matricies into objects and then merge
them into a single merged object. Each analysis workflow (Seurat,
Scater, Scranpy, etc) has its own way of storing data. We will add
dataset labels as cell.ids just in case you have overlapping barcodes
between the datasets. After that we add a column Chemistry
in the metadata for plotting later on.
sce <- SingleCellExperiment(assays = list(counts = cbind(cov.1, cov.15, cov.17, ctrl.5,
ctrl.13, ctrl.14)))
dim(sce)## [1] 33538 9000# Adding metadata
sce@colData$sample <- unlist(sapply(c("cov.1", "cov.15", "cov.17", "ctrl.5", "ctrl.13",
"ctrl.14"), function(x) rep(x, ncol(get(x)))))
sce@colData$type <- ifelse(grepl("cov", sce@colData$sample), "Covid", "Control")Once you have created the merged Seurat object, the count matrices and individual count matrices and objects are not needed anymore. It is a good idea to remove them and run garbage collect to free up some memory.
# remove all objects that will not be used.
rm(cov.15, cov.1, cov.17, ctrl.5, ctrl.13, ctrl.14)
# run garbage collect to free up memory
gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 9046114 483.2 16839967 899.4 12029455 642.5
## Vcells 35108062 267.9 97347381 742.8 72827206 555.7Here it is how the count matrix and the metatada look like for every cell.
head(counts(sce)[, 1:10])
head(sce@colData, 10)## 6 x 10 sparse Matrix of class "dgCMatrix"
##
## MIR1302-2HG . . . . . . . . . .
## FAM138A . . . . . . . . . .
## OR4F5 . . . . . . . . . .
## AL627309.1 . . . . . . . . . .
## AL627309.3 . . . . . . . . . .
## AL627309.2 . . . . . . . . . .
## DataFrame with 10 rows and 2 columns
## sample type
## <character> <character>
## AGGGTCCCATGACCCG-1 cov.1 Covid
## TACCCACAGCGGGTTA-1 cov.1 Covid
## CCCAACTTCATATGGC-1 cov.1 Covid
## TCAAGTGTCCGAACGC-1 cov.1 Covid
## ATTCCTAGTGACTGTT-1 cov.1 Covid
## GTGTTCCGTGGGCTCT-1 cov.1 Covid
## CCTAAGACAGATTAAG-1 cov.1 Covid
## AATAGAGAGGGTTAGC-1 cov.1 Covid
## GGGTCACTCACCTACC-1 cov.1 Covid
## TCCTCTTGTACAGTCT-1 cov.1 CovidHaving the data in a suitable format, we can start calculating some quality metrics. We can for example calculate the percentage of mitocondrial and ribosomal genes per cell and add to the metadata. This will be helpfull to visualize them across different metadata parameteres (i.e. datasetID and chemistry version). There are several ways of doing this, and here manually calculate the proportion of mitochondrial reads and add to the metadata table.
Citing from “Simple Single Cell” workflows (Lun, McCarthy & Marioni, 2017): “High proportions are indicative of poor-quality cells (Islam et al. 2014; Ilicic et al. 2016), possibly because of loss of cytoplasmic RNA from perforated cells. The reasoning is that mitochondria are larger than individual transcript molecules and less likely to escape through tears in the cell membrane.”
First, let Scran calculate some general qc-stats for genes and cells
with the function perCellQCMetrics. It can also calculate
proportion of counts for specific gene subsets, so first we need to
define which genes are mitochondrial, ribosomal and hemoglogin.
# Mitochondrial genes
mito_genes <- rownames(sce)[grep("^MT-", rownames(sce))]
# Ribosomal genes
ribo_genes <- rownames(sce)[grep("^RP[SL]", rownames(sce))]
# Hemoglobin genes - includes all genes starting with HB except HBP.
hb_genes <- rownames(sce)[grep("^HB[^(P)]", rownames(sce))]sce <- addPerCellQC(sce, flatten = T, subsets = list(mt = mito_genes, hb = hb_genes,
ribo = ribo_genes))
head(colData(sce))## DataFrame with 6 rows and 14 columns
## sample type sum detected subsets_mt_sum
## <character> <character> <numeric> <integer> <numeric>
## AGGGTCCCATGACCCG-1 cov.1 Covid 7698 2140 525
## TACCCACAGCGGGTTA-1 cov.1 Covid 13416 3391 952
## CCCAACTTCATATGGC-1 cov.1 Covid 16498 3654 1253
## TCAAGTGTCCGAACGC-1 cov.1 Covid 1425 608 141
## ATTCCTAGTGACTGTT-1 cov.1 Covid 7535 1808 470
## GTGTTCCGTGGGCTCT-1 cov.1 Covid 4378 1345 352
## subsets_mt_detected subsets_mt_percent subsets_hb_sum
## <integer> <numeric> <numeric>
## AGGGTCCCATGACCCG-1 11 6.81995 2
## TACCCACAGCGGGTTA-1 11 7.09600 6
## CCCAACTTCATATGGC-1 12 7.59486 1
## TCAAGTGTCCGAACGC-1 10 9.89474 1
## ATTCCTAGTGACTGTT-1 11 6.23756 4
## GTGTTCCGTGGGCTCT-1 10 8.04020 1
## subsets_hb_detected subsets_hb_percent subsets_ribo_sum
## <integer> <numeric> <numeric>
## AGGGTCCCATGACCCG-1 1 0.02598077 2564
## TACCCACAGCGGGTTA-1 2 0.04472272 2264
## CCCAACTTCATATGGC-1 1 0.00606134 2723
## TCAAGTGTCCGAACGC-1 1 0.07017544 444
## ATTCCTAGTGACTGTT-1 3 0.05308560 3397
## GTGTTCCGTGGGCTCT-1 1 0.02284148 1588
## subsets_ribo_detected subsets_ribo_percent total
## <integer> <numeric> <numeric>
## AGGGTCCCATGACCCG-1 82 33.3074 7698
## TACCCACAGCGGGTTA-1 85 16.8754 13416
## CCCAACTTCATATGGC-1 87 16.5050 16498
## TCAAGTGTCCGAACGC-1 68 31.1579 1425
## ATTCCTAGTGACTGTT-1 81 45.0829 7535
## GTGTTCCGTGGGCTCT-1 79 36.2723 4378Here is an example on how to calculate proportion mitochondria in another way:
# Way2: Doing it manually
sce@colData$percent_mito <- Matrix::colSums(counts(sce)[mito_genes, ])/sce@colData$totalNow we can plot some of the QC-features as violin plots.
# total is total UMIs per cell detected is number of detected genes. the
# different gene subset percentages are listed as subsets_mt_percent etc.
plot_grid(plotColData(sce, y = "detected", x = "sample", colour_by = "sample"), plotColData(sce,
y = "total", x = "sample", colour_by = "sample"), plotColData(sce, y = "subsets_mt_percent",
x = "sample", colour_by = "sample"), plotColData(sce, y = "subsets_ribo_percent",
x = "sample", colour_by = "sample"), plotColData(sce, y = "subsets_hb_percent",
x = "sample", colour_by = "sample"), ncol = 3)