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

參考文獻(xiàn)：Yang Z. 2014. Molecular evolution: a statistical approach. Oxford
(England): Oxford University Press. Chapter 9.
以及其他網(wǎng)絡(luò)資源、wiki

以下所有內(nèi)容建立在locus/site內(nèi)部沒有重組，并且任何兩個(gè)loci之間自由重組的情況下掷酗。

Fisher-Wright模型是群體遺傳學(xué)中一個(gè)理想化的模型调违，它假設(shè)群體的數(shù)量在世代時(shí)間是穩(wěn)定不變的，世代之間沒有重疊泻轰、隨機(jī)交配技肩、中性進(jìn)化。在一些一年生草本植物中確實(shí)如此浮声。

a中展現(xiàn)的是Fisher-Wright模型的基本概念虚婿。每一列代表一個(gè)世代，每一行代表20個(gè)基因（N=10泳挥，10個(gè)個(gè)體然痊，都是二倍體），每個(gè)時(shí)代的個(gè)體數(shù)和對應(yīng)的基因數(shù)都是相同的屉符。右圖是選出的一個(gè)n=5的樹剧浸，T表示的是每一個(gè)lineage所經(jīng)歷的世代數(shù)（coalescent waiting time）。

這個(gè)模型可能存在一些偏差矗钟，因?yàn)闆]有考慮到遺傳漂變唆香。如果使用有效群體大小的概念 $N_e$ 可以一定程度上進(jìn)行矯正。 $Ne = 4NmNf /(Nm + Nf )$
其中 $N_m和N_f$ 分別為male和female的個(gè)體數(shù)吨艇。理解一下：在Nm和Nf組成的群體中隨機(jī)取兩條序列躬它，它們coalescent的waiting time(世代數(shù)*世代間隔)大約相當(dāng)于Ne/2個(gè)male個(gè)體和Ne/2個(gè)female個(gè)體組成的群體中取樣的coalescent waiting time。也就是說东涡，你兩個(gè)性別對半分冯吓，那最好，說明能夠自由交配疮跑，有效群體大小等于群體大小组贺。如果不能對半分，那我就要矯正一下祸挪，因?yàn)槟惝a(chǎn)生后代的數(shù)量肯定沒有對半分的群體高锣披。一般來說N指的就是Ne。

兩個(gè)基因的溯祖(coalescent)

假設(shè)我們現(xiàn)在說的是二倍體生物贿条，個(gè)體數(shù)為N，基因池里基因的數(shù)量為 $2N$ （后面我們說的都是基因的傳代增热，不說個(gè)體整以，這個(gè)思維要轉(zhuǎn)化一下），那么現(xiàn)存的這2N個(gè)基因中隨機(jī)挑選兩個(gè)峻仇，它們在上一世代來源于同一parents的概率就是 $1/2N$ 公黑。來源于上一世代不同parents的概率就是 $1-1/2N$ 。所以，兩個(gè)基因在 $i$ 個(gè)世代內(nèi)都沒有溯祖的概率就為： $Pr\{{T'_2>i\}=(1- \frac{1}{2N})^i}$
剛好在 $i$ 世代發(fā)生溯祖的概率為(因?yàn)檫@里只涉及到兩個(gè)基因凡蚜，所以第一次溯祖就是所有基因溯祖)： $Pr \{{T'_2=i\}=(1- \frac{1}{2N})^{i-1}} \times \frac{1}{2N}$
這里有人會問了：有性繁殖的生物（也就是群體大小為N）的基因來自父母雙親人断，在上一時(shí)代不可能溯祖。這里的模型確實(shí)更加適用于無性繁殖的二倍體生物（也就是群體大小為2N）朝蜘，但是在經(jīng)過很多世代后恶迈，這一點(diǎn)的影響微乎其微，可以不考慮谱醇。

現(xiàn)在我們重新scale一下時(shí)間變量暇仲。令 $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}}$

最后那個(gè)約等號的推導(dǎo)過程

這里副渴， $\frac {i}{2N}$ 就成了一個(gè)總體的變量奈附，所以 $T_2$ 就符合了指數(shù)分布，其均值為1煮剧，方差為1: $f(T_2)=e^{-T_2}$

復(fù)習(xí)：指數(shù)分布 $Pr\{X> \chi\} = e^{-\lambda X}$ 的均值為 $1/\lambda$

現(xiàn)在我們要將“世代時(shí)間” $T$ 轉(zhuǎn)換成核酸替換的數(shù)量（假設(shè)我們可以通過核算替換數(shù)量來評估分化的時(shí)間）斥滤，那么：
1） $t_2=T'_2\mu$ （ $\mu$ 為核酸替換式數(shù)量/每個(gè)位點(diǎn)/每世代）。

2）現(xiàn)在再引入一個(gè)變量 $\theta$ 勉盅，為種群內(nèi)部差異參數(shù)(也叫群體大小參數(shù))中跌， $\theta =4N\mu$ ，它是指在有效群體大小為N的群體中菇篡，隨機(jī)抽取兩條序列漩符，它們的平均序列差異。舉個(gè)例子驱还，人類群體中的theta大約為0.0006嗜暴，含義是這個(gè)群體中隨機(jī)抽取兩條基因組序列，其差異平均為0.6每kb（這是一個(gè)群體遺傳中非常常用的參數(shù)）议蟆。

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

所以以核酸替換數(shù)量單位的時(shí)間衡量標(biāo)準(zhǔn)的速率就為 $2/\theta$ 闷沥。

小結(jié)一下：
在目前的文獻(xiàn)中中，主要有三種coalescent waiting time 的衡量標(biāo)準(zhǔn)：
1） $T'_2$ 衡量的是generations咐容，均值是2N舆逃。
2） $T_2=T'_2/(2N)$ 衡量的是2N*generation（這里的2N算一個(gè)常量），均值是1.
3） $t_2 = T′_2 \mu= T_2 × 2N \mu = T_2 × \frac{ \theta}{2}$ 衡量的是每個(gè)位點(diǎn)上的核酸替換數(shù)量戳粒，均值是 $\theta/2$ .

image.png

關(guān)于參數(shù) $\theta$ ：同樣地路狮，我們也可以用 $\theta$ 來計(jì)算有效群體大小。比如說通過實(shí)驗(yàn)抽樣得到人類群體中的 $\theta$ 大約為0.0006蔚约，mutation rate $\mu$ 大約為 $\mu ≈ 2.4 × 10^{–8}$ 奄妨，則 $N=\theta/(4 \mu) \approx 6250$ 。這個(gè)數(shù)字非常有趣苹祟，因?yàn)樗c人類現(xiàn)存的七十多億現(xiàn)實(shí)群體大小差別巨大砸抛。所以遺傳學(xué)家認(rèn)為人類曾經(jīng)經(jīng)歷過瓶頸時(shí)期评雌，現(xiàn)在我們之間的遺傳差異比較小，有效群體大小只有6250左右直焙。

n個(gè)基因的溯祖(coalescent)

假設(shè)現(xiàn)在有n個(gè)基因景东。
和上文類似的，經(jīng)過一個(gè)世代奔誓，n個(gè)基因都沒有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}$ 斤吐，表示隨機(jī)取兩個(gè)組成一對，有多少種可能丝里。上面這個(gè)公式怎么理解呢曲初。假設(shè)第一個(gè)基因確定了它的上一世代的親本，第二個(gè)基因取到相同親本（溯祖了）的概率為 $(1- \frac{1}{2N})$ 杯聚，第三個(gè)又相同的概率為 $(1- \frac{2}{2N})$ 臼婆，依次類推。

上面這個(gè)約等號是這樣推導(dǎo)的：這個(gè)連乘展開來是多次的幌绍，包含1/(2N)^2颁褂， 1/(2N)^3這些高次的項(xiàng)，但是由于我們默認(rèn)n是遠(yuǎn)小于N的傀广，因此幾乎不可能有超過兩個(gè)gene在同一世代coalesce颁独，所以這些高次的項(xiàng)都可以被刪除，只剩下一次項(xiàng)加一起（我大聲高呼“妙啊”!）伪冰。

和兩基因一樣的誓酒，將上面多基因的式子擴(kuò)展到每一個(gè)generation，在 $i$ 世代剛好發(fā)生第一次基因coalesce的概率為：

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

同樣的這也是一個(gè)幾何分布(形似 (1-k)^i x k,則期望E(p)= 1/k)贮聂，其期望（均值）為 $E(T'_n)=2N/ {n \choose 2}$ 靠柑，每一對基因溯祖的概率為1/2N 每世代。

所以吓懈，假設(shè)溯祖過程中任意一段存在 $j$ 個(gè)基因歼冰，那么到下一次coalescent的世代時(shí)間平均期望為 $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$

這里在引入一個(gè)概念：labelled history. 簡單來說就是在溯祖過程中存在多少種歷史可能（祖先節(jié)點(diǎn)的先后順序）甘穿，就有多少種labelled history腮恩。比如說((a,b),(c,d))這顆樹有兩種labelled histoy，因?yàn)閍b祖先和cd祖先分化的先后不確定扒磁；而(((a,b),c),d)這棵樹只有一種labelled histoy庆揪。

大家應(yīng)該發(fā)現(xiàn)了，我們正在構(gòu)建的這顆樹是隨機(jī)coalesce的妨托，這種樹稱為genealogical 樹缸榛，存在很多種labelled history；與此對應(yīng)的rooted tree只有一種history兰伤。所以這棵樹存在的history有這么多個(gè)：
$H_n= {n \choose 2}·{n-1 \choose 2}...{2 \choose 2}= \frac {n!(n-1)!}{2^{n-1}}$

怎么理解上面這個(gè)公式：在最底層世代内颗，有n個(gè)基因，隨機(jī)挑選兩個(gè)coalesce的組合有 ${n \choose 2}$ 種（也就是說在最底層世代有這么多種labdelled history）敦腔，倒數(shù)第二層就只剩下n-1個(gè)lineage了均澳，以此類推。

所以現(xiàn)在有 $H_n$ 種history符衔，每一種的概率都是一樣的 $1/H_n$ 找前。對于任意一種歷史 $G$ ， $Tn, T_{n–1}, . . . , T{2}$ 的分布都是相對獨(dú)立的指數(shù)分布判族，其概率密度為：
$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=幾代進(jìn)去就行)

聯(liá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前面那個(gè)系數(shù)（累乘）正好和Hn這個(gè)分布消掉了）

現(xiàn)在計(jì)算樹高的期望以及方差：

$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比較大時(shí) $T_{MRCA} \approx 2$ 形帮，又因?yàn)?img class="math-inline" src="https://math.jianshu.com/math?formula=E(T_2)%3D1" alt="E(T_2)=1" mathimg="1">槽惫，所以所有基因溯祖的時(shí)間大約為最后兩條lineage溯祖的時(shí)間的兩倍（n足夠大時(shí)）。當(dāng)n很大時(shí)辩撑，
$V(T_{MRCA}) \approx \frac {8{ \pi}^2}{6} -12 \approx 1.16$
又因?yàn)?img class="math-inline" src="https://math.jianshu.com/math?formula=V(T_2)%3D1" alt="V(T_2)=1" mathimg="1">界斜，因此 $T_{MRCA}$ 的方差主要來源于 $T_2$ 。

現(xiàn)在計(jì)算樹長：

樹長(tree length)的定義為gelealogical樹內(nèi)所有的支長(branch length)的總和(這里的支長聽起來像是個(gè)長度合冀，但其實(shí)可以用時(shí)間累加各薇，因?yàn)橹徊贿^差一個(gè)替換速率的系數(shù))：
$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}$

當(dāng)n很大時(shí)， $E(T_{total}) \approx 2( \gamma + \log n)$ 君躺，此處歐拉常數(shù) $\gamma = \lim_{n \to + \infty} \sum_{j=1}^n \frac{1}{j} -\log(n) \approx 0.577$ （這里太難了...無能為力峭判，就當(dāng)作看個(gè)結(jié)果）
$V(T_{total}) ≈ 2 {\pi}^2/3 ≈ 6.579$
因此當(dāng)n很大并且逐漸變大時(shí)，樹長 $T_{total}$ 的變化很小晰洒，而且方差幾乎不變朝抖。
現(xiàn)在如果考慮樹高M(jìn)RCA的話，當(dāng)我對一個(gè)群體進(jìn)行取樣谍珊，得到n個(gè)基因治宣，那么n個(gè)基因計(jì)算得到的MRCA正好為整個(gè)群體的MRCA的概率為 $\frac {n-1}{n+1}$ (有論文支持)，所以即便是很小的采樣（比如我只踩了三四個(gè)個(gè)體砌滞，換算成二倍體就是六八個(gè)基因）侮邀，也有很高的概率能夠coalesce到root節(jié)點(diǎn)。

這三棵樹的葉很短而主干很長贝润，由于coalescent rate等于??

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

绊茧，其主要是由

n^2

決定的，所以n越小coalescent rate越低打掘，就需要花越長時(shí)間(世代)來coalesce华畏。

知道這個(gè)樹長有什么用呢鹏秋？可以用來評估建樹能預(yù)測多久以前的歷史。比如說我們已經(jīng)知道人類的有效群體大小評估大約為 $N\approx 1e4$ 亡笑，我們假設(shè)generation gap為 $g=20年(假設(shè)平均二十歲生子)$ 侣夷，那么MRCA一般不超過期望正負(fù)兩個(gè)標(biāo)準(zhǔn)差： $E(T_{MRCA})+2 \sqrt{V(T_{MRCA})} \approx 2+2 \sqrt{1.16} \approx4.15(也就是說4.15*2N generations的時(shí)間)$
所以我們用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的原理部分暫時(shí)寫道這里，后面希望寫一篇文章通過編程語言來重現(xiàn)這一過程晰甚。
具體各模型之間的關(guān)系衙传，誰優(yōu)誰劣我還分不太清楚，但總有一天會匯總一下討論厕九，說不定寫篇中文綜述蓖捶。這個(gè)Fissher-Wright模型使用的是Ne不變的模式，所以不太靠譜止剖，應(yīng)該會有更靠譜的方法腺阳，后面去查查。

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

能看耐心完的人應(yīng)該寥寥無幾吧???♂?歡迎討論

coalescent theory的補(bǔ)充資料
https://en.wikipedia.org/wiki/Coalescent_theory

延伸的還有一個(gè)模型叫平均雜合度(mean heterozygosity)： $\overline{H}$

theta就是之前說的那個(gè)theta亭引，很重要，所以theta越大平均雜合度越高皮获，和有效群體大小和突變速率都有關(guān)

另外補(bǔ)充一點(diǎn)：之所以稱為中性溯祖焙蚓，是因?yàn)檫@個(gè)coalescent模擬過程完全是隨機(jī)發(fā)生的，認(rèn)為物種形成過程是由于隨機(jī)突變以及基因漂移導(dǎo)致的洒宝，在此基礎(chǔ)之上反過來溯祖购公。它取決于世代數(shù)、有效群體大小等雁歌，但是和自然選擇完全無關(guān)宏浩，完全不考慮。這一點(diǎn)值得深思靠瞎。