【基于R語言群體遺傳學】-16-中性檢驗Tajima‘s D及連鎖不平衡 linkage disequilibrium (LD)

Tajima's D Test

已經開發了幾種中性檢驗，用于識別模型假設的潛在偏差。在這里，我們將說明一種有影響力的中性檢驗，即Tajima's D（Tajima 1989）。Tajima's D通過比較數據集中的兩個𝜃 = 4N𝜇估計值來工作。我們已經推導出了𝜃is，等于平均成對雜合性（average pairwise heterozygosity），當我們討論共祖時（也稱為Tajima的估計器）。當考慮DNA序列集合中的等位基因或SNP總數以及將它們聯合到共同祖先的共祖樹內包含的預期世代數時，可以得出另一種推導。這被稱為Watterson的𝜃估計器或𝜃W（Yong 2019）。正如我們在共祖中所展示的，n個譜系的預期共祖時間是：

一組n個初始譜系到一個單一祖先的所有共祖時間的總和是：

在每次共祖事件之間的步驟中，有i + 1個譜系可能發生突變。因此，當考慮一段時間內所有譜系可能產生的等位基因數量時，我們乘以i + 1：

我們將樹上所有譜系的總時間乘以每代突變率𝜇，得到我們期望在n個DNA序列樣本中的等位基因總數，而4N是一個常數，所以我們可以將其放在求和之外

其中S是序列集合中SNP的數量。這可以重新排列以從SNP的數量估計𝜃W=4N𝜇：

請注意，Watterson的𝜃估計器需要了解一組譜系的共祖，但這是在1975年發表的，當時還沒有發表超過兩個譜系的共祖（Kingman 1982）。Tajima's D是平均成對雜合性𝜃估計值與從樣本中SNP數量估計的𝜃之間的差異，除以該差異的預期方差的平方根：

而：

這看起來有點亂，但在各種位置反復出現兩個不同的n求和，只需要計算一次然后填入。這可以用以下代碼計算（平均成對差異和S也可以從數據集中計算，但為了簡潔起見，我們在這里省略了）。?

Tajima's D是用來評估一個種群中中性突變（即沒有自然選擇影響的突變）的假設是否成立。它通過比較兩個不同的估計器來衡量種群的遺傳多樣性和種群規模的變化，我們通過R語言實現：
?

# Calculates Tajima’s D# 平均成對差異，用于估計theta_IS
theta_IS <- 2.8# 數據集中的SNP數量
S <- 16# 采樣的等位基因拷貝數
n <- 20# 初始化求和變量
i1_sum <- 0.0# 循環計算i1的和，這是Watterson's theta的一部分
for(i in 2:n-1){i1_sum <- i1_sum + 1/i
}# 計算Watterson's theta，它是基于序列多態性的一個種群規模的估計器
theta_W <- S / i1_sum# 初始化第二個求和變量
i2_sum <- 0.0# 循環計算i2的和，用于后續計算
for(i in 2:n-1){i2_sum <- i2_sum + 1/i^2
}# 計算期望值e1，它是Tajima's D公式中的項
e1 <- ((n+1)/(3*(n-1)) - 1/i1_sum) / i1_sum# 計算期望值e2，它也是Tajima's D公式中的項
e2 <- (2*(n^2+n+3)/(9*n*(n-1)) - (n+2)/(n*i1_sum) + i2_sum/i1_sum^2) / (i1_sum^2 + i2_sum)# 計算Tajima's D值，它衡量的是theta_IS和theta_W之間的標準化差異
(D <- (theta_IS - theta_W) / sqrt(e1*S + e2*S*(S-1)))

前三個變量將根據您的數據進行調整。在這個例子中，返回的D = -1.409。大于或小于2的D被認為是顯著的；然而，實際的p值是通過模擬確定的。D的正值表示中間頻率等位基因過多，這可能是由于人口減少或平衡選擇，因為這兩種情況都會延長人口歷史較老部分的共祖事件時間。在更大的種群中，有更多的祖先可供選擇，共祖是一種罕見的事件，并會膨脹𝜃IS，因為相對于𝜃W，較老的譜系在后代中以更高的頻率共享。這種負D表示稀有頻率等位基因過多（共祖的最近尖端被放大；它們對S和𝜃W的貢獻比對平均成對差異的貢獻更多，因為它們很稀有），這表明人口擴張、選擇性清除或對有害等位基因的低效凈化選擇。關鍵是要計算多個位點的D值，并尋找異常值以標記假定選擇候選者。人口統計學效應，如人口規模的變化，應該影響基因組中的所有位點，而選擇（通常被認為）在其影響上是位點特異性的。還有許多其他的中性檢驗，如HKA檢驗（Hudson等人，1987）、McDonald-Kreitman檢驗（McDonald和Kreitman，1991）、Fay和Wu的H（Fay和Wu，2000）以及dN/dS比率（Yang和Bielawski，2000）。其中許多也利用了物種之間發生的遺傳變化，它們都有各自的優點和缺點。

linkage disequilibrium (LD)

群體遺傳學的一個獨特性質，在進化博弈論等類似領域中并未發現，即不同位點甚至不同染色體上的等位基因可以“鏈接”（盡管并不總是與經典遺傳學中的重組圖譜同義），并且比隨機預期更頻繁地一起遺傳。連鎖不平衡（LD）的程度由𝒟量化，不要與Tajima的D混淆。考慮兩個位點：一個具有A/a多態性，另一個具有B/b多態性。我們對于跨位點一起遺傳的等位基因之間的關聯感到好奇。使用概率的乘法規則，我們期望AB單倍型（pAB）的頻率是兩個等位基因頻率pApB的乘積，如果它們是獨立遺傳的話。這兩者之間的差異由𝒟量化，作為連鎖不平衡的一種度量。

根據AB單倍型是過量還是不足，𝒟可以是正數或負數（或者如果你從ab單倍型任意計算𝒟，符號會改變）。𝒟也可以從所有單倍型頻率計算得出。

為了說明，假設我們有一個包含兩個SNP的單倍型頻率的小型數據集。一個是A/G多態性，另一個具有C/T等位基因：

讓我們關注A-C單倍型。A等位基因的頻率是0.57，C等位基因的頻率是0.65

遺傳漂變、種群結構和強選擇是推動𝒟偏離零的力量。在大種群中，預測𝒟隨時間呈指數衰減軌跡返回到零，就像在小種群中遺傳漂變下的雜合性一樣：

其中r是感興趣的位點對之間預期的重組分數。這可以用來估計單倍型的年齡。最后，即使對于在不同染色體上獨立分配的位點，𝒟也需要時間衰減。哈代-溫伯格基因型可以在一代中恢復，但過去的事件對LD有持續影響，這可以用來推斷更遠的過去的過程，如種群中不再存在的種群結構。當從實際數據集計算𝒟時，雙重雜合子是不明確的。假設我們有一個個體的C/T，A/G SNP集合。C等位基因與第二個位置的A還是G等位基因相關聯？通常我們不知道。但是，不明確的單倍型頻率為我們提供了關于解決雙重雜合子可能方式的信息。如果C-G單倍型非常常見，而C-A單倍型很少見，那么這表明C/T，A/G個體可能具有C-G/T-A單倍型。使用這種方法計算𝒟太繁瑣，無法手工完成。幸運的是，這正是EM算法發揮作用的地方。Kalinowski和Hedrick（2001）使用大角羊（Ovis canadensis）數據集（Boyce等人，1997）來估計LD。這個物種很罕見，樣本量很小，所以我們需要從可用的數據中獲得盡可能多的信息。以下R代碼實現了Kalinowski和Hedrick（2001）給出的方程式。它從猜測相等的單倍型頻率和𝒟 = 0開始。然后它更新這個猜測，并迅速達到最大似然解𝒟 ≈0.0779和na-B單倍型頻率基本為零。

# 定義各個復合基因型的頻率
AABB <- 2  # 兩個位點都是純合子AABB的個體數量
AaBB <- 0  # 一個位點是雜合子，另一個是純合子AaBB的個體數量
aaBB <- 0  # 一個位點是純合子aa，另一個是純合子BB的個體數量
AABb <- 0  # 第一個位點是純合子AA，第二個位點是雜合子Bb的個體數量
AaBb <- 1  # 兩個位點都是雜合子AaBb的個體數量（雙雜合子）
aaBb <- 0  # 第一個位點是純合子aa，第二個位點是雜合子Bb的個體數量
AAbb <- 1  # 第一個位點是純合子AA，第二個位點是純合子bb的個體數量
Aabb <- 0  # 第一個位點是雜合子Aa，第二個位點是純合子bb的個體數量
aabb <- 0  # 兩個位點都是純合子aabb的個體數量# 使用上述輸入運行函數‘Dcalc’
Dcalc(AABB, AaBB, aaBB, AABb, AaBb, aaBb, AAbb, Aabb, aabb)