RPT方法在Hardy-Weinberg平衡檢驗中的應用

朱 凱 李 悅

Hardy-Weinberg平衡定律(Hardy-Weinberg equilibrium ,HWE)是群體遺傳學的基本法則之一。該法則指出在一個沒有發生人為選擇、外來物種遷入、基因突變等情況的隨機婚配的群體中，經若干代繁殖后各基因型頻率將由各等位基因的頻率決定。設指定的基因位點具有r種不同的等位基因類型，分別為A1,A2,…,Ar,對應的各等位基因頻率分別為p1,p2,…,pr,則該樣本中各基因型可表示為AiAj，其頻率分別為xij，其中1≤j≤i≤r,當i

即純合子基因型頻率為其基因頻率的平方，雜合子基因型頻率為構成雜合子兩基因頻率乘積的2倍[1]。為檢驗隨機樣本中基因型分布是否符合Hardy - Weinberg平衡定律，多種基于χ2分布的擬合優度檢驗方法被提出，常用的為 Pearsonχ2檢驗和似然比χ2檢驗，這些基于漸近檢驗理論的方法在樣本量較大的情形下通常較為可靠。然而，隨著分子生物學技術及微衛星標記技術的發展，具有 10 個以上等位基因的遺傳位點已屢見不鮮，在樣本量一定的情況下隨著遺傳位點的增多，必然造成各基因型頻數的減少，使得傳統的漸近算法的可靠性受到質疑[2]。此時，針對小樣本確切檢驗方法的研究與應用成為必然，Louis等提出了Hardy - Weinberg平衡條件似然函數方法，但存在著如算法運行效率低和過度排列等問題[3]。為此，筆者將RPT(randomized permutation test)方法與似然函數方法相結合，并在matlab軟件中實現這一算法。

等位基因排列的似然函數

設從某人群中隨機抽取n人并測定其某一等位基因，該基因位點具有r種不同的等位基因類型，其數據可排列成以下基因型頻數矩陣f

A1f11A2f12f22?………Arfr1fr2…frrA1A2…Ar

其中fij(1≤j≤i≤r)是基因型AiAj的觀察頻數，記

(1)

(2)

Hardy-Weinberg平衡指出在基因型概率經若干代隨機繁殖后將達到平衡，故樣本數據越接近平衡概率，則這些形成該基因型頻數矩陣的概率就會越大。對于上述基因型頻數矩陣f陣，其是否滿足Hardy-Weinberg平衡取決于基因型頻數矩陣f確切概率的大小，p=p(L(p|f)>L(p|g))(其中g是與f各等位基因數目均相同的隨機矩陣)。然而對于大樣本、多等位基因的數據，與各等位基因數目與f均相同的矩陣g數目相當龐大，無法一一枚舉進行計算。

RPT方法的算法改進

Permutation檢驗是1935年R.A.Fisher提出的思想，其思路為：設已知各基因型滿足Hardy-Weinberg平衡，求出各基因所有可能的組合情況，將各種情況對應的似然函數Li與樣本對應的似然函數Lobs進行比較，然后求出Li小于或等于臨界值Lobs的頻率，即為p值。該方法稱為exact permutation test(簡稱EPT)；當組合數據太多難以導出確切的理論分布時，可采用隨機抽取部分樣本估計近似分布。將部分樣本的似然函數Li與Lobs進行比較，做出推斷性結論，該方法稱之為random permutation test(簡稱RPT)[6]。因此采用本方法求模擬求解首要任務是對n個人的2n個基因隨機地重排，并根據重排后的基因型構建新的頻數矩陣，并求其相應的概率Li，其主要步驟如下：

1.計算現有樣本對應的似然函數Lobs(本文中Lobs=L(p|f)),稱Lobs為似然函數的臨界值。

2.隨機地對這n個人基因位點進行標號，分別記為1號和2號(可將1、2分別視為父、母體遺傳)，將這n個人對應的等位基因標號記入M1,M2兩個n維向量；

3.將M1向量隨機重排(在matlab軟件中可用permrnd命令)，重排后按順序與M2向量進行配對(可將該過程視為隨機婚配)；

4.根據以上配對重新構建基因型頻數矩陣，利用公式(2)計算出現該隨機配對對應的似然函數，將該值記為Li(本文中Li=L(p|g))；

5.重復步驟3、4(重復次數可根據實際情況設定)，得到檢驗統計量Li的經驗抽樣分布；

6.計算概率p，在符合Hardy-Weinberg平衡這一假設成立的條件下

再根據p值做出相應的推斷。

實例應用與結果分析

本文就多等位基因Hardy-Weinberg平衡條件的RPT似然函數算法，分別使用本研究算法和Pearsonχ2檢驗和似然比χ2檢驗，并對三種方法進行對比。表1是1997年Le Coutre等對德國27名高歇氏病患者的七個致病基因數據[7]。

表1 27名高歇氏病患者致病基因位點各基因型數據

使用本文所介紹的RPT似然函數法在不同的模擬次數n=500000時對其進行Hardy-Weinberg平衡檢驗，其結果P=0.00232。本文對該數據進行了Peaersonχ2檢驗和似然比χ2檢驗，三種方法所得結果的對比見表2。

表2 27名高歇氏病患者致病基因位點各基因型數據三種方法結果

由于該樣本各基因型理論頻數和實際頻數均較小，不滿足χ2檢驗的條件，在本例中兩種χ2檢驗之間出現了截然相反的結論(α=0.05)，但本文提出的RPT似然函數法P=0.00232，較兩種經典方法的P值小得多，這充分說了似然函數方法的靈敏性。因此，此次檢驗的最終結論應為不能認為該數據滿足Hardy-Weinberg平衡。

討 論

在等位基因類型較多，而樣本例數較小時，可使用RPT似然函數法對數據是否滿足Hardy-Weinberg平衡條件進行統計推斷[4-6]。在應用RPT方法進行Hardy-Weinberg平衡檢驗時需注意如下幾個問題:

1.確切概率法作為一種非參數檢驗，由于不依賴樣本數據量及統計量的分布類型，從而避免了復雜的理論推導，較經典方法適用范圍更加廣泛；

2.在求似然函數的過程中，會反復調用階乘函數，若數據量太大有可能出現超出計算機識別范圍，從而出現溢值現象，故本方法僅適用于樣本量不是太大的情況(fij≤100(1≤j≤i≤r))；

3.本方法所使用的程序中模擬次數、不同基因型人數、基因位點種類這些參數均可以改變，一定程度上保證了該方法的應用范圍；

4.模擬誤差來源于Monte Carlo模擬抽樣。理論上，無限次的模擬將會完全消除模擬誤差，但顯然這是不可能也沒有必要的。因此確定模擬誤差并且計算可行的RPT抽樣次數是有效控制模擬誤差的必要步驟[7]；

5.由于程序運行中要用到隨機數，最終的結果會有一定的誤差，因此建議反復運行程序，必要時增加模擬次數，待結果相對穩定時再下結論[7-8]。

由于本研究采用的是基于確切概率的抽樣方法，故本方法所得結果應優于近似分布的擬合優度檢驗方法。同時，RPT作為EPT的一種近似方法，具有使用方便、誤差小、執行效率高的優點，在樣本量較小的情況下處理能力較強，是一種有效提高檢驗效率的好方法。此外，對于多等位基因，由于組合數目十分巨大，本方法可大幅度提高計算效率。

(歡迎生物統計愛好者來信交流統計方法或索要該方法計算機程序)

參 考 文 獻

1.李照海，覃紅，張洪.遺傳學中的統計方法.科學出版社,2006,4-11.

2.Cazeneuve C,Sarkisian T.MEFV-gene analysis in American patients with familial Mediterranean fever:Diagnostic value and unfavorable renal prognosis of the M694V homozygous genetype-Genetic and therapeutic implications .American Journal of Human Genetics,1999,65:88-97.

3.Louis EJ,Dempster ER.An exact test for Hardy-Weinberg.Biometrics,1976,32:183-186．

4.韓宏，王彤.Hardy-Weinberg平衡條件確切檢驗的完全排列算法研究與應用.中國衛生統計,2007,24(5):472-475.

5.Sun Wei Guo,Elizabeth A.thompson.Performing the exact test of Hardy-Weinberg proportion for multiple alleles.Biometrics,48,361-372.

6.荀鵬程,趙楊,柏建嶺等.Permutation Test在假設檢驗中的應用.數理統計與管理，2006,26(5):616.

7.Le Coutre P,Demina A.Molecular analysis of Gaucer disease:distribuion of eight mutations and the complete gene deletion in 27 patients from Germany.Hum Genet,1997,99:816-821.

8.朱凱,李悅.RPT對秩和檢驗的改進及Matlab實現.中國衛生統計,2012,29(4):597，602.