suppressPackageStartupMessages({
library(scater)
library(scran)
library(dplyr)
library(patchwork)
library(ggplot2)
library(pheatmap)
library(scPred)
library(scmap)
})
Note
Code chunks run R commands unless otherwise specified.
Celltype prediction can either be performed on indiviudal cells where each cell gets a predicted celltype label, or on the level of clusters. All methods are based on similarity to other datasets, single cell or sorted bulk RNAseq, or uses known marker genes for each cell type.
We will select one sample from the Covid data, ctrl_13 and predict celltype by cell on that sample.
Some methods will predict a celltype to each cell based on what it is most similar to, even if that celltype is not included in the reference. Other methods include an uncertainty so that cells with low similarity scores will be unclassified.
There are multiple different methods to predict celltypes, here we will just cover a few of those.
We will use a reference PBMC dataset from the scPred package which is provided as a Seurat object with counts. And we will test classification based on the scPred and scMap methods. Finally we will use gene set enrichment predict celltype based on the DEGs of each cluster.
1 Read data
First, lets load required libraries
Let’s read in the saved Covid-19 data object from the clustering step.
# download pre-computed data if missing or long compute
fetch_data <- TRUE
# url for source and intermediate data
path_data <- "https://export.uppmax.uu.se/naiss2023-23-3/workshops/workshop-scrnaseq"
path_file <- "data/covid/results/bioc_covid_qc_dr_int_cl.rds"
if (!dir.exists(dirname(path_file))) dir.create(dirname(path_file), recursive = TRUE)
if (fetch_data && !file.exists(path_file)) download.file(url = file.path(path_data, "covid/results/bioc_covid_qc_dr_int_cl.rds"), destfile = path_file)
alldata <- readRDS(path_file)Let’s read in the saved Covid-19 data object from the clustering step.
ctrl.sce <- alldata[, alldata@colData$sample == "ctrl.13"]
# remove all old dimensionality reductions as they will mess up the analysis further down
reducedDims(ctrl.sce) <- NULL2 Reference data
Load the reference dataset with annotated labels.
reference <- scPred::pbmc_1
referenceAn object of class Seurat
32838 features across 3500 samples within 1 assay
Active assay: RNA (32838 features, 0 variable features)Convert to a SCE object.
ref.sce <- Seurat::as.SingleCellExperiment(reference)Rerun analysis pipeline. Run normalization, feature selection and dimensionality reduction
# Normalize
ref.sce <- computeSumFactors(ref.sce)
ref.sce <- logNormCounts(ref.sce)
# Variable genes
var.out <- modelGeneVar(ref.sce, method = "loess")
hvg.ref <- getTopHVGs(var.out, n = 1000)
# Dim reduction
ref.sce <- runPCA(ref.sce,
exprs_values = "logcounts", scale = T,
ncomponents = 30, subset_row = hvg.ref
)
ref.sce <- runUMAP(ref.sce, dimred = "PCA")plotReducedDim(ref.sce, dimred = "UMAP", colour_by = "cell_type")
Run all steps of the analysis for the ctrl sample as well. Use the clustering from the integration lab with resolution 0.5.
# Normalize
ctrl.sce <- computeSumFactors(ctrl.sce)
ctrl.sce <- logNormCounts(ctrl.sce)
# Variable genes
var.out <- modelGeneVar(ctrl.sce, method = "loess")
hvg.ctrl <- getTopHVGs(var.out, n = 1000)
# Dim reduction
ctrl.sce <- runPCA(ctrl.sce, exprs_values = "logcounts", scale = T, ncomponents = 30, subset_row = hvg.ctrl)
ctrl.sce <- runUMAP(ctrl.sce, dimred = "PCA")plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "louvain_SNNk15")
3 scMap
The scMap package is one method for projecting cells from a scRNA-seq experiment on to the cell-types or individual cells identified in a different experiment. It can be run on different levels, either projecting by cluster or by single cell, here we will try out both.
For scmap cell type labels must be stored in the cell_type1 column of the colData slots, and gene ids that are consistent across both datasets must be stored in the feature_symbol column of the rowData slots.
3.1 scMap cluster
# add in slot cell_type1
ref.sce@colData$cell_type1 <- ref.sce@colData$cell_type
# create a rowData slot with feature_symbol
rd <- data.frame(feature_symbol = rownames(ref.sce))
rownames(rd) <- rownames(ref.sce)
rowData(ref.sce) <- rd
# same for the ctrl dataset
# create a rowData slot with feature_symbol
rd <- data.frame(feature_symbol = rownames(ctrl.sce))
rownames(rd) <- rownames(ctrl.sce)
rowData(ctrl.sce) <- rdThen we can select variable features in both datasets.
# select features
counts(ctrl.sce) <- as.matrix(counts(ctrl.sce))
logcounts(ctrl.sce) <- as.matrix(logcounts(ctrl.sce))
ctrl.sce <- selectFeatures(ctrl.sce, suppress_plot = TRUE)
counts(ref.sce) <- as.matrix(counts(ref.sce))
logcounts(ref.sce) <- as.matrix(logcounts(ref.sce))
ref.sce <- selectFeatures(ref.sce, suppress_plot = TRUE)Then we need to index the reference dataset by cluster, default is the clusters in cell_type1.
ref.sce <- indexCluster(ref.sce)Now we project the Covid-19 dataset onto that index.
project_cluster <- scmapCluster(
projection = ctrl.sce,
index_list = list(
ref = metadata(ref.sce)$scmap_cluster_index
)
)
# projected labels
table(project_cluster$scmap_cluster_labs)
B cell CD4 T cell CD8 T cell cDC cMono ncMono
70 104 122 38 213 161
NK cell pDC Plasma cell unassigned
287 2 1 175 Then add the predictions to metadata and plot UMAP.
# add in predictions
ctrl.sce@colData$scmap_cluster <- project_cluster$scmap_cluster_labs
plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "scmap_cluster")
4 scMap cell
We can instead index the refernce data based on each single cell and project our data onto the closest neighbor in that dataset.
ref.sce <- indexCell(ref.sce)Again we need to index the reference dataset.
project_cell <- scmapCell(
projection = ctrl.sce,
index_list = list(
ref = metadata(ref.sce)$scmap_cell_index
)
)We now get a table with index for the 5 nearest neigbors in the reference dataset for each cell in our dataset. We will select the celltype of the closest neighbor and assign it to the data.
cell_type_pred <- colData(ref.sce)$cell_type1[project_cell$ref[[1]][1, ]]
table(cell_type_pred)cell_type_pred
B cell CD4 T cell CD8 T cell cDC cMono ncMono NK cell
101 170 292 52 212 176 169
pDC
1 Then add the predictions to metadata and plot umap.
# add in predictions
ctrl.sce@colData$scmap_cell <- cell_type_pred
plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "scmap_cell")
Plot both:
wrap_plots(
plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "scmap_cluster"),
plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "scmap_cell"),
ncol = 2
)
5 scPred
scPred will train a classifier based on all principal components. First, getFeatureSpace() will create a scPred object stored in the @misc slot where it extracts the PCs that best separates the different celltypes. Then trainModel() will do the actual training for each celltype.
scPred works with Seurat objects, so we will convert both objects to seurat objects. You may see a lot of warnings about renaming things, but as long as you do not see an Error, you should be fine.
suppressPackageStartupMessages(library(Seurat))
reference <- Seurat::as.Seurat(ref.sce)
ctrl <- Seurat::as.Seurat(ctrl.sce)The loadings matrix is lost when converted to Seurat object, and scPred needs that information. So we need to rerun PCA with Seurat and the same hvgs.
VariableFeatures(reference) <- hvg.ref
reference <- reference %>%
ScaleData(verbose = F) %>%
RunPCA(verbose = F)
VariableFeatures(ctrl) <- hvg.ctrl
ctrl <- ctrl %>%
ScaleData(verbose = F) %>%
RunPCA(verbose = F)reference <- getFeatureSpace(reference, "cell_type")● Extracting feature space for each cell type...
DONE!reference <- trainModel(reference)● Training models for each cell type...
DONE!scPred will train a classifier based on all principal components. First, getFeatureSpace() will create a scPred object stored in the @misc slot where it extracts the PCs that best separates the different celltypes. Then trainModel() will do the actual training for each celltype.
get_scpred(reference)'scPred' object
✔ Prediction variable = cell_type
✔ Discriminant features per cell type
✔ Training model(s)
Summary
|Cell type | n| Features|Method | ROC| Sens| Spec|
|:-----------|----:|--------:|:---------|-----:|-----:|-----:|
|B cell | 280| 50|svmRadial | 1.000| 1.000| 1.000|
|CD4 T cell | 1620| 50|svmRadial | 0.994| 0.972| 0.963|
|CD8 T cell | 945| 50|svmRadial | 0.973| 0.859| 0.971|
|cDC | 26| 50|svmRadial | 0.994| 0.727| 0.999|
|cMono | 212| 50|svmRadial | 1.000| 0.957| 0.997|
|ncMono | 79| 50|svmRadial | 1.000| 0.962| 0.999|
|NK cell | 312| 50|svmRadial | 0.998| 0.926| 0.995|
|pDC | 20| 50|svmRadial | 1.000| 0.950| 1.000|
|Plasma cell | 6| 50|svmRadial | 1.000| 1.000| 1.000|You can optimize parameters for each dataset by chaining parameters and testing different types of models, see more at: https://powellgenomicslab.github.io/scPred/articles/introduction.html. But for now, we will continue with this model. Now, let’s predict celltypes on our data, where scPred will align the two datasets with Harmony and then perform classification.
ctrl <- scPredict(ctrl, reference)● Matching reference with new dataset...
─ 1000 features present in reference loadings
─ 938 features shared between reference and new dataset
─ 93.8% of features in the reference are present in new dataset
● Aligning new data to reference...
● Classifying cells...
DONE!DimPlot(ctrl, group.by = "scpred_prediction", label = T, repel = T) + NoAxes()
Now plot how many cells of each celltypes can be found in each cluster.
ggplot(ctrl@meta.data, aes(x = louvain_SNNk15, fill = scpred_prediction)) +
geom_bar() +
theme_classic()
Add the predictions into the SCE object
ctrl.sce@colData$scpred_prediction <- ctrl$scpred_prediction6 Compare results
Now we will compare the output of the two methods using the convenient function in scPred crossTab() that prints the overlap between two metadata slots.
crossTab(ctrl, "scmap_cell", "scpred_prediction")7 GSEA with celltype markers
Another option, where celltype can be classified on cluster level is to use gene set enrichment among the DEGs with known markers for different celltypes. Similar to how we did functional enrichment for the DEGs in the differential expression exercise. There are some resources for celltype gene sets that can be used. Such as CellMarker, PanglaoDB or celltype gene sets at MSigDB. We can also look at overlap between DEGs in a reference dataset and the dataset you are analyzing.
7.1 DEG overlap
First, lets extract top DEGs for our Covid-19 dataset and the reference dataset. When we run differential expression for our dataset, we want to report as many genes as possible, hence we set the cutoffs quite lenient.
# run differential expression in our dataset, using clustering at resolution 0.3
DGE_list <- scran::findMarkers(
x = alldata,
groups = as.character(alldata@colData$louvain_SNNk15),
pval.type = "all",
min.prop = 0
)# Compute differential gene expression in reference dataset (that has cell annotation)
ref_DGE <- scran::findMarkers(
x = ref.sce,
groups = as.character(ref.sce@colData$cell_type),
pval.type = "all",
direction = "up"
)
# Identify the top cell marker genes in reference dataset
# select top 50 with hihgest foldchange among top 100 signifcant genes.
ref_list <- lapply(ref_DGE, function(x) {
x$logFC <- rowSums(as.matrix(x[, grep("logFC", colnames(x))]))
x %>%
as.data.frame() %>%
filter(p.value < 0.01) %>%
top_n(-100, p.value) %>%
top_n(50, logFC) %>%
rownames()
})
unlist(lapply(ref_list, length)) B cell CD4 T cell CD8 T cell cDC cMono ncMono
50 50 19 17 50 50
NK cell pDC Plasma cell
50 50 24 Now we can run GSEA for the DEGs from our dataset and check for enrichment of top DEGs in the reference dataset.
suppressPackageStartupMessages(library(fgsea))
# run fgsea for each of the clusters in the list
res <- lapply(DGE_list, function(x) {
x$logFC <- rowSums(as.matrix(x[, grep("logFC", colnames(x))]))
gene_rank <- setNames(x$logFC, rownames(x))
fgseaRes <- fgsea(pathways = ref_list, stats = gene_rank, nperm = 10000)
return(fgseaRes)
})
names(res) <- names(DGE_list)
# You can filter and resort the table based on ES, NES or pvalue
res <- lapply(res, function(x) {
x[x$pval < 0.1, ]
})
res <- lapply(res, function(x) {
x[x$size > 2, ]
})
res <- lapply(res, function(x) {
x[order(x$NES, decreasing = T), ]
})
res$`1`
pathway pval padj ES NES nMoreExtreme size
1: cMono 0.0001102050 0.0004959224 0.9563653 1.702528 0 47
2: ncMono 0.0001096852 0.0004959224 0.8757519 1.564034 0 49
3: NK cell 0.0011299435 0.0017341040 -0.7183156 -1.792347 0 49
4: CD8 T cell 0.0005042864 0.0015128593 -0.9226336 -1.877822 0 18
5: CD4 T cell 0.0011560694 0.0017341040 -0.8765216 -2.191726 0 50
6: B cell 0.0010775862 0.0017341040 -0.9056237 -2.249271 0 47
leadingEdge
1: S100A8,S100A9,LYZ,S100A12,VCAN,FCN1,...
2: AIF1,S100A11,S100A4,SERPINA1,PSAP,LST1,...
3: GNLY,NKG7,B2M,CTSW,GZMA,GZMM,...
4: IL32,CCL5,GZMH,CD3D,CD2,CD8A,...
5: RPL3,RPS4X,RPS27A,PIK3IP1,RPS3,EEF1A1,...
6: CXCR4,RPS5,CD52,CD79B,MS4A1,RPL18A,...
$`10`
pathway pval padj ES NES nMoreExtreme size
1: cMono 0.0031746032 0.0057142857 0.8297227 1.992805 0 47
2: CD8 T cell 0.0309092834 0.0397405072 -0.7878850 -1.304310 291 18
3: cDC 0.0139771283 0.0209656925 -0.8182753 -1.347905 131 17
4: NK cell 0.0019601774 0.0044103993 -0.7630056 -1.350465 18 49
5: B cell 0.0008258491 0.0037163208 -0.7773761 -1.372465 7 47
6: Plasma cell 0.0012603718 0.0037811154 -0.8420529 -1.421526 11 24
7: CD4 T cell 0.0001030822 0.0009277394 -0.9105099 -1.613758 0 50
leadingEdge
1: S100A9,S100A8,LYZ,VCAN,S100A12,CD14,...
2: CCL5,IL32,GZMH,CD3D,CD2,LYAR,...
3: HLA-DPA1,HLA-DPB1,HLA-DQB1,HLA-DRA,HLA-DRB5,HLA-DMA
4: BIN2,GNLY,HCST,CST7,JAK1,CTSW,...
5: RPL23A,CD52,RPS23,RPS5,RPS11,RPL13A,...
6: RPL36AL,COX7A2,SUB1,DAD1,CYCS,PEBP1,...
7: RPS29,RPL35A,RPL21,RPL38,RPS25,RPS16,...
$`2`
pathway pval padj ES NES nMoreExtreme size
1: B cell 0.0001516990 0.0006656805 0.9641111 1.943584 0 47
2: CD4 T cell 0.0001508523 0.0006656805 0.8943277 1.817384 0 50
3: cDC 0.0009815148 0.0012619476 0.9382779 1.633196 5 17
4: CD8 T cell 0.0007702182 0.0011553273 -0.9083063 -1.716291 2 18
5: cMono 0.0005865103 0.0010557185 -0.8175352 -1.824417 1 47
6: ncMono 0.0002958580 0.0006656805 -0.8882121 -1.991298 0 49
7: NK cell 0.0002958580 0.0006656805 -0.9006980 -2.019290 0 49
leadingEdge
1: MS4A1,CD37,TNFRSF13C,CXCR4,CD79B,BANK1,...
2: RPS29,RPS6,RPL32,RPL36,RPL5,RPS25,...
3: HLA-DRA,HLA-DPB1,HLA-DQB1,HLA-DRB1,HLA-DPA1,HLA-DMA,...
4: CCL5,IL32,GZMH,CD3D,CD2,LYAR,...
5: S100A9,LYZ,S100A6,S100A8,TYROBP,FCN1,...
6: S100A4,FCER1G,IFITM3,S100A11,AIF1,TIMP1,...
7: HCST,NKG7,ITGB2,MYO1F,GNLY,CST7,...
$`3`
pathway pval padj ES NES nMoreExtreme size
1: CD4 T cell 0.0001574555 0.0008207934 0.9797631 2.007371 0 50
2: NK cell 0.0662106703 0.0993160055 -0.6135878 -1.368627 241 49
3: cDC 0.0004725898 0.0010633270 -0.9266903 -1.727259 1 17
4: pDC 0.0008190008 0.0014742015 -0.7974487 -1.763510 2 47
5: cMono 0.0002730003 0.0008207934 -0.9207862 -2.036264 0 47
6: ncMono 0.0002735978 0.0008207934 -0.9490592 -2.116907 0 49
leadingEdge
1: IL7R,LDHB,PIK3IP1,RPS29,RPS12,RPL3,...
2: NKG7,GNLY,ITGB2,MYO1F,FGFBP2,CST7,...
3: HLA-DRA,HLA-DRB1,HLA-DPA1,HLA-DPB1,HLA-DQB1,HLA-DMA,...
4: PLEK,NPC2,CTSB,PTPRE,IRF8,PLAC8,...
5: S100A9,LYZ,S100A8,TYROBP,FCN1,APLP2,...
6: FCER1G,IFITM3,SAT1,PSAP,FTH1,LYN,...
$`4`
pathway pval padj ES NES nMoreExtreme size
1: B cell 0.0001577287 0.0006158057 0.9602781 1.951884 0 47
2: CD4 T cell 0.0001569859 0.0006158057 0.9032855 1.851547 0 50
3: cDC 0.0003421143 0.0006158057 0.9472752 1.664561 1 17
4: CD8 T cell 0.0004867364 0.0007021924 -0.9106557 -1.715908 1 18
5: cMono 0.0005461496 0.0007021924 -0.7951206 -1.769638 1 47
6: ncMono 0.0002739726 0.0006158057 -0.9007143 -2.018766 0 49
7: NK cell 0.0002739726 0.0006158057 -0.9029475 -2.023771 0 49
leadingEdge
1: MS4A1,CD37,CXCR4,BANK1,CD79B,TNFRSF13C,...
2: RPS29,RPS6,RPL32,RPS3A,RPL3,RPS25,...
3: HLA-DRA,HLA-DQB1,HLA-DPB1,HLA-DRB1,HLA-DPA1,HLA-DMA,...
4: CCL5,IL32,GZMH,CD3D,CD2,LYAR,...
5: S100A9,S100A6,LYZ,S100A8,TYROBP,FCN1,...
6: S100A4,FCER1G,IFITM3,S100A11,AIF1,PSAP,...
7: ITGB2,HCST,NKG7,GNLY,MYO1F,CST7,...
$`5`
pathway pval padj ES NES nMoreExtreme size
1: NK cell 0.0001380072 0.0004555809 0.9389534 1.873438 0 49
2: CD4 T cell 0.0001374948 0.0004555809 0.9146247 1.828931 0 50
3: CD8 T cell 0.0001518603 0.0004555809 0.9579395 1.658312 0 18
4: Plasma cell 0.0204893868 0.0263434973 0.8141737 1.475513 138 24
5: cDC 0.0055072464 0.0082608696 -0.8728717 -1.658161 18 17
6: ncMono 0.0003628447 0.0006531205 -0.8943373 -2.039555 0 49
7: cMono 0.0003610108 0.0006531205 -0.9056652 -2.049700 0 47
leadingEdge
1: NKG7,GNLY,CST7,GZMA,CTSW,GZMM,...
2: IL7R,RPS29,RPS3,RPL3,RPS6,RPL34,...
3: CCL5,IL32,GZMH,CD3D,LYAR,CD8A,...
4: RPL36AL,FKBP11,PPIB,SUB1,PEBP1,CYCS,...
5: HLA-DRA,HLA-DMA,HLA-DRB1,BASP1,HLA-DQB1,HLA-DRB5,...
6: FCER1G,IFITM3,FTH1,AIF1,COTL1,LST1,...
7: S100A9,S100A8,LYZ,TYROBP,FCN1,VCAN,...
$`6`
pathway pval padj ES NES nMoreExtreme size
1: NK cell 0.0001635590 0.0007672634 0.9817766 2.014939 0 49
2: CD8 T cell 0.0053735483 0.0080603224 0.8999420 1.595264 30 18
3: Plasma cell 0.0399521940 0.0513671065 0.7995328 1.477481 233 24
4: cDC 0.0021246459 0.0038243626 -0.9062204 -1.672658 8 17
5: ncMono 0.0010288066 0.0023148148 -0.7814587 -1.741362 3 49
6: B cell 0.0002557545 0.0007672634 -0.8690717 -1.925379 0 47
7: cMono 0.0002557545 0.0007672634 -0.8771297 -1.943231 0 47
leadingEdge
1: GNLY,NKG7,FGFBP2,CST7,CTSW,PRF1,...
2: CCL5,GZMH,IL32,LYAR,CD2,LINC01871,...
3: PPIB,FKBP11,SDF2L1,SPCS2,MANF,SUB1,...
4: HLA-DRA,HLA-DRB1,HLA-DQB1,HLA-DPA1,HLA-DMA,HLA-DRB5,...
5: COTL1,FTH1,SAT1,AIF1,LST1,IFITM3,...
6: MS4A1,CD79B,BANK1,CD37,TNFRSF13C,LINC00926,...
7: S100A9,S100A8,LYZ,FCN1,TKT,VCAN,...
$`7`
pathway pval padj ES NES nMoreExtreme size
1: NK cell 0.0737128601 0.132683148 -0.6699334 -1.262411 649 49
2: cDC 0.0027180336 0.006115576 -0.8740506 -1.462567 22 17
3: B cell 0.0007953642 0.002386092 -0.7912947 -1.484294 6 47
4: cMono 0.0004544938 0.002045222 -0.8032182 -1.506660 3 47
5: CD4 T cell 0.0001132375 0.001019137 -0.8955614 -1.691480 0 50
leadingEdge
1: ITGB2,NKG7,MYO1F,GNLY,IFITM1,JAK1,...
2: HLA-DRB1,HLA-DRA,HLA-DPB1,HLA-DPA1,HLA-DQB1,HLA-DMA,...
3: RPS23,RPL18A,CD52,RPL12,FAU,RPL23A,...
4: JUND,S100A6,NFKBIA,TYROBP,LYZ,TKT,...
5: RPL34,RPS3,RPL32,RPL13,RPL3,EEF1A1,...
$`8`
pathway pval padj ES NES nMoreExtreme size
1: ncMono 0.0001044605 0.000471797 0.9270402 1.646546 0 49
2: cMono 0.0001048438 0.000471797 0.9124335 1.616877 0 47
3: cDC 0.0026725540 0.004008831 0.9056296 1.478728 22 17
4: NK cell 0.0023310023 0.004008831 -0.7132923 -1.887183 0 49
5: CD8 T cell 0.0007446016 0.002233805 -0.9433388 -2.029006 0 18
6: B cell 0.0021551724 0.004008831 -0.8488506 -2.244658 0 47
leadingEdge
1: AIF1,S100A11,PSAP,LST1,SERPINA1,FCER1G,...
2: LYZ,FCN1,TYROBP,S100A9,S100A8,S100A6,...
3: HLA-DRA,HLA-DRB1,HLA-DPA1,HLA-DPB1,HLA-DQB1,HLA-DMA,...
4: NKG7,CST7,GZMM,CTSW,GZMA,CD247,...
5: CCL5,IL32,CD3D,GZMH,CD2,CD8A,...
6: CXCR4,MS4A1,TNFRSF13C,CD79B,BANK1,LINC00926,...
$`9`
pathway pval padj ES NES nMoreExtreme size
1: ncMono 0.0001061008 0.0009549072 0.9743962 1.718131 0 49
2: cDC 0.0276117619 0.0414176428 0.8495606 1.373365 230 17
3: B cell 0.0046801872 0.0084243370 -0.6804063 -1.722364 2 47
4: NK cell 0.0017331023 0.0039823009 -0.7548175 -1.902989 0 49
5: CD8 T cell 0.0006439150 0.0028976175 -0.9315008 -1.982888 0 18
6: CD4 T cell 0.0017699115 0.0039823009 -0.7961247 -2.021346 0 50
leadingEdge
1: LST1,AIF1,FCGR3A,COTL1,FCER1G,SERPINA1,...
2: HLA-DPA1,HLA-DRA,HLA-DRB1,HLA-DPB1,HLA-DRB5,MTMR14,...
3: CXCR4,MS4A1,BANK1,TNFRSF13C,LINC00926,RPL13A,...
4: NKG7,GNLY,CST7,CTSW,GZMA,GZMM,...
5: CCL5,IL32,GZMH,CD3D,CD2,KLRG1,...
6: LDHB,IL7R,RPL3,RPS27A,MGAT4A,RPL13,...Selecting top significant overlap per cluster, we can now rename the clusters according to the predicted labels. OBS! Be aware that if you have some clusters that have non-significant p-values for all the gene sets, the cluster label will not be very reliable. Also, the gene sets you are using may not cover all the celltypes you have in your dataset and hence predictions may just be the most similar celltype. Also, some of the clusters have very similar p-values to multiple celltypes, for instance the ncMono and cMono celltypes are equally good for some clusters.
new.cluster.ids <- unlist(lapply(res, function(x) {
as.data.frame(x)[1, 1]
}))
alldata@colData$ref_gsea <- new.cluster.ids[as.character(alldata@colData$louvain_SNNk15)]
wrap_plots(
plotReducedDim(alldata, dimred = "UMAP", colour_by = "louvain_SNNk15"),
plotReducedDim(alldata, dimred = "UMAP", colour_by = "ref_gsea"),
ncol = 2
)
Compare the results with the other celltype prediction methods in the ctrl_13 sample.
ctrl.sce@colData$ref_gsea <- alldata@colData$ref_gsea[alldata@colData$sample == "ctrl.13"]
wrap_plots(
plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "ref_gsea"),
plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "scmap_cell"),
plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "scpred_prediction"),
ncol = 3
)
7.2 With annotated gene sets
We have downloaded the celltype gene lists from http://bio-bigdata.hrbmu.edu.cn/CellMarker/CellMarker_download.html and converted the excel file to a csv for you. Read in the gene lists and do some filtering.
path_file <- file.path("data/human_cell_markers.txt")
if (!file.exists(path_file)) download.file(file.path(path_data, "human_cell_markers.txt"), destfile = path_file)markers <- read.delim("data/human_cell_markers.txt")
markers <- markers[markers$speciesType == "Human", ]
markers <- markers[markers$cancerType == "Normal", ]
# Filter by tissue (to reduce computational time and have tissue-specific classification)
# sort(unique(markers$tissueType))
# grep("blood",unique(markers$tissueType),value = T)
# markers <- markers [ markers$tissueType %in% c("Blood","Venous blood",
# "Serum","Plasma",
# "Spleen","Bone marrow","Lymph node"), ]
# remove strange characters etc.
celltype_list <- lapply(unique(markers$cellName), function(x) {
x <- paste(markers$geneSymbol[markers$cellName == x], sep = ",")
x <- gsub("[[]|[]]| |-", ",", x)
x <- unlist(strsplit(x, split = ","))
x <- unique(x[!x %in% c("", "NA", "family")])
x <- casefold(x, upper = T)
})
names(celltype_list) <- unique(markers$cellName)
# celltype_list <- lapply(celltype_list , function(x) {x[1:min(length(x),50)]} )
celltype_list <- celltype_list[unlist(lapply(celltype_list, length)) < 100]
celltype_list <- celltype_list[unlist(lapply(celltype_list, length)) > 5]# run fgsea for each of the clusters in the list
res <- lapply(DGE_list, function(x) {
x$logFC <- rowSums(as.matrix(x[, grep("logFC", colnames(x))]))
gene_rank <- setNames(x$logFC, rownames(x))
fgseaRes <- fgsea(pathways = celltype_list, stats = gene_rank, nperm = 10000)
return(fgseaRes)
})
names(res) <- names(DGE_list)
# You can filter and resort the table based on ES, NES or pvalue
res <- lapply(res, function(x) {
x[x$pval < 0.01, ]
})
res <- lapply(res, function(x) {
x[x$size > 5, ]
})
res <- lapply(res, function(x) {
x[order(x$NES, decreasing = T), ]
})
# show top 3 for each cluster.
lapply(res, head, 3)$`1`
pathway pval padj ES NES
1: Neutrophil 0.0001045697 0.01027098 0.9075102 1.678442
2: CD1C+_B dendritic cell 0.0001092657 0.01027098 0.9190331 1.658198
3: Adipose-derived stem cell 0.0005841121 0.02745327 0.8796298 1.518646
nMoreExtreme size leadingEdge
1: 0 82 S100A8,S100A9,S100A12,MNDA,CD14,S100A11,...
2: 0 54 S100A8,S100A9,LYZ,S100A12,VCAN,FCN1,...
3: 4 32 CD14,CD44,ANPEP,KLF4,PECAM1,ICAM1,...
$`10`
pathway pval padj ES NES
1: CD1C+_B dendritic cell 0.003496503 0.07571486 0.8559383 2.050422
2: Monocyte derived dendritic cell 0.003690037 0.07571486 0.8907285 1.867305
3: DCLK1+ progenitor cell 0.004149378 0.07571486 0.7092887 1.782269
nMoreExtreme size leadingEdge
1: 0 54 S100A9,S100A8,LYZ,VCAN,S100A12,CD163,...
2: 1 18 S100A9,S100A8,CST3,CD14,SIRPA,ITGAM
3: 0 68 VCAN,IFI6,BASP1,GBP1,TPST1,CYBRD1,...
$`2`
pathway pval padj ES NES
1: Follicular B cell 0.007232084 0.07085634 0.8625481 1.559597
2: Lake et al.Science.Ex5 0.009654204 0.07085634 0.9605640 1.444209
3: Pyramidal cell 0.004878049 0.06113821 -0.9678714 -1.526969
nMoreExtreme size leadingEdge
1: 43 22 MS4A1,CD69,CD22,FCER2,CD40,PAX5,...
2: 54 6 RCSD1,RNF182
3: 20 6 NRGN,CD3E
$`3`
pathway pval padj ES NES nMoreExtreme
1: Naive CD8+ T cell 0.0001491647 0.009256524 0.8683195 1.928582 0
2: Naive CD4+ T cell 0.0001632387 0.009256524 0.9276762 1.794550 0
3: CD4+ T cell 0.0001680108 0.009256524 0.9353639 1.734658 0
size leadingEdge
1: 91 LDHB,PIK3IP1,NPM1,TCF7,NOSIP,RPS8,...
2: 34 IL7R,EEF1B2,TCF7,NOSIP,RPS5,MAL,...
3: 25 IL7R,LTB,CD3E,CD3D,CD3G,CD2,...
$`4`
pathway pval padj ES NES
1: Follicular B cell 0.007506255 0.07055880 0.8583276 1.573944
2: Endothelial progenitor cell 0.009442653 0.07718343 -0.9271852 -1.526440
3: CD4+CD25+ regulatory T cell 0.002736602 0.03674866 -0.9710324 -1.526959
nMoreExtreme size leadingEdge
1: 44 22 MS4A1,CD69,CD22,CD40,FCER2,PAX5,...
2: 40 8 PECAM1,PTPRC
3: 11 6 CD3E,CD3D,CD3G,CD4,PTPRC
$`5`
pathway pval padj ES NES
1: CD4+ cytotoxic T cell 0.0001291990 0.009386858 0.8796627 1.875976
2: Natural killer cell 0.0001296849 0.009386858 0.8005797 1.702411
3: CD8+ T cell 0.0001497903 0.009386858 0.9520965 1.672226
nMoreExtreme size leadingEdge
1: 0 86 NKG7,CCL5,GNLY,GZMH,CST7,GZMA,...
2: 0 84 NKG7,GNLY,CD3D,GZMA,CD3E,CD3G,...
3: 0 19 NKG7,CD3D,CD3E,CD3G,CD8A,GZMK,...
$`6`
pathway pval padj ES
1: CD4+ cytotoxic T cell 0.0001536098 0.009706733 0.9416682
2: Effector CD8+ memory T (Tem) cell 0.0001548947 0.009706733 0.9037846
3: Natural killer cell 0.0001540357 0.009706733 0.8487983
NES nMoreExtreme size leadingEdge
1: 2.070368 0 86 GNLY,NKG7,GZMB,CCL5,FGFBP2,CST7,...
2: 1.967557 0 79 GNLY,GZMB,FGFBP2,GZMH,KLRD1,ARL4C,...
3: 1.862851 0 84 GNLY,NKG7,GZMB,GZMA,KLRD1,CD247,...
$`7`
pathway pval padj ES NES
1: Megakaryocyte 0.0014306152 0.09059150 0.8651089 1.821106
2: Radial glial cell 0.0014456090 0.09059150 -0.9730759 -1.484748
3: Morula cell (Blastomere) 0.0001082485 0.02035073 -0.7508912 -1.493162
nMoreExtreme size leadingEdge
1: 1 26 PPBP,PF4,GP9,ITGA2B,CD9,RASGRP2,...
2: 11 7 VIM
3: 0 88 RPL34,RPL11,RPL23,RPL39,RPL7,RPL22,...
$`8`
pathway pval padj ES NES
1: CD1C+_B dendritic cell 0.0001035947 0.006690153 0.8760942 1.561950
2: Neutrophil 0.0001013377 0.006690153 0.8542972 1.550760
3: Stromal cell 0.0001067578 0.006690153 0.8710348 1.522084
nMoreExtreme size leadingEdge
1: 0 54 LYZ,FCN1,S100A9,S100A8,CSTA,MNDA,...
2: 0 82 S100A9,S100A11,LST1,S100A8,MNDA,CD14,...
3: 0 38 VIM,CD44,TIMP2,TIMP1,PECAM1,ICAM1,...
$`9`
pathway pval padj ES NES
1: Mesenchymal cell 0.0002090083 0.01964678 0.8438168 1.502770
2: Hemangioblast 0.0001342462 0.01964678 0.9905418 1.475874
3: Endothelial progenitor cell 0.0010739697 0.06730210 0.9756572 1.453697
nMoreExtreme size leadingEdge
1: 1 61 COTL1,S100A4,CTSC,HES4,ZEB2,VIM,...
2: 0 8 PECAM1,CD34
3: 7 8 PECAM1,PTPRC#CT_GSEA8:
new.cluster.ids <- unlist(lapply(res, function(x) {
as.data.frame(x)[1, 1]
}))
alldata@colData$cellmarker_gsea <- new.cluster.ids[as.character(alldata@colData$louvain_SNNk15)]
wrap_plots(
plotReducedDim(alldata, dimred = "UMAP", colour_by = "cellmarker_gsea"),
plotReducedDim(alldata, dimred = "UMAP", colour_by = "ref_gsea"),
ncol = 2
)
Discuss
Do you think that the methods overlap well? Where do you see the most inconsistencies?
In this case we do not have any ground truth, and we cannot say which method performs best. You should keep in mind, that any celltype classification method is just a prediction, and you still need to use your common sense and knowledge of the biological system to judge if the results make sense.
Finally, lets save the data with predictions.
saveRDS(ctrl.sce, "data/covid/results/bioc_covid_qc_dr_int_cl_ct-ctrl13.rds")8 Session info
Click here
sessionInfo()R version 4.3.0 (2023-04-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: Etc/UTC
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] fgsea_1.28.0 caret_6.0-94
[3] lattice_0.21-8 SeuratObject_4.1.3
[5] Seurat_4.3.0 scmap_1.24.0
[7] scPred_1.9.2 pheatmap_1.0.12
[9] patchwork_1.1.2 dplyr_1.1.2
[11] scran_1.30.0 scater_1.30.1
[13] ggplot2_3.4.2 scuttle_1.12.0
[15] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
[17] Biobase_2.62.0 GenomicRanges_1.54.1
[19] GenomeInfoDb_1.38.5 IRanges_2.36.0
[21] S4Vectors_0.40.2 BiocGenerics_0.48.1
[23] MatrixGenerics_1.14.0 matrixStats_1.0.0
loaded via a namespace (and not attached):
[1] spatstat.sparse_3.0-1 bitops_1.0-7
[3] lubridate_1.9.2 httr_1.4.6
[5] RColorBrewer_1.1-3 tools_4.3.0
[7] sctransform_0.3.5 utf8_1.2.3
[9] R6_2.5.1 lazyeval_0.2.2
[11] uwot_0.1.14 withr_2.5.0
[13] sp_1.6-1 gridExtra_2.3
[15] progressr_0.13.0 cli_3.6.1
[17] spatstat.explore_3.2-1 labeling_0.4.2
[19] spatstat.data_3.0-1 randomForest_4.7-1.1
[21] proxy_0.4-27 ggridges_0.5.4
[23] pbapply_1.7-0 harmony_1.2.0
[25] parallelly_1.36.0 limma_3.58.1
[27] rstudioapi_0.14 FNN_1.1.3.2
[29] generics_0.1.3 ica_1.0-3
[31] spatstat.random_3.1-5 Matrix_1.5-4
[33] ggbeeswarm_0.7.2 fansi_1.0.4
[35] abind_1.4-5 lifecycle_1.0.3
[37] yaml_2.3.7 edgeR_4.0.7
[39] recipes_1.0.6 SparseArray_1.2.3
[41] Rtsne_0.16 grid_4.3.0
[43] promises_1.2.0.1 dqrng_0.3.0
[45] crayon_1.5.2 miniUI_0.1.1.1
[47] beachmat_2.18.0 cowplot_1.1.1
[49] pillar_1.9.0 knitr_1.43
[51] metapod_1.10.1 future.apply_1.11.0
[53] codetools_0.2-19 fastmatch_1.1-3
[55] leiden_0.4.3 googleVis_0.7.1
[57] glue_1.6.2 data.table_1.14.8
[59] vctrs_0.6.2 png_0.1-8
[61] gtable_0.3.3 kernlab_0.9-32
[63] gower_1.0.1 xfun_0.39
[65] S4Arrays_1.2.0 mime_0.12
[67] prodlim_2023.03.31 survival_3.5-5
[69] timeDate_4022.108 iterators_1.0.14
[71] hardhat_1.3.0 lava_1.7.2.1
[73] statmod_1.5.0 bluster_1.12.0
[75] ellipsis_0.3.2 fitdistrplus_1.1-11
[77] ROCR_1.0-11 ipred_0.9-14
[79] nlme_3.1-162 RcppAnnoy_0.0.20
[81] irlba_2.3.5.1 vipor_0.4.5
[83] KernSmooth_2.23-20 rpart_4.1.19
[85] colorspace_2.1-0 nnet_7.3-18
[87] tidyselect_1.2.0 compiler_4.3.0
[89] BiocNeighbors_1.20.2 DelayedArray_0.28.0
[91] plotly_4.10.2 scales_1.2.1
[93] lmtest_0.9-40 stringr_1.5.0
[95] digest_0.6.31 goftest_1.2-3
[97] spatstat.utils_3.0-3 rmarkdown_2.22
[99] XVector_0.42.0 RhpcBLASctl_0.23-42
[101] htmltools_0.5.5 pkgconfig_2.0.3
[103] sparseMatrixStats_1.14.0 fastmap_1.1.1
[105] rlang_1.1.1 htmlwidgets_1.6.2
[107] shiny_1.7.4 DelayedMatrixStats_1.24.0
[109] farver_2.1.1 zoo_1.8-12
[111] jsonlite_1.8.5 BiocParallel_1.36.0
[113] ModelMetrics_1.2.2.2 BiocSingular_1.18.0
[115] RCurl_1.98-1.12 magrittr_2.0.3
[117] GenomeInfoDbData_1.2.11 munsell_0.5.0
[119] Rcpp_1.0.10 viridis_0.6.3
[121] reticulate_1.30 stringi_1.7.12
[123] pROC_1.18.2 zlibbioc_1.48.0
[125] MASS_7.3-58.4 plyr_1.8.8
[127] parallel_4.3.0 listenv_0.9.0
[129] ggrepel_0.9.3 deldir_1.0-9
[131] splines_4.3.0 tensor_1.5
[133] locfit_1.5-9.8 igraph_1.4.3
[135] spatstat.geom_3.2-1 reshape2_1.4.4
[137] ScaledMatrix_1.10.0 evaluate_0.21
[139] foreach_1.5.2 httpuv_1.6.11
[141] RANN_2.6.1 tidyr_1.3.0
[143] purrr_1.0.1 polyclip_1.10-4
[145] future_1.32.0 scattermore_1.2
[147] rsvd_1.0.5 xtable_1.8-4
[149] e1071_1.7-13 later_1.3.1
[151] viridisLite_0.4.2 class_7.3-21
[153] tibble_3.2.1 beeswarm_0.4.0
[155] cluster_2.1.4 timechange_0.2.0
[157] globals_0.16.2