家養(yǎng)動物數(shù)量性狀全基因組關(guān)聯(lián)研究進(jìn)展
——單標(biāo)記回歸分析策略

郭家中,王小龍,劉小林*

(1.四川農(nóng)業(yè)大學(xué) 動物科技學(xué)院,四川 雅安 625014; 2.西北農(nóng)林科技大學(xué) 動物科技學(xué)院,陜西 楊凌 712100)

全基因組關(guān)聯(lián)分析(genome-wide association study，GWAS)是一種旨在充分利用群體水平的連鎖不平衡，以SNP芯片技術(shù)為基礎(chǔ)，在全基因組范圍內(nèi)定位影響表型性狀的遺傳因素的統(tǒng)計遺傳學(xué)分析方法。該方法以候選基因關(guān)聯(lián)分析為理論原型，由人類遺傳學(xué)家為開展復(fù)雜疾病的遺傳研究而首先提出[1]，其理論依據(jù)是“常見疾病，常見變異”(Common disease-common variants，CDCV)的遺傳學(xué)假說[2]，即從遺傳學(xué)的角度，人類常見疾病主要由少數(shù)等位基因頻率(MAF)高于0.01(未知)常見遺傳變異位點(diǎn)所引起[3]。

2005年，Klein等[4]采用GWAS分析的方法，成功地鑒定到影響年齡相關(guān)性視網(wǎng)膜黃斑變異的重要遺傳因子補(bǔ)體因子H，標(biāo)志著全基因組關(guān)聯(lián)分析方法真正開始應(yīng)用到復(fù)雜疾病和數(shù)量性狀的遺傳分析研究中。受人類復(fù)雜疾病的GWAS研究所取得大量進(jìn)展的影響，家養(yǎng)動物的GWAS研究也陸續(xù)開展起來[5-8]，尤其是對于質(zhì)量性狀遺傳機(jī)制的研究取得了巨大的成功[9-11]。

當(dāng)前GWAS研究中主要采用的分析策略是基于群體水平樣本數(shù)據(jù)的單標(biāo)記統(tǒng)計分析，但此類分析策略卻面臨著各種混雜因素(confounding factors)導(dǎo)致的高強(qiáng)度假陽性率的巨大挑戰(zhàn)。而混雜因素主要分為兩大類：群體分層(population stratification)和家系親緣關(guān)系(familial relatedness)。因此，在過去的幾年中，群體水平GWAS統(tǒng)計方法的發(fā)展重點(diǎn)和熱點(diǎn)就是保持較高的統(tǒng)計分析功效的同時，如何清除或校正群體混雜因素的影響。本文主要綜述了群體分層的主要檢測方法、基于單標(biāo)記線性混合模型策略的GWAS分析方法的研究進(jìn)展及全基因組關(guān)聯(lián)分析樣本量的估計，為定位影響家養(yǎng)動物數(shù)量性狀的主要QTL提供理論依據(jù)。

1 群體混雜程度的度量和檢測

1.1 群體混雜程度的度量

關(guān)聯(lián)研究的基本假設(shè)之一是遺傳標(biāo)記與表型之間相互獨(dú)立[12,13]，即樣本個體來自同一群體。然而，在實(shí)際GWAS研究中，上述假設(shè)經(jīng)常會由于樣本個體來自多個群體遺傳特征差異較大的群體而無法成立；與此同時，群體分層對關(guān)聯(lián)分析的影響與樣本容量成正比[14]。因此，在實(shí)施群體水平的關(guān)聯(lián)分析時，忽略混雜因素的干擾將會引起偽關(guān)聯(lián)。為清除群體分層對統(tǒng)計推斷的影響，Devlin 和Roeder[12]提出基因組控制(genomic control，GC)，同時采用基因組通脹因子(genomic inflation factor)度量混雜因素對關(guān)聯(lián)分析的干擾程度。GC方法的基本思想是：比較關(guān)聯(lián)分析中大量(假定)非關(guān)聯(lián)位點(diǎn)的卡方檢驗統(tǒng)計量與原假設(shè)下自由度為1的卡方分布對應(yīng)的統(tǒng)計量之間的大小。為獲得穩(wěn)健的比值關(guān)系，一般使用兩種分布的中位數(shù)進(jìn)行比較，其公式如下：

(1)

其中，λ定義為通脹因子，反映了混雜因素對關(guān)聯(lián)分析的干擾程度。在GWAS分析中，理論上僅有非常小的一部分位點(diǎn)與目標(biāo)性狀表現(xiàn)出真實(shí)的關(guān)聯(lián)，所以在實(shí)踐中通常所有的待檢測位點(diǎn)都被用于計算λ。而群體分層因素對關(guān)聯(lián)分析的干擾與樣本容量成正比，因而通脹因子λ也是樣本大小的函數(shù)。據(jù)此，不同GWAS研究的通脹程度的比較可以采用λ1000。λ1000 被定義為當(dāng)樣本大小為1000時GWAS分析的通脹大小。

一般認(rèn)為，如果λ1000 < 1.05則群體混雜因素的干擾可以被忽略；如果λ1000 > 1.05則需要對關(guān)聯(lián)分析進(jìn)行校正, 即原始的檢驗統(tǒng)計量除以λ。目前，λ或λ1000是GWAS研究中度量關(guān)聯(lián)結(jié)果通脹大小的標(biāo)準(zhǔn)方法。總體來說，GC在群體分層程度不嚴(yán)重的GWAS研究中是一種較為合適的校正方法；但是，當(dāng)多個遺傳異質(zhì)群體混雜或?qū)?shù)量性狀進(jìn)行GWAS分析時，GC則顯示出過于保守，簡單地使用該方法會損失統(tǒng)計功效[15,16]。另外，Yang等[17]發(fā)現(xiàn)λ與目標(biāo)性狀的遺傳力呈正比，對于遺傳力較高的數(shù)量性狀往往λ的估計值理論上較大。因此，對于大樣本容量的GWAS研究需要新的策略或方法校正群體混雜因素。

1.2 群體分層的主要檢測方法

一般而言，在人類GWAS研究中，群體分層是最主要的混雜干擾因素。針對樣本中可能存在的異質(zhì)性群體混雜，Prichard等[18,19]首先提出基于病例-對照設(shè)計的二值性狀結(jié)構(gòu)群體關(guān)聯(lián)檢驗(structured population association test，STRAT)方法。該方法的基本設(shè)計思路是：建立在亞群水平上的標(biāo)記位點(diǎn)與表型的關(guān)聯(lián)分析檢驗?zāi)芘懦后w分層的干擾。其實(shí)施過程包括兩個階段：首先根據(jù)全基因組的SNPs的基因型信息，確定樣本群體可能的群體結(jié)構(gòu)[18]；然后以亞群內(nèi)標(biāo)記位點(diǎn)與表型相互獨(dú)立為原假設(shè)，構(gòu)造似然比(likelihood ratio)統(tǒng)計量，檢驗標(biāo)記位點(diǎn)與表型是否關(guān)聯(lián)[19]。在STRAT方法基礎(chǔ)上，Thornberry等[20]又發(fā)展出基于Logistic回歸模型的數(shù)量性狀的結(jié)構(gòu)關(guān)聯(lián)分析。理論上講，樣本個體的歸類對未知的亞群數(shù)高度敏感，且亞群數(shù)目的推斷依賴于模型的選擇[21]?？傮w來說，由于首先需要推斷樣本群體的群體結(jié)構(gòu)，STRAT分析通常需要龐大的計算量。因此，該方法在實(shí)踐中僅適合于樣本量較小的GWAS研究，對于樣本量較大(幾千個個體)的GWAS研究，該方法在計算效率方面已無法滿足要求。

當(dāng)前，由Patterson等[21]發(fā)展的主成分分析(principal component analysis，PCA)是另一種被大量使用的異質(zhì)性群體混雜的檢測方法。該方法的基本原理是：個體遺傳相似矩陣，即親緣關(guān)系矩陣，包含了樣本群體的群體遺傳特征。從線性代數(shù)的角度，親緣關(guān)系矩陣是實(shí)對稱矩陣，對該矩陣進(jìn)行特征值和特征向量分析較為容易；而特征值大小反映了與該特征值對應(yīng)的特征向量對原矩陣信息的負(fù)載量。因此，原始矩陣的信息可以采用少數(shù)幾個特征值較大的特征向量近似代替。與STRAT不同，PCA直接目標(biāo)不是將樣本個體歸類到不同的亞群而是個體在變異軸上的坐標(biāo)。而該方法另一個顯著特征是：PCA方法提供了群體是否分層的統(tǒng)計檢驗方法，使得群體混雜的推斷更加可靠；更為重要的是，PCA有著極高的計算效率。Price等[15]又進(jìn)一步將PCA應(yīng)用到病例-對照設(shè)計的全基因組關(guān)聯(lián)分析中，PCA應(yīng)用于GWAS的基本思想是：根據(jù)PCA分析得到特征向量對基因型和表型值分別進(jìn)行校正，然后利用校正后的數(shù)據(jù)進(jìn)行關(guān)聯(lián)分析和統(tǒng)計推斷。總體來說，PCA是一種快速、高效且穩(wěn)定的群體遺傳結(jié)構(gòu)的分析方法。目前，在群體水平的GWAS研究中被廣泛采用，除了EIGENSOFT軟件外，該方法也被整合到其他程序中[22-23]。

在人類復(fù)雜疾病的GWAS研究中，主要的混雜因素是群體分層，因此上述三種校正方法單一或者綜合地使用能有效的降低混雜因素引起的假陽性[15]。而在家養(yǎng)動植物和模式生物中，除了群體分層，復(fù)雜的家系關(guān)系往往是更為重要的混雜因素。在這種情形下，主要針對群體分層的上述方法就不能很好地校正這些混雜因素對關(guān)聯(lián)分析的干擾。

2 數(shù)量性狀GWAS研究的線性混合模型

在家養(yǎng)動物中，諸如產(chǎn)奶量、乳脂量、乳蛋白量等奶牛重要經(jīng)濟(jì)性狀屬于受大量基因和諸多環(huán)境因素共同影響的數(shù)量性狀。對于數(shù)量性狀的遺傳關(guān)聯(lián)分析，基本的統(tǒng)計推斷策略仍然是：采用基于線性回歸模型或廣義線性模型分析的t檢驗或F檢驗等。而數(shù)量性狀的遺傳理論基礎(chǔ)微效多基因假說，但這些效應(yīng)微小的遺傳位點(diǎn)通常很難單獨(dú)檢測，所以只能將這些未知位點(diǎn)的效應(yīng)當(dāng)作一個整體進(jìn)行分析，即微效多基因加性效應(yīng)和。據(jù)此，根植于動物育種的線性混合模型方法被建議用于數(shù)量性狀的GWAS分析，并得到巨大的發(fā)展。

針對近交系雜交的玉米樣本群體，Yu等[16]首先提出全基因組關(guān)聯(lián)分析的線性混合模型方法。在該方法中，群體分層和微效多基因效應(yīng)分別被當(dāng)作固定效應(yīng)和隨機(jī)效應(yīng)包括在模型中；同時，應(yīng)用標(biāo)記信息估計微效多基因效應(yīng)的相關(guān)系數(shù)矩陣，而群體分層效應(yīng)的關(guān)聯(lián)矩陣則需要利用STRUCTURE估計，并采用似然率進(jìn)行統(tǒng)計檢驗。總體來看，Yu等[16]開創(chuàng)了GWAS線性混合模型分析的先河，并通過具體數(shù)據(jù)說明線性混合模型在降低假陽性率的同時還能保持較高的功效；但是，上述模型的方差組分參數(shù)較多。因此，在計算效率上該方法只能適用于樣本量較小的分析，對于樣本量較大的人類或家養(yǎng)動物數(shù)據(jù)，應(yīng)用該方法進(jìn)行分析是不現(xiàn)實(shí)的。故關(guān)于GWAS線性混合模型方法后續(xù)發(fā)展的重點(diǎn)是如何提高計算效率。

為提高計算效率，Aulchenko等[22,24]又提出一種線性混合模型的近似方法(genome-wide rapid association using mixed model and regression，GRAMMAR)。GRAMMAR方法在應(yīng)用上包括兩大步驟：首先預(yù)測每個個體的微效多基因效應(yīng)的實(shí)現(xiàn)值，得到表型值的剩余值；然后將剩余值當(dāng)作標(biāo)準(zhǔn)的表型值，采用簡單線性回歸模型進(jìn)行關(guān)聯(lián)分析。相對于Yu等[16]的方法，GRAMMAR僅需進(jìn)行一次方差組分參數(shù)的估計，極大地減少計算任務(wù)。值得注意的是，對于遺傳力低的性狀該方法有著較好的近似效果；而對于遺傳力較高的數(shù)量性狀，應(yīng)用該方法會降低功效甚至?xí)饦?biāo)記效應(yīng)估計的偏差。

除計算效率低，在Yu等[16]的方法中，采用SPAGeDi軟件估計的親緣關(guān)系矩陣有時候并不是正定矩陣，此時無法進(jìn)一步估計方差組分參數(shù)。圍繞上述問題，Kang等[25]提出了高效混合模型關(guān)聯(lián)分析算法EMMA(efficient mixed-model association)，并開發(fā)出基于R語言[26]的程序包EMMA。該方法的最大優(yōu)點(diǎn)是其關(guān)于標(biāo)記與性狀關(guān)聯(lián)檢驗的統(tǒng)計量是精確統(tǒng)計量而不是近似值；該方法的另一個突出特點(diǎn)是該算法僅包含兩個方差組分，非常適合家養(yǎng)動物數(shù)據(jù)。盡管如此，方差組分的估計問題使得EMMA的計算效率理依舊不是太高，同時R本身的設(shè)計原理也決定了EMMA軟件不擅長處理大樣本的GWAS數(shù)據(jù)。為提高EMMA的計算效率，Kang等[27]進(jìn)一步發(fā)展出更為高效的EMMAX算法和軟件(EMMA eXpedited)。EMMAX算法的基本思想是：絕大多數(shù)作用于數(shù)量性狀的遺傳位點(diǎn)的效應(yīng)往往很小，故對于同一個數(shù)據(jù)集內(nèi)所有SNP的關(guān)聯(lián)分析和檢驗，不需要每次重新估計方差組分，而只需基于零假設(shè)下的模型一次性估計出方差組分(或遺傳力)，并應(yīng)用在后續(xù)的關(guān)聯(lián)分析中。根據(jù)相同的思想，Zhang等[28]提出了壓縮線性混合模型方法——compressed MLM(compressed mixed linear model)及P3D算法(population parameters previously determined)，并將該方法實(shí)現(xiàn)在TASSEL軟件中。類似地，Lippert等[29]又提出FaST-LMM(factored spectrally transformed linear mixed model)算法及軟件。最近，Svishcheva等[30]基于GRAMMAR[24]和FASTA方法[31]又發(fā)展了更為計算效率更高但效應(yīng)估計值為無偏估計的GRMMAR-Gamma方法。

如前所述，EMMAX、TASSEL和FaST-LMM三種算法及軟件擁有非常高的計算效率，應(yīng)用相關(guān)軟件可以在普通臺式計算機(jī)上進(jìn)行上千個體規(guī)模的GWAS分析。但需要強(qiáng)調(diào)的是，從統(tǒng)計推斷的角度，上述三種方法都是近似算法。毫無疑問，精確檢驗統(tǒng)計量的實(shí)現(xiàn)仍然是線性混合模型方法在理論上追求的目標(biāo)，而計算效率卻是最大的障礙。最近，Zhou和Stephens[30]開發(fā)出實(shí)現(xiàn)精確檢驗統(tǒng)計量的高效的GWAS分析程序GEMMA(genome-wide efficient mixed-model association)，該程序的運(yùn)行速度大約是EMMA的N倍。

總體而論，基于單標(biāo)記策略的GWAS分析方法發(fā)展的已經(jīng)相當(dāng)成熟，尤其是線型混合模型方法的提出和應(yīng)用。對于不同的樣本群體，研究人員可根據(jù)數(shù)據(jù)的規(guī)模和特征選擇合適的方法和軟件。

3 GWAS研究的樣本量估計

從統(tǒng)計學(xué)的角度，一項研究采用的樣本量從根本上決定統(tǒng)計分析的功效；從數(shù)量遺傳學(xué)的角度，由于大多數(shù)遺傳位點(diǎn)對數(shù)量性狀效應(yīng)太小。因此，必須增大樣本，減弱隨機(jī)誤差的干擾，才有可能定位到QTL。但在當(dāng)前的情形下，商業(yè)SNP芯片的價格依舊非常昂貴，樣本量的增大所帶來的試驗費(fèi)用的增加往往是非常巨大的。所以，在正式開展一項全基因關(guān)聯(lián)分析研究之前，很有必要對數(shù)量性狀的全基因關(guān)聯(lián)分析所需要的樣本量進(jìn)行簡單地估計。而針對家養(yǎng)動物全基因組關(guān)聯(lián)分析，Goddard和Hayes[32]提出了一個計算樣本量的簡單近似公式：

rm*t=rm*q×rq*g×rg*t

(2)

(3)

4 結(jié) 語

以QTL檢測和定位為前提的數(shù)量性狀遺傳機(jī)制的探索是一項在幾乎沒有任何先驗信息的情況下，從零開始的研究。由于數(shù)量性狀遺傳機(jī)制的復(fù)雜性，對于數(shù)量性狀遺傳機(jī)制的研究不管是從試驗設(shè)計、統(tǒng)計分析方法還是后續(xù)的分子生物學(xué)工作都要比質(zhì)量性狀的遺傳定位研究更加復(fù)雜和困難。而充分利用群體水平所積累的歷史重組事件的全基因組關(guān)聯(lián)分析已經(jīng)被大量研究證明是研究家養(yǎng)動物重要經(jīng)濟(jì)性狀遺傳因素的一種高效方法。

參考文獻(xiàn)：

[1] Risch N, Merikangas K. The future of genetic studies of complex human diseases[J]. Science, 1996, 273(5281): 1 516-1 517.

[2] Reich D E, Lander E S. On the allelic spectrum of human diseases[J]. Trends Genet, 2001, 17(9): 502-510.

[3] Wang W Y S, Barratt B J, Clayton D G, et al. Genome-wide association studies: theoretical and practical concerns[J]. Nat Rev Genet, 2005, 6(2): 109-118.

[4] Klein R J, Zeiss C, Chew E Y, et al. 2005. Complement factor H polymorphism in age-related Macular degeneration[J]. Science, 308(5720):385-389.

[5] Abasht B, Lamont S J. Genome-wide association analysis reveals cryptic alleles as an important factor in heterosis for fatness in chicken F2 populations[J]. Anim Genet, 2007, 38:491-498.

[6] Daetwyler H D, Schenkel F S, Sargolzaei M, et al. A genome scan to detect quantitative trait loci for economically important traits in Holstein cattle using two methods and a dense single nucleotide polymorphism map[J]. J Dairy Sci, 2008, 91(8):3 225-3 236.

[7] Cole J B, Wiggans G R, Ma L, et al. Genome-wide association analysis of thirty one production, health, reproduction and body conformation traits in contemporary U.S. Holstein Cow[J]. BMC Genomics, 2011, 12:408.

[8] Gregersen V R, Conley L N, Sorensen K K, et al. Genome-wide association scan and phased haplotype construction for quantitative traits loci affecting boar taint in three pig breeds [J]. BMC Genomics 2012, 13:22.

[9] Karlsson E K, Baranowska I, Wade C M, et al. Efficient mapping of mendelian traits in dogs through genome-wide association[J]. Nat Genet, 2007, 39: 1 321-1 328.

[10] Dorshorst B, Molin A M, Rubin C J, et al. A complex genomic rearrangement involving the Endothelin 3 locus cuases dermal hyperpigmentation in the chicken[J]. PLoS Genet, 2011, 7(12):e1002412.

[11] Andersson L S, Larhammar M, Memic F, et al. Mutations in DMRT3 affect locomotion in horses and spinal circuit function in mice[J]. Nature. 2012. 488: 642-646.

[12] Devlin B, Roeder K. Genomic control for association studies[J]. Biometrics, 1999, 55(4):997-1 004.

[13] Reich D, Goldstein D. Detecting association in a case-control study while allowing for population stratification[J]. Genet Epidemiol, 2001, 20(1): 4-16.

[14] Marchini J, Cardon L R, Phillips M S, et al. The effects of population structure on large genetic association studies[J]. Nat Genet, 2004, 36(5): 512-517.

[15] Price A L, Patterson N J, Plenge R M, et al. Principal component analysis corrects for stratification in genome-wide association studies[J]. Nat Genet, 2006, 38(8):904-909.

[16] Yu J, Pressoir G, Briggs W, Vroh B I, et al. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness[J]. Nat Genet, 2006, 38(2):203-208.

[17] Yang J, Weedon M N, Purcell S, et al. Genomic inflation factors under polygenic inheritance[J]. Europ J Hum Genet, 2011, 19(7):1-6.

[18] Prichard J K, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data[J]. Genetics, 2000, 155(2): 945-959.

[19] Prichard J K, Stephens M, Rosenberg N A, et al. Association mapping in structured populations[J]. Am J Hum Genet, 2000. 67(1): 170-181.

[20] Thornsberry J M, Goodman M M, Doebley J, et al. Dwarf8 polymorphisms associate with variation in flowering time[J]. Nat Genet, 2001, 28(3):286-289.

[21] Patterson N, Price A L, Reich D. Population structure and eigen analysis[J]. PLoS Genet, 2006, 2: e190.

[22] Aulchenko Y S, Ripke S, Isaacs A, et al. GenABEL: an R library for genome-wide association analysis[J]. Bioinformatics, 2007, 23(10):1 294-1 296.

[23] Yang J, Lee S H, Goddard M E, et al. GCTA: a tool for Genome-wide complex trait analysis[J]. Am J Hum Genet, 2011, 88(1):76-78.

[24] Aulchenko Y S, de Konning D J, Haley C. Genome-wide rapid association using mixed model and regression: a fast and simple method for genome-wide pedigree-based quantitative trait loci association analysis[J]. Genetics, 2007, 177(1): 577-585.

[25] Kang H M, Zaitlen N A, Wade C M, et al. Efficient control for population structure in model organism association mapping[J]. Genetics, 2008, 178(3):1 709-1 723.

[26] R development Core Team. R: A language and environment for statistical computing[M]// R Foundation for Statistical Computing. Vienna, Austria: 2005.

[27] Kang H M, Sul J H, Service S K, et al. Variance component model to account for sample structure in genome-wide association studies[J]. Nat Genet, 2010, 42(4): 348-353.

[28] Zhang Z W, Ersoz E, Lai C Q, et al. 2010. Mixed linear model approach adapted for genome-wide association studies[J]. Nat Genet, 42(4): 355-360.

[29] Lippert C, Listgarten J, Liu Y, et al. FaST linear mixed models for genome-wide association studies[J]. Nat Methods, 2011,8: 833-835.

[30] Svishcheva G R, Axenovich T I, Belonogova N M, et al. Rapid variance component-based method for whole-genome association analysis [J]. Nat Genet, 2012, 44(10): 1166-1170.

[31] Chen W M, Abecasis G R. Family-based association tests for genome wide association scans[J]. Am J Hum Genet, 2007, 81(5): 913-926.

[32] Goddard M E, Hayes B J. Mapping genes for complex traits in domestic animals and their use in breeding programmes[J]. Nat Rev Genet, 2009,10(6): 381-391.

[33] Pausch H, Flisikowski K, Jung S, et al. Genome-wide association study identifies two major loci affecting calving ease and growth-related traits in cattle[J]. Genetics, 2011,187(1):289-297.