将bulk数据特征整合进单细胞数据中：R包scAB

伴随着单细胞测序技术的发展和价格的下降，越来越多的课题组开始进行单细胞测序以进行更好的研究。尽管，单细胞测序可以为我们探究疾病的发生和发展带来更为精确的结果，但与bulk测序相比更为昂贵的价格使得单细胞测序的样本数量不会太大（大佬课题组请无视😄），这就导致很难将临床特征与疾病联系起来，特别是肿瘤纷繁复杂的临床分级和病理分级。因此，将大数据量的bulk测序得到的某些临床特征映射到单细胞数据中是解决这个问题很好的手段。之前 凑齐六个字老师 已经分享过了使用Scissor包将bulk数据映射到单细胞数据中，本文将介绍另一种scAB算法。

推文：Scissor算法-从含有表型的bulkRNA数据中提取信息进而鉴别单细胞亚群

开发这个包的作者是Jin，也是R包CellChat的第一作者（不愧是大佬）。（Zhang Q, Jin S, Zou X. scAB detects multiresolution cell states with clinical significance by integrating single-cell genomics and bulk sequencing data. Nucleic Acids Res. 2022;50(21):12112-12130. doi:10.1093/nar/gkac1109）

简而言之，scAB可以将具有临床特征（分期，分级，生存等）的bulk RNA数据投射到单细胞数据的具体类型的细胞上，为我们研究疾病的临床特征与某一类细胞的功能或状态提供依据。#需要注意的是，目前这个包能处理的临床特征有生存数据，或是二元分组数据，如响应与无响应、野生型与突变型、正常与疾病等，如果是多个分组可能没办法很好的使用。

scAB分析的步骤

数据准备：

需要三类数据：单细胞表达数据；bulk RNA数据；bulk RNA的临床特征数据

scAB R包github在 https://github.com/jinworks/scAB.

1、导入数据和加载R包

devtools::install_github("jinworks/scAB")
library(Seurat)
library(preprocessCore)
library(scAB)

scRNA <- readRDS("scRNA_AAA.rds") #这里用自己处理的单细胞数据
dataset <- readRDS("Allmerged.rds")#这里也是自己处理的bulkRNA数据
pd <- read.table("cluster.txt", header = TRUE, row.names = 1, check.names = FALSE, sep = "\t")#bulk数据的临床分类数据

# 确保 dataset 的列名与 pd 的行名一致
colnames(dataset) <- rownames(pd)

sc_dataset <- scRNA
#dim(sc_dataset)
#[1] 41617 18189

bulk_dataset <- dataset

#GSM6204824 GSM6204823 GSM6204822 GSM6204821 GSM6204820 GSM6204819 GSM6204818 GSM6204817 GSM6204816
#A1CF      5.221065   4.873666   5.366132   4.957653   5.511200   5.297416   5.240153   4.934747   5.095085
#A2M      11.108126  13.428196  12.875538  12.607416  10.259988  11.939849  11.491157  10.451503  11.617009

phenotype <- pd
colnames(phenotype)
#             Group1 Group2 Group3 Group4
#GSM6204824      0      0      1      0
#GSM6204823      0      0      1      0
#GSM6204822      0      0      0      1
#GSM6204821      0      0      0      1
#GSM6204820      0      0      1      0
#phenotype <- phenotype[,"Group1"]
#如果有多个表型分组，可以取其中一个，运行上面的代码

2、单细胞数据处理和分群

#建议使用自己已经处理好的单细胞数据
UMAP_celltype <- DimPlot(sc_dataset, reduction ="umap",
                         group.by="cell_type",label = T)
UMAP_celltype

默认颜色有点丑。。。忽视它！！！

3、创建scAB对象和分析

#按流程创建scAB对象：

scAB_data <- create_scAB(sc_dataset,bulk_dataset,phenotype)

#敲黑板！！！，这里如果用自己合并的单细胞数据可能会报错！因为scAB包默认使用单细胞数据中graphs的RNA_snn作为获取的参数。
#所以，需要看一下我们单细胞数据的相关名称
print(names(sc_dataset@graphs))
#[1] "integrated_nn"  "integrated_snn"

#果然没有，所以我们需要创建一个
sc_dataset@graphs$RNA_snn <- sc_dataset@graphs$integrated_snn

#再看看数据
print(names(sc_dataset@graphs))
#[1] "integrated_nn"  "integrated_snn" "RNA_snn"

#这回有了，可以按流程走了。需要注意的是，这一步及其耗费时间，我的单细胞数据包含1.8w个细胞，跑了4个半小时。。。

scAB_data <- create_scAB(sc_dataset,bulk_dataset,phenotype,method = c("binary"))

#如果临床特征数据是二分类变量用"binary"，如果是生存数据用"survival"

#自动选择合适的K值
K <- select_K(scAB_data)
K
#[1] 2

scAB_result <- scAB(Object=scAB_data, K=K)
sc_dataset <- findSubset(sc_dataset, 
                         scAB_Object = scAB_result, 
                         tred = 2)
table(sc_dataset$scAB_select)
#Other cells scAB+ cells 
#      16888        1301 
#得到1301个scAB+ cells

#可视化一下看看
UMAP_scAB <- DimPlot(sc_dataset,group.by="scAB_select",
+                      cols=c("#80b1d3","red"),
+                      pt.size=0.001,
+                      order=c("scAB+ cells","Other cells"))
> patchwork::wrap_plots(plots = list(UMAP_celltype,UMAP_scAB), ncol = 2)

成功出图～

除了这些基础功能外，scAB还利用findSubset函数对所有的scAB+细胞进行亚群细分。。。

UMAP_subset3 <- DimPlot(sc_dataset,group.by="scAB_Subset3",
                        cols=c("#80b1d3","red"),
                        pt.size=0.001,
                        order=c("scAB+ cells","Other cells"))
UMAP_subset5 <- DimPlot(sc_dataset,
                        group.by="scAB_Subset4",
                        cols=c("#80b1d3","red"),
                        pt.size=0.001,
                        order=c("scAB+ cells","Other cells"))
UMAP_subset <- patchwork::wrap_plots(plots = list(UMAP_subset3,
                                                  UMAP_subset5), 
                                     ncol = 2)
UMAP_subset

大家可以自己尝试一下。。。

最后，可以通过FindMarkersl来识别差异表达的特征基因，这有助于选择潜在的生物标志物。截止值可以自己设定：

markers <- FindMarkers(sc_dataset, ident.1 = "scAB+ cells", 
                        group.by = 'scAB_select', 
                        logfc.threshold = 1)
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s  
 markers <- markers[which(markers$p_val_adj<0.05),]
head(markers)
              p_val avg_log2FC pct.1 pct.2    p_val_adj
KRT19  4.443975e-33  -1.751544 0.028 0.148 1.849449e-28
AGR2   5.202878e-30  -1.831747 0.023 0.130 2.165282e-25
IGLL5  6.451323e-26  -1.271237 0.161 0.076 2.684847e-21
SPINK1 1.758940e-25  -2.368104 0.022 0.114 7.320182e-21
TFF1   1.078211e-22  -2.451591 0.020 0.103 4.487192e-18
KRT8   2.286126e-21  -1.627475 0.049 0.138 9.514170e-17

总的来说，用起来还是很方便的，需要准备的文件就三个，只是这个算法使用的是NMF非负矩阵和相关系数矩阵，会耗费很多时间，对计算机资源也有一定的要求。但相比于Scissor包，scAB能得到更多的阳性细胞。大家可以根据需求选择哪种算法作为自己研究的使用工具。

注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟