ixxmu
commented
4 months ago 大海哥带你解码bulk RNA-seq的反卷积新工具—— BayesPrism by 生信滩
引言
小伙伴们对CIBERSORTx去卷积的算法应该耳熟能详了吧！下面大海哥给大家介绍一下BayesPrism，一种基于贝叶斯模型的前沿技术，它巧妙地将单细胞RNA测序（scRNA-seq）作为参考谱融入Bulk RNA测序（RNA-seq）数据中，从而实现细胞类型评分和基因表达的后验分布推断。
介绍
BayesPrism 利用从匹配或相似组织类型采集的样本scRNA-seq，对Bulk RNA-seq和空间转录组学进行细胞类型和基因表达去卷积。它将 scRNA-seq 作为先验信息处理，并估计细胞类型分数和细胞类型特异性基因表达在每个批次中的联合后验分布。
大海哥送福利：生信人必备神器——服务器
平时生信分析学习中有要的小伙伴可以联系大海哥租赁，粉丝福利都是市场超低价格，赶快找大海哥领取免费的试用账号吧！
定制生信分析
云服务器租赁
（加微信备注99领取使用）
BayesPrism 由去卷积模块和嵌入学习模块组成。解卷积模块根据scRNA-seq的细胞类型特异性表达谱建立先验模型，以联合估计肿瘤（或非肿瘤）样本的大量RNA-seq表达的细胞类型组成和细胞类型特异性基因表达的后验分布。嵌入学习模块使用期望最大化（EM）技术，利用恶性基因程序的线性组合近似计算肿瘤表达量，同时以推断表达量和解卷积模块估计的非恶性细胞比例为条件。
示例
R包安装
#方法一library("devtools");install_github("Danko-Lab/BayesPrism/BayesPrism")#如果方法一报错可以下载R包到本地自行安装（https://github.com/Danko-Lab/BayesPrism）devtools::install_local('E:/BayesPrism/BayesPrism.zip')##选择自己的路径哦！！
加载BayesPrism包
suppressWarnings(library(BayesPrism))加载数据集
rdata 文件存储在 "your_git_repository/tutorial.dat/tutorial.gbm.rdata"
地址：https://github.com/Danko-Lab/BayesPrism
rdata 文件包含运行 BayesPrism 所需的对象。
1.bk.dat：批量 RNA-seq 表达的按基因样本的原始计数矩阵。rownames是样本 ID，而 colnames是基因名称ID。
2.sc.dat：批量 RNA-seq 表达的逐个基因的原始计数矩阵。rownames是细胞ID，而 colnames是基因名称ID。这里大海提醒大家，如果你的sc.dat矩阵是稀疏矩阵，你应该将其转换为密集矩阵，即sc.dat <- as.matrix(sc.dat)
3.cell.type.labels：是一个字符向量，长度与nrow(sc.dat)相同，表示sc.dat中每个细胞的类型。
4.cell.state.labels：表示sc.dat中每个细胞的亚类型。
#加载数据load("../tutorial.dat/tutorial.gbm.rdata")ls()#查看数据1. bk.datdim(bk.dat)head(rownames(bk.dat))head(colnames(bk.dat))#查看数据sc.datdim(sc.dat)head(rownames(sc.dat))head(colnames(sc.dat))sort(table(cell.state.labels))table(cbind.data.frame(cell.state.labels, cell.type.labels))
• 请注意，BayesPrism 将对mixture和 scRNA-seq 参考之间共享的基因执行反卷积，即通过交叉 colnames(mixture) 和 colnames(reference)。因此，请确保使用一致的基因注释。
•建议使用未标准化和未转换的计数数据。当原始计数不可用时，线性归一化（例如 TPM、RPM、RPKM、FPKM）也是可以接受的，因为 BayesPrism 对于参考和mixture之间的线性乘法差异具有鲁棒性。理想情况下，如果使用标准化数据，最好提供参考并通过相同方法进行批量转换。应避免对数转换。
•请确保所有细胞状态包含合理数量的细胞，例如>20或>50。
细胞状态之间和细胞类型之间的相关性分析
我们建议首先绘制细胞状态之间和细胞类型之间的成对相关矩阵。这将使我们对数据质量有所了解。在细胞类型/状态没有用足够的信息量表示的情况下（细胞计数低和/或库大小低），低质量的细胞类型/状态往往聚集在一起。用户可以以更高的置信度重新聚类数据，或者将这些单元类型/状态与最相似的单元类型/状态合并，或者删除它们（如果重新聚类和合并不合适）。
plot.cor.phi (input=sc.dat,input.labels=cell.state.labels,title="cell state correlation",pdf.prefix="gbm.cor.cs",cexRow=0.2, cexCol=0.2,margins=c(2,2))
plot.cor.phi (input=sc.dat,input.labels=cell.type.labels,title="cell type correlation",#specify pdf.prefix if need to output to pdf#pdf.prefix="gbm.cor.ct",cexRow=0.5, cexCol=0.5,)
过滤异常基因
高强度表达的基因，例如核糖体蛋白基因和线粒体基因，可能主导分布并使推论产生偏差。这些基因通常无法提供区分细胞类型的信息，并且可能是大量虚假方差的来源。因此，它们可能不利于反卷积。我们建议去除这些基因。
用户可以从 scRNA-seq 参考中可视化异常基因的分布。我们计算所有细胞类型中每个基因的平均表达和它们的细胞类型特异性得分,可视化并确定 scRNA-seq 数据中的异常基因
sc.stat <- plot.scRNA.outlier(input=sc.dat, #make sure the colnames are gene symbol or ENSMEBL IDcell.type.labels=cell.type.labels,species="hs",return.raw=TRUE #return the data used for plotting.)
如图所示，核糖体蛋白基因通常表现出高平均表达和低细胞类型特异性得分。
如果需要，用户还可以根据函数输出的统计数据对 sc.dat 中的基因进行子集化。
head(sc.stat)sc.stat显示每个基因的归一化平均表达（x 轴）和最大特异性（y 轴）的对数，以及每个基因是否属于潜在的异常值类别。other_Rb是一组主要由核糖体假基因组成的基因。
同样，我们也可以可视化来自Bulk RNA 测序的异常基因。我们计算所有细胞类型中每个基因的平均表达值的和。由于我们没有来自Bulk RNA 测序数据的细胞类型水平信息，因此我们根据 scRNA-seq 计算细胞类型特异性得分，与上述相同。
bk.stat <- plot.bulk.outlier(bulk.input=bk.dat, #确保列名为基因符号或 ENSMEBL IDsc.input=sc.dat, #确保列名为基因符号或 ENSMEBL IDcell.type.labels=cell.type.labels,species="hs", #currently only human(hs) and mouse(mm) annotations are supportedreturn.raw=TRUE)
head(bk.stat)从 scRNA-seq 数据中过滤异常基因
接下来，我们从选定的组中删除基因。请注意，当mixture和混合物之间的性别不相同时，我们建议排除 chrX 和 chrY 中的基因。我们还删除了低转录基因，因为这些基因的转录测量往往容易受到噪音影响。去除这些基因也可以加快计算速度。
sc.dat.filtered <- cleanup.genes (input=sc.dat,input.type="count.matrix",species="hs",gene.group=c( "Rb","Mrp","other_Rb","chrM","MALAT1","chrX","chrY") ,exp.cells=5)dim(sc.dat.filtered)
接下来，我们检查不同类型基因的基因表达的一致性。由于批量和单细胞数据通常是通过不同的实验方案收集的，因此它们对不同类型的基因可能具有不同的敏感性。
plot.bulk.vs.sc (sc.input = sc.dat.filtered,bulk.input = bk.dat)
我们观察到蛋白质编码基因是两个测定之间最一致的组。为了减少批次效应并加快计算速度，我们对蛋白质编码基因进行反卷积。
##蛋白质编码基因的子集sc.dat.filtered.pc <- select.gene.type (sc.dat.filtered,gene.type = "protein_coding")
我们提供了一种通过使用来自不同细胞类型的细胞状态之间的成对 t 检验执行差异表达来选择基因的函数。也可以使用其他差异表达分析。
diff.exp.stat <- get.exp.stat(sc.dat=sc.dat[,colSums(sc.dat>0)>3],cell.type.labels=cell.type.labels,cell.state.labels=cell.state.labels,pseudo.count=0.1,n.cores=1)
理想情况下，我们希望为每种细胞类型选择足够数量的基因（> 50）。如果没有，用户可以降低截止值。
##要根据特征基因对计数矩阵进行子集化，请执行以下操作sc.dat.filtered.pc.sig <- select.marker (sc.dat=sc.dat.filtered.pc,stat=diff.exp.stat,pval.max=0.01,lfc.min=0.1)dim(sc.dat.filtered.pc.sig)
构造Prism对象
Prism 对象包含运行 BayesPrism 所需的所有数据，即 scRNA-seq 参考矩阵、每行参考的细胞类型和细胞状态标签以及批量 RNA-seq 的混合矩阵。
当使用 scRNA-seq 计数矩阵作为输入（推荐）时，用户需要指定input.type = "count.matrix"。另一个选项input.type是"GEP"（基因表达谱），它是基因矩阵的细胞状态。当使用从其他测定得出的参考（例如排序的批量数据）时，使用此选项。
cell.type.labels 参数键是与恶性细胞类型相对应的字符。NULL如果没有恶性细胞或参考和混合物之间的恶性细胞来自匹配样本，则设置为，在这种情况下，所有细胞类型将被同等对待。
myPrism <- new.prism(reference=sc.dat.filtered.pc,mixture=bk.dat,input.type="count.matrix",cell.type.labels = cell.type.labels,cell.state.labels = cell.state.labels,key="tumor",outlier.cut=0.01,outlier.fraction=0.1,)
运行BayesPrism
gibbs.control可以使用和指定控制吉布斯采样和优化的参数opt.control。我们建议使用默认参数。
bp.res <- run.prism(prism = myPrism, n.cores=50)#> Run Gibbs sampling...#> Current time: 2022-06-21 17:18:56#> Estimated time to complete: 46mins#> Estimated finishing time: 2022-06-21 18:04:15#> Start run...#> Warning in searchCommandline(parallel, cpus = cpus, type = type, socketHosts = socketHosts, : Unknown option on commandline: rmarkdown::render('/home/chut/tmp/vignette/vignettes/tutorial_deconvolution.Rmd',~+~~+~encoding~+~#> R Version: R version 4.1.2 (2021-11-01)#> snowfall 1.84-6.1 initialized (using snow 0.4-4): parallel execution on 50 CPUs.#>#> Stopping cluster#> Update the reference matrix ...#> Warning in searchCommandline(parallel, cpus = cpus, type = type, socketHosts = socketHosts, : Unknown option on commandline: rmarkdown::render('/home/chut/tmp/vignette/vignettes/tutorial_deconvolution.Rmd',~+~~+~encoding~+~#> snowfall 1.84-6.1 initialized (using snow 0.4-4): parallel execution on 50 CPUs.#>#>#> Stopping cluster#> Run Gibbs sampling using updated reference ...#> Current time: 2022-06-21 17:53:16#> Estimated time to complete: 17mins#> Estimated finishing time: 2022-06-21 18:10:09#> Start run...#> Warning in searchCommandline(parallel, cpus = cpus, type = type, socketHosts = socketHosts, : Unknown option on commandline: rmarkdown::render('/home/chut/tmp/vignette/vignettes/tutorial_deconvolution.Rmd',~+~~+~encoding~+~#> snowfall 1.84-6.1 initialized (using snow 0.4-4): parallel execution on 50 CPUs.#>#>#> Stopping cluster
提取结果
BayesPrism 预计会生成以下结果
slotNames(bp.res) 我们提供了从 run.prism 输出中提取解卷积细胞类型分数和表达的实用功能。
theta <- get.fraction (bp=bp.res,which.theta="final",state.or.type="type")head(theta)
#> tumor myeloid pericyte endothelial tcell oligo#> TCGA-06-2563-01A-01R-1849-01 0.8392297 0.04329259 2.999022e-02 0.07528272 6.474488e-04 0.0115573149#> TCGA-06-0749-01A-01R-1849-01 0.7090654 0.17001073 8.995526e-07 0.01275526 1.179331e-06 0.1081665709#> TCGA-06-5418-01A-01R-1849-01 0.8625322 0.09839143 9.729268e-03 0.02416954 7.039913e-07 0.0051768589#> TCGA-06-0211-01B-01R-1849-01 0.8893449 0.04482991 1.131622e-02 0.05435490 2.508238e-06 0.0001515524#> TCGA-19-2625-01A-01R-1850-01 0.9406438 0.03546026 1.932740e-03 0.01309753 4.535897e-06 0.0088610997#> TCGA-19-4065-02A-11R-2005-01 0.6763166 0.08374439 1.849921e-01 0.01918132 3.541126e-04 0.0354114398
theta.cv <- bp.res@posterior.theta_f@theta.cvhead(theta.cv)
#> tumor myeloid pericyte endothelial tcell oligo#> TCGA-06-2563-01A-01R-1849-01 0.0001722829 0.0016025331 0.0026568433 0.001853843 0.05452368 0.005606057#> TCGA-06-0749-01A-01R-1849-01 0.0002333853 0.0006859332 0.8109050326 0.005683516 0.74729658 0.001276784#> TCGA-06-5418-01A-01R-1849-01 0.0001601128 0.0009389782 0.0070872942 0.004131925 0.83670176 0.011362855#> TCGA-06-0211-01B-01R-1849-01 0.0001175529 0.0014412303 0.0055064706 0.001673105 0.60905434 0.280609474#> TCGA-19-2625-01A-01R-1850-01 0.0001327408 0.0023661566 0.0335184780 0.006286823 0.73453327 0.006138184#> TCGA-19-4065-02A-11R-2005-01 0.0002862600 0.0012227981 0.0009100677 0.005660069 0.06371147 0.002243368
Z.tumor <- get.exp (bp=bp.res,state.or.type="type",cell.name="tumor")head(t(Z.tumor[1:5,]))
#> TCGA-06-2563-01A-01R-1849-01 TCGA-06-0749-01A-01R-1849-01 TCGA-06-5418-01A-01R-1849-01 TCGA-06-0211-01B-01R-1849-01 TCGA-19-2625-01A-01R-1850-01#> ENSG00000130876.10 55.980 444.980 9.000 56.996 15.996#> ENSG00000165244.6 206.344 38.096 291.872 530.732 507.788#> ENSG00000173597.7 17.100 6.588 21.804 28.744 10.800#> ENSG00000158022.6 0.000 0.476 2.616 0.960 7.824#> ENSG00000167220.10 2273.824 1060.976 2775.180 2368.016 2084.620#> ENSG00000126106.12 678.468 685.948 481.000 698.888 1296.080
保存结果
save(bp.res, file="bp.res.rdata") 紧接着可以使用上述结果进行潜在的下游分析，例如：
结语
在这里，大海哥给大家介绍了BayesPrism。基于贝叶斯模型，对单细胞参考数据和批量数据之间的基因表达差异进行明确建模和边际化。BayesPrism 在肿瘤和非肿瘤环境下的细胞类型分数推断方面都大大优于目前的其它的方法。BayesPrism的出现为整合分析Bulk和 scRNA-seq 数据引入了一个强大的新工具。
大海哥还提供思路设计、定制生信分析、文献思路复现；有需要的小伙伴欢迎直接扫码咨询~
定制生信分析
服务器租赁
扫码咨询大海哥
往期回顾
01 |
02 |
03 |
04 |
https://mp.weixin.qq.com/s/I_wIRu4v3KhdbkmF26MFXA