在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）。

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)

一句话：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")

如图所示，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))

跑出来长这样：

这张图表明了随着时间间隔的延长，两次计算之间的自相关性下降。下降得越快越好。再根据效率函数计算出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%。比如下图的情况：

作者建议在调整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等是不是相似。

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

# 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老师们的教程，花了一整天，终于搞懂了