ixxmu
commented
1 year ago 整合单细胞数据和Bulk数据的多种方法（一）：R包scAB by 生信菜鸟团
最近刷到NAR最新的一篇算法文章《scAB detects multiresolution cell states with clinical significance by integrating single-cell genomics and bulk sequencing data》，DOI10.1093/nar/gkac1109。通过整合单细胞基因组学（包括scRNA-seq和sc-ATAC-seq数据）和Bulk RNA 测序数据及表型信息，在这里作者提出了一种称为 scAB 的新算法。
值得一提的是，这篇文章的通讯作者Jin，也是R包CellChat的第一作者。前面我已多次分享了CellChat的经验帖：
实际上，随着单细胞测序技术的蓬勃发展，目前已经提出了层出不穷的方法用于衔接单细胞数据和Bulk测序数据，例如Scissor， scPrognosis 和DEGAS 等工具。其他工具我们后面再做介绍。本文将介绍：
scAB算法的工作流程； 使用R实现方法（github: https://github.com/jinworks/scAB）。
一. scAB算法的工作流程
具体的算法作者在method部分写的很清楚，但是我不是学数学的，这部分对我来说确实有点晦涩难懂。简单来说的话，
（A）输入数据为
单细胞数据，可以是scATAC-seq或者scRNA-seq数据； 其次是对应组织的Bulk RNA-seq数据，以及表型数据。
（B）然后，scAB基于单细胞和Bulk RNA-seq数据计算单个细胞和Bulk样本之间的成对的 Pearson 相关性，生成相似性矩阵。通过对单个细胞和bulk数据的每个样本之间的相似性矩阵进行知识和图形引导的矩阵分解，揭示了一组细胞分层模式。每个模式通常对应于与特定表型相关的已知生物过程/信号，是从模型输出的细胞加载矩阵的每一行中获得的。加载值（即权重）表示每个细胞在每个模式中的贡献，具有高加载值的细胞被定义为表型相关细胞。
（C）scAB 能够从推断的细胞加载矩阵 H 中同时检测粗粒和细粒表型相关细胞状态。每个细粒细胞状态由每个学习模式中的表型相关细胞组成，粗粒表型相关细胞状态由所有细粒细胞状态的表型相关细胞的联合定义。不同的细粒细胞亚群可能与某些共有和特定的肿瘤标志物相关。
（D） 使用Bulk RNA-seq 数据集评估识别的表型相关细胞状态和特征的临床意义，揭示了临床相关的基因特征，这些特征能够进行各种预后分析，例如生存分析、免疫疗法疗效预测和协同治疗预测。
二. scAB的代码实现
scAB R包github在 https://github.com/jinworks/scAB.
安装：
devtools::install_github("jinworks/scAB")
Installation of other dependencies if not automatically installed
Install Seurat using install.packages('Seurat').Install diptest using install.packages('diptest').Install multimode using install.packages('multimode').Install preprocessCore using if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager"); BiocManager::install("preprocessCore").
Step1. R包加载：
library(Seurat)
library(preprocessCore)
library(scAB)
Step2. 加载数据
如上所述，该算法需要三种类型的数据输入:
单细胞RNA-seq数据是Seurat对象， Bulk RNA-seq表达矩阵， 表型数据可以是一个有两列的矩阵，时间和状态，也可以是一个向量。
data("data_survival")
dim(sc_dataset)
sc_dataset
# An object of class Seurat
# 21287 features across 8853 samples within 1 assay
# Active assay: RNA (21287 features, 0 variable features)
head(bulk_dataset[,1:10])
head(phenotype)
# time status
# TCGA-2Y-A9GS-01 724 1
# TCGA-2Y-A9GT-01 1624 1
# TCGA-2Y-A9GU-01 1939 0
# TCGA-2Y-A9GV-01 2532 1
# TCGA-2Y-A9GW-01 1271 1
# TCGA-2Y-A9GX-01 2442 0
#样本名称需保持一致
colnames(bulk_dataset) == row.names(phenotype)
这里的表型数据是“生存与否”。
Step3. 处理单细胞数据
示例的单细胞数据sc_dataset是Seurat对象，可以看到并没有进行过降维聚类等处理。
作者提供了一个run_seurat函数，进行常规处理：
sc_dataset <- run_seurat(sc_dataset,verbose = FALSE)
sc_dataset
An object of class Seurat
# 21287 features across 8853 samples within 1 assay
# Active assay: RNA (21287 features, 3000 variable features)
# 2 dimensional reductions calculated: pca, umap
UMAP_celltype <- DimPlot(sc_dataset, reduction ="umap",
group.by="celltype",label = T)
UMAP_celltype
run_seurat本质上就是把标准化、高变基因特征选择、归一化、PCA、降维聚类等基本流程打包为一个函数，最后用预先注释好的信息进行可视化。
Step4. 创建scAB对象
使用“create_scAB”函数计算样本细胞相似度矩阵和样本表型得分。
scAB_data <- create_scAB(sc_dataset,bulk_dataset,phenotype)
Step5. 选择K值
函数'select_K'可以帮助选择k。当重构误差下降趋势减缓时，k是一个合适的值。
K <- select_K(scAB_data)
K
# [1] 4
Step5和Step6的运行速度比较慢。
Step6. 运行scAB模型
使用默认参数运行scAB模型。正则化参数α1和α2也可以通过'select_alpha'函数进行交叉验证。
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
# 7938 915
识别到915个scAB+细胞。
Step7. 可视化与表型相关的细胞
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+细胞以外，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
作者使用logfc.threshold = 1和p_val_adj < 0.05作为截断点来识别差异表达最多的特征基因，这有助于选择潜在的生物标志物。截止值可能取决于特定的数据集，用户可以调整该参数值。较低的截断值导致更多的候选特征基因：
markers <- FindMarkers(sc_dataset, ident.1 = "scAB+ cells",
group.by = 'scAB_select',
logfc.threshold = 1)
markers <- markers[which(markers$p_val_adj<0.05),]
head(markers)
# p_val avg_log2FC pct.1 pct.2 p_val_adj
# CFHR1 0 1.860090 0.314 0.014 0
# APOB 0 1.385898 0.439 0.045 0
# CPS1 0 1.182203 0.306 0.016 0
# TF 0 3.084592 0.487 0.044 0
# AHSG 0 3.219864 0.405 0.046 0
# KNG1 0 2.078271 0.374 0.030 0
总体来说，这个R包还是很不错的。一个小问题是，因为这个算法依赖NMF非负矩阵和相关系数矩阵，导致运行速度比较慢，测试数据只有8000多个细胞，但运行整个流程需要花好数个小时，因此如果有python版本的话，会不会更适合当前的大样本单细胞数据？
- END -
https://mp.weixin.qq.com/s/2cLFOgaOecDrDypmv7mBpA