溯祖理论(coalescent theory)：Fisher-Wright模型以及中性溯祖

参考文献：Yang Z. 2014. Molecular evolution: a statistical approach. Oxford
(England): Oxford University Press. Chapter 9.
以及其他网络资源、wiki

以下所有内容建立在locus/site内部没有重组，并且任何两个loci之间自由重组的情况下。

Fisher-Wright模型是群体遗传学中一个理想化的模型，它假设群体的数量在世代时间是稳定不变的，世代之间没有重叠、随机交配、中性进化。在一些一年生草本植物中确实如此。

a中展现的是Fisher-Wright模型的基本概念。每一列代表一个世代，每一行代表20个基因（N=10，10个个体，都是二倍体），每个时代的个体数和对应的基因数都是相同的。右图是选出的一个n=5的树，T表示的是每一个lineage所经历的世代数（coalescent waiting time）。

这个模型可能存在一些偏差，因为没有考虑到遗传漂变。如果使用有效群体大小的概念 $N_e$ 可以一定程度上进行矫正。 $Ne = 4NmNf /(Nm + Nf )$
其中 $N_m和N_f$ 分别为male和female的个体数。理解一下：在Nm和Nf组成的群体中随机取两条序列，它们coalescent的waiting time(世代数*世代间隔)大约相当于Ne/2个male个体和Ne/2个female个体组成的群体中取样的coalescent waiting time。也就是说，你两个性别对半分，那最好，说明能够自由交配，有效群体大小等于群体大小。如果不能对半分，那我就要矫正一下，因为你产生后代的数量肯定没有对半分的群体高。一般来说N指的就是Ne。

两个基因的溯祖(coalescent)

假设我们现在说的是二倍体生物，个体数为N，基因池里基因的数量为 $2N$ （后面我们说的都是基因的传代，不说个体，这个思维要转化一下），那么现存的这2N个基因中随机挑选两个，它们在上一世代来源于同一parents的概率就是 $1/2N$ 。来源于上一世代不同parents的概率就是 $1-1/2N$ 。所以，两个基因在 $i$ 个世代内都没有溯祖的概率就为： $Pr\{{T'_2>i\}=(1- \frac{1}{2N})^i}$
刚好在 $i$ 世代发生溯祖的概率为(因为这里只涉及到两个基因，所以第一次溯祖就是所有基因溯祖)： $Pr \{{T'_2=i\}=(1- \frac{1}{2N})^{i-1}} \times \frac{1}{2N}$
这里有人会问了：有性繁殖的生物（也就是群体大小为N）的基因来自父母双亲，在上一时代不可能溯祖。这里的模型确实更加适用于无性繁殖的二倍体生物（也就是群体大小为2N），但是在经过很多世代后，这一点的影响微乎其微，可以不考虑。

现在我们重新scale一下时间变量。令 $T_2=T'_2/(2N)$ ，那么 $Pr \{{T_2> \frac{i}{2N}}\}= Pr \{{T'_2=i\}=(1- \frac{1}{2N})^{i}} \approx e^{- \frac{i}{2N}}$

最后那个约等号的推导过程

这里， $\frac {i}{2N}$ 就成了一个总体的变量，所以 $T_2$ 就符合了指数分布，其均值为1，方差为1: $f(T_2)=e^{-T_2}$

复习：指数分布 $Pr\{X> \chi\} = e^{-\lambda X}$ 的均值为 $1/\lambda$

现在我们要将“世代时间” $T$ 转换成核酸替换的数量（假设我们可以通过核算替换数量来评估分化的时间），那么：
1） $t_2=T'_2\mu$ （ $\mu$ 为核酸替换式数量/每个位点/每世代）。

2）现在再引入一个变量 $\theta$ ，为种群内部差异参数(也叫群体大小参数)， $\theta =4N\mu$ ，它是指在有效群体大小为N的群体中，随机抽取两条序列，它们的平均序列差异。举个例子，人类群体中的theta大约为0.0006，含义是这个群体中随机抽取两条基因组序列，其差异平均为0.6每kb（这是一个群体遗传中非常常用的参数）。

因此 $t_2$ 的概率分布密度就为： $f(t_2)= \frac{2}{ \theta} e^{- \frac{2}{ \theta} t_2}$

所以以核酸替换数量单位的时间衡量标准的速率就为 $2/\theta$ 。

小结一下：
在目前的文献中中，主要有三种coalescent waiting time 的衡量标准：
1） $T'_2$ 衡量的是generations，均值是2N。
2） $T_2=T'_2/(2N)$ 衡量的是2N*generation（这里的2N算一个常量），均值是1.
3） $t_2 = T′_2 \mu= T_2 × 2N \mu = T_2 × \frac{ \theta}{2}$ 衡量的是每个位点上的核酸替换数量，均值是 $\theta/2$ .

image.png

关于参数 $\theta$ ：同样地，我们也可以用 $\theta$ 来计算有效群体大小。比如说通过实验抽样得到人类群体中的 $\theta$ 大约为0.0006，mutation rate $\mu$ 大约为 $\mu ≈ 2.4 × 10^{–8}$ ，则 $N=\theta/(4 \mu) \approx 6250$ 。这个数字非常有趣，因为它与人类现存的七十多亿现实群体大小差别巨大。所以遗传学家认为人类曾经经历过瓶颈时期，现在我们之间的遗传差异比较小，有效群体大小只有6250左右。

n个基因的溯祖(coalescent)

假设现在有n个基因。
和上文类似的，经过一个世代，n个基因都没有coalesce的概率为：
$(1- \frac{1}{2N})(1- \frac{2}{2N})(1- \frac{3}{2N})...(1- \frac{ {n-1}}{ 2N}) \approx 1- \frac {1+2+3...+n-1}{2N}=1- {n \choose 2} \frac{1}{2N}$

其中 ${n \choose 2}=\frac {n(n-1)}{2}$ ，表示随机取两个组成一对，有多少种可能。上面这个公式怎么理解呢。假设第一个基因确定了它的上一世代的亲本，第二个基因取到相同亲本（溯祖了）的概率为 $(1- \frac{1}{2N})$ ，第三个又相同的概率为 $(1- \frac{2}{2N})$ ，依次类推。

上面这个约等号是这样推导的：这个连乘展开来是多次的，包含1/(2N)^2， 1/(2N)^3这些高次的项，但是由于我们默认n是远小于N的，因此几乎不可能有超过两个gene在同一世代coalesce，所以这些高次的项都可以被删除，只剩下一次项加一起（我大声高呼“妙啊”!）。

和两基因一样的，将上面多基因的式子扩展到每一个generation，在 $i$ 世代刚好发生第一次基因coalesce的概率为：

$Pr\{ T'_n=i\}=\Bigg[1- {n \choose 2} \frac{1}{2N} \Bigg]^{i-1} \times {n \choose 2} \frac{1}{2N}$

同样的这也是一个几何分布(形似 (1-k)^i x k,则期望E(p)= 1/k)，其期望（均值）为 $E(T'_n)=2N/ {n \choose 2}$ ，每一对基因溯祖的概率为1/2N 每世代。

所以，假设溯祖过程中任意一段存在 $j$ 个基因，那么到下一次coalescent的世代时间平均期望为 $E(T'_j)=2N/ {j \choose 2}$ 。

同样地，令 $T_j=T'_j/(2N)$ ， $T_j$ 的概率分布为：

$f(T_j)= \frac{ j ( j-1)}{2} exp \bigg\{ -\frac{j(j-1)}{2} T_j \bigg \}$
其均值（期望）为 $2/j(j-1)$ ，方差为 $[2/j(j-1)]^2$

这里在引入一个概念：labelled history. 简单来说就是在溯祖过程中存在多少种历史可能（祖先节点的先后顺序），就有多少种labelled history。比如说((a,b),(c,d))这颗树有两种labelled histoy，因为ab祖先和cd祖先分化的先后不确定；而(((a,b),c),d)这棵树只有一种labelled histoy。

大家应该发现了，我们正在构建的这颗树是随机coalesce的，这种树称为genealogical 树，存在很多种labelled history；与此对应的rooted tree只有一种history。所以这棵树存在的history有这么多个：
$H_n= {n \choose 2}·{n-1 \choose 2}...{2 \choose 2}= \frac {n!(n-1)!}{2^{n-1}}$

怎么理解上面这个公式：在最底层世代，有n个基因，随机挑选两个coalesce的组合有 ${n \choose 2}$ 种（也就是说在最底层世代有这么多种labdelled history），倒数第二层就只剩下n-1个lineage了，以此类推。

所以现在有 $H_n$ 种history，每一种的概率都是一样的 $1/H_n$ 。对于任意一种历史 $G$ ， $Tn, T_{n–1}, . . . , T{2}$ 的分布都是相对独立的指数分布，其概率密度为：
$f(T_n,T_{n-1}...T_2|G)= \prod_{j=2}^n \frac{ j(j-1)}{2} exp \Bigg \{- \frac{j(j-1)}{2} T_j \Bigg \}$
(这里的意思是，要求第几段的coalescent waiting time，把n=几代进去就行)

联合概率分布（把历史G的分布和在历史G的情况下Tn的分布乘起来）：

$f(G,T_n,T_{n-1}...T_2)=f(G)f(T_n,T_{n-1}...T_2|G)=\prod_{j=2}^nexp \Bigg\{-\frac{j(j-1)}{2} T_j\Bigg\}$ （exp前面那个系数（累乘）正好和Hn这个分布消掉了）

现在计算树高的期望以及方差：

$T_{MRCA}=T_n+T_{n-1}+...+T_2$
$E(T_{MRCA})=E(T_n+T_{n-1}+...+T_2)= \sum_{j=2}^n \frac {2}{j(j-1)}=2 \sum_{j=2}^n (\frac {1}{j-1} - \frac {1}{j})=2(1- \frac{1}{ n})$
$V(T_{MRCA})=\sum_{j=2}^n V(T_j) = \sum_{j=2}^n \bigg[\frac {2}{j(j-1)}\bigg]^{2}=8 \sum_{j=1}^{n-1} \frac {1}{j^2}-4\bigg( 3-\frac{2}{n}-\frac{1}{n^2}\bigg)$

n比较大时 $T_{MRCA} \approx 2$ ，又因为 $E(T_2)=1$ ，所以所有基因溯祖的时间大约为最后两条lineage溯祖的时间的两倍（n足够大时）。当n很大时，
$V(T_{MRCA}) \approx \frac {8{ \pi}^2}{6} -12 \approx 1.16$
又因为 $V(T_2)=1$ ，因此 $T_{MRCA}$ 的方差主要来源于 $T_2$ 。

现在计算树长：

树长(tree length)的定义为gelealogical树内所有的支长(branch length)的总和(这里的支长听起来像是个长度，但其实可以用时间累加，因为只不过差一个替换速率的系数)：
$T_{total}=\sum_{j=2}^n jT_j$
$V_{total}= \sum_{j=2}^n j^2V(T_j)= \sum_{j=2}^nj^2( \frac{2}{j (j-1)})^2=4 \sum_{j=1}^{n-1} \frac{1}{j^2}$

当n很大时， $E(T_{total}) \approx 2( \gamma + \log n)$ ，此处欧拉常数 $\gamma = \lim_{n \to + \infty} \sum_{j=1}^n \frac{1}{j} -\log(n) \approx 0.577$ （这里太难了...无能为力，就当作看个结果）
$V(T_{total}) ≈ 2 {\pi}^2/3 ≈ 6.579$
因此当n很大并且逐渐变大时，树长 $T_{total}$ 的变化很小，而且方差几乎不变。
现在如果考虑树高MRCA的话，当我对一个群体进行取样，得到n个基因，那么n个基因计算得到的MRCA正好为整个群体的MRCA的概率为 $\frac {n-1}{n+1}$ (有论文支持)，所以即便是很小的采样（比如我只踩了三四个个体，换算成二倍体就是六八个基因），也有很高的概率能够coalesce到root节点。

这三棵树的叶很短而主干很长，由于coalescent rate等于􏰇

{n \choose 2}\frac{2}{\theta}=\frac{n(n-1)}{\theta}

，其主要是由

n^2

决定的，所以n越小coalescent rate越低，就需要花越长时间(世代)来coalesce。

知道这个树长有什么用呢？可以用来评估建树能预测多久以前的历史。比如说我们已经知道人类的有效群体大小评估大约为 $N\approx 1e4$ ，我们假设generation gap为 $g=20年(假设平均二十岁生子)$ ，那么MRCA一般不超过期望正负两个标准差： $E(T_{MRCA})+2 \sqrt{V(T_{MRCA})} \approx 2+2 \sqrt{1.16} \approx4.15(也就是说4.15*2N generations的时间)$
所以我们用DNA取样建树来重建人类的历史，信息的有效性不超过 $4.15 \times 2N \times g=4.15 \times 2*1e4 \times 20 \approx 1.7$ million years（别忘了这里的 $T$ 是用 $T'/2N$ 得来的，而 $T'$ 又是以generations为单位的）

好了，这一篇coalescent的原理部分暂时写道这里，后面希望写一篇文章通过编程语言来重现这一过程。
具体各模型之间的关系，谁优谁劣我还分不太清楚，但总有一天会汇总一下讨论，说不定写篇中文综述。这个Fissher-Wright模型使用的是Ne不变的模式，所以不太靠谱，应该会有更靠谱的方法，后面去查查。

BTW，MRCA的意思就是most recent common ancestor.

能看耐心完的人应该寥寥无几吧🤦‍♂️欢迎讨论

coalescent theory的补充资料
https://en.wikipedia.org/wiki/Coalescent_theory

延伸的还有一个模型叫平均杂合度(mean heterozygosity)： $\overline{H}$

theta就是之前说的那个theta，很重要，所以theta越大平均杂合度越高，和有效群体大小和突变速率都有关

另外补充一点：之所以称为中性溯祖，是因为这个coalescent模拟过程完全是随机发生的，认为物种形成过程是由于随机突变以及基因漂移导致的，在此基础之上反过来溯祖。它取决于世代数、有效群体大小等，但是和自然选择完全无关，完全不考虑。这一点值得深思。