阅读 87 SEO

在R语言中实现贝叶斯-MCMC(核酸替换速率和分子距离)

主要参考资料:来自帝国理工学院的副研究员 Fabricia Nascimento的教程网站
https://thednainus.wordpress.com/2017/03/03/tutorial-bayesian-mcmc-phylogenetics-using-r/
用到的R代码都在:https://github.com/thednainus/Bayesian_tutorial
文章中有一些原数据来源于杨子恒老师的Molecular Evolution: A statistic approach一书,可以去经管之家下载。

MCMC已经被广泛应用于phylogeny、divergence time和物种界定等领域。
这篇教程的主要目标是使用贝叶斯MCMC来评估两个物种之间的分子距离(记为参数d),K80模型下核酸转换/颠换速率的比值(transition/transversion rasio; 记为参数k)。 核酸数据通过三项分布来计算。对于参数d和k施加gamma分布。

不了解贝叶斯公式和MCMC的朋友请先自行百度大致了解一下


Part1:计算似然值likelihood、先验prior和后验分布posterior distribution(un-normalized非常态化的)(为避免歧义,下面的部分术语还是用英文)

假设现在有一个人类和大猩猩的pairwise alignment,一共948bp,两者之间存在 84 个 transitions 和 6个 transversions。(为了避免数字上存在问题,下面都使用对数化来计算)

n <- 948 # length of alignment in bp
ns <- 84 # total number of transitions (23+30+10+21)
nv <- 6  # total number of transversions (1+0+2+0+2+1+0+0)

首先,需要对K80模型有个基本的认识:K80模型只有v和s两个参数,分别代表颠换率和转换率,如下表所示,应该有4种转换和8中颠换(下图(c)中右上角打错了,AG之间应该为s,而不是v)。


进化基因组学的统计理论与方法,谷迅 et al.。图(c)中右上角打错了,AG之间应该为s,而不是v。

log-likelihood函数可以表示为: (在已知距离d和转换/颠换速率比值k的情况下,观测到数据集D的概率函数)
那么对于数据 ,在K80模型中,有如下表示

其中:


别问是怎么得来的,这就是K80模型的内容,要看具体推导可以看1980年kimura的论文

与此对应的,写成R的代码应该是:

# Kimura's (1980) likelihood for two sequences
k80.lnL <- function(d, k, n=948, ns=84, nv=6){
 
  p0 <- .25 + .25 * exp(-4*d/(k+2)) + 
        .5 * exp(-2*d*(k+1)/(k+2))
 
  p1 <- .25 + .25 * exp(-4*d/(k+2)) - 
        .5 * exp(-2*d*(k+1)/(k+2))
 
  p2 <- .25 - .25 * exp(-4*d/(k+2))
 
return ((n - ns - nv) * log(p0/4) +
        ns * log(p1/4) + nv * log(p2/4))
}

# 这里return的值就是log似然值了,这样看是不是明白许多

那么进一步,根据贝叶斯公式可以推导出后验概率(在观测到D数据集的情况下,参数d和k的概率分布,也就是我们最终想求的):

这里的z是常态化常数(normalized constant),
而 和 分别是对参数d和k施加的的marginal priors。也就是分别定义d和k概率分布的函数。在这里,我们对其施加gamma分布: , ,这两个参数的prior mean分别为2/20=0.1和2/0.1=20.

常数z是多重积分,基本上涉及多少个参数就是几重积分,所以很难解。其实以上各步骤都没有涉及到MCMC,都是贝叶斯的内容,而正是因为z这个边缘似然值很难解、错误率高,所以才引入了MCMC链来避免计算z值。

所以我们暂时放弃计算z值,那么根据(1)式,未常态化(un-normalized)的后验概率就应该是:

现在就可以在R中把likelihood、prior和posterior画出来了,先设定一个参数网格:

dim <- 100  # dimension for the plot
d.v <- seq(from=0, to=0.3, len=dim) # vector of d values
##从0到0.3,长度为100。假设这两个物种最大分子距离不超过0.3,也就是序列差异不超过0.3,0-0.3都有可能。
k.v <- seq(from=0, to=100, len=dim) # vector of k values
##从0到100,长度为100。同理,假设这两个物种最大转换/颠换速率比不超过100,0-100都有可能。
dk <- expand.grid(d=d.v, k=k.v)
##画一个填满的网格
 
par(mfrow=c(1, 3))
#分面板,1列,每列三个面板

分别绘制三个面板的参数图:

# Prior surface, f(D)f(k)
Pri <- matrix(dgamma(dk$d, 2, 20) * dgamma(dk$k, 2, .1),
       ncol=dim)
##dgamma 密度gamma函数。两个gamma分布相乘,这两个分布分别就是f(d)和f(k),这个相乘的结果就是f(d)f(k),也就是先验了(prior)。
 
image(d.v, k.v, -Pri, las=1, col=heat.colors(50),
      main="Prior", xlab="distance, d",
      ylab="kappa, k", cex.main=2.0,
      cex.lab=1.5, cex.axis=1.5)
 
contour(d.v, k.v, Pri, nlev=10, drawlab=FALSE, add=TRUE)
 
# Likelihood surface, f(D|d,k)
lnL <- matrix(k80.lnL(d=dk$d, k=dk$k), ncol=dim)
##把之前设定的k80模型拿来用,输入拟定的d和k的向量

# for numerical reasons, we scale the likelihood to be 1
# at the maximum, i.e. we subtract max(lnL)
L <- exp(lnL - max(lnL))  
##把log-likelihood转化为likelihood。这里作者因为数字有点问题的原因,要把L控制在0-1之间。

 
image(d.v, k.v, -L, las=1, col=heat.colors(50),
      main="Likelihood", xlab="distance, d",
      ylab="kappa, k", cex.main=2.0,
      cex.lab=1.5, cex.axis=1.5)
 
contour(d.v, k.v, L, nlev=10,
        drawlab=FALSE, add=TRUE) # creates a contour plot
 
# Unscaled posterior surface, f(d)f(k)f(D|d,k)
Pos <- Pri * L
##f(d)f(k)就是prior,d(D|d,k)就是likelihood,两者都是之前算好了的,直接乘一起就是没有均一化/或者说常态化的后验概率(posterior)
 
image(d.v, k.v, -Pos, las=1, col=heat.colors(50),
      main="Posterior", xlab="distance, d",
      ylab="kappa, k", cex.main=2.0,
      cex.lab=1.5, cex.axis=1.5)
 
contour(d.v, k.v, Pos, nlev=10,
        drawlab=FALSE, add=TRUE)
image.png

一句话:unscaled posteriors=likelihood*priors。likelihood是模型的基础,prior影响likelihood,posterior是影响后的结果。


Part2 马尔可夫链蒙特卡罗方法 Markov chain Monte Carlo (MCMC)

尽管我们得到了posterior distribution,但是因为没有计算边缘似然常数z,所以要通过MCMC抽样来获得更加准确的posterior distribution。

MCMC算法的基本逻辑:
1、对参数d和k设置初始状态。
2、提出一个新的状态 d*提案(从一个合适的提案密度(proposal density)中得出)。
3、对这个新的状态的提案进行概率评估: 如果接受了提案( ),设定d=d*,否则d=d
4、现在d进行了一次状态检验,再对k进行一次,重复2-3.
5、 储存目前状态的d和k。
6、回到第二步继续进行状态的提案。

Gamma分布公式 :
由于 是个常数,因此在第3步评估状态的时候会被消掉,所以在下面计算prior的时候也不必考虑了。

把part1中的代码重新梳理一遍:

# function returning the logarithm of the unscaled posterior:
# f(d) * f(k) * f(D|d,k)
# by default we set the priors as: 
# f(d) = Gamma(d | 2, 20) and f(k) = Gamma(k | 2, .1)
ulnPf <- function(d, k, n=948, ns=84, nv=6, a.d=2, b.d=20, 
        a.k=2, b.k=.1) {
  # The normalizing constant in the prior densities can 
  # be ignored:
  lnpriord <- (a.d - 1)*log(d) - b.d * d
  lnpriork <- (a.k - 1)*log(k) - b.k * k
  ##分别计算d和k的log-prior
  ##这里作者写的是log,但严格按照上文gamma分布公式对应应该是ln。其实无所谓,就是前面多一个常数比例,最后MCMC链里都会消掉的。
  ## ulnPf:unnormalized posterior 

  # log-Likelihood (K80 model):
  expd1 <- exp(-4*d/(k+2));
  expd2 <- exp(-2*d*(k+1)/(k+2))
  ##这里作者解释说,因为多了一步exp函数,所以运算会比part1快一些

  p0 <- .25 + .25 * expd1 + .5 * expd2
  p1 <- .25 + .25 * expd1 - .5 * expd2
  p2 <- .25 - .25 * expd1
  lnL <- ((n - ns - nv) * log(p0/4) + 
         ns * log(p1/4) + nv * log(p2/4))
   
  # Return unnormalised posterior:
  return (lnpriord + lnpriork + lnL)
  ## 这里就是加,不是乘,因为是ln
}



终于到了写MCMC的时候了:

mcmcf <- function(init.d, init.k, N, w.d, w.k) {
  # init.d and init.k :初始d和k的状态
  # N :MCMC重复次数
  # w.d :d的sliding-window 提案的宽度
  # w.k :k的sliding-window 提案的宽度
   
  # We keep the visited states (d, k) in sample.d and sample.k 
  # for easy plotting. In practical MCMC applications, these 
  # are usually written into a file.
  sample.d <- sample.k <- numeric(N+1)
  ##这就是个计数的作用,看看d和k经历了多少个循环,先设置这么多个空,后面再填上状态
 
  # STEP 1: Initialise the MCMC chain
  ## 第一步,启动MCMC链
  d <- init.d;  sample.d[1] <- init.d
  k <- init.k;  sample.k[1] <- init.k
  ulnP <- ulnPf(d, k)
  ## 这个是上一个block写好的ulnPf函数,输入b和k,输出log-posterior
  acc.d <- 0;  acc.k <- 0 # number of acceptances
  ## 对接受状态的次数进行计数
   
  for (i in 1:N) {
    # STEP 2: Propose new state d* 
    # we use a uniform sliding window of width w with reflection
    # to propose new values d* and k* 
    # propose d* and accept/reject the proposal
    ## 第二步: 提出一个新的状态d*,下面是用dprop表示的
    ## 使用sliding-window来进行状态的评估
    dprop <- d + runif(1, -w.d/2, w.d/2)
    ## 通过runif随机生成一个数值,这个数值位于设定的滑窗之内,所以不会太大也不会太小。使新的d*状态增加这个数值。
    # reflect if dprop is negative
    if (dprop < 0) dprop <- -dprop
    ## 如果不小心小于0了,取绝对值
     
    ulnPprop <- ulnPf(dprop, k)
    ## 计算新的d*提案下的log-posterior数值
    lnalpha <- ulnPprop - ulnP
    ## 判断是否接受新状态。ulnPprop > ulnP才接受状态提案,所以lnalpha>0接受提案。
    
 
    # STEP 3: Accept or reject the proposal
    # if ru < alpha accept proposed d*:
    ## 第三部:判断是否接受状态提案
    if (lnalpha > 0 || runif(1) < exp(lnalpha)){
    ## 这个runif什么作用我没看懂,难道是允许存在一定的干扰和偏差?
      d <- dprop;  ulnP <- ulnPprop; 
      acc.d  <- acc.d + 1
    }
    # else reject and stay where we are (so that nothing needs 
    # to be done).
     
    # STEP 4: Repeat steps 2-3 for k
    # propose k* and accept/reject the proposal
    ## 第四步:对参数 k 重复同样的动作
    kprop <- k + runif(1, -w.k/2, w.k/2)
    # reflect if kprop is negative
    if (kprop < 0) kprop <- -kprop
     
    ulnPprop <- ulnPf(d, kprop)
    lnalpha <- ulnPprop - ulnP
    # if ru < alpha accept proposed k*:
    if (lnalpha > 0 || runif(1) < exp(lnalpha)){
      k <- kprop;  ulnP <- ulnPprop
      acc.k  <- acc.k + 1
    }
    # else reject and stay where we are.
     
    # STEP 5: Save chain state
    ## 第五步: 储存链的状态
    sample.d[i+1] <- d;  sample.k[i+1] <- k
  }
   
  # print out the proportion of times
  # the proposals were accepted
  print("Acceptance proportions (d, k):")
  print(c(acc.d/N, acc.k/N))
   
  # return vector of d and k visited during MCMC
   
  return (list(d=sample.d, k=sample.k))
  ##输出状态列表
}



现在来测试一下mcmc链,设置init.d=0.2(初始分子距离为0.2), init.k=20(初始转换/颠换速率比kappa为20), N=1e4(重复10000次), w.d=0.12(d的window宽为0.12), w.k(k的window宽为180):

# Test run-time:
## 可以在自己电脑上测试一下重复一万次用了多少时间,我这台Mac pro A1708大概0.13秒。
system.time(mcmcf(0.2, 20, 1e4, .12, 180)) 

# Run again and save MCMC output:
## 储存结果
dk.mcmc <- mcmcf(0.2, 20, 1e4, .12, 180)

现在进行trace plot,看看参数状态随着MCMC链的变化:

par(mfrow=c(1,3))
 
# trace plot of d
plot(dk.mcmc$d, ty='l', xlab="Iteration", ylab="d", 
     main="Trace of d")
 
# trace plot of k
plot(dk.mcmc$k, ty='l', xlab="Iteration", ylab="k", 
     main="Trace of k")
 
# We can also plot the joint sample of d and k
# (points sampled from posterior surface)
plot(dk.mcmc$d, dk.mcmc$k, pch='.', xlab="d", ylab="k", 
     main="Joint of d and k")
参数d、k和联合dk的状态变化图(trace plot)。如果再记录一个log-posterior

如图所示,MCMC链后期的参数混合比较好(图中看起来变化比较大是因为比例尺和窗口大小设置的问题)。其实可以发现,d和k的联合作图做出来的和之前我们画的prior图基本一致,从公式上来说都是 ,观测值比较符合我们预期的概率分布。右边log-likelihood也可以看出最佳的后验概率。


Part3 评估MCMC链的效率

作者语:
MCMC链是自相关的(autocorrelated),因为任何一个状态要么是前一个状态变化而来的,要么就是前一个状态没变。从直觉上我们会想,一条MCMC链的效率与其自相关程度是密不可分的,它越是自相关,效率就越低,就要花更多时间去使其收敛以得到较好的后验概率分布。

现在定义一条MCMC链的效率为:

为滞后 时间之后,两者之间的自相关系数。总而言之,这是一个时间序列的自相关函数,两次观测时间相隔越近(也是就是 越小),自相关应该越高,相隔越大(也就是 越大),自相关越低,是一个衰减的序列。

关于自相关可以来这里补课:
https://blog.csdn.net/yuting_sunshine/article/details/95317735

以下使用R语言的acf函数来计算自相关系数并绘制自相关图:

# run a very long chain (1e6 generations take about
# 40s in my MacBook Air) to calculate efficiency
dk.mcmc2 <- mcmcf(0.2, 20, 1e7, .12, 180)

# R's acf function (for AutoCorrelation Function) 
par(mfrow=c(1,2))
acf(dk.mcmc2$d)
acf(dk.mcmc2$k)
## 使用acf函数计算一组数据的自相关

# Define efficiency function
eff <- function(acf) 1 / (1 + 2 * sum(acf$acf[-1]))
## 定义efficiency函数

# the efficiencies are roughly 22% for d and 
# 20% for k:
# [1] 0.2255753 # mcmcf(0.2, 20, 1e7, .12, 180)
eff(acf(dk.mcmc2$d))
# [1] 0.2015054 # mcmcf(0.2, 20, 1e7, .12, 180)
eff(acf(dk.mcmc2$k))

跑出来长这样:


ACF自相关图

这张图表明了随着时间间隔的延长,两次计算之间的自相关性下降。下降得越快越好。再根据效率函数计算出k和d在重复1e7次后的链的效率在0.2左右。这意味着这条链的抽样模拟效果几乎相当于0.2*N个独立样本拟合后的计算结果(真的有这么靠谱吗......)。作者语:20%效率的链就算非常高效了。

这就是计算有效样本大小(effective sample size, EFF)的公式:

另一方面,由于MCMC链的自相关度很高,在处理结果的时候可以间隔取值,以去除高自相关的iteration,这样又可以节省磁盘空间又可以保留信息量。‘

效率低的链

效率低的链可能是由于状态提案密度较低,或者是由于过度设置复杂参数导致后验概率复杂。
什么情况下会出现效率低的链呢?作者聚了两个例子,我简单描述一下。
第一,比如说 和 的步长设置得不合适。比如说d作为分子距离,是在0-1之间的,而设置3就会使链堵住。k作为转换/颠换速率比,应该是比较大的,如果设置每次变动5以内,就会效率太低(作者原话叫baby-sized steps)。如此计算得到这两条链的eff大概为1.5%和0.3%。比如下图的情况:

右侧两图d与k步长设置不合理
自相关总是不为0,只有在几个点才是0。这是糟糕的链。

作者建议在调整sliding-windows步长的时候,使最后状态的接受率在30%左右。



Part4 如何判断MCMC链已经收敛

这里涉及到一个burn-in的概念,其实很简单。就是MCMC链在收敛的过程中,最初给定的初始值一般是不合适的,因此要把初始一定范围内的接受状态去除掉。

举个例子,这是分别设置较高和较低初始值的结果:

dk.mcmc.l <- mcmcf(0.01, 20, 1e4, .12, 180)
dk.mcmc.h <- mcmcf(0.4, 20, 1e4, .12, 180)

plot(dk.mcmc.l$d, xlim = c(1,200), ylim = c(0,0.4), ty = "l")
lines(dk.mcmc.h$d, col="red")

# We use the low chain to calculate the mean
# and sd of d. We could have used the high chain
# as well.
mean.d <- mean(dk.mcmc.l$d)
sd.d <- sd(dk.mcmc.l$d)
abline(h = mean.d + 2 * c(-sd.d, sd.d), lty = 2)
不同初始值的收敛情况

作者语:判断MCMC链有没有收敛的最简单方法就是多跑几条链,然后看他们的95% CI、后验均值和histogram等是不是相似。

高效率链(左)和低效率链(右)的参数k密度分布图。高效的链收敛快、两次run的重合度高,而低效链恰恰相反。

也可以计算后验均值和标准误:

# posterior means (similar for efficient chains):
mean(dk.mcmc$d); mean(dk.mcmc.b$d)
mean(dk.mcmc$k); mean(dk.mcmc.b$k)

# posterior means (not so similar for the inefficient chains):
mean(dk.mcmc3$d); mean(dk.mcmc3.b$d)
mean(dk.mcmc3$k); mean(dk.mcmc3.b$k)
## 计算后验均值

# efficient chain, standard error of the means
sqrt(1/1e4 * var(dk.mcmc$d) / 0.23) # roughly 2.5e-4
sqrt(1/1e4 * var(dk.mcmc$k) / 0.20) # roughly 0.2

# inefficient chain, standard error of the means
sqrt(1/1e4 * var(dk.mcmc3$d) / 0.015) # roughly 9.7e-4
sqrt(1/1e4 * var(dk.mcmc3$k) / 0.003) # roughly 1.6
##计算标准误,其中var()函数计算的是variance,即MCMC样本之间的差值除以样本量(抽样次数)

最后:如果模型太过复杂,一定要多跑几条链并且通过差别比较大的初始值来判断是否收敛。R语言coda包中的gelman.diag可以诊断MCMC链的收敛。

感谢Fabricia Nascimento、Mario dos Reis和Ziheng Yang老师们的教程,花了一整天,终于搞懂了

作者:杨康yang

原文链接:https://www.jianshu.com/p/b19d6f1b639b

文章分类
后端
版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容, 请发送邮件至 gxwowoo@163.com 举报,一经查实,本站将立刻删除。
相关推荐