之前的三篇博客，我們對(duì)于哈代溫伯格遺傳比例有了一個(gè)全面的認(rèn)識(shí)，沒有看的朋友可以先看一下前面的博客：

群體遺傳學(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ā)生的機(jī)會(huì)。

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

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

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

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

（6）Autogamic: Self-fertilizing.

自花授粉：自我受精。

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

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

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

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

Genotyoe	AA	SS	AS
number	366	12	123

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

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

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

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

為了確定我們所看到的基因型數(shù)量是否真的符合我們的預(yù)期，我們將使用內(nèi)置的R函數(shù)pchisq()來計(jì)算來自χ2分布的概率值（P值）。在pchisq()函數(shù)中，我們希望將參數(shù)lower.tail設(shè)置為FALSE，因?yàn)槲覀兿肟吹轿覀兊?span style="color:#fe2c24;">χ2值高于實(shí)際值的概率。隨著我們的觀察和預(yù)期差異越來越大，我們的χ2值應(yīng)該增加，粗略地說，得到一個(gè)非常大的χ2值的概率應(yīng)該越來越小。

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

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，則認(rèn)為一個(gè)值與預(yù)期顯著不同）。因此，這個(gè)觀察數(shù)據(jù)似乎完全符合我們從哈代-溫伯格預(yù)測中所期望的結(jié)果。 χ2檢驗(yàn)實(shí)際上是對(duì)似然比檢驗(yàn)的一種便捷近似，這種檢驗(yàn)被稱為G檢驗(yàn)或擬合優(yōu)度檢驗(yàn)，也常用于評(píng)估模型預(yù)測與實(shí)際現(xiàn)實(shí)世界數(shù)據(jù)之間的一致性：

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

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檢驗(yàn)得出的P值（0.668），與χ2檢驗(yàn)（0.664）非常相似，因此我們?cè)俅蜗喈?dāng)確信我們的觀察數(shù)據(jù)與我們的預(yù)期沒有太大差異。

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

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

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

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基因中的至少一個(gè)變體超出哈代-溫伯格條件下預(yù)期的東西；也就是說，我們沒有在每個(gè)新世代中隨機(jī)結(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"))

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

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

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