vignettes/dataset.Rmd
dataset.RmdIn this workshop, we use example data from the TENxPBMCData package. This package provides an R / Bioconductor resource for representing and manipulating different single-cell RNA-seq data sets profiling peripheral blood mononuclear cells (PBMC) generated by 10x Genomics (https://support.10xgenomics.com/single-cell-gene-expression/datasets).
library(TENxPBMCData)The man page for the TENxPBMCData() function gives an idea of the datasets that are available from this package. It can be opened with the following command.
help(TENxPBMCData)Here, we use the pbmc3k dataset, which contains gene expression profiles for 2,700 single peripheral blood mononuclear cells. The first time this dataset is loaded, this command downloads the dataset to a local cache, which takes some time, depending on the speed of your internet connection. Subsequent times, the same command loads the dataset directly from the local cache.
sce <- TENxPBMCData(dataset = "pbmc3k")At this point we can inspect the dataset in the console.
sce
#> class: SingleCellExperiment
#> dim: 32738 2700
#> metadata(0):
#> assays(1): counts
#> rownames(32738): ENSG00000243485 ENSG00000237613 ... ENSG00000215616
#> ENSG00000215611
#> rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
#> colnames: NULL
#> colData names(11): Sample Barcode ... Individual Date_published
#> reducedDimNames(0):
#> altExpNames(0):The dataset is provided as an object of the SingleCellExperiment class. In particular, this summary view indicates that the following pieces of information are available:
counts
NULL
Barcode) and the donor identifier (Individual).Note that a SingleCellExperiment object (or, more generally, any SummarizedExperiment object) like this one already contains sufficient information to launch an interactive application instance to visualize the available data and metadata, using the iSEE() function.
For the purpose of this workshop, we first apply some preprocessing to the SingleCellExperiment object, in order to populate it with more information that can be visualized with iSEE.
We start by adding column names to the object, and use gene symbols instead of Ensembl IDs as row names. In the case where multiple Ensembl identifiers correspond to the same gene symbol, the scater::uniquifyFeatureNames function concatenates the Ensembl ID and the gene symbol in order to generate unique feature names.
library(scater)
colnames(sce) <- paste0("Cell", seq_len(ncol(sce)))
rownames(sce) <- scater::uniquifyFeatureNames(
ID = rowData(sce)$ENSEMBL_ID,
names = rowData(sce)$Symbol_TENx
)
head(rownames(sce))
#> [1] "MIR1302-10" "FAM138A" "OR4F5" "RP11-34P13.7" "RP11-34P13.8"
#> [6] "AL627309.1"Next, we use the scater package to calculate gene- and cell-level quality metrics. These metrics are added as columns to the rowData and colData slots of the SingleCellExperiment object, respectively. We also add some additional metrics that are not automatically computed by the scater package.
MT <- rownames(sce)[grep("^MT-", rownames(sce))]
sce <- scater::addPerCellQC(sce, subsets = list(MT = MT))
sce <- scater::addPerFeatureQC(sce)
sce$log10_total <- log10(sce$total)
rowData(sce)$n_cells <- as.integer(rowData(sce)$detected * ncol(sce))
rowData(sce)$log10_total <- log10(rowSums(assay(sce, "counts")) + 1)
sce
#> class: SingleCellExperiment
#> dim: 32738 2700
#> metadata(0):
#> assays(1): counts
#> rownames(32738): MIR1302-10 FAM138A ... AC002321.2 AC002321.1
#> rowData names(7): ENSEMBL_ID Symbol_TENx ... n_cells log10_total
#> colnames(2700): Cell1 Cell2 ... Cell2699 Cell2700
#> colData names(18): Sample Barcode ... total log10_total
#> reducedDimNames(0):
#> altExpNames(0):We filter out a few cells with a large fraction of the counts coming from mitochondrial genes, since these may be damaged cells. Notice the reduced number of columns in the dataset below.
(sce <- sce[, sce$subsets_MT_percent < 5])
#> class: SingleCellExperiment
#> dim: 32738 2643
#> metadata(0):
#> assays(1): counts
#> rownames(32738): MIR1302-10 FAM138A ... AC002321.2 AC002321.1
#> rowData names(7): ENSEMBL_ID Symbol_TENx ... n_cells log10_total
#> colnames(2643): Cell1 Cell2 ... Cell2699 Cell2700
#> colData names(18): Sample Barcode ... total log10_total
#> reducedDimNames(0):
#> altExpNames(0):Next, we calculate size factors and normalized and log-transformed expression values, using the scran and scater packages. Note that it is typically recommended to pre-cluster the cells before computing the size factors, as follows:
# set.seed(1000)
# clusters <- scran::quickCluster(sce, BSPARAM = IrlbaParam())
# sce <- scran::computeSumFactors(sce, cluster = clusters, min.mean = 0.1)However, for time reasons, we will skip the pre-clustering step in this workshop.
library(scran)
assay(sce, "counts") <- as(assay(sce, "counts"), "sparseMatrix") # dirty fix for DockerHub
sce <- scran::computeSumFactors(sce, min.mean = 0.1)
summary(sizeFactors(sce))
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.1164 0.7309 0.9004 1.0000 1.1635 9.8580
sce <- scater::logNormCounts(sce)In order to extract the most informative genes, we first model the mean-variance trend and decompose the variance into biological and technical components.
dec <- scran::modelGeneVar(sce)
top.dec <- dec[order(dec$bio, decreasing = TRUE), ]
head(top.dec)
#> DataFrame with 6 rows and 6 columns
#> mean total tech bio p.value FDR
#> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
#> LYZ 1.63950 3.92209 0.884955 3.03713 4.59557e-117 2.54472e-113
#> S100A9 1.07323 3.36899 0.778205 2.59079 2.72222e-110 9.04431e-107
#> HLA-DRA 1.56937 3.16577 0.880146 2.28562 5.73028e-68 7.93262e-65
#> FTL 3.67846 2.88499 0.708226 2.17676 2.46521e-94 5.85029e-91
#> CD74 2.23653 2.75839 0.854474 1.90392 1.34728e-50 1.01732e-47
#> CST3 1.10995 2.59172 0.790231 1.80149 7.24671e-53 6.01912e-50Next, we apply Principal Components Analysis (PCA) and t-distributed Stochastic Neighbor Embedding (t-SNE) to generate low-dimensional representations of the cells in our data set. These low-dimensional representations are added to the reducedDim slot of the SingleCellExperiment object.
library(BiocSingular)
set.seed(1000)
sce <- scran::denoisePCA(sce, technical = dec)
ncol(reducedDim(sce, "PCA"))
#> [1] 5
set.seed(1000)
sce <- scater::runTSNE(sce, dimred = "PCA", perplexity = 30)
sce <- scater::runUMAP(sce, dimred = "PCA")
sce
#> class: SingleCellExperiment
#> dim: 32738 2643
#> metadata(0):
#> assays(2): counts logcounts
#> rownames(32738): MIR1302-10 FAM138A ... AC002321.2 AC002321.1
#> rowData names(7): ENSEMBL_ID Symbol_TENx ... n_cells log10_total
#> colnames(2643): Cell1 Cell2 ... Cell2699 Cell2700
#> colData names(19): Sample Barcode ... log10_total sizeFactor
#> reducedDimNames(3): PCA TSNE UMAP
#> altExpNames(0):After this, we cluster the cells using a graph-based algorithm, and find ‘marker genes’ for each cluster as the genes that are significantly upregulated in the cluster compared to each of the other inferred clusters. The adjusted p-values from this test, for each cluster, are added to the rowData slot of the object.
snn.gr <- scran::buildSNNGraph(sce, use.dimred = "PCA")
clusters <- igraph::cluster_walktrap(snn.gr)
sce$Cluster <- factor(clusters$membership)
table(sce$Cluster)
#>
#> 1 2 3 4 5 6 7 8 9 10 11 12 13 14
#> 143 330 212 362 224 389 338 144 142 129 79 115 25 11
markers <- scran::findMarkers(sce, groups = sce$Cluster,
test.type = "t",
direction = "up", pval.type = "all")
for (i in names(markers)) {
rowData(sce)[, paste0("FDR_cluster", i)] <-
markers[[i]]$FDR[match(rownames(sce),
rownames(markers[[i]]))]
}
sce
#> class: SingleCellExperiment
#> dim: 32738 2643
#> metadata(0):
#> assays(2): counts logcounts
#> rownames(32738): MIR1302-10 FAM138A ... AC002321.2 AC002321.1
#> rowData names(21): ENSEMBL_ID Symbol_TENx ... FDR_cluster13
#> FDR_cluster14
#> colnames(2643): Cell1 Cell2 ... Cell2699 Cell2700
#> colData names(20): Sample Barcode ... sizeFactor Cluster
#> reducedDimNames(3): PCA TSNE UMAP
#> altExpNames(0):Finally, we assign a label to each cell, based on their individual transcriptome profiles, using the SingleR method and the Monaco immune data (https://doi.org/10.1016/j.celrep.2019.01.041) as a reference.
For each prediction, we assign the labels values to a specific colData element
library(SingleR)
library(celldex)
ref_monaco <- MonacoImmuneData()Here we assign the cell type according to the main classification scheme (this includes B cells, Basophils, CD4+ T cells, CD8+ T cells, Dendritic cells, Monocytes, Neutrophils, NK cells, Progenitors, and T cells)
pred_monaco_main <- SingleR(test = sce, ref = ref_monaco, labels = ref_monaco$label.main)
table(pred_monaco_main$labels)
#>
#> B cells CD4+ T cells CD8+ T cells Dendritic cells Monocytes
#> 348 894 294 42 632
#> NK cells Progenitors T cells
#> 160 15 258
sce$labels_main <- pred_monaco_main$labelsWe do something similar with a more fine-grained classification, this time including the cell subtypes (e.g., for B cells, the subtypes would include Exhausted B cells, Naive B cells, Non-switched memory B cells, Plasmablasts, and Switched memory B cells)
pred_monaco_fine <- SingleR(test = sce, ref = ref_monaco, labels = ref_monaco$label.fine)
table(pred_monaco_fine$labels)
#>
#> Central memory CD8 T cells Classical monocytes
#> 63 320
#> Effector memory CD8 T cells Exhausted B cells
#> 40 41
#> Follicular helper T cells Intermediate monocytes
#> 130 243
#> MAIT cells Myeloid dendritic cells
#> 79 42
#> Naive B cells Naive CD4 T cells
#> 202 250
#> Naive CD8 T cells Natural killer cells
#> 73 150
#> Non classical monocytes Non-switched memory B cells
#> 64 70
#> Non-Vd2 gd T cells Plasmablasts
#> 16 6
#> Plasmacytoid dendritic cells Progenitor cells
#> 4 14
#> Switched memory B cells T regulatory cells
#> 29 95
#> Terminal effector CD4 T cells Terminal effector CD8 T cells
#> 4 65
#> Th1 cells Th1/Th17 cells
#> 144 125
#> Th17 cells Th2 cells
#> 65 179
#> Vd2 gd T cells
#> 130
sce$labels_fine <- pred_monaco_fine$labelsSimilarly, we use the information contained in the cell ontology labels.
pred_monaco_ont <- SingleR(test = sce, ref = ref_monaco, labels = ref_monaco$label.ont)
table(pred_monaco_ont$labels)
#>
#> CL:0000236 CL:0000545 CL:0000546 CL:0000623 CL:0000782 CL:0000784 CL:0000788
#> 43 130 185 154 42 4 201
#> CL:0000798 CL:0000815 CL:0000860 CL:0000875 CL:0000895 CL:0000899 CL:0000900
#> 108 99 319 64 235 58 80
#> CL:0000907 CL:0000912 CL:0000913 CL:0000940 CL:0000970 CL:0000972 CL:0000980
#> 70 131 62 86 69 28 6
#> CL:0001044 CL:0001062 CL:0002038 CL:0002043 CL:0002393
#> 3 63 143 13 247
sce$labels_ont <- pred_monaco_ont$labelsThe next table shows the relationship between the coarse and fine grained assignments in the data at hand.
table(sce$labels_main,
sce$labels_fine)
#>
#> Central memory CD8 T cells Classical monocytes
#> B cells 0 0
#> CD4+ T cells 5 0
#> CD8+ T cells 49 0
#> Dendritic cells 0 4
#> Monocytes 0 315
#> NK cells 0 0
#> Progenitors 0 1
#> T cells 9 0
#>
#> Effector memory CD8 T cells Exhausted B cells
#> B cells 0 41
#> CD4+ T cells 0 0
#> CD8+ T cells 29 0
#> Dendritic cells 0 0
#> Monocytes 0 0
#> NK cells 0 0
#> Progenitors 0 0
#> T cells 11 0
#>
#> Follicular helper T cells Intermediate monocytes MAIT cells
#> B cells 0 0 0
#> CD4+ T cells 125 0 1
#> CD8+ T cells 4 0 2
#> Dendritic cells 0 1 0
#> Monocytes 0 242 0
#> NK cells 0 0 0
#> Progenitors 0 0 0
#> T cells 1 0 76
#>
#> Myeloid dendritic cells Naive B cells Naive CD4 T cells
#> B cells 0 200 0
#> CD4+ T cells 0 0 200
#> CD8+ T cells 0 0 46
#> Dendritic cells 33 0 0
#> Monocytes 9 1 0
#> NK cells 0 0 0
#> Progenitors 0 0 0
#> T cells 0 1 4
#>
#> Naive CD8 T cells Natural killer cells
#> B cells 0 0
#> CD4+ T cells 9 0
#> CD8+ T cells 63 0
#> Dendritic cells 0 0
#> Monocytes 0 0
#> NK cells 0 148
#> Progenitors 0 0
#> T cells 1 2
#>
#> Non classical monocytes Non-switched memory B cells
#> B cells 0 70
#> CD4+ T cells 0 0
#> CD8+ T cells 0 0
#> Dendritic cells 0 0
#> Monocytes 64 0
#> NK cells 0 0
#> Progenitors 0 0
#> T cells 0 0
#>
#> Non-Vd2 gd T cells Plasmablasts Plasmacytoid dendritic cells
#> B cells 0 6 0
#> CD4+ T cells 0 0 0
#> CD8+ T cells 3 0 0
#> Dendritic cells 0 0 4
#> Monocytes 0 0 0
#> NK cells 4 0 0
#> Progenitors 0 0 0
#> T cells 9 0 0
#>
#> Progenitor cells Switched memory B cells T regulatory cells
#> B cells 0 29 0
#> CD4+ T cells 0 0 90
#> CD8+ T cells 0 0 4
#> Dendritic cells 0 0 0
#> Monocytes 0 0 0
#> NK cells 0 0 0
#> Progenitors 14 0 0
#> T cells 0 0 1
#>
#> Terminal effector CD4 T cells Terminal effector CD8 T cells
#> B cells 0 0
#> CD4+ T cells 0 0
#> CD8+ T cells 2 37
#> Dendritic cells 0 0
#> Monocytes 0 0
#> NK cells 0 5
#> Progenitors 0 0
#> T cells 2 23
#>
#> Th1 cells Th1/Th17 cells Th17 cells Th2 cells Vd2 gd T cells
#> B cells 0 1 0 0 1
#> CD4+ T cells 122 104 62 172 4
#> CD8+ T cells 15 3 1 7 29
#> Dendritic cells 0 0 0 0 0
#> Monocytes 0 0 0 0 1
#> NK cells 0 0 0 0 3
#> Progenitors 0 0 0 0 0
#> T cells 7 17 2 0 92This concludes the preparation of the data. We have now a SingleCellExperiment object that contains different types of abundance values, representations in reduced dimensions, as well as a range of row (feature) and column (cell) metadata.
We use a function implemented in this workshop package, to cache the precomputed object in an file that we will use as a starting point in the recipes vignette.
library(iSEEWorkshopEuroBioc2020)
cache_demo_sce(sce, "pbmc3k")
#> Saving object to /github/home/.cache/iSEEWorkshopEuroBioc2020/pbmc3k.rds