suppressPackageStartupMessages({
library(scater)
library(scran)
library(dplyr)
library(patchwork)
library(ggplot2)
library(pheatmap)
library(scPred)
library(scmap)
library(SingleR)
})
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.
Ideally celltype predictions should be run on each sample separately and not using the integrated data. In this case 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
<- TRUE
fetch_data
# url for source and intermediate data
<- "https://nextcloud.dc.scilifelab.se/public.php/webdav"
path_data <- "-u zbC5fr2LbEZ9rSE:scRNAseq2025"
curl_upass
<- "data/covid/results/bioc_covid_qc_dr_int_cl.rds"
path_file 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/bioc_covid_qc_dr_int_cl.rds"), destfile = path_file, method = "curl", extra = curl_upass)
<- readRDS(path_file) alldata
Let’s read in the saved Covid-19 data object from the clustering step.
<- alldata[, alldata$sample == "ctrl.13"]
ctrl.sce
# remove all old dimensionality reductions as they will mess up the analysis further down
reducedDims(ctrl.sce) <- NULL
2 Reference data
Load the reference dataset with annotated labels that is provided by the scPred
package, it is a subsampled set of cells from human PBMCs.
<- scPred::pbmc_1
reference reference
An object of class Seurat
32838 features across 3500 samples within 1 assay
Active assay: RNA (32838 features, 0 variable features)
2 layers present: counts, data
Convert to a SCE object.
<- Seurat::as.SingleCellExperiment(reference) ref.sce
Rerun analysis pipeline. Run normalization, feature selection and dimensionality reduction
# Normalize
<- computeSumFactors(ref.sce)
ref.sce <- logNormCounts(ref.sce)
ref.sce
# Variable genes
<- modelGeneVar(ref.sce, method = "loess")
var.out <- getTopHVGs(var.out, n = 1000)
hvg.ref
# Dim reduction
<- runPCA(ref.sce,
ref.sce exprs_values = "logcounts", scale = T,
ncomponents = 30, subset_row = hvg.ref
)<- runUMAP(ref.sce, dimred = "PCA") ref.sce
Run all steps of the analysis for the ctrl sample as well. Use the clustering from the integration lab with resolution 0.5.
# Normalize
<- computeSumFactors(ctrl.sce)
ctrl.sce <- logNormCounts(ctrl.sce)
ctrl.sce
# Variable genes
<- modelGeneVar(ctrl.sce, method = "loess")
var.out <- getTopHVGs(var.out, n = 1000)
hvg.ctrl
# Dim reduction
<- runPCA(ctrl.sce, exprs_values = "logcounts", scale = T, ncomponents = 30, subset_row = hvg.ctrl)
ctrl.sce <- runUMAP(ctrl.sce, dimred = "PCA") ctrl.sce
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
$cell_type1 <- ref.sce$cell_type
ref.sce# create a rowData slot with feature_symbol
<- data.frame(feature_symbol = rownames(ref.sce))
rd rownames(rd) <- rownames(ref.sce)
rowData(ref.sce) <- rd
# same for the ctrl dataset
# create a rowData slot with feature_symbol
<- data.frame(feature_symbol = rownames(ctrl.sce))
rd rownames(rd) <- rownames(ctrl.sce)
rowData(ctrl.sce) <- rd
Then 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))
<- selectFeatures(ctrl.sce, suppress_plot = TRUE)
ctrl.sce
counts(ref.sce) <- as.matrix(counts(ref.sce))
logcounts(ref.sce) <- as.matrix(logcounts(ref.sce))
<- selectFeatures(ref.sce, suppress_plot = TRUE) ref.sce
Then we need to index the reference dataset by cluster, default is the clusters in cell_type1
.
<- indexCluster(ref.sce) ref.sce
Now we project the Covid-19 dataset onto that index.
<- scmapCluster(
project_cluster 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
74 109 113 31 203 138
NK cell pDC Plasma cell unassigned
279 2 1 152
Then add the predictions to metadata and plot UMAP.
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.
<- indexCell(ref.sce) ref.sce
Again we need to index the reference dataset.
<- scmapCell(
project_cell 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.
<- colData(ref.sce)$cell_type1[project_cell$ref[[1]][1, ]]
cell_type_pred table(cell_type_pred)
cell_type_pred
B cell CD4 T cell CD8 T cell cDC cMono ncMono
97 181 249 49 170 173
NK cell pDC Plasma cell
180 2 1
Then add the predictions to metadata and plot umap.
# add in predictions
$scmap_cell <- cell_type_pred
ctrl.sce
plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "scmap_cell")
Plot both:
5 SinlgeR
SingleR is performs unbiased cell type recognition from single-cell RNA sequencing data, by leveraging reference transcriptomic datasets of pure cell types to infer the cell of origin of each single cell independently.
There are multiple datasets included in the celldex
package that can be used for celltype prediction, here we will test two different ones, the DatabaseImmuneCellExpressionData
and the HumanPrimaryCellAtlasData
. In addition we will use the same reference dataset that we used for label transfer above but using SingleR instead.
5.1 Immune cell reference
= celldex::DatabaseImmuneCellExpressionData()
immune <- SingleR(test = ctrl.sce, ref = immune, assay.type.test=1,
singler.immune labels = immune$label.main)
head(singler.immune)
DataFrame with 6 rows and 4 columns
scores labels delta.next
<matrix> <character> <numeric>
AGGTCATGTGCGAACA-13 0.0679680:0.0760411:0.186694:... T cells, CD4+ 0.1375513
CCTATCGGTCCCTCAT-13 0.0027079:0.0960641:0.386088:... NK cells 0.1490740
TCCTCCCTCGTTCATT-13 0.0361115:0.1067465:0.394579:... NK cells 0.1220681
TACGGTATCGGATTAC-13 -0.0131813:0.0717678:0.283882:... NK cells 0.0620657
AATAGAGAGTTCGGTT-13 0.0841091:0.1367749:0.273738:... T cells, CD4+ 0.0660296
TGACGCGTCGCAATTG-13 -0.0160540:0.0912025:0.396588:... NK cells 0.2057787
pruned.labels
<character>
AGGTCATGTGCGAACA-13 T cells, CD4+
CCTATCGGTCCCTCAT-13 NK cells
TCCTCCCTCGTTCATT-13 NK cells
TACGGTATCGGATTAC-13 NK cells
AATAGAGAGTTCGGTT-13 T cells, CD4+
TGACGCGTCGCAATTG-13 NK cells
5.2 HPCA reference
<- HumanPrimaryCellAtlasData()
hpca <- SingleR(test = ctrl.sce, ref = hpca, assay.type.test=1,
singler.hpca labels = hpca$label.main)
head(singler.hpca)
DataFrame with 6 rows and 4 columns
scores labels delta.next
<matrix> <character> <numeric>
AGGTCATGTGCGAACA-13 0.141378:0.310009:0.275987:... T_cells 0.4918992
CCTATCGGTCCCTCAT-13 0.145926:0.300045:0.277827:... NK_cell 0.3241970
TCCTCCCTCGTTCATT-13 0.132119:0.311754:0.274127:... NK_cell 0.0640608
TACGGTATCGGATTAC-13 0.125120:0.283118:0.250322:... T_cells 0.1545913
AATAGAGAGTTCGGTT-13 0.191441:0.374422:0.329988:... T_cells 0.5063484
TGACGCGTCGCAATTG-13 0.121131:0.279338:0.249160:... NK_cell 0.3018100
pruned.labels
<character>
AGGTCATGTGCGAACA-13 T_cells
CCTATCGGTCCCTCAT-13 NK_cell
TCCTCCCTCGTTCATT-13 NK_cell
TACGGTATCGGATTAC-13 T_cells
AATAGAGAGTTCGGTT-13 T_cells
TGACGCGTCGCAATTG-13 NK_cell
5.3 With own reference data
<- SingleR(test=ctrl.sce, ref=ref.sce, labels=ref.sce$cell_type, de.method="wilcox")
singler.ref head(singler.ref)
DataFrame with 6 rows and 4 columns
scores labels delta.next
<matrix> <character> <numeric>
AGGTCATGTGCGAACA-13 0.741719:0.840093:0.805977:... CD4 T cell 0.0423204
CCTATCGGTCCCTCAT-13 0.649491:0.736753:0.815987:... NK cell 0.0451715
TCCTCCCTCGTTCATT-13 0.669603:0.731356:0.823308:... NK cell 0.0865526
TACGGTATCGGATTAC-13 0.708827:0.776244:0.808044:... CD8 T cell 0.0905218
AATAGAGAGTTCGGTT-13 0.729010:0.847462:0.816299:... CD4 T cell 0.0409309
TGACGCGTCGCAATTG-13 0.625165:0.681015:0.794166:... NK cell 0.0597978
pruned.labels
<character>
AGGTCATGTGCGAACA-13 CD4 T cell
CCTATCGGTCCCTCAT-13 NK cell
TCCTCCCTCGTTCATT-13 NK cell
TACGGTATCGGATTAC-13 CD8 T cell
AATAGAGAGTTCGGTT-13 CD4 T cell
TGACGCGTCGCAATTG-13 NK cell
Compare results:
$singler.immune = singler.immune$pruned.labels
ctrl.sce$singler.hpca = singler.hpca$pruned.labels
ctrl.sce$singler.ref = singler.ref$pruned.labels
ctrl.sce
wrap_plots(
plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "singler.immune"),
plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "singler.hpca"),
plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "singler.ref"),
ncol = 3
)
6 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.
table(ctrl.sce$scmap_cell, ctrl.sce$singler.hpca)
B_cell Macrophage Monocyte NK_cell Platelets T_cells
B cell 96 0 0 0 0 0
CD4 T cell 0 0 0 2 1 177
CD8 T cell 0 0 0 136 0 105
cDC 1 0 46 0 0 0
cMono 0 1 157 0 0 0
ncMono 0 1 170 0 0 0
NK cell 0 0 0 168 0 7
pDC 1 0 0 0 0 0
Plasma cell 1 0 0 0 0 0
Or plot onto umap:
wrap_plots(
plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "scmap_cluster"),
plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "scmap_cell"),
plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "singler.immune"),
plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "singler.hpca"),
plotReducedDim(ctrl.sce, dimred = "UMAP", colour_by = "singler.ref"),
ncol = 3
)
As you can see, the methods using the same reference all have similar results. While for instance singleR with different references give quite different predictions. This really shows that a relevant reference is the key in having reliable celltype predictions rather than the method used.
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
<- scran::findMarkers(
DGE_list x = alldata,
groups = as.character(alldata$leiden_k20),
pval.type = "all",
min.prop = 0
)
# Compute differential gene expression in reference dataset (that has cell annotation)
<- scran::findMarkers(
ref_DGE x = ref.sce,
groups = as.character(ref.sce$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.
<- lapply(ref_DGE, function(x) {
ref_list $logFC <- rowSums(as.matrix(x[, grep("logFC", colnames(x))]))
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
<- lapply(DGE_list, function(x) {
res $logFC <- rowSums(as.matrix(x[, grep("logFC", colnames(x))]))
x<- setNames(x$logFC, rownames(x))
gene_rank <- fgsea(pathways = ref_list, stats = gene_rank, nperm = 10000)
fgseaRes return(fgseaRes)
})names(res) <- names(DGE_list)
# You can filter and resort the table based on ES, NES or pvalue
<- lapply(res, function(x) {
res $pval < 0.1, ]
x[x
})<- lapply(res, function(x) {
res $size > 2, ]
x[x
})<- lapply(res, function(x) {
res order(x$NES, decreasing = T), ]
x[
}) res
$`1`
pathway pval padj ES NES nMoreExtreme
<char> <num> <num> <num> <num> <num>
1: cMono 0.0051546392 0.0067669173 0.9222327 2.294590 0
2: ncMono 0.0052631579 0.0067669173 0.6362900 1.573839 0
3: B cell 0.0005097879 0.0011470228 -0.7928592 -1.374113 4
4: CD8 T cell 0.0047966632 0.0067669173 -0.8495939 -1.386710 45
5: Plasma cell 0.0004142931 0.0011470228 -0.8672090 -1.445820 3
6: NK cell 0.0001019160 0.0004586221 -0.8480240 -1.473063 0
7: CD4 T cell 0.0001018537 0.0004586221 -0.9096249 -1.582175 0
size leadingEdge
<int> <list>
1: 47 S100A8, ....
2: 49 AIF1, S1....
3: 47 RPS5, RP....
4: 18 CCL5, IL....
5: 24 RPL36AL,....
6: 49 B2M, NKG....
7: 50 RPL3, RP....
$`10`
pathway pval padj ES NES nMoreExtreme size
<char> <num> <num> <num> <num> <num> <int>
1: B cell 0.0095038015 0.0171068427 -0.7417620 -1.227837 94 47
2: CD8 T cell 0.0339512392 0.0436515932 -0.7979091 -1.271178 336 18
3: cMono 0.0021008403 0.0047268908 -0.7724947 -1.278708 20 47
4: cDC 0.0309475806 0.0436515932 -0.8059082 -1.279949 306 17
5: NK cell 0.0002001001 0.0006003002 -0.8092281 -1.340929 1 49
6: ncMono 0.0001000500 0.0004502251 -0.8579812 -1.421715 0 49
7: CD4 T cell 0.0001000400 0.0004502251 -0.8976778 -1.488726 0 50
leadingEdge
<list>
1: FAU, RPL....
2: CCL5, IL....
3: JUND, S1....
4: HLA-DRB1....
5: B2M, HCS....
6: FTH1, S1....
7: RPS29, R....
$`11`
pathway pval padj ES NES nMoreExtreme
<char> <num> <num> <num> <num> <num>
1: ncMono 0.0002187227 0.0004921260 0.9815952 2.117847 0
2: cDC 0.0086166744 0.0110785814 0.8857459 1.604149 36
3: Plasma cell 0.0017934003 0.0026901004 -0.8394435 -1.625605 9
4: CD8 T cell 0.0005265929 0.0009478673 -0.9084665 -1.675881 2
5: B cell 0.0001838235 0.0004921260 -0.8089817 -1.763016 0
6: NK cell 0.0001841621 0.0004921260 -0.8476554 -1.858533 0
7: CD4 T cell 0.0001844678 0.0004921260 -0.9165981 -2.015372 0
size leadingEdge
<int> <list>
1: 49 LST1, AI....
2: 17 HLA-DPA1....
3: 24 ISG20, R....
4: 18 CCL5, IL....
5: 47 CXCR4, R....
6: 49 NKG7, GN....
7: 50 RPL31, R....
$`12`
pathway pval padj ES NES nMoreExtreme
<char> <num> <num> <num> <num> <num>
1: CD4 T cell 0.0001001101 0.0009009911 0.9518024 1.571553 0
2: Plasma cell 0.0023352625 0.0070057874 0.8574356 1.376587 22
3: B cell 0.0016033671 0.0070057874 0.7964622 1.311655 15
4: CD8 T cell 0.0339905525 0.0559006211 0.8136559 1.286572 330
5: NK cell 0.0395593390 0.0559006211 0.7258037 1.197098 394
6: ncMono 0.0588235294 0.0661764706 -0.4386277 -1.325709 0
7: cMono 0.0434782609 0.0559006211 -0.6808314 -2.152243 0
8: cDC 0.0033333333 0.0075000000 -0.8964181 -2.228026 0
size leadingEdge
<int> <list>
1: 50 IL7R, RP....
2: 24 RPL36AL,....
3: 47 RPS5, RP....
4: 18 IL32, CD....
5: 49 IFITM1, ....
6: 49 FCER1G, ....
7: 47 S100A8, ....
8: 17 HLA-DRA,....
$`13`
pathway pval padj ES NES nMoreExtreme size
<char> <num> <num> <num> <num> <num> <int>
1: cMono 0.0002861230 0.000516203 0.9237564 2.137192 0 47
2: ncMono 0.0002867795 0.000516203 0.8495391 1.979065 0 49
3: cDC 0.0554525795 0.062384152 -0.7508028 -1.380756 358 17
4: CD8 T cell 0.0010830883 0.001624633 -0.8666052 -1.605641 6 18
5: Plasma cell 0.0013903909 0.001787645 -0.8330442 -1.618347 8 24
6: NK cell 0.0001534919 0.000461042 -0.8068779 -1.761721 0 49
7: B cell 0.0001536807 0.000461042 -0.8446545 -1.831407 0 47
8: CD4 T cell 0.0001536098 0.000461042 -0.9229947 -2.021512 0 50
leadingEdge
<list>
1: S100A8, ....
2: AIF1, S1....
3: HLA-DPB1....
4: IL32, GZ....
5: ISG20, R....
6: GNLY, NK....
7: RPS5, RP....
8: RPL14, R....
$`14`
pathway pval padj ES NES nMoreExtreme size
<char> <num> <num> <num> <num> <num> <int>
1: cMono 0.000100040 0.0004502251 0.9178528 1.452446 0 47
2: ncMono 0.000100050 0.0004502251 0.9130436 1.446765 0 49
3: cDC 0.015832741 0.0284989343 0.8411903 1.275312 155 17
4: CD4 T cell 0.009904952 0.0222861431 0.7618107 1.207520 98 50
5: pDC 0.052621048 0.0789315726 0.7272444 1.150820 525 47
6: CD8 T cell 0.007575758 0.0222861431 -0.9739951 -2.428217 0 18
leadingEdge
<list>
1: S100A9, ....
2: AIF1, S1....
3: HLA-DRA,....
4: RPS13, R....
5: CTSB, NP....
6: IL32, CC....
$`15`
pathway pval padj ES NES nMoreExtreme
<char> <num> <num> <num> <num> <num>
1: NK cell 0.0181818182 0.0181818182 0.8291322 2.166520 0
2: CD8 T cell 0.0067114094 0.0100671141 0.9883234 2.113793 0
3: pDC 0.0141794047 0.0159518302 -0.7349782 -1.219553 140
4: B cell 0.0025140788 0.0049250101 -0.7703560 -1.278256 24
5: cDC 0.0083282551 0.0107077566 -0.8413649 -1.333409 81
6: Plasma cell 0.0027361167 0.0049250101 -0.8330300 -1.340330 26
7: cMono 0.0001005632 0.0003016895 -0.8458496 -1.403523 0
8: CD4 T cell 0.0001004823 0.0003016895 -0.8841618 -1.470774 0
9: ncMono 0.0001005328 0.0003016895 -0.8851646 -1.470914 0
size leadingEdge
<int> <list>
1: 49 NKG7, GN....
2: 18 CCL5, GZ....
3: 47 NPC2, YP....
4: 47 RPS11, C....
5: 17 HLA-DRA,....
6: 24 SUB1, RP....
7: 47 S100A9, ....
8: 50 TPT1, RP....
9: 49 FTH1, CO....
$`2`
pathway pval padj ES NES nMoreExtreme size
<char> <num> <num> <num> <num> <num> <int>
1: B cell 0.0008802817 0.0019806338 0.9759274 2.286607 0 47
2: cDC 0.0017094017 0.0030769231 0.9364487 1.819993 2 17
3: CD4 T cell 0.0102803738 0.0132176235 0.6651708 1.567700 10 50
4: CD8 T cell 0.0047044632 0.0070566948 -0.8797950 -1.455393 38 18
5: cMono 0.0001127904 0.0003383713 -0.8465424 -1.546511 0 47
6: ncMono 0.0001121202 0.0003383713 -0.8885809 -1.631038 0 49
7: NK cell 0.0001121202 0.0003383713 -0.9126109 -1.675146 0 49
leadingEdge
<list>
1: MS4A1, C....
2: HLA-DRA,....
3: RPS6, RP....
4: CCL5, IL....
5: S100A6, ....
6: S100A4, ....
7: HCST, NK....
$`3`
pathway pval padj ES NES nMoreExtreme
<char> <num> <num> <num> <num> <num>
1: cMono 0.0002200220 0.0003960396 0.9539785 2.056921 0
2: ncMono 0.0002193463 0.0003960396 0.8821112 1.915328 0
3: Plasma cell 0.0029955947 0.0038514789 -0.8344716 -1.600177 16
4: CD8 T cell 0.0005242006 0.0007863009 -0.9104821 -1.656777 2
5: NK cell 0.0001837222 0.0003960396 -0.8350365 -1.809856 0
6: B cell 0.0001832509 0.0003960396 -0.8791693 -1.892953 0
7: CD4 T cell 0.0001828822 0.0003960396 -0.9284623 -2.023550 0
size leadingEdge
<int> <list>
1: 47 S100A9, ....
2: 49 AIF1, S1....
3: 24 ISG20, P....
4: 18 CCL5, IL....
5: 49 NKG7, GN....
6: 47 CXCR4, R....
7: 50 RPL3, RP....
$`4`
pathway pval padj ES NES nMoreExtreme size
<char> <num> <num> <num> <num> <num> <int>
1: CD4 T cell 0.0012106538 0.0021791768 0.9835702 2.454592 0 50
2: cDC 0.0018957346 0.0028436019 -0.8984204 -1.480202 15 17
3: NK cell 0.0003279047 0.0007377856 -0.8152480 -1.492383 2 49
4: pDC 0.0002192261 0.0006576784 -0.8331320 -1.519475 1 47
5: cMono 0.0001096131 0.0004932588 -0.9043158 -1.649300 0 47
6: ncMono 0.0001093016 0.0004932588 -0.9322838 -1.706628 0 49
leadingEdge
<list>
1: IL7R, LD....
2: HLA-DRA,....
3: NKG7, GN....
4: PLEK, PL....
5: S100A9, ....
6: FCER1G, ....
$`5`
pathway pval padj ES NES nMoreExtreme size
<char> <num> <num> <num> <num> <num> <int>
1: NK cell 0.0007974482 0.0017942584 0.9427049 2.171396 0 49
2: CD8 T cell 0.0005299417 0.0015898251 0.9757739 1.934063 0 18
3: CD4 T cell 0.0160384924 0.0258044554 0.6841476 1.584800 19 50
4: pDC 0.0828932262 0.1065770051 -0.6791309 -1.275989 721 47
5: cDC 0.0172029703 0.0258044554 -0.8349892 -1.400588 138 17
6: ncMono 0.0001143118 0.0005166475 -0.9090886 -1.714372 0 49
7: cMono 0.0001148106 0.0005166475 -0.9257589 -1.739368 0 47
leadingEdge
<list>
1: NKG7, GN....
2: CCL5, GZ....
3: IL7R, RP....
4: NPC2, CT....
5: HLA-DRA,....
6: FCER1G, ....
7: S100A9, ....
$`6`
pathway pval padj ES NES nMoreExtreme
<char> <num> <num> <num> <num> <num>
1: B cell 0.0001007557 0.0004534005 0.9362193 1.593408 0
2: CD4 T cell 0.0001005935 0.0004534005 0.9245421 1.577215 0
3: cDC 0.0005278159 0.0015834477 0.9262262 1.485196 4
4: Plasma cell 0.0065257924 0.0117464264 0.8261552 1.356850 62
5: pDC 0.0129974811 0.0149253731 0.7484645 1.273857 128
6: ncMono 0.0149253731 0.0149253731 -0.5308483 -1.572343 0
7: cMono 0.0129870130 0.0149253731 -0.6492382 -1.903165 0
8: CD8 T cell 0.0020449898 0.0046012270 -0.9208483 -2.178479 0
9: NK cell 0.0149253731 0.0149253731 -0.7908323 -2.342401 0
size leadingEdge
<int> <list>
1: 47 CD37, MS....
2: 50 RPS6, RP....
3: 17 HLA-DRA,....
4: 24 ISG20, R....
5: 47 IRF8, BC....
6: 49 S100A4, ....
7: 47 TYROBP, ....
8: 18 CCL5, IL....
9: 49 NKG7, HC....
$`7`
pathway pval padj ES NES nMoreExtreme
<char> <num> <num> <num> <num> <num>
1: NK cell 0.0001021450 0.0004596527 0.9546415 1.672175 0
2: CD4 T cell 0.0001020200 0.0004596527 0.8733240 1.531929 0
3: Plasma cell 0.0004268032 0.0012804097 0.8984626 1.506395 3
4: CD8 T cell 0.0010966115 0.0024673758 0.9003202 1.473335 9
5: cDC 0.0998937301 0.1284347958 -0.6073699 -1.379830 93
6: ncMono 0.0141509434 0.0212264151 -0.5229824 -1.472616 2
7: cMono 0.0045248869 0.0081447964 -0.7893952 -2.210234 0
size leadingEdge
<int> <list>
1: 49 NKG7, GN....
2: 50 RPL3, RP....
3: 24 PPIB, RP....
4: 18 CCL5, GZ....
5: 17 HLA-DRA,....
6: 49 COTL1, A....
7: 47 S100A9, ....
$`8`
pathway pval padj ES NES nMoreExtreme
<char> <num> <num> <num> <num> <num>
1: Plasma cell 0.0185744657 0.0278616986 -0.7754252 -1.341161 165
2: NK cell 0.0014034330 0.0031577243 -0.7580622 -1.409313 12
3: B cell 0.0009742368 0.0029227105 -0.7613088 -1.409610 8
4: cDC 0.0062471604 0.0112448887 -0.8528729 -1.423184 54
5: cMono 0.0004329942 0.0019484737 -0.7894528 -1.461720 3
6: CD4 T cell 0.0001079214 0.0009712929 -0.8908780 -1.658522 0
size leadingEdge
<int> <list>
1: 24 PPIB, RP....
2: 49 ITGB2, N....
3: 47 RPS23, R....
4: 17 HLA-DRB1....
5: 47 JUND, S1....
6: 50 RPL34, R....
$`9`
pathway pval padj ES NES nMoreExtreme size
<char> <num> <num> <num> <num> <num> <int>
1: NK cell 0.0010504202 0.0018907563 0.9859981 2.203758 0 49
2: CD8 T cell 0.0232052212 0.0298352844 0.8506292 1.698669 31 18
3: cDC 0.0058119261 0.0087178891 -0.8722269 -1.428073 49 17
4: ncMono 0.0004419890 0.0009944751 -0.8169229 -1.466723 3 49
5: B cell 0.0001107420 0.0003322259 -0.8439880 -1.509311 0 47
6: cMono 0.0001107420 0.0003322259 -0.8862629 -1.584912 0 47
7: CD4 T cell 0.0001103266 0.0003322259 -0.8893830 -1.599761 0 50
leadingEdge
<list>
1: NKG7, GN....
2: CCL5, GZ....
3: HLA-DRA,....
4: COTL1, F....
5: RPL18A, ....
6: S100A9, ....
7: RPS13, R....
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.
<- unlist(lapply(res, function(x) {
new.cluster.ids as.data.frame(x)[1, 1]
}))
$ref_gsea <- new.cluster.ids[as.character(alldata$leiden_k20)]
alldata
wrap_plots(
plotReducedDim(alldata, dimred = "UMAP", colour_by = "leiden_k20"),
plotReducedDim(alldata, dimred = "UMAP", colour_by = "ref_gsea"),
ncol = 2
)
Compare the results with the other celltype prediction methods in the ctrl_13 sample.
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.
<- file.path("data/human_cell_markers.txt")
path_file if (!file.exists(path_file)) download.file(file.path(path_data, "misc/cell_marker_human.csv"), destfile = path_file, method = "curl", extra = curl_upass)
<- read.delim("data/human_cell_markers.txt")
markers <- markers[markers$speciesType == "Human", ]
markers <- markers[markers$cancerType == "Normal", ]
markers
# 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.
<- lapply(unique(markers$cellName), function(x) {
celltype_list <- 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)
x
})names(celltype_list) <- unique(markers$cellName)
# celltype_list <- lapply(celltype_list , function(x) {x[1:min(length(x),50)]} )
<- celltype_list[unlist(lapply(celltype_list, length)) < 100]
celltype_list <- celltype_list[unlist(lapply(celltype_list, length)) > 5] celltype_list
# run fgsea for each of the clusters in the list
<- lapply(DGE_list, function(x) {
res $logFC <- rowSums(as.matrix(x[, grep("logFC", colnames(x))]))
x<- setNames(x$logFC, rownames(x))
gene_rank <- fgsea(pathways = celltype_list, stats = gene_rank, nperm = 10000)
fgseaRes return(fgseaRes)
})names(res) <- names(DGE_list)
# You can filter and resort the table based on ES, NES or pvalue
<- lapply(res, function(x) {
res $pval < 0.01, ]
x[x
})<- lapply(res, function(x) {
res $size > 5, ]
x[x
})<- lapply(res, function(x) {
res order(x$NES, decreasing = T), ]
x[
})
# show top 3 for each cluster.
lapply(res, head, 3)
$`1`
pathway pval padj ES NES
<char> <num> <num> <num> <num>
1: CD1C+_B dendritic cell 0.005714286 0.09215686 0.9240174 2.278437
2: Neutrophil 0.009803922 0.10112794 0.8330121 2.225166
3: Monocyte derived dendritic cell 0.002100840 0.07405411 0.9258710 1.966663
nMoreExtreme size leadingEdge
<num> <int> <list>
1: 0 54 S100A8, ....
2: 0 82 S100A8, ....
3: 0 18 S100A8, ....
$`10`
pathway pval padj ES NES
<char> <num> <num> <num> <num>
1: Megakaryocyte erythroid cell 0.00779922 0.08625020 -0.7071694 -1.187934
2: Natural killer cell 0.00380000 0.07937778 -0.7167473 -1.204219
3: Neutrophil 0.00219978 0.05169483 -0.7256130 -1.218451
nMoreExtreme size leadingEdge
<num> <int> <list>
1: 77 83 PTPRC, C....
2: 37 84 PTPRC, N....
3: 21 82 PTPRC, I....
$`11`
pathway pval padj ES NES nMoreExtreme
<char> <num> <num> <num> <num> <num>
1: Stromal cell 0.000646134 0.04049106 0.8411401 1.728948 2
2: Mesenchymal cell 0.003513123 0.09435243 0.7389426 1.637908 16
3: Neutrophil 0.001625356 0.06111337 0.7063061 1.630526 7
size leadingEdge
<int> <list>
1: 38 PECAM1, ....
2: 61 COTL1, S....
3: 82 LST1, FC....
$`12`
pathway pval padj ES NES nMoreExtreme
<char> <num> <num> <num> <num> <num>
1: T helper9 (Th9) cell 0.002138351 0.03382502 0.9414120 1.428838 19
2: Naive T cell 0.002313354 0.03382502 0.9275264 1.428112 21
3: CD8+ T cell 0.001742160 0.03382502 0.8930096 1.415301 16
size leadingEdge
<int> <list>
1: 10 CD3E, CD....
2: 12 IL7R, CD....
3: 19 IL7R, CD....
$`13`
pathway pval padj ES NES
<char> <num> <num> <num> <num>
1: Neutrophil 0.0002906977 0.01481015 0.8834234 2.234207
2: CD1C+_B dendritic cell 0.0002794857 0.01481015 0.9185908 2.172957
3: Megakaryocyte 0.0002732240 0.01481015 0.8960407 1.880966
nMoreExtreme size leadingEdge
<num> <int> <list>
1: 0 82 S100A8, ....
2: 0 54 S100A8, ....
3: 0 26 PPBP, PF....
$`14`
pathway pval padj ES NES
<char> <num> <num> <num> <num>
1: Monocyte derived dendritic cell 0.000303859 0.010767623 0.9303951 1.420317
2: Neutrophil 0.000099990 0.009403762 0.8820943 1.416893
3: Hemangioblast 0.001068262 0.015273941 0.9721348 1.412983
nMoreExtreme size leadingEdge
<num> <int> <list>
1: 2 18 S100A9, ....
2: 0 82 S100A9, ....
3: 9 8 PECAM1, CDH1
$`15`
pathway pval padj ES
<char> <num> <num> <num>
1: Naive CD8+ T cell 0.006802721 0.07856247 -0.7007477
2: Stem cell 0.005329847 0.07157223 -0.7492465
3: Specialist antigen presenting cell 0.006542526 0.07856247 -0.7565157
NES nMoreExtreme size leadingEdge
<num> <num> <int> <list>
1: -1.189770 67 91 RPS8, AI....
2: -1.249298 52 52 VIM, CD4....
3: -1.255061 64 46 CD48, CD....
$`2`
pathway pval padj ES NES
<char> <num> <num> <num> <num>
1: CD4+CD25+ regulatory T cell 0.007720492 0.05805810 -0.9655285 -1.415802
2: T helper17 (Th17) cell 0.008414164 0.05858751 -0.8162433 -1.418422
3: Granulocyte 0.009432821 0.06333466 -0.8610190 -1.425316
nMoreExtreme size leadingEdge
<num> <int> <list>
1: 58 6 CD3E, PT....
2: 71 28 CD3E, CD....
3: 77 18 CD63, PT....
$`3`
pathway pval padj ES NES
<char> <num> <num> <num> <num>
1: Neutrophil 0.0002157032 0.02053298 0.9126227 2.148260
2: CD1C+_B dendritic cell 0.0002184360 0.02053298 0.9239943 2.048965
3: Stromal cell 0.0024763620 0.08186371 0.8242241 1.722839
nMoreExtreme size leadingEdge
<num> <int> <list>
1: 0 82 S100A9, ....
2: 0 54 S100A9, ....
3: 10 38 VIM, TIM....
$`4`
pathway pval padj ES NES nMoreExtreme
<char> <num> <num> <num> <num> <num>
1: Naive CD8+ T cell 0.0028328612 0.03672018 0.8335019 2.296178 0
2: Naive CD4+ T cell 0.0009813543 0.02049940 0.9020382 2.098437 0
3: CD4+ T cell 0.0007936508 0.02049940 0.9168569 2.016731 0
size leadingEdge
<int> <list>
1: 91 LDHB, PI....
2: 34 IL7R, TC....
3: 25 LTB, IL7....
$`5`
pathway pval padj ES NES
<char> <num> <num> <num> <num>
1: CD4+ cytotoxic T cell 0.0010660981 0.02863235 0.8810115 2.234672
2: Natural killer cell 0.0010559662 0.02863235 0.7810959 1.971997
3: CD8+ T cell 0.0005668934 0.02664399 0.9445287 1.914053
nMoreExtreme size leadingEdge
<num> <int> <list>
1: 0 86 NKG7, CC....
2: 0 84 NKG7, GN....
3: 0 19 NKG7, CD....
$`6`
pathway pval padj ES NES
<char> <num> <num> <num> <num>
1: Follicular B cell 0.001248959 0.03772828 0.8849481 1.448385
2: Lake et al.Science.Ex5 0.005774216 0.06308120 0.9631247 1.413750
3: Germinal center B cell 0.009589652 0.07511894 0.9175945 1.395307
nMoreExtreme size leadingEdge
<num> <int> <list>
1: 11 22 MS4A1, C....
2: 48 6 RCSD1, A....
3: 85 9 IRF8, PA....
$`7`
pathway pval padj ES
<char> <num> <num> <num>
1: CD4+ cytotoxic T cell 0.0001004218 0.006301324 0.9001302
2: Effector CD8+ memory T (Tem) cell 0.0001005530 0.006301324 0.8819448
3: Cytotoxic T cell 0.0003219921 0.015133627 0.8840783
NES nMoreExtreme size leadingEdge
<num> <num> <int> <list>
1: 1.613471 0 86 NKG7, GN....
2: 1.577777 0 79 GNLY, FG....
3: 1.478784 2 23 GZMB, GZ....
$`8`
pathway pval padj ES NES nMoreExtreme
<char> <num> <num> <num> <num> <num>
1: Megakaryocyte 0.001955034 0.1115692 0.8280536 1.809771 1
2: Platelet 0.005188067 0.1464878 0.6629861 1.640075 3
3: Natural killer cell 0.009988555 0.1877848 -0.6753616 -1.313413 95
size leadingEdge
<int> <list>
1: 26 PPBP, PF....
2: 45 GP9, ITG....
3: 84 PTPRC, N....
$`9`
pathway pval padj ES NES
<char> <num> <num> <num> <num>
1: CD4+ cytotoxic T cell 0.001464129 0.03937997 0.9473821 2.311910
2: Effector CD8+ memory T (Tem) cell 0.001388889 0.03937997 0.8677227 2.094062
3: Natural killer cell 0.001466276 0.03937997 0.7976773 1.935329
nMoreExtreme size leadingEdge
<num> <int> <list>
1: 0 86 NKG7, GN....
2: 0 79 GNLY, GZ....
3: 0 84 NKG7, GN....
#CT_GSEA8:
<- unlist(lapply(res, function(x) {
new.cluster.ids as.data.frame(x)[1, 1]
}))$cellmarker_gsea <- new.cluster.ids[as.character(alldata$leiden_k20)]
alldata
wrap_plots(
plotReducedDim(alldata, dimred = "UMAP", colour_by = "cellmarker_gsea"),
plotReducedDim(alldata, dimred = "UMAP", colour_by = "ref_gsea"),
ncol = 2
)
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.3 (2024-02-29)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS Big Sur ... 10.16
Matrix products: default
BLAS/LAPACK: /Users/asabjor/miniconda3/envs/seurat5_u/lib/libopenblasp-r0.3.28.dylib; LAPACK version 3.12.0
locale:
[1] sv_SE.UTF-8/sv_SE.UTF-8/sv_SE.UTF-8/C/sv_SE.UTF-8/sv_SE.UTF-8
time zone: Europe/Stockholm
tzcode source: system (macOS)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] fgsea_1.28.0 celldex_1.12.0
[3] Seurat_5.1.0 SeuratObject_5.0.2
[5] sp_2.1-4 SingleR_2.4.0
[7] scmap_1.24.0 scPred_1.9.2
[9] pheatmap_1.0.12 patchwork_1.3.0
[11] dplyr_1.1.4 scran_1.30.0
[13] scater_1.30.1 ggplot2_3.5.1
[15] scuttle_1.12.0 SingleCellExperiment_1.24.0
[17] SummarizedExperiment_1.32.0 Biobase_2.62.0
[19] GenomicRanges_1.54.1 GenomeInfoDb_1.38.1
[21] IRanges_2.36.0 S4Vectors_0.40.2
[23] BiocGenerics_0.48.1 MatrixGenerics_1.14.0
[25] matrixStats_1.5.0
loaded via a namespace (and not attached):
[1] spatstat.sparse_3.1-0 bitops_1.0-9
[3] lubridate_1.9.4 httr_1.4.7
[5] RColorBrewer_1.1-3 tools_4.3.3
[7] sctransform_0.4.1 R6_2.5.1
[9] lazyeval_0.2.2 uwot_0.1.16
[11] withr_3.0.2 gridExtra_2.3
[13] progressr_0.15.1 cli_3.6.3
[15] spatstat.explore_3.3-4 fastDummies_1.7.4
[17] labeling_0.4.3 spatstat.data_3.1-4
[19] randomForest_4.7-1.2 proxy_0.4-27
[21] ggridges_0.5.6 pbapply_1.7-2
[23] harmony_1.2.3 parallelly_1.41.0
[25] limma_3.58.1 RSQLite_2.3.9
[27] FNN_1.1.4.1 generics_0.1.3
[29] ica_1.0-3 spatstat.random_3.3-2
[31] Matrix_1.6-5 ggbeeswarm_0.7.2
[33] abind_1.4-5 lifecycle_1.0.4
[35] yaml_2.3.10 edgeR_4.0.16
[37] BiocFileCache_2.10.1 recipes_1.1.0
[39] SparseArray_1.2.2 Rtsne_0.17
[41] blob_1.2.4 grid_4.3.3
[43] promises_1.3.2 dqrng_0.3.2
[45] ExperimentHub_2.10.0 crayon_1.5.3
[47] miniUI_0.1.1.1 lattice_0.22-6
[49] beachmat_2.18.0 cowplot_1.1.3
[51] KEGGREST_1.42.0 pillar_1.10.1
[53] knitr_1.49 metapod_1.10.0
[55] future.apply_1.11.2 codetools_0.2-20
[57] fastmatch_1.1-6 leiden_0.4.3.1
[59] googleVis_0.7.3 glue_1.8.0
[61] spatstat.univar_3.1-1 data.table_1.15.4
[63] vctrs_0.6.5 png_0.1-8
[65] spam_2.11-0 gtable_0.3.6
[67] cachem_1.1.0 gower_1.0.1
[69] xfun_0.50 S4Arrays_1.2.0
[71] mime_0.12 prodlim_2024.06.25
[73] survival_3.8-3 timeDate_4041.110
[75] iterators_1.0.14 hardhat_1.4.0
[77] lava_1.8.0 statmod_1.5.0
[79] bluster_1.12.0 interactiveDisplayBase_1.40.0
[81] fitdistrplus_1.2-2 ROCR_1.0-11
[83] ipred_0.9-15 nlme_3.1-165
[85] bit64_4.5.2 filelock_1.0.3
[87] RcppAnnoy_0.0.22 irlba_2.3.5.1
[89] vipor_0.4.7 KernSmooth_2.23-26
[91] rpart_4.1.24 DBI_1.2.3
[93] colorspace_2.1-1 nnet_7.3-20
[95] tidyselect_1.2.1 curl_6.0.1
[97] bit_4.5.0.1 compiler_4.3.3
[99] BiocNeighbors_1.20.0 DelayedArray_0.28.0
[101] plotly_4.10.4 scales_1.3.0
[103] lmtest_0.9-40 rappdirs_0.3.3
[105] stringr_1.5.1 digest_0.6.37
[107] goftest_1.2-3 spatstat.utils_3.1-2
[109] rmarkdown_2.29 XVector_0.42.0
[111] htmltools_0.5.8.1 pkgconfig_2.0.3
[113] sparseMatrixStats_1.14.0 dbplyr_2.5.0
[115] fastmap_1.2.0 rlang_1.1.4
[117] htmlwidgets_1.6.4 shiny_1.10.0
[119] DelayedMatrixStats_1.24.0 farver_2.1.2
[121] zoo_1.8-12 jsonlite_1.8.9
[123] BiocParallel_1.36.0 ModelMetrics_1.2.2.2
[125] BiocSingular_1.18.0 RCurl_1.98-1.16
[127] magrittr_2.0.3 GenomeInfoDbData_1.2.11
[129] dotCall64_1.2 munsell_0.5.1
[131] Rcpp_1.0.13-1 viridis_0.6.5
[133] reticulate_1.40.0 stringi_1.8.4
[135] pROC_1.18.5 zlibbioc_1.48.0
[137] MASS_7.3-60.0.1 AnnotationHub_3.10.0
[139] plyr_1.8.9 parallel_4.3.3
[141] listenv_0.9.1 ggrepel_0.9.6
[143] deldir_2.0-4 Biostrings_2.70.1
[145] splines_4.3.3 tensor_1.5
[147] locfit_1.5-9.10 igraph_2.0.3
[149] spatstat.geom_3.3-4 RcppHNSW_0.6.0
[151] reshape2_1.4.4 ScaledMatrix_1.10.0
[153] BiocVersion_3.18.1 evaluate_1.0.1
[155] BiocManager_1.30.25 foreach_1.5.2
[157] httpuv_1.6.15 RANN_2.6.2
[159] tidyr_1.3.1 purrr_1.0.2
[161] polyclip_1.10-7 future_1.34.0
[163] scattermore_1.2 rsvd_1.0.5
[165] xtable_1.8-4 e1071_1.7-16
[167] RSpectra_0.16-2 later_1.4.1
[169] viridisLite_0.4.2 class_7.3-23
[171] tibble_3.2.1 memoise_2.0.1
[173] AnnotationDbi_1.64.1 beeswarm_0.4.0
[175] cluster_2.1.8 timechange_0.3.0
[177] globals_0.16.3 caret_6.0-94