簡(jiǎn)介
在單細(xì)胞數(shù)據(jù)分析過(guò)程中渐溶，我們經(jīng)常會(huì)遇到不同樣本之間整合的批次效應(yīng)和細(xì)胞類(lèi)型注釋的困難胖眷，projectLSI包利用term frequency–inverse document frequency (TF-IDF) transformation and latent semantic indexing (LSI)算法進(jìn)行數(shù)據(jù)降維轉(zhuǎn)換饶号，可以將query的單細(xì)胞或bulk轉(zhuǎn)錄組數(shù)據(jù)集映射到reference參考數(shù)據(jù)集中以消除潛在的批次效應(yīng)，同時(shí)也可以利用bulk轉(zhuǎn)錄組數(shù)據(jù)驗(yàn)證單細(xì)胞注釋分群的結(jié)果。
R包安裝
devtools::install_github("sajuukLyu/projectLSI")
實(shí)例演示
接下來(lái)撑毛，我們將使用兩個(gè)單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)集pbmc3k和pbmc4k牙躺，以及一個(gè)bulk轉(zhuǎn)錄組數(shù)據(jù)集bulk.data進(jìn)行實(shí)例演示projectLSI包的使用流程愁憔。我們將以pbmc3k數(shù)據(jù)集作為參考數(shù)據(jù)集，使用projectLSI程序?qū)bmc4k和bulk.data數(shù)據(jù)集映射到參考數(shù)據(jù)集中孽拷。
加載所需R包和示例數(shù)據(jù)集
pbmc3k and pbmc4k datasets are from package TENxPBMCData, and bulk.data is part of GSE74246.
library(Seurat)
library(projectLSI)
library(patchwork)
data(pbmc3k)
data(pbmc4k)
data(bulk.data)
pbmc3k
## An object of class Seurat
## 32738 features across 2700 samples within 1 assay
## Active assay: RNA (32738 features)
pbmc4k
## An object of class Seurat
## 33694 features across 4340 samples within 1 assay
## Active assay: RNA (33694 features)
bulk.data[1:5,1:5]
# CD4T_1 CD4T_2 CD4T_3 CD4T_4 CD8T_1
#A1BG 0 3 7 4 0
#A1BG-AS1 3 1 1 3 0
#A1CF 10 15 3 0 1
#A2M 141 273 870 92 351
#A2M-AS1 14 23 154 18 31
dim(bulk.data)
## [1] 25498 20
names(bulk.data)
# [1] "CD4T_1" "CD4T_2" "CD4T_3" "CD4T_4" "CD8T_1" "CD8T_2" "CD8T_3" "CD8T_4"
# [9] "NK_1" "NK_2" "NK_3" "NK_4" "B_1" "B_2" "B_3" "B_4"
# [17] "Mono_1" "Mono_2" "Mono_3" "Mono_4"
單細(xì)胞數(shù)據(jù)預(yù)處理
# for pbmc3k
# 計(jì)算線粒體含量
pbmc3k$pct.mt <- PercentageFeatureSet(pbmc3k, pattern = "^MT-")
# 數(shù)據(jù)質(zhì)控
FeatureScatter(pbmc3k, feature1 = "nCount_RNA", feature2 = "pct.mt") +
FeatureScatter(pbmc3k, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
image.png
# 質(zhì)控過(guò)濾
pbmc3k <- subset(pbmc3k, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & pct.mt < 5)
# 數(shù)據(jù)標(biāo)準(zhǔn)化
pbmc3k <- NormalizeData(pbmc3k)
## Performing log-normalization
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
# 篩選高變異基因
pbmc3k <- FindVariableFeatures(pbmc3k, nfeatures = 2000)
## Calculating gene variances
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## Calculating feature variances of standardized and clipped values
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
# for pbmc4k
pbmc4k$pct.mt <- PercentageFeatureSet(pbmc4k, pattern = "^MT-")
FeatureScatter(pbmc4k, feature1 = "nCount_RNA", feature2 = "pct.mt") +
FeatureScatter(pbmc4k, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
image.png
pbmc4k <- subset(pbmc4k, subset = nFeature_RNA > 200 & nFeature_RNA < 3500 & pct.mt < 8)
pbmc4k <- NormalizeData(pbmc4k)
## Performing log-normalization
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
pbmc4k <- FindVariableFeatures(pbmc4k, nfeatures = 2000)
## Calculating gene variances
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## Calculating feature variances of standardized and clipped values
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
進(jìn)行線性降維處理
這里吨掌，我們將使用projectLSI包對(duì)標(biāo)準(zhǔn)化的數(shù)據(jù)進(jìn)行TF-IDF和LSI線性降維處理。注意脓恕，我們使用篩選出的高變異基因作為輸入膜宋。
# 使用calcLSI函數(shù)進(jìn)行線性降維
pbmc3k.lsi <- calcLSI(pbmc3k[["RNA"]]@data[VariableFeatures(pbmc3k), ])
pbmc3k.lsi$matSVD[1:5,1:5]
# PC_1 PC_2 PC_3 PC_4
#3k_AAACATACAACCAC-1 -0.03787607 -1.038326e-02 -0.0009909208 0.005607445
#3k_AAACATTGAGCTAC-1 -0.03799807 -4.920159e-05 0.0066228430 -0.018501420
#3k_AAACATTGATCAGC-1 -0.03885814 -8.351669e-03 -0.0021662553 0.005449835
# 將LSI降維信息添加到seurat對(duì)象中
pbmc3k[["pca"]] <- CreateDimReducObject(
embeddings = pbmc3k.lsi$matSVD,
loadings = pbmc3k.lsi$fLoad,
assay = "RNA",
stdev = pbmc3k.lsi$sdev,
key = "PC_")
ElbowPlot(pbmc3k)
image.png
細(xì)胞聚類(lèi)分群
pbmc3k <- FindNeighbors(pbmc3k, dims = 1:10)
## Computing nearest neighbor graph
## Computing SNN
pbmc3k <- FindClusters(pbmc3k, resolution = 0.6)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 2638
## Number of edges: 97177
##
## Running Louvain algorithm...
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## Maximum modularity in 10 random starts: 0.8542
## Number of communities: 9
## Elapsed time: 0 seconds
UMAP降維可視化
Notice that ret_model parameter should be TRUE for later projection.
set.seed(42)
umap.pbmc3k <- uwot::umap(pbmc3k.lsi$matSVD[, 1:10],
n_neighbors = 30,
min_dist = 0.5,
metric = "euclidean",
ret_model = T,
verbose = T)
## 00:58:06 UMAP embedding parameters a = 0.583 b = 1.334
## 00:58:06 Read 2638 rows and found 10 numeric columns
## 00:58:06 Using Annoy for neighbor search, n_neighbors = 30
## 00:58:06 Building Annoy index with metric = euclidean, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 00:58:06 Writing NN index file to temp file /tmp/RtmpbWcqgH/file17b95ca52051
## 00:58:06 Searching Annoy index using 8 threads, search_k = 3000
## 00:58:06 Annoy recall = 100%
## 00:58:07 Commencing smooth kNN distance calibration using 8 threads
## 00:58:07 Initializing from normalized Laplacian + noise
## Spectral initialization failed to converge, using random initialization instead
## 00:58:07 Commencing optimization for 500 epochs, with 110674 positive edges
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 00:58:12 Optimization finished
# 提取UMAP降維信息
umap.pbmc3k.emb <- umap.pbmc3k$embedding
rownames(umap.pbmc3k.emb) <- colnames(pbmc3k)
colnames(umap.pbmc3k.emb) <- paste0("UMAP_", seq_len(ncol(umap.pbmc3k.emb)))
head(umap.pbmc3k.emb)
# UMAP_1 UMAP_2
#3k_AAACATACAACCAC-1 4.462377 1.675100
#3k_AAACATTGAGCTAC-1 -1.583569 -9.764045
#3k_AAACATTGATCAGC-1 7.612832 1.880275
# 將UMAP降維信息添加到seurat對(duì)象中
pbmc3k[["umap"]] <- CreateDimReducObject(
embeddings = umap.pbmc3k.emb,
assay = "RNA",
key = "UMAP_")
DimPlot(pbmc3k, label = T)
image.png
# 查看marker基因的表達(dá)情況
FeaturePlot(pbmc3k, c("MS4A1", "GNLY", "CD3E",
"CD14", "FCER1A", "FCGR3A",
"LYZ", "PPBP", "CD8A"), order = T)
image.png
細(xì)胞類(lèi)型注釋
image.png
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono",
"B", "NK", "FCGR3A+ Mono",
"CD8 T", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc3k)
pbmc3k <- RenameIdents(pbmc3k, new.cluster.ids)
DimPlot(pbmc3k, label = T) + NoLegend()
image.png
將query單細(xì)胞數(shù)據(jù)集映射到參考數(shù)據(jù)集中
接下來(lái)，我們將使用projectLSI函數(shù)將pbmc4k查詢數(shù)據(jù)集映射到參考數(shù)據(jù)集pbmc3k中进肯。
matSVD.pbmc4k <- projectLSI(pbmc4k[["RNA"]]@data, pbmc3k.lsi)
head(matSVD.pbmc4k)
# PC_1 PC_2 PC_3 PC_4
#4k_AAACCTGAGAAGGCCT-1 -0.03410563 0.02377144 0.0002275122 0.001904910
#4k_AAACCTGAGACAGACC-1 -0.03532231 0.02358280 0.0023570705 0.001242694
#4k_AAACCTGAGATAGTCA-1 -0.03497681 0.02380977 0.0012423880 0.002034311
pbmc4k[["pca"]] <- CreateDimReducObject(
embeddings = matSVD.pbmc4k,
loadings = pbmc3k.lsi$fLoad,
assay = "RNA",
key = "PC_")
# cluster cells using projected LSI
pbmc4k <- FindNeighbors(pbmc4k, dims = 1:10)
## Computing nearest neighbor graph
## Computing SNN
pbmc4k <- FindClusters(pbmc4k, resolution = 0.6)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 4284
## Number of edges: 154662
##
## Running Louvain algorithm...
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## Maximum modularity in 10 random starts: 0.8670
## Number of communities: 10
## Elapsed time: 0 seconds
# perform UMAP using first 10 PCs, just like pbmc3k
umap.pbmc4k.proj <- uwot::umap_transform(matSVD.pbmc4k[, 1:10], umap.pbmc3k, verbose = T)
## 01:37:48 Read 4284 rows and found 10 numeric columns
## 01:37:48 Processing block 1 of 1
## 01:37:48 Writing NN index file to temp file /tmp/RtmpbWcqgH/file17b933607747
## 01:37:48 Searching Annoy index using 8 threads, search_k = 3000
## 01:37:48 Commencing smooth kNN distance calibration using 8 threads
## 01:37:48 Initializing by weighted average of neighbor coordinates using 8 threads
## 01:37:48 Commencing optimization for 167 epochs, with 128520 positive edges
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 01:37:50 Finished
head(umap.pbmc4k.proj)
# [,1] [,2]
#[1,] -11.620971 1.3273895
#[2,] -12.708588 -0.8110412
#[3,] -12.193903 -0.9287748
rownames(umap.pbmc4k.proj) <- colnames(pbmc4k)
colnames(umap.pbmc4k.proj) <- paste0("UMAP_", seq_len(ncol(umap.pbmc4k.proj)))
pbmc4k[["umap"]] <- CreateDimReducObject(
embeddings = umap.pbmc4k.proj,
assay = "RNA",
key = "UMAP_")
DimPlot(pbmc4k, label = T)
image.png
FeaturePlot(pbmc4k, c("MS4A1", "GNLY", "CD3E",
"CD14", "FCER1A", "FCGR3A",
"LYZ", "PPBP", "CD8A"), order = T)
image.png
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "B", "Memory CD4 T",
"CD8 T", "CD14+ Mono", "NK", "FCGR3A+ Mono",
"DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc4k)
pbmc4k <- RenameIdents(pbmc4k, new.cluster.ids)
DimPlot(pbmc4k, label = T) + NoLegend()
image.png
接下來(lái)激蹲，我們將映射好的pbmc4k和pbmc3k數(shù)據(jù)集合并到一起進(jìn)行展示
pbmc7k <- merge(pbmc3k, pbmc4k)
pbmc7k[["umap"]] <- CreateDimReducObject(
embeddings = rbind(pbmc3k[["umap"]]@cell.embeddings,
pbmc4k[["umap"]]@cell.embeddings),
assay = "RNA", key = "UMAP_")
DimPlot(pbmc7k, label = T) + NoLegend()
image.png
pbmc7k$celltype <- Idents(pbmc7k)
Idents(pbmc7k) <- pbmc7k$orig.ident
DimPlot(pbmc7k)
image.png
可以看到，pbmc4k和參考數(shù)據(jù)集pbmc3k的細(xì)胞類(lèi)型重疊的很好江掩，不存在明顯的批次效應(yīng)学辱。
將bulk轉(zhuǎn)錄組數(shù)據(jù)映射到參考數(shù)據(jù)集中
接下來(lái)，我們將使用projectLSI包將bulk轉(zhuǎn)錄組數(shù)據(jù)映射到單細(xì)胞參考數(shù)據(jù)集中环形。首先策泣，我們將使用psudoSC函數(shù)將bulk轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行down-sampleing重取樣構(gòu)建psudo-single-cell數(shù)據(jù)。
# 設(shè)置n=100每個(gè)樣本重取樣100次
psudo.all <- psudoSC(bulk.data, n = 100, depth = 3000)
## downsampling counts...
## merging all samples...
dim(psudo.all)
## [1] 25498 2000
psudo.all[1:5,1:5]
#5 x 5 sparse Matrix of class "dgCMatrix"
# CD4T_1_1 CD4T_1_2 CD4T_1_3 CD4T_1_4 CD4T_1_5
#A1BG . . . . .
#A1BG-AS1 . . . . .
#A1CF . . . . .
# 構(gòu)建seurat對(duì)象
psudo.so <- CreateSeuratObject(psudo.all, project = "bulk")
# 數(shù)據(jù)標(biāo)準(zhǔn)化
psudo.so <- NormalizeData(psudo.so)
## Performing log-normalization
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
Idents(psudo.so) <- rep(c("CD4T.bulk", "CD8T.bulk", "NK.bulk", "B.bulk", "Mono.bulk"), rep(400, 5))
table(Idents(psudo.so))
#CD4T.bulk CD8T.bulk NK.bulk B.bulk Mono.bulk
# 400 400 400 400 400
# 使用projectLSI函數(shù)進(jìn)行數(shù)據(jù)映射
bulk.matSVD <- projectLSI(psudo.so[["RNA"]]@data, pbmc3k.lsi)
head(bulk.matSVD)
# PC_1 PC_2 PC_3 PC_4 PC_5
#CD4T_1_1 -0.03980147 -0.01410879 -0.0014396458 5.488960e-04 0.011402060
#CD4T_1_2 -0.04032496 -0.01581425 0.0003935543 1.469935e-03 0.009831236
#CD4T_1_3 -0.04006459 -0.01520021 0.0001330849 4.762865e-04 0.009996544
umap.bulk.proj <- uwot::umap_transform(bulk.matSVD[, 1:10], umap.pbmc3k, verbose = T)
## 13:16:17 Read 2000 rows and found 10 numeric columns
## 13:16:17 Processing block 1 of 1
## 13:16:17 Writing NN index file to temp file /tmp/RtmplGf1gl/file1a6a91b3753
## 13:16:17 Searching Annoy index using 8 threads, search_k = 3000
## 13:16:17 Commencing smooth kNN distance calibration using 8 threads
## 13:16:17 Initializing by weighted average of neighbor coordinates using 8 threads
## 13:16:17 Commencing optimization for 167 epochs, with 60000 positive edges
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 13:16:18 Finished
rownames(umap.bulk.proj) <- colnames(psudo.so)
colnames(umap.bulk.proj) <- paste0("UMAP_", seq_len(ncol(umap.bulk.proj)))
psudo.so[["umap"]] <- CreateDimReducObject(
embeddings = umap.bulk.proj,
assay = "RNA",
key = "UMAP_")
DimPlot(psudo.so, label = T)
image.png
接下來(lái)抬吟，我們將映射好的bulk轉(zhuǎn)錄組數(shù)據(jù)和pbmc3k參考數(shù)據(jù)集合并到一起進(jìn)行展示
pbmc.mix <- merge(pbmc3k, psudo.so)
pbmc.mix[["umap"]] <- CreateDimReducObject(
embeddings = rbind(pbmc3k[["umap"]]@cell.embeddings,
psudo.so[["umap"]]@cell.embeddings),
assay = "RNA", key = "UMAP_")
DimPlot(pbmc.mix, label = T) + NoLegend()
image.png
可以看到萨咕，bulk轉(zhuǎn)錄數(shù)據(jù)與單細(xì)胞參考數(shù)據(jù)集的細(xì)胞類(lèi)型可以很好的重疊在一起，雖然bulk轉(zhuǎn)錄組數(shù)據(jù)中的CD8 T細(xì)胞與單細(xì)胞數(shù)據(jù)中存在一定的偏差火本，但其他的細(xì)胞類(lèi)型注釋的很好危队。
