1简介

R包希尔达是在贝叶斯框架下开发的，允许在统计上测试两组之间突变签名的突变负担是否有变化。突变特征的定义基于Shiraishi等人提出的独立模型。

1．1纸

  • 白石等人。一种简单的基于模型的方法来推断和可视化癌症突变特征，bioRxiv, doi:http://dx.doi.org/10.1101/019901．

  • 杨直， Priyatama Pandey, Darryl Shibata, David V. Conti, Paul Marjoram, Kimberly D. Siegmund。HiLDA:研究突变特征差异的统计方法，bioRxiv, doi:https://doi.org/10.1101/577452

2安装和加载包

2.1安装

2.1.1Bioconductor

希尔达需要安装几个CRAN和Bioconductor R包。当使用以下命令安装包时，依赖关系通常会自动处理:

install.packages(“BiocManager”)BiocManager::安装(“希尔达”)

[注意:如果您已经安装了BiocManager．]

您也可以从GitHub下载最新版本devtools：

devtools: install_github(“USCbiostats /希尔达”)

2.1.2只是另一个吉布斯采样器(JAGS)

为了运行HiLDA，还需要安装一个名为Just Another Gibbs Sampler (JAGS)的外部程序，从这个网站下载http://mcmc-jags.sourceforge.net/．有关详细信息，请按照INSTALL文件安装该程序。

3.输入数据

希尔达一个包是否建立在一些基本的功能上pmsignature包括如何读取输入数据。这里有一个例子pmsignature在输入数据上，突变特性是否有用于对突变进行分类的元素，例如:

  • 6个替换(C>A, C>G, C>T, T>A, T>C和T>G)
  • 2个侧翼基地(A, C, G, T)
  • 转录方向。

3．1突变位置格式

样品1 chr1 100a C样品1 chr1 200a T样品1 chr2 100g T样品2 chr1 300t C样品3 chr3 400t C样品
  • 第一列显示样本的名称
  • 第二列是染色体的名称
  • 第三列是染色体的坐标
  • 第4列显示参考基(A、C、G或T)。
  • 第五列表示交替基底(A、C、G或T)。

4工作流

4．1获取输入数据

在这里,inputFile输入文件的路径。numBases包括中心碱基在内的侧翼碱基的数量是多少(如果你想考虑两个5 '和3 '碱基，那么设置5)。另外，你可以使用trDir．numSig设置从输入数据估计的突变签名数。您将看到一条关于某些突变正在被删除的警告消息。

库(希尔达)
##加载所需的包:ggplot2
inputFile <- system.file("extdata/食管.mp.txt.gz"， package="HiLDA") G <- hildaReadMPFile(inputFile, numBases=5, trDir=TRUE)
在hildaReadMPFile(inputFile, numBases = 5, trDir = TRUE)中警告:在24861个突变中，我们可以获得24728个突变的转录方向信息。其他突变被移除。

此外，我们还提供了一个包含10个突变目录的小型模拟数据集，并用它来演示HiLDA的关键功能。我们首先加载存储为extdata/sample.rdata的示例数据集G。

负载(执行(“extdata /样品。rdata"， package = "HiLDA"))类(G)
# #[1]“MutationFeatureData”# # attr(“包”)# #[1]“希尔达”

如果你想在手稿中使用南加州大学的数据，请从OSF主页下载数据https://osf.io/a8dzx/

4.2运行测试希尔达

4.2.1执行全局测试和本地测试

读入示例数据G之后，我们可以从HiLDA运行本地和全局测试。在这里，我们指定inputG作为G，突变特征数为3，参照组的指数为1:4，迭代次数为1000。localTest被假意味着全局测试在运行时被调用真正的意味着将调用局部测试。

set.seed(123) hildaGlobal <- hildaTest(inputG=G, numSig=3, localTest=FALSE, refGroup=1:4, nIter=1000)
模块GLM已加载
##编译模型图##解析未声明变量##分配节点##图信息:##观察到的随机节点:6370 ##未观察到的随机节点:1304 ##图总大小:14890 ## ##初始化模型
hildaLocal <- hildaTest(inputG=G, numSig=3, localTest=TRUE, refGroup=1:4, nIter=1000)
##编译模型图##解析未声明变量##分配节点##图信息:##观察到的随机节点:6370 ##未观察到的随机节点:1305 ##图总大小:14878 ## ##初始化模型

4.３获得签名pmsignature

该节点用于为运行MCMC采样提供一个初始值，以减少EM算法的运行时间pmsignature由Shiraishi等人开发的包装。

Param <- pmgetSignature(G, K = 3)
## #试用:1;#迭代:64;时间(s): 0.08;融合:真实;Loglikelihood: -8202.3184 ## #trial: 2;#迭代:27个;时间(s): 0.05;融合:真实;Loglikelihood: -8202.3184 ## #trial: 3;迭代:# 39; time(s): 0.04; convergence: TRUE; loglikelihood: -8202.3200 ## #trial: 4; #iteration: 21; time(s): 0.03; convergence: TRUE; loglikelihood: -8202.3187 ## #trial: 5; #iteration: 63; time(s): 0.08; convergence: TRUE; loglikelihood: -8202.3185 ## #trial: 6; #iteration: 12; time(s): 0.03; convergence: TRUE; loglikelihood: -8202.3184 ## #trial: 7; #iteration: 22; time(s): 0.03; convergence: TRUE; loglikelihood: -8202.3334 ## #trial: 8; #iteration: 22; time(s): 0.03; convergence: TRUE; loglikelihood: -8202.3184 ## #trial: 9; #iteration: 12; time(s): 0.03; convergence: TRUE; loglikelihood: -8202.3184 ## #trial: 10; #iteration: 20; time(s): 0.03; convergence: TRUE; loglikelihood: -8202.3187

4.3.1使用初始值执行全局测试和本地测试

以与运行HiLDA测试非常相似的方式，只需要指定useInits是参数由前一个函数返回，以允许在MCMC采样中使用初始值。

set.seed(123) hildaGlobal <- hildaTest(inputG=G, numSig=3, useInits = Param, localTest=TRUE, refGroup=1:4, nIter=1000)
##编译模型图##解析未声明变量##分配节点##图信息:##观察到的随机节点:6370 ##未观察到的随机节点:1305 ##图总大小:14878 ## ##初始化模型
hildaLocal <- hildaTest(inputG=G, numSig=3, useInits = Param, localTest=TRUE, refGroup=1:4, nIter=1000)
##编译模型图##解析未声明变量##分配节点##图信息:##观察到的随机节点:6370 ##未观察到的随机节点:1305 ##图总大小:14878 ## ##初始化模型

4.3.2评估MCMC链的收敛性

在MCMC采样结束后，我们可以计算潜在的尺度缩减统计量来检验两个链的收敛性。通常建议小于1.10。如果没有，可以通过增加数量来实现硝酸钠．

hildaRhat (hildaGlobal)
## [1] 1.055494
hildaRhat (hildaLocal)
## [1] 1.050305

4.4可视化pmsignature的突变签名

为了允许用户比较来自pmsignature和HiLDA的突变签名，这个函数用于绘制来自pmsignature的结果。

pmPlots <- pmBarplot(G, Param, refGroup=1:4, sigOrder=c(1,3,2))
##警告:' gather_() '在tidyr 1.2.0中已弃用。请使用“gather()”代替。此警告每8小时显示一次。调用' lifecycle::last_lifecycle_warnings() '查看此警告是在哪里生成的。
cowplot::plot_grid(pmPlots$sigPlot, pmPlots$propPlot, rel_width = c(1,3))

4．5可视化来自HiLDA的突变签名

相反，下面的函数用于绘制来自HiLDA的结果。

hildaPlots <- hildaBarplot(G, hildaLocal, refGroup=1:4, sigOrder=c(1,3,2)) cowplot::plot_grid(pmPlots$sigPlot, pmPlots$propPlot, rel_width =c(1,3))

4.6输出曝光平均差值的后验分布

为了可视化暴露的平均差异的95%可信区间，下面的函数将差异与突变特征一起绘制出来。

hildaDiffPlots <- hildaDiffPlot(G, hildaLocal, sigOrder=c(1,3,2)) cowplot::plot_grid(hildaDiffPlots$sigPlot, hildaDiffPlots$diffPlot, rel_width =c(1,3))