Home Bioinformatics Evolution

# PCMs (Phylogenetic comparative methods)

系统发育比较方法（PCM）旨在分析系统发育框架中物种性状的数据集。它们用于各种目的，包括检测选择，测量性状进化的速率，估计祖先性状，控制多物种数据的系统发育相互依赖性，估计协变量对特征的影响。

# 手写 MCMC

MCMC 用于 https://lukejharmon.github.io/pcm/chapter2_stats / 中的 2.4b 参数推断（为什么手写呢，我也觉得很奇怪）

Metropolis-Hastings 采样算法

  • 初始化时间 t=1
  • 设置𝑢的值，并初始化初始状态 u𝜃(𝑡)=𝑢
  • 重复以下的过程
    • 令 t=t+1
    • 从已知分布𝑞(𝜃∣𝜃(𝑡−1)) 中生成一个候选状态 θ(∗)
    • 计算接受概率 $$\alpha = min (1, \fracp)()\)t)h)e)t)a)(p\theta【{t-1】))}\fracq(\theta【{t-1】|\theta^*)}q|\theta^{t-1(|\theta^{t-1\|\theta^{t-1t|\theta^{t-1h|\theta^{t-1e|\theta^{t-1t|\theta^{t-1a|\theta^{t-1)})$$
    • 从均匀分布 Uniform (0,1) 生成一个随机值 a
    • 如果 a⩽α𝑎⩽𝛼，接受新生成的值：θ(t)=θ(∗)𝜃(𝑡)=𝜃(∗)；否则：θ(t)=θ(t−1)
  • 直到 t=T

基于扔硬币的角度（使用 2.4b 参数推断的例子）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#!/usr/bin/env /public/home/daijichen/anaconda3/envs/base-analysis/bin/python

import random
import math
from scipy.stats import binom
import matplotlib.pyplot as plt

T = 10000 # MCMC times
# 初始化
i = 0
sigma = 0.01 # proposal distribution
p_list = [0.0] * (T + 1) # 储存链
p0 = random.uniform(0, 1) # 均一分布抽样

while i < T:
    i += 1
    if i == 0:
        pi = p0
    pj = random.uniform(pi - sigma, pi + sigma) # 均已分布再次抽样
    prior_odds_ration = 1/1 # 均匀分布的性质，P曲线下面积为1
    proposal_density_ration = 1 # 提议分布的对称性
    likelihood_ration = binom.pmf(k=63, n=100, p=pi)/binom.pmf(k=63, n=100, p=pj) # 似然值
    accept = prior_odds_ration * proposal_density_ration * likelihood_ration

    alpha = random.uniform(0, 1) # 随机采样
    if accept > alpha:
        pass # reject
    else:
        pi = pj
    p_list[i] = pi # chain

# Brownian motion (BM)

布朗运动是 “随机游走” 模型的一个例子，因为特征值在任何时间间隔内，在方向和距离上都是随机变化的。最开始用于描述粒子在空间中的随机运动。

布朗运动作为一种模型占据主导地位的主要原因可能是它具有一些非常方便的统计特性，可以对树进行相对简单的分析和计算。

令群体中性状（trait）的平均值为 $$z$$, 或为随时间 $$t$$ 变化的函数 $$z (t)$$。可以用布朗运动过程来模拟平均性状值随时间的变化。布朗运动可以描述为，总体平均特征起始值 $$z (0)$$，在任何性状变化发生之前，在祖先种群中看到的平均性状值。演化速率 $$\sigma2$$**，决定了性状在时间中随机行走的速度。布朗运动符合正态分布，因此具有均值（Mean）0 与方差 $$variance=\sigma2t$$，即性状值在任意时间间隔内的变化都是由均值为 0 的正态分布得出的，方差与进化速率和时间长度成正 ** 比。

布朗运动的三个性质，

1.E[z(t)]=z(0)2.“行走”的每个连续间隔都是独立的3.z(t):N(z(0),σ2t)1. E[z(t)] = z(0)\\ 2. “行走”的每个连续间隔都是独立的\\ 3. z(t):N(z(0), \sigma^2t)

  1. 说明性状在任何时间 t 的期望值等于该性状在时间 0 的值。这里的期望值是指在 $$z (t)$$ 许多重复上的平均值，直观含义是布朗运动没有 “趋势”，在正负方向上都是相等的。
  2. 布朗运动是一个连续时间的过程，所以时间没有离散的 “步骤”。
  3. z(t)$$是从均值$$z(0)$$和方差$$\sigma^2$$t的正态分布中得出的，过程的总体方差是该速率乘以经过的时间。

# 遗传漂变下的布朗运动

当中性进化时，性状变化仅取决于遗传漂变。

我们假设 **（1）一个性状受到许多基因的影响，每个基因的影响都很小，并且性状的值不影响适应度（明显与抗病不符）。（2）突变是随机的，对字符的影响很小 **。考虑性状平均值为 $$z$$，群体有效种群大小为 $$Ne$$，由于没有选择，表型性状只会由于突变和遗传漂变而改变。

考虑最简单的 **“infinite alleles” model**，在这个模型下，突变是随机发生的，具有随机的表型效应。假设突变率符合均值为 0，方差为 $$\sigma_m^2$$ 的正态分布随机而得，$$mutations:N (0, \sigma_m^2)$$。模型假设等位基因的数量是如此之大，以至于同一等位基因实际上不可能发生多次突变。随着时间的推移，由于遗传漂变，种群中的等位基因的频率会发生变化。漂变和突变一起决定了平均性状随时间的动态变化。在此模型下，模拟多次后，多个种群之间具有相同的性状均值（性质一），表型不影响生存或生殖，且突变率被假设为随机与对称。

紧接着考虑模型的表型方差，$$\sigma_B2$$，在一段时间内，在许多独立的进化变化 “运行” 中，平均性状值的方差。由于模型的假设，我们可以仅仅聚焦于加性遗传方差在时间 t 的每个群体内，$$\sigma_a2$$。加性遗传方差衡量加性作用的遗传变异的总量 (即每个等位基因的贡献加在一起来预测最终表型，每个等位位点仅一次的突变与遗传漂变)。排除了显性效应与上位性效应。由于遗传漂变 (倾向于降低 $$\sigma_a2$$) 和突变 (倾向于增加 $$\sigma_a2$$)，种群中的加性遗传方差会随着时间而变化。我们可以将 $$\sigma_a^2$$ 从一代到下一代的期望值建模为：

E[σa2(t+1)]=(1−12Ne)E[σa2(t)]+σm2E[\sigma_a^2(t+1)] = (1 - \frac{1}{2N_e})E[\sigma_a^2(t)]+\sigma_m^2

其中 t 为每代的时间间隔，$$N_e$$ 为有效种群大小，$$\sigma_m^2$$ 为突变方差。$$(1 - \frac {1}{2N_e}) E [\sigma_a2 (t)]$$ 表示由于遗传漂变导致每一代 ** 加性遗传变异的减少 **，第二部分描述了由于每一代新的突变 ($$\sigma_a2$$)，加性遗传方差是如何增加的。

当假设 0 时间点的其实值为 $$\sigma_{st}^2$$，则任意时间点的期望加性方差为：

E[σa2(t)]=(1−12Ne)t[σst2−2Neσm2]+2Neσm2E[\sigma_a^2(t)] = (1 - \frac{1}{2N_e})^t[\sigma_{st}^2-2N_e\sigma_m^2]+2N_e\sigma_m^2

t 非常大时 $$(1 - \frac {1}{2N_e})^t$$ 趋近于 0，这意味着在进化的种群中，累加性遗传变异最终会在遗传漂变和新突变之间达到平衡，这样累加性遗传变异就会从一代到下一代开始停止变化，则有

lim⁡t→∞E[σa2(t)]=2Neσm2\lim\limits_{t\rightarrow\infty}E[\sigma_a^2(t)]=2N_e\sigma_m^2

因此，平衡遗传变异取决于种群大小和自发突变。

我们现在可以推导出时间 t 的群体间表型方差，$$\sigma_B2 (t)$$。假定加性方差 $$\sigma_a2 (t)$$ 不变且处于稳定的状态。在独立进化的种群中，平均性状值会彼此偏离。在时间 t 后经过一段时间后，总体间方差的期望值为：

σB2(t)=tσa2Ne\sigma_B^2(t)=t\frac{\sigma_a^2}{N_e}

带入上式后有：

σB2(t)=tσa2Ne=t2Neσm2Ne=2tσm2\sigma_B^2(t)=t\frac{\sigma_a^2}{N_e}=t\frac{2N_e\sigma_m^2}{N_e}=2t\sigma_m^2

该方程表明，两个分化种群之间的变异取决于它们分化后的两倍时间和突变输入的速率。注意，对于这个模型，种群之间的变异量与种群的起始状态和有效种群大小无关。因此，这个模型预测，长期的进化速度是由种群中新突变的供应所决定的。

尽管我们必须为这个推导做出特定的假设，Lynch 和 Hill (1986) 表明，方程 3.6 是一个普遍的结果，适用于一系列模型，甚至包括优势、连锁、非随机交配和其他过程。尽管方程 3.6 有些用处，但我们不能经常测量任何自然种群的突变方差 $$\sigma_m^2$$(但见 Turelli 1984)。相比之下，我们有时确实知道某一特定性状的遗传力。遗传力描述了一个群体中由于加性遗传效应 $$\sigma_a2 (t)$$ 引起的 ** 群体内 ** 总表型变异 $$\sigma_b2$$ 的比例（这里指狭义遗传力）。

带入则有：

lim⁡t→∞E[σa2(t)]=2Neσm2h2=σa2(t)σb2=2Neσm2σb2\lim\limits_{t\rightarrow\infty}E[\sigma_a^2(t)]=2N_e\sigma_m^2\\ h^2=\frac{\sigma_a^2(t)}{\sigma_b^2}=\frac{2N_e\sigma_m^2}{\sigma_b^2}

则有：

σm2=h2σb22Ne\sigma_m^2=\frac{h^2\sigma_b^2}{2N_e}

此处 h2 为遗传力，$$\sigma_b^2$$ 群体内总表型变异（包括种群内的所有变异来源，包括非加性遗传效应和环境效应）

则有：

σB2(t)=2tσm2=t∗h2σb2Ne\sigma_B^2(t)=2t\sigma_m^2=\frac{t*h^2\sigma_b^2}{N_e}

在经过时间 t 后，种群平均表型的期望值等于初始值（布朗运动的性质），方差与时间、遗传力和性状方差呈正相关，与有效种群大小呈负相关。

为了得出这个结果，我们不得不对新突变的正常性做出特殊的假设，这似乎是不现实的。值得注意的是，如果表型受到足够多突变的影响，中心极限定理保证种群内表型的分布将是正态的 —— 无论这些突变的潜在分布是什么。我们还必须假设性状是中性的，这是一个更可疑的假设。

# 选择下的布朗运动

Hansen-Martins models 三个模型

模型都要求选择的强度相对较弱，否则随着时间的推移，性状的遗传变异将被选择耗尽，性状进化的动态将发生变化。

模型一假设种群的进化是由于直接选择，但选择的强度和方向在一代到下一代之间随机变化。将每代选择的模型建模为均值为 0，方差为 $$\sigma_s^2$$ 的正态分布。

此时的布朗运动模型为，

σB2(t)=t(h2σb2Ne+σs2)\sigma_B^2(t)=t(\frac{h^2\sigma_b^2}{N_e} +\sigma_s^2)

在特殊情况下，选择的变异比漂移的变异大得多:

σB2(t)=σs2\sigma_B^2(t)=\sigma_s^2

也就是说，当选择 (平均而言) 比漂变强得多时，进化的速度完全由选择决定。这并不是那么牵强，因为许多研究表明，自然选择比漂移更强大，而且一代一代地在方向和幅度上都发生了变化。

# 参考

Bayesian Analyses of Comparative Data with the Ornstein–Uhlenbeck Model: Potential Pitfalls

https://lukejharmon.github.io/

Edited on