上一期内容提到了BayesPrism算法用于单细胞数据的反卷积，BayesPrism算法在实际应用中非常占用计算资源以及消耗使用者的时间。那么是否有较好的替代包呢？
曾老师告诉了我一个R包—InstaPrism，他希望我将其和BayesPrism算法做个对比。
开发者给InstaPrism的直接定义就是“用于快速实现BayesPrism”。
简单来说这个包是基于BayesPrism分析框架分析并把里边非常耗时的Gibbs采样步骤替换成了fixed-point algorithm，从而减少了运行时间。
接下来分别运行这两种算法并对结果做一个粗糙的对比。
BayesPrism分析步骤详见推文：https://mp.weixin.qq.com/s/Rj8PL6wJqwWhsmLG8tx89w。
BayesPrism
1、输入数据情况
dim(sc.dat)
#[1] 5042 20124dim(bk.dat)
#[1] 424 60533
sc.dat: 5042个细胞和20124个基因的单细胞数据
bk.dat：424个样本和60533个基因的bulkRNA数据(TCGA-LIHC的count数据）
2、运行BayesPrism
关键步骤 bp.res <- run.prism(prism = myPrism, n.cores=20)，设置了20线程
#核心代码
start_time <- Sys.time() # 记录初始时间
bp.res <- run.prism(prism = myPrism, n.cores=20)
end_time <- Sys.time() # 记录终止时间
end_time - start_time # 计算时间差# Time difference of 2.865344 hours
需要将近3个小时，如果能设置50线程的话时间可以进一步缩短。建议采用较高配置的电脑或者服务器运行。
InstaPrism-code from github
1、安装及加载R包
if (!require("BiocManager", quietly = TRUE))install.packages("BiocManager")
if (!require("Biobase", quietly = TRUE))BiocManager::install("Biobase")
if (!require("devtools", quietly = TRUE))install.packages("devtools")
devtools::install_github("humengying0907/InstaPrism")library(InstaPrism)
2、数据处理
开发者建议使用TPM或者count的bulkRNA数据，并且从例子来看TPM是更优选。但为了前后分析一致，本次运算仍使用count数据。
# 处理单细胞数据
a <- as.data.frame(GetAssayData(sc_dataset, layer = "counts")) # 行为基因，列为细胞名/样本名cell_meta <- sc_dataset@meta.data # 处理好的单细胞注释文件
a_ref <- refPrepare(sc_Expr = a,cell.type.labels = cell_meta$celltype, # 细胞类型cell.state.labels = cell_meta$state_type, # 细胞状态)
head(a_ref)[1:3,1:3]# bulkRNA数据
b_bulk <- exp # 行为基因，列为细胞名/样本名# 运行函数
start_time <- Sys.time() # 记录初始时间
deconv_res <- InstaPrism(bulk_Expr = exp, # bulkRNA数据refPhi_cs = a_ref) # 单细胞数据
end_time <- Sys.time() # 记录终止时间
end_time - start_time # 计算时间差# Time difference of 1.086436 mins
真的非常快哦。本来需要将近3小时，现在只需要1分钟。
3、提取数据
# 用@Post.ini.ct@theta可获得细胞的theta估计值
estimated_frac <- t(deconv_res@Post.ini.ct@theta)# 提取之前保存好的BayesPrism分析结果进行比较
theta <- as.data.frame(theta) #BayesPrism分析结果
estimated_frac <- as.data.frame(estimated_frac) #IntraPrism分析结果
names(theta) <- c("BayesPrism_1","BayesPrism_2","BayesPrism_3")
head(theta)[1:5,1:3]
# BayesPrism_1 BayesPrism_2 BayesPrism_3
#TCGA-ZP-A9CY-01A-11R-A38B-07 4.718213e-07 0.9999993 2.523228e-07
#TCGA-DD-AAVZ-01A-11R-A41C-07 6.814688e-02 0.9314891 3.640295e-04
#TCGA-DD-A1EC-01A-21R-A131-07 5.272019e-01 0.4400359 3.276221e-02
#TCGA-DD-A1EC-11A-11R-A131-07 2.532253e-07 0.9999996 1.089799e-07
#TCGA-DD-AADN-01A-11R-A41C-07 3.301156e-01 0.5707665 9.911796e-02names(estimated_frac) <- c("InstaPrism_1","InstaPrism_2","InstaPrism_3")
head(estimated_frac)[1:5,1:3]
# InstaPrism_1 InstaPrism_2 InstaPrism_3
#TCGA-ZP-A9CY-01A-11R-A38B-07 0.01309068 0.9869093 1.698626e-12
#TCGA-DD-AAVZ-01A-11R-A41C-07 0.05282472 0.9471751 2.133425e-07
#TCGA-DD-A1EC-01A-21R-A131-07 0.18103636 0.8130102 5.953452e-03
#TCGA-DD-A1EC-11A-11R-A131-07 0.01004673 0.9899533 8.360430e-13
#TCGA-DD-AADN-01A-11R-A41C-07 0.45022208 0.4644255 8.535239e-02identical(rownames(estimated_frac),rownames(theta))
# [1] TRUE
dat <- cbind(estimated_frac,theta)
其实从这里粗看就可以发现数据结果的差异还是非常大的。
4、相关性分析(粗糙的比较一下)
# 懒的写循环了，就三张图，徒手修改代码出图 hhh
library(ggstatsplot)
colnames(dat)
ggscatterstats(data = dat,x = InstaPrism_1, y = BayesPrism_1, xlab = "InstaPrism_1",ylab = "BayesPrism_1",bf.message = FALSE)
相关性的值分别在0.84，0.87，0.67。说实话，这个相关性的数据，应该不能让使用者满意吧。
为什么会出现这样的情况？目前能想到的有两个可能的主要原因，1、没有进行异常基因过滤；2、算法上的差异。
因此笔者重新回溯了开发者提供在github中的代码，发现并没有提供相应的过滤函数(英语不好，如果理解没有误的话)。同时笔者也尝试采用了BayesPrism算法中提供的过滤函数过滤了数据之后再跑IntraPrism，相关性更低了....
接下来了解了Gibbs 采样和Fixed-Point算法。Gibbs 采样通常用于贝叶斯推断，特别是当直接采样较难实现时。它在处理高维分布和复杂模型时具有优势，而定点算法常用于确定性优化问题，例如在确定的条件下求解系统的稳定状态或最优解。Fixed-Point算法在某些特定的情况下可以替代Gibbs采样：1）当所使用的模型相对简单、确定且具有良好的收敛性质时，定点算法更适用；2）在计算资源有限或需要快速结果的情况下，定点算法的高效性显得尤为重要；3）在某些高维度数据分析中，如果定点算法能够被证明在这些维度上具有稳定的收敛性，它可能是一个更好的选择；4）当有较强的先验信息或假设可以指导模型的构建和优化时，定点算法可以利用这些信息进行快速收敛；5）在对模型的鲁棒性要求较低的情况下，定点算法可以作为替代。总而言之，要使用 Fixed-Point算法需要满足一定的条件。
那么是否InstaPrism能平替BayesPrism呢? 恐怕大多数情况下还是不行的，还是老老实实的用BayesPrism进行分析吧~
参考资料：
1、github：https://github.com/humengying0907/InstaPrism?tab=readme-ov-file
2、生信碱移：https://mp.weixin.qq.com/s/NTOxq_soYmPVbSItXVKGPw
致谢：感谢曾老师，小洁老师以及生信技能树团队全体成员。
注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟
- END -
