【群體遺傳學(xué)】 π （pi）的計(jì)算

雜合度 heterozygosity

某個(gè)位點(diǎn)的第 $i$ 個(gè)等位基因的樣本頻率為 $p_i$ 赁豆，那么該位點(diǎn)所有等位基因的頻率和應(yīng)該是1。先考慮二倍體的雙等位基因相嵌，那就是 $p_1 + p_2 = 1$ 咨油。衡量單個(gè)多態(tài)位點(diǎn)變異（variation）的一個(gè)方法是計(jì)算樣本雜合度（heterozygosity）囚巴，公式如下：

$h = \frac{n}{n-1}(1-\sum{p_i^2})$

在公式中原在， $n$ 代表的是樣本中序列的數(shù)量。

$\pi$

上面這個(gè)公式是針對(duì)一個(gè)位點(diǎn)的彤叉，如果是正對(duì)一條序列的話庶柿，那其實(shí)就就是將整條序列的雜合度加起來(lái)即可。

$\pi = \sum_{j=1}^{S}h_j$

其中 $S$ 表示的是分離位點(diǎn)的數(shù)量秽浇， $h_j$ 表示的是第 $j$ 個(gè)分離位點(diǎn)的雜合度浮庐。在Wright-Fisher模型（無(wú)限位點(diǎn)的二倍體）下， $E(\pi) = \theta$ 柬焕，因此有時(shí)這個(gè)統(tǒng)計(jì)量也叫 $\theta_{\pi}$ 审残。我們需要注意的是在單態(tài)位點(diǎn)（monomorphic site）時(shí)雜合度是0。

先看這樣一個(gè)例子：

假設(shè)現(xiàn)在有4個(gè)樣本斑举，15個(gè)位點(diǎn)搅轿，但是只有6個(gè)位點(diǎn)是分離位點(diǎn)，我們先計(jì)算每個(gè)分離位點(diǎn)的雜合度：

根據(jù)公式可知富玷，對(duì)分離位點(diǎn)1（圖中的第二列序列）璧坟，有兩個(gè)等為位點(diǎn)，分別是T和C赎懦，其中T有3個(gè)雀鹃，C有1個(gè)，那么對(duì)T來(lái)說(shuō)励两，它的頻率就是0.75黎茎，對(duì)C來(lái)說(shuō)它的頻率就是0.25。根據(jù)公式可得：

$h = \frac{4}{3}\sum[1-(0.75^2 + 0.25^2) ]= 0.50$

我們以此計(jì)算就能得到其他5個(gè)分離位點(diǎn)的雜合度分別為：0.667伐蒋，0.5工三，0.667，0.5先鱼，0.5。

那么就能計(jì)算 $\pi$ 值了：

$\pi = 0.5 + 0.667 + 0.5 +0.667 + 0.5 + 0.5 = 3.33$

但是我們通常關(guān)注的是每個(gè)位點(diǎn) $\pi$ 的均值：

$\pi = \frac{3.33}{15} = 0.222$

我們將 $\pi$ 的計(jì)算進(jìn)行推廣就能得到下面這個(gè)公式：

$\pi = \frac{\sum_{i < j}k_{ij}}{n(n-1)/2}$

其中 $k_{ij}$ 表示的是第 $i$ 條序列和第 $j$ 條序列之間不同核苷酸的數(shù)量奸鬓，分母表示的是 $n$ 個(gè)序列之間進(jìn)行比較的唯一次數(shù)（非重復(fù)比較）”号希現(xiàn)在我們將這個(gè)公式應(yīng)用到上面的序列中。

現(xiàn)在是有4條序列串远，所以 $n = 4$ . 然后以此進(jìn)行比較:

第一條VS第二條:3個(gè)不同的核苷酸

第一條VS第三條:4個(gè)不同的核苷酸

第一條VS第四條:3個(gè)不同的核苷酸

第二條VS第三條:5個(gè)不同的核苷酸

第二條VS第四條:0個(gè)不同的核苷酸

第三條VS第四條:5個(gè)不同的核苷酸

所以, $\pi = \frac{3+4+3+5+0+5}{4(4-1)/2} = 3.33$

需要注意的是當(dāng)數(shù)據(jù)量很大的時(shí)候,使用公式 $\pi = \sum_{j=1}^{S}h_j$ 計(jì)算更快宏多。

正如前面說(shuō)到的儿惫，我們?cè)谟?jì)算序列之間的差異時(shí)通常是省略indel將其變成缺失值進(jìn)行處理的。當(dāng)使用公式 $\pi = \sum_{j=1}^{S}h_j$ 并且將indel變成缺失值時(shí)伸但，針對(duì)不同位點(diǎn) $n$ 是不同的肾请。使用公式 $\pi = \frac{\sum_{i < j}k_{ij}}{n(n-1)/2}$ 的話，通常會(huì)省略gap位置更胖。

比如這個(gè)例子：

如果用第一個(gè)公式铛铁，那么 $\pi = 3.49$ ，但是如果用第二個(gè)公式的話却妨， $\pi = 2.83$ 饵逐。原因是第一個(gè)公式將indel當(dāng)作缺失值進(jìn)行處理，而第二個(gè)公式將indel當(dāng)作gap直接省略了這些位點(diǎn)（哪怕是在這些位點(diǎn)并不是分離位點(diǎn)）彪标。不同的公式給出的結(jié)果也不一樣倍权，尤其是正對(duì)平均的每個(gè)位點(diǎn)時(shí)。因此捞烟，在處理基因組這種大數(shù)據(jù)時(shí)薄声，通常使用 $\pi = \sum_{j=1}^{S}h_j$ 這個(gè)公式。

我們可以把 $\pi$ 的期望方差表示成參數(shù)為 $\theta$ 的函數(shù)题画。雖然在中性進(jìn)化模型下奸柬，這個(gè)參數(shù)沒(méi)啥用??。

如果沒(méi)有重組發(fā)生的話：

$Var(\pi) = \frac{n + 1}{3(n + 1)}\theta + \frac{n(n^2 + n + 3)}{9n(n - 1)}\theta^2$

從公式可以看出婴程，和 $\pi$ 相關(guān)的方差很大廓奕，即使樣本很大時(shí)，方差也不接近于0档叔。

$\theta_W$

$\theta_W$ 通常叫 $Watterson‘s \theta桌粉。用$ \pi $表示核苷酸變異的另外一種方法是利用樣品中所有分離位點(diǎn)的數(shù)量$ S $進(jìn)行衡量，但是需要注意的是樣本量太大時(shí)會(huì)得到很大的$ S $衙四，因此需要對(duì)$ S$進(jìn)行校正：

$\theta_W = \frac{S}{a}$

$a = \sum_{i = 1}^{n-1}\frac{1}{i}$

對(duì)類似于Wright-Fisher模型處于平衡狀態(tài)且有無(wú)限突變位點(diǎn)的群體铃肯， $\theta_W$ 也是 $\theta$ 的估計(jì)量。

那么綜上：

$\theta:E(S) = \theta{a}$

將這個(gè)公式應(yīng)用到這個(gè)例子上：

$\theta_W = \frac{6}{1/1+1/2+1/3} = 3.28$

可以看到這個(gè)公式得到的結(jié)果和前面公式計(jì)算得到的3.33很接近传蹈。

還是和前面說(shuō)的一樣押逼，遇到indel不同的處理方式得到的結(jié)果不一樣：

如果將indel當(dāng)作缺失值進(jìn)行處理，那 $\theta_W = 5/(1/1+1/2+1/3) = 2.73$
如果將indel當(dāng)作gap進(jìn)行處理惦界，那 $\theta_W = 1/(1/1+1/2) = 0.667$

將這兩種不同方法得到的結(jié)果相加：

$\theta_W = 2.73 + 0.667 = 3.40$

同樣挑格，我們可以用參數(shù)為 $\theta$ 的函數(shù)來(lái)表示 $S$ 的期望方差（Wright-Fisher模型，沒(méi)有重組發(fā)生）：

$Var(s) = \sum_{i=1}^{n-1}\frac{1}{i}\theta + \sum_{i=1}^{n-1}\frac{1}{i^2}\theta^{2}$

如果是自由重組的話沾歪，就只是前半部分漂彤。

還可以從這個(gè)公式推斷出：

$Var(\theta_W) = \frac{Var(S)}{a^2}$

我們通常會(huì)看到關(guān)于 $\theta$ 的兩種估計(jì)值： $\pi$ 和 $\theta_W$ ，測(cè)序錯(cuò)誤等會(huì)造成不同的影響，因此通常需要兩個(gè)值都看挫望，還有更多的統(tǒng)計(jì)參數(shù)可以使用（如Tajima's D）立润。

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者

人面猴
序言：七十年代末，一起剝皮案震驚了整個(gè)濱河市媳板，隨后出現(xiàn)的幾起案子桑腮，更是在濱河造成了極大的恐慌，老刑警劉巖蛉幸，帶你破解...
沈念sama閱讀 218,682評(píng)論 6贊 507
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件破讨，死亡現(xiàn)場(chǎng)離奇詭異，居然都是意外死亡巨缘，警方通過(guò)查閱死者的電腦和手機(jī)添忘，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 93,277評(píng)論 3贊 395
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進(jìn)店門(mén)，熙熙樓的掌柜王于貴愁眉苦臉地迎上來(lái)若锁，“玉大人搁骑，你說(shuō)我怎么就攤上這事∮止蹋” “怎么了仲器？”我有些...
開(kāi)封第一講書(shū)人閱讀 165,083評(píng)論 0贊 355
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵，是天一觀的道長(zhǎng)仰冠。經(jīng)常有香客問(wèn)我乏冀，道長(zhǎng)，這世上最難降的妖魔是什么洋只？我笑而不...
開(kāi)封第一講書(shū)人閱讀 58,763評(píng)論 1贊 295
?港島之戀（遺憾婚禮）
正文為了忘掉前任辆沦，我火速辦了婚禮，結(jié)果婚禮上识虚，老公的妹妹穿的比我還像新娘肢扯。我一直安慰自己，他們只是感情好担锤，可當(dāng)我...
茶點(diǎn)故事閱讀 67,785評(píng)論 6贊 392
惡毒庶女頂嫁案：這布局不是一般人想出來(lái)的
文/花漫我一把揭開(kāi)白布蔚晨。她就那樣靜靜地躺著，像睡著了一般肛循。火紅的嫁衣襯著肌膚如雪铭腕。梳的紋絲不亂的頭發(fā)上，一...
開(kāi)封第一講書(shū)人閱讀 51,624評(píng)論 1贊 305
城市分裂傳說(shuō)
那天多糠，我揣著相機(jī)與錄音累舷，去河邊找鬼。笑死熬丧，一個(gè)胖子當(dāng)著我的面吹牛笋粟，可吹牛的內(nèi)容都是我干的怀挠。我是一名探鬼主播析蝴，決...
沈念sama閱讀 40,358評(píng)論 3贊 418
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開(kāi)眼害捕，長(zhǎng)吁一口氣：“原來(lái)是場(chǎng)噩夢(mèng)啊……” “哼！你這毒婦竟也來(lái)了闷畸？” 一聲冷哼從身側(cè)響起尝盼，我...
開(kāi)封第一講書(shū)人閱讀 39,261評(píng)論 0贊 276
萬(wàn)榮殺人案實(shí)錄
序言：老撾萬(wàn)榮一對(duì)情侶失蹤，失蹤者是張志新（化名）和其女友劉穎佑菩，沒(méi)想到半個(gè)月后盾沫，有當(dāng)?shù)厝嗽跇?shù)林里發(fā)現(xiàn)了一具尸體，經(jīng)...
沈念sama閱讀 45,722評(píng)論 1贊 315
?護(hù)林員之死
正文獨(dú)居荒郊野嶺守林人離奇死亡殿漠，尸身上長(zhǎng)有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點(diǎn)故事閱讀 37,900評(píng)論 3贊 336
?白月光啟示錄
正文我和宋清朗相戀三年赴精，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了。大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片绞幌。...
茶點(diǎn)故事閱讀 40,030評(píng)論 1贊 350
活死人
序言：一個(gè)原本活蹦亂跳的男人離奇死亡蕾哟，死狀恐怖，靈堂內(nèi)的尸體忽然破棺而出莲蜘，到底是詐尸還是另有隱情谭确，我是刑警寧澤，帶...
沈念sama閱讀 35,737評(píng)論 5贊 346
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布票渠，位于F島的核電站逐哈，受9級(jí)特大地震影響，放射性物質(zhì)發(fā)生泄漏问顷。R本人自食惡果不足惜昂秃，卻給世界環(huán)境...
茶點(diǎn)故事閱讀 41,360評(píng)論 3贊 330
男人毒藥：我在死后第九天來(lái)索命
文/蒙蒙一、第九天我趴在偏房一處隱蔽的房頂上張望杜窄。院中可真熱鬧肠骆，春花似錦、人聲如沸羞芍。這莊子的主人今日做“春日...
開(kāi)封第一講書(shū)人閱讀 31,941評(píng)論 0贊 22
一樁弒父案，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽(yáng)荷科。三九已至唯咬，卻和暖如春，著一層夾襖步出監(jiān)牢的瞬間畏浆，已是汗流浹背胆胰。一陣腳步聲響...
開(kāi)封第一講書(shū)人閱讀 33,057評(píng)論 1贊 270
情欲美人皮
我被黑心中介騙來(lái)泰國(guó)打工，沒(méi)想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留刻获，地道東北人蜀涨。一個(gè)月前我還...
沈念sama閱讀 48,237評(píng)論 3贊 371
代替公主和親
正文我出身青樓，卻偏偏與公主長(zhǎng)得像，于是被迫代替她去往敵國(guó)和親厚柳。傳聞我的和親對(duì)象是個(gè)殘疾皇子氧枣，可洞房花燭夜當(dāng)晚...
茶點(diǎn)故事閱讀 44,976評(píng)論 2贊 355