之前的三篇博客，我們對于哈代溫伯格遺傳比例有了一個全面的認識，沒有看的朋友可以先看一下前面的博客：

群體遺傳學(xué)_tRNA做科研的博客-CSDN博客

1.一些新名詞

（1）Algorithm: A series of operations executed in a specific order.

算法：按照特定順序執(zhí)行的一系列操作。

（2）Probability: The chance of an occurrence given repeated attempts.

概率：在重復(fù)嘗試中發(fā)生的可能性。

（3）Likelihood: The chance of an occurrence given a model assumption.

可能性：在給定模型假設(shè)下發(fā)生的機會。

（4）Machine Learning: A process where computational results are validated to improve accuracy.

機器學(xué)習(xí)：驗證計算結(jié)果以提高準(zhǔn)確性的過程。

（5）Parthenogenesis: Development of an embryo without fertilization.

單性生殖：胚胎未經(jīng)受精而發(fā)育。

（6）Autogamic: Self-fertilizing.

自花授粉：自我受精。

2.期望的偏差（Deviation from exception）

在這一點上，我們已經(jīng)花費了大量時間來驗證我們對哈代-溫伯格假設(shè)下的等位基因頻率的期望是否合理。我們可以做出的最有趣的觀察之一是，某個種群正在違背我們的期望，當(dāng)這種情況發(fā)生時，我們就可以開始探索其他可能性。在這個探索中，我們可以使用的一個特定的統(tǒng)計工具叫做χ2（卡方）統(tǒng)計檢驗。我們可以使用這個檢驗來看我們的觀察到的基因型頻率是否真的偏離了基于哈代-溫伯格預(yù)測的期望。關(guān)于R語言統(tǒng)計相關(guān)的知識，可以看我寫的博客：

【R語言從0到精通】-3-R統(tǒng)計分析（列聯(lián)表、獨立性檢驗、相關(guān)性檢驗、t檢驗）_r 列聯(lián)表分析-CSDN博客

我們通過一個真實的數(shù)據(jù)集，該數(shù)據(jù)集包含了來自尼日利亞拉各斯501人樣本的基因型計數(shù)（Taiwo等人，2011年）。這些是產(chǎn)生血紅蛋白的基因的基因型，該基因與鐮狀細胞貧血癥（血紅蛋白S）相關(guān)。首先，我們計算等位基因和預(yù)期基因型頻率，然后我們可以對這些數(shù)據(jù)進行χ2檢驗。 首先，將每個觀察到的基因型計數(shù)保存為它們自己的變量。我們將使用AA表示純合子非鐮刀基因型，SS表示純合子鐮刀等位基因基因型，AS表示雜合子，三者的和就是總?cè)藬?shù)N。

Genotyoe	AA	SS	AS
number	366	12	123

我們計算鐮刀型等位基因S的等位基因頻率：

AA <- 366
AS <- 123
SS <- 12
N <- AA + AS + SS
p <- (SS + (AS/2))/N
p

根據(jù)觀察到的等位基因頻率p，我們現(xiàn)在可以計算預(yù)期的基因型。因為我們想要追蹤兩種不同的等位基因（S和A），因此有兩種不同的純合性，我們將SS純合子定義為p2，而AA純合子定義為（1-p）2。這里的含義是，只有兩種可能的等位基因：S的頻率為p，A的頻率為非p的所有部分。 現(xiàn)在，通過將我們計算出的基因型頻率乘以實際抽樣個體的數(shù)量，我們可以得到我們預(yù)期的個體基因型數(shù)量：

ExpAA <- N*(1-p)^2
ExpAS <- N*2*p*(1-p)
ExpSS <- N*p^2

為了確定我們所看到的基因型數(shù)量是否真的符合我們的預(yù)期，我們將使用內(nèi)置的R函數(shù)pchisq()來計算來自χ2分布的概率值（P值）。在pchisq()函數(shù)中，我們希望將參數(shù)lower.tail設(shè)置為FALSE，因為我們想看到我們的χ2值高于實際值的概率。隨著我們的觀察和預(yù)期差異越來越大，我們的χ2值應(yīng)該增加，粗略地說，得到一個非常大的χ2值的概率應(yīng)該越來越小。

其中E是預(yù)期的計數(shù)數(shù)量，O是觀察到的計數(shù)數(shù)量，這會在所有類別上進行求和。我們希望找到這個χ2統(tǒng)計量在分布中的位置，但為了有一個合適的分布，我們必須告訴函數(shù)考慮多少自由度（df）來進行測試。一般來說，我們在計算自由度時，從數(shù)據(jù)的類別數(shù)減一開始，所以在這個例子中有三個類別（ExpAA、ExpAS和ExpSS）減一。然而，我們還必須從觀測數(shù)據(jù)中估計一個參數(shù)p，以生成每個類別的預(yù)期值。這意味著我們又失去了一個自由度。因此，df = 3 - 2 = 1（通過從觀察數(shù)據(jù)中估計參數(shù)，預(yù)期數(shù)值“擬合”觀察數(shù)據(jù)更緊密，所以這是有代價的）

chi2 <- (ExpAA-AA)^2/ExpAA+(ExpAS-AS)^2/ExpAS+(ExpSS-SS)^2/ExpSS
pvalue <- pchisq(chi2, df = 1, lower.tail = FALSE)
chi2
pvalue

結(jié)果得到的P值（0.664>0.5），這表明我們的觀察值與預(yù)期值相當(dāng)一致（如果P值小于0.05，則認為一個值與預(yù)期顯著不同）。因此，這個觀察數(shù)據(jù)似乎完全符合我們從哈代-溫伯格預(yù)測中所期望的結(jié)果。 χ2檢驗實際上是對似然比檢驗的一種便捷近似，這種檢驗被稱為G檢驗或擬合優(yōu)度檢驗，也常用于評估模型預(yù)測與實際現(xiàn)實世界數(shù)據(jù)之間的一致性：

這種方法，顧名思義，關(guān)注的是我們的觀察值與預(yù)期值的似然比。這種G檢驗方法使用與χ2檢驗相同的分布，并且表現(xiàn)相似。χ2檢驗通常被教授而不是G檢驗，因為它不需要你計算對數(shù)值；我們進行稍微更簡化的G檢驗統(tǒng)計量的計算：

geno <- c(AA, AS, SS)
expe <- c(ExpAA, ExpAS, ExpSS)
G <- 2 * sum(geno * log(geno/expe))
pvalue <- pchisq(G, df = 1, lower.tail = FALSE)
G
pvalue

G檢驗得出的P值（0.668），與χ2檢驗（0.664）非常相似，因此我們再次相當(dāng)確信我們的觀察數(shù)據(jù)與我們的預(yù)期沒有太大差異。

如果你還記得前一章的我們說哈代-溫伯格預(yù)測的必要條件之一是沒有任何遺傳變異受到自然選擇的影響。在這里，我們處理的等位基因?qū)δ骋槐硇陀兄卮笥绊?#xff0c;例如在純合子時導(dǎo)致鐮刀型貧血，在雜合子時賦予抗瘧疾能力，這明顯違反了這一假設(shè)（Luzzatto 2012）。但是，正如我們從剛才分析的血紅蛋白S數(shù)據(jù)中看到的，這些假設(shè)經(jīng)常被違反，然而與哈代-溫伯格預(yù)期的偏離可能看起來非常小。我們看一個違法的例子：

哈代-溫伯格假設(shè)之一是每一代配子的有效隨機結(jié)合，無論潛在的等位基因頻率如何。這在克隆物種中嚴重破壞，其中一個親本產(chǎn)生一個與自己基因相同的后代。水蚤就是這樣一種物種，雌性通常通過孤雌生殖（未受精的卵發(fā)育成胚胎）繁殖，一些種群甚至必須進行孤雌生殖（Paland等人，2005）

讓我們來看一個例子，采集的118只水蚤個體，關(guān)于磷酸葡萄糖異構(gòu)酶（PGI）的兩個等位基因，我們再次稱之為“A”和“S”，以便我們可以重用之前的代碼（Hebert和Crease 1983）。發(fā)現(xiàn)了100個AS雜合子和34個AA純合子，而SS純合子在樣本中完全缺失。我們再次進行卡方檢驗：?

AA <- 34
AS <- 100
SS <- 0
N <- AA + AS + SS
p <- (SS + (AS/2))/N
p
ExpAA <- N*(1-p)^2
ExpAS <- N*2*p*(1-p)
ExpSS <- N*p^2chi2 <- (ExpAA-AA)^2/ExpAA+(ExpAS-AS)^2/ExpAS+(ExpSS-SS)^2/ExpSS
pvalue <- pchisq(chi2, df = 1, lower.tail = FALSE)
chi2
pvalue

我們得到新的p值：

我們可以看到，這與我們的預(yù)期有很大的偏差，P值為5.56×10?12，我們可以得出結(jié)論，PGI基因中的至少一個變體超出哈代-溫伯格條件下預(yù)期的東西；也就是說，我們沒有在每個新世代中隨機結(jié)合配子。 我們可以使用R函數(shù)barplot()將這些數(shù)據(jù)可視化為條形圖。

dat <- matrix(c(geno,expe), nrow = 2, byrow = T)
barplot(dat,beside=T,col=c("turquoise4", "sienna1"),names.arg=c("AA", "SA", "SS"))
legend(x="topright", legend=c("Observed","Expected"),pch=15, col=c("turquoise4","sienna1"))

在處理較小的樣本量時，考慮使用替代的檢驗方法可能更為合適，例如“精確檢驗”（exact test）。在精確檢驗中，會使用所有可能的等位基因基因型配置來為觀察到的配置分配一個P值。關(guān)于這一背景下的精確檢驗的進一步討論，可以參考Guo和Thompson（1992年）、Wigginton等人（2005年）、Engels（2009年）以及其中的參考文獻。然而，一般來說，如果樣本量足夠大，能夠檢驗感興趣效應(yīng)的大小，并且不過分關(guān)注接近顯著性截斷邊界的結(jié)果（例如，Johnson 1999年），這些不同的統(tǒng)計方法在最終解釋上將會是一致的。

原書內(nèi)容寫的有點不清晰，很多地方重復(fù)冗余，我進行提煉總結(jié)，許多R語言的錯誤我也進行了糾正，如果有什么問題，歡迎大家進行討論。

下一篇博客我們將不只討論兩個等位基因的情況，而是進行一些拓展，下個博客見！