上一期内容提到了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 -