果叔带你解码bulk RNA-seq的反卷积新工具—— BayesPrism

#方法一library("devtools");install_github("Danko-Lab/BayesPrism/BayesPrism")#如果方法一报错可以下载R包到本地自行安装（https://github.com/Danko-Lab/BayesPrism）devtools::install_local('E:/BayesPrism/BayesPrism.zip')##选择自己的路径哦！！

suppressWarnings(library(BayesPrism))

#加载数据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))

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,                         )

sc.stat <- plot.scRNA.outlier(  input=sc.dat, #make sure the colnames are gene symbol or ENSMEBL ID  cell.type.labels=cell.type.labels,  species="hs", return.raw=TRUE #return the data used for plotting.)

head(sc.stat)

bk.stat <- plot.bulk.outlier(  bulk.input=bk.dat, #确保列名为基因符号或 ENSMEBL ID        sc.input=sc.dat, #确保列名为基因符号或 ENSMEBL ID  cell.type.labels=cell.type.labels,  species="hs", #currently only human(hs) and mouse(mm) annotations are supportedreturn.raw=TRUE)

head(bk.stat)

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")     

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)

##要根据特征基因对计数矩阵进行子集化，请执行以下操作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)          

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,)          

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          

slotNames(bp.res)          

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.cv
head(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")