projectLSI：将你的单细胞或bulk转录组数据映射到参考数据集中

原创 Davey1220 bioinfomics

简介

在单细胞数据分析过程中，我们经常会遇到不同样本之间整合的批次效应和细胞类型注释的困难，projectLSI包利用term frequency–inverse document frequency (TF-IDF) transformation and latent semantic indexing (LSI)算法进行数据降维转换，可以将query的单细胞或bulk转录组数据集映射到reference参考数据集中以消除潜在的批次效应，同时也可以利用bulk转录组数据验证单细胞注释分群的结果。

R包安装

devtools::install_github("sajuukLyu/projectLSI")

实例演示

接下来，我们将使用两个单细胞转录组数据集pbmc3k和pbmc4k，以及一个bulk转录组数据集bulk.data进行实例演示projectLSI包的使用流程。我们将以pbmc3k数据集作为参考数据集，使用projectLSI程序将pbmc4k和bulk.data数据集映射到参考数据集中。

加载所需R包和示例数据集

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"

单细胞数据预处理

# for pbmc3k
# 计算线粒体含量
pbmc3k$pct.mt <- PercentageFeatureSet(pbmc3k, pattern = "^MT-")
# 数据质控
FeatureScatter(pbmc3k, feature1 = "nCount_RNA", feature2 = "pct.mt") +
FeatureScatter(pbmc3k, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")

# 质控过滤
pbmc3k <- subset(pbmc3k, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & pct.mt < 5)
# 数据标准化
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")

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%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|

进行线性降维处理

这里，我们将使用projectLSI包对标准化的数据进行TF-IDF和LSI线性降维处理。注意，我们使用筛选出的高变异基因作为输入。

# 使用calcLSI函数进行线性降维
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对象中
pbmc3k[["pca"]] <- CreateDimReducObject(
  embeddings = pbmc3k.lsi$matSVD,
  loadings = pbmc3k.lsi$fLoad,
  assay = "RNA",
  stdev = pbmc3k.lsi$sdev,
  key = "PC_")

ElbowPlot(pbmc3k)

细胞聚类分群

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对象中
pbmc3k[["umap"]] <- CreateDimReducObject(
  embeddings = umap.pbmc3k.emb,
  assay = "RNA",
  key = "UMAP_")

DimPlot(pbmc3k, label = T)

# 查看marker基因的表达情况
FeaturePlot(pbmc3k, c("MS4A1", "GNLY", "CD3E",
                      "CD14", "FCER1A", "FCGR3A",
                      "LYZ", "PPBP", "CD8A"), order = T)

细胞类型注释

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()

将query单细胞数据集映射到参考数据集中

接下来，我们将使用projectLSI函数将pbmc4k查询数据集映射到参考数据集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)

FeaturePlot(pbmc4k, c("MS4A1", "GNLY", "CD3E",
                      "CD14", "FCER1A", "FCGR3A",
                      "LYZ", "PPBP", "CD8A"), order = T)

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()

接下来，我们将映射好的pbmc4k和pbmc3k数据集合并到一起进行展示

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()

pbmc7k$celltype <- Idents(pbmc7k)
Idents(pbmc7k) <- pbmc7k$orig.ident
DimPlot(pbmc7k)

可以看到，pbmc4k和参考数据集pbmc3k的细胞类型重叠的很好，不存在明显的批次效应。

将bulk转录组数据映射到参考数据集中

接下来，我们将使用projectLSI包将bulk转录组数据映射到单细胞参考数据集中。首先，我们将使用psudoSC函数将bulk转录组数据进行down-sampleing重取样构建psudo-single-cell数据。

# 设置n=100每个样本重取样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            .        .        .        .        .

# 构建seurat对象
psudo.so <- CreateSeuratObject(psudo.all, project = "bulk")
# 数据标准化
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函数进行数据映射
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)

接下来，我们将映射好的bulk转录组数据和pbmc3k参考数据集合并到一起进行展示

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()