全基因組關聯分析中混合模型的原理、優化與應用

譚力治，趙毅強

全基因組關聯分析中混合模型的原理、優化與應用

中國農業大學生物學院，北京 100193

全基因組關聯分析（genome-wide association study，GWAS）是定位基因組中與性狀顯著關聯的變異位點的有效方法。隨著表型記錄的完善、高通量基因型分型技術的發展，以及統計方法的改進，全基因組關聯分析在人類疾病、動物植物遺傳等領域得到了廣泛的應用。假陽性是影響全基因組關聯分析結果可靠性的重要因素之一。為了控制假陽性，除了校正值，GWAS模型從最簡單的方差分析（或用于質量性狀的卡方檢驗）到加入固定效應協變量的普通線性模型（general linear model，GLM），再到加入隨機效應的混合線性模型（mixed linear model，MLM）持續改進，控制了多種混雜因素導致的假陽性。將個體的遺傳效應擬合為由基因組親緣關系矩陣（genomic relationships matrix，GRM）定義的隨機效應是目前常用的方法。由于MLM的參數估計大量消耗計算資源，研究人員不斷嘗試模型求解優化和GRM的構建優化（GRM的構建優化同時也提高了計算效率），最終將基于MLM計算的時間復雜度由O（MN3）逐步改進到O（MN），實現了計算速度與統計功效的飛躍。針對質量性狀病例對照比失衡帶來的假陽性問題，研究人員進一步對廣義混合線性模型（generalized linear mixed model，GLMM）進行了校正。本文較全面地介紹了GWAS的基本原理和發展，著重闡述了GWAS中MLM模型的改進和優化細節，同時，列舉了GWAS在農業中的應用，包括在植物、動物和微生物方面的研究成果，以及基于單倍型的GWAS應用。最后，從進一步提高GWAS統計功效和GWAS試驗設計2個角度對GWAS未來的發展進行了展望。

全基因組關聯分析；復雜性狀；隨機效應；基因組親緣關系矩陣；混合線性模型

1 概述

1.1 GWAS的基礎

20世紀以來，作為一種有效的候選基因定位方法，連鎖分析廣泛用于定位孟德爾性狀和常見疾病的基因和變異[1]。連鎖分析在家系內尋找性狀與標記等位基因的共分離，對于單基因性狀的定位具有較高的精度。對于人類遺傳病等復雜性狀的遺傳位點，由于單個變異的邊際效應過小，其方法具有較大的局限[2]。人類基因組計劃（human genome project，HGP）在2001年發表了人類基因組草圖，成為基因組研究的一個重大進步。人類基因組單體型圖計劃（the international hapmap project，HapMap）于2002年啟動，旨在建立人類全基因組遺傳變異圖譜。HapMap計劃基于2個重要的遺傳概念：其一是遺傳變異和表型變異存在關聯；其二是標記之間的連鎖不平衡（linkage disequilibrium，LD）。在LD區間內，變異之間的信息冗余，LD的程度決定需要多少個遺傳變異來對全基因組進行“標記”。這兩點也成為全基因組關聯分析（genome-wide association study，GWAS）重要的理論基礎。

在人類基因組計劃完成后的十年，以Illumina Infinium和Affymetrix原位光刻為代表的高密度芯片分型技術，以及以Solexa為代表的下一代測序技術（next generation sequencing，NGS）飛速發展，大大提高了基因分型的通量并降低了分型成本。眾多物種基因組圖譜的繪制和高通量基因型分型技術的進步，為GWAS研究提供了豐富的標記信息，促進了人類和動植物疾病以及復雜性狀的遺傳定位。另一方面，GWAS的統計模型從非參數卡方檢驗到普通線性模型（general linear model，GLM）再到混合線性模型（mixed linear model，MLM）持續改進，MLM成為當前GWAS的首選方法。與其他方法相比，MLM可以同時捕獲由于群體分層、家系結構和潛在關聯而產生的混雜效應，實現了更高的統計功效[2-4]。盡管與簡單的模型相比，基于MLM的GWAS分析計算量非常大，但隨著研究者的不斷努力，目前大規模GWAS的計算已經被優化到人們可以接受的程度，其算法的時間復雜度實現了O（MN3）—O（MN2）—O（MN1.5）—O（MN）的巨大進步[5-7]。

1.2 GWAS的基本原理

盡管連鎖不平衡也是GWAS的理論基礎之一，GWAS并不依賴于家系而是利用群體的歷史重組信息。常見疾病-常見變異假設（common disease common variant，CDCV）認為常見疾病由多個在群體中普遍存在的變異共同導致，每個變異都對疾病發生產生貢獻。基于CDCV假設，GWAS使用豐富的單核苷酸多態（single nucleotide polymorphism，SNP）為標記，掃描表型在整個基因組中的關聯信號。

GWAS中各位點通常獨立進行檢驗。對于二元質量性狀，數據可以表示為一個列聯表，表中每個元素為具有特定基因型-表型組合類別的個體數量，使用卡方檢驗判斷類別之間的獨立性。遺傳模型的類型決定了列表的形式，例如，在顯性或隱性模型的情況下是一個2×2的列聯表，而加性模型則用一個2×3的列聯表表示。由于加性模型認為基因型與表型具有有序關系，在趨勢或有序的假設下，也可以使用Cochran- Armitage趨勢檢驗來捕捉這種關系。

在實際應用中，研究人員可能希望加入協變量來控制混雜因素，例如，某種疾病的患病概率往往隨著年齡的增長而增加，或者與性別相關。使用線性模型可以方便地將年齡和性別作為協變量加入，對模型進行調整。對于質量性狀，最常用的線性模型是logistic回歸模型[8]。而對于數量性狀，則使用GLM。

群體分層是另一類主要的混雜因素。群體分層指亞群體間等位基因頻率的系統性差異。從群體遺傳學的角度看，群體分層可能由于選擇壓力或者遺傳漂變所致。群體分層可能造成基因型和表型的虛假關聯。如圖1所示，子群體1和子群體2之間存在明顯的等位基因頻率差異，單獨來看，2個子群體的OR值（odds ratio）均為1。將2個子群體結合后，合并后群體的OR值達到2.87，提示群體分層造成了假陽性結果。為了識別真正的關聯信號，控制群體分層很有必要[9]。

圖1 群體分層示例

控制群體分層的第一種方法是使用單一群體確保群體同質性。第二種方法是基于家系的設計，從直系子女和全同胞個體中挑選不同表型的個體進行分析。第三種方法是使用Devlin等[10]提出的基因組控制方法，使用Cochran-Armitage趨勢檢驗來計算膨脹系數，用來校正卡方統計量的膨脹。然而這些方法存在一些局限性，比如難以獲得足夠符合要求的樣本，方法的適用性有限，以及使用統一的調整忽略了個體間基因型的差異。為了解決這些局限，人們開發了另一種替代方法，即用基因型主成分來捕捉群體結構[11]。主成分分析抽取數據中的主要變異，被證明能夠準確反映種群之間整體的遺傳差異，人們將基因型前幾個主成分用作協變量來校正群體分層產生的影響。由于其計算簡單，效果較好，基于主成分的校正是目前最常用的控制群體分層的方法。

由于目前的全基因組分型芯片或全基因組測序提供的基因座數量可達上百萬，因此，GWAS中對每個基因座單獨檢驗后的顯著性校正尤為重要。Bonferroni校正是多重檢驗校正的經典方法，然而，使用Bonferroni校正GWAS閾值可能由于LD的存在而變得過于保守[12]。人中GWAS普遍使用的閾值為5×10-8，該閾值基于假設獨立的SNP數目大約為100萬個，其閾值隨SNP數目的增加而逐漸嚴格。相比于Bonferroni校正，FDR（false discovery rate）校正是一種“相對溫和”的校正方法。該方法不追求完全避免假陽性的結果，而是將假陽性結果和真陽性的比例控制在一定范圍內，一般選擇信號的前5%作為FDR校正的標準。此外，使用置換檢驗（permutation test）來獲得調整后的顯著性值被認為是最好的校正方法，但這種方法的計算量大，盡管使用其近似方法[13]可以提高計算效率，但在實際應用中幾乎仍是不可行的。

2 混合線性模型在GWAS中的發展

2.1 從普通線性模型到混合線性模型

在線性模型中，通常使用樣本來自亞群的比例或全部基因型計算所得的主成分（principal component，PC）表示群體結構，稱其為Q矩陣[10, 14]。Q矩陣中的協變量被用作固定效應進行擬合。該模型表示為=++，其中，和分別是表型和單個遺傳標記（SNP），是殘差。這個GLM也被稱為Q模型。

除了群體分層，人們意識到來自個體的效應（或等價的多基因效應）也是產生結果偏倚的因素。個體效應的自相關結構可以通過個體之間的親緣關系矩陣來指示[15]，早期在動物中多使用基于系譜的親緣關系矩陣，由于植物系譜往往未知，因此，無法將在動物中使用的策略直接用于植物[16]。而使用全基因組遺傳標記可方便地計算基因組親緣關系矩陣（genomic relationships matrix，GRM）（此處也稱為K），此時個體的遺傳效應被擬合為由K定義的隨機效應。同時具有Q和K的模型是MLM，表示為=+++，也稱為Q+K模型[4]。研究表明，Q和Q+K模型較好地控制了假陽性，而且Q+K模型比Q模型或單獨的K模型表現更好[4, 15]。

在最初的MLM模型中[4]，Q+K方法的混合模型方程表示為：

=++++

式中，為表型向量；是除SNP或種群結構以外的固定效應；是SNP固定效應；是群體結構固定效應；是多基因隨機效應；是殘差；、、、是分別與、、、對應的設計矩陣。隨機效應的方差()=KV，是個體親緣關系n階方陣，V為多基因遺傳方差；表型的方差=KV+RV，是n階單位矩陣，V為殘差方陣。通過求解混合模型方程，可獲得、、（固定效應）的最佳線性無偏估計值（best linear unbiased estimate，BLUE）和（隨機效應）的最佳線性無偏預測值（best linear unbiased prediction，BLUP）。MLM進一步避免了個體相關導致的假陽性結果[17]。

2.2 GWAS中混合線性模型算法的優化

最大似然法（maximum likelihood，ML）或約束最大似然法（restricted maximum likelihood，REML）常用于MLM的方差組分估計。與GLM相比，MLM的計算量非常龐大，研究人員從不同角度對MLM用于GWAS進行計算效率和統計功效上的優化，表1總結了GWAS中MLM的優化模型。Kang等[18]提出了高效混合模型關聯（efficient mixed model association，EMMA）法。在似然估計中，EMMA通過把遺傳方差和殘差方差2個組分的優化簡化為對兩者比值的優化，并通過特征分解簡化參數估計中的迭代運算，顯著提高了MLM的求解速度。

GWAS中使用成百上千的個體和成千上萬的標記，對每個標記進行檢測時都估計一次隨機效應的方差使得全基因組分析效率低下。Kang等[19]又提出了改進的高效混合模型關聯（efficient mixed-model association expedited，EMMAX）法，該方法認為由位點多基因效應捕獲的個體隨機效應在模型中的貢獻很小。基于這樣的假設，EMMAX法把對隨機效應方差的多次估計改為單次估計并將其在模型中固定。EMMAX法可以視作EMMA法的簡化近似方法，在隨機效應方差組分估計完成后，采用GLM單獨估計每個標記的效應，計算速度相對于EMMA獲得了大幅提升。

壓縮MLM法（compressed mixed linear model，CMLM）[20]通過聚類把相近個體劃分為不同組別，通過組間GRM來替代個體間GRM實現對隨機效應的壓縮。由于隨機效應的計算和個體數的三次方成正比，此方法大大節省了計算時間。Li等[21]對CMLM算法進行改進，從8種聚類算法與3種組間親緣關系算法的24種組合中計算最優組合，稱其為增強CMLM（enriched CMLM，ECMLM），把CMLM的檢測功效提高了10%左右。在提出CMLM的同時，Zhang等[20]也提出了提前確定模型的群體參數（population parameters previously determined，P3D）的兩步法優化策略，第一步通過沒有標記效應的簡化模型估計總的遺傳方差、殘差以及聚類數等群體參數，這些參數作為先驗信息在第二步的模型中固定使用。第二步依然使用MLM，把原始表型和已估計的參數用于模型中估計標記效應。CMLM和P3D可單獨或聯合使用來優化計算和提供統計功效，相對于常規MLM獲得數千倍的效率提升。

Lippert等[5]提出因式譜變換線性混合模型（factored spectrally transformed linear mixed models，Fast-LMM）。該方法的核心是對GRM進行一次特征分解，轉化為多個不相關的矩陣后使其能使用GLM進行高效求解。由于方法并不要求每個標記具有相同的效應，Fast-LMM方法得到的是各標記效應的精確估計。作者將Fast-LMM與EMMAX進行比較，在（wellcome trust case control consortium，WTCCC）克羅恩病（Crohn’s disease）數據集的4 000與8 000個標記中，Fast-LMM的計算速度分別為EMMAX的11和5倍。

Aulchenko等[22]首次提出了簡單快速的GRAMMAR法（genome wide rapid association using mixed model and regression）進行GWAS。GRAMMAR分為兩步：第一步使用GRM作為隨機效應和除標記效應外的協變量做固定效應建模。第一步模型的殘差作為第二步模型的表型，第二步使用GLM僅對每個標記單獨建模。GRAMMAR法第一步中MLM相對耗時，而剝離了隨機效應的第二步GLM非常高效。為了改進GRAMMAR在標記效應估計上的偏差，作者團隊提出了改進的二步法GRAMMAR-Gamma[23]。該方法在第一步中通過特征分解加速矩陣運算，并構造了不考慮標記相關性的簡化得分檢驗（score test）統計量和GRAMMAR-Gamma校正因子。在第二步中仍然使用GLM對每個標記單獨建模，但對檢驗統計量除以GRAMMAR-Gamma因子進行校正以獲得更精確的標記效應的估計。該方法取得了與標準似然檢驗幾乎相同的功效，但極大降低了運算時間。作者使用人類和擬南芥的數據發現GRAMMAR-Gamma的運行時間遠小于EMMAX和Fast-LMM，在運行速度上分別為二者的38和10倍，再次提高了計算效率。

全基因組高效混合模型關聯（genome-wide efficient mixed-model association，GEMMA）法[24]在對GRM進行特征分解之前，對矩陣求一階和二階導數。為了避免EMMA中對每一個標記復雜的特征分解步驟，使用矩陣向量乘法替代特征分解，轉化為只涉及標量的遞歸乘法優化求解。GEMMA法通過優化大型矩陣運算，顯著提高了GWAS的運行速度，并獲得和EMMA法一致的精確解。使用雜交小鼠的高密度脂蛋白膽固醇（high-density lipoprotein cholesterol，HDL-C）數據，作者報道GEMMA的運行速度為Fast-LMM的13.6倍，而使用WTCCC的克羅恩病數據，GEMMA的運行速度為Fast-LMM的1.87倍。作者同時指出，EMMAX與GRAMMAR等近似方法可能導致假陰性，造成GWAS檢測的功效降低。

Fast-LMM的作者進一步提出使用全部標記的子集構建GRM來降低計算開銷，稱為FaST-LMM- Select[25]。在操作中，作者首先使用GLM對標記進行檢驗并對值進行升序排序。選取達到基因組控制因子λ時對應的標記集合，剔除其中強連鎖的標記，使用剩余標記構建性狀特異的GRM。該方法進一步降低了計算成本，且作者發現相較于使用全部標記，該方法能顯著降低假陽性和假陰性率。Wang等[26]提出的SUPER（settlement of mlm under progressively exclusive relationship）方法運用了類似的概念來優化GRM的構建。該方法對標記預先獲得的值或效應值進行排序，在劃分好的若干染色體片段中，選取每個片段中值最低的標記并剔除與待測標記連鎖的標記后，同樣使用剩余標記構建性狀特異的GRM。經作者比較，該方法比FaST-LMM-Select的假陽性更低，并提高了對遺傳力的估計。之后發表的BOLT-LMM[6]法包括兩部分，基礎部分和主流方法一樣采用微效多基因假設下的混合線性模型。但是在估計方差組分時使用共軛梯度法實現近似計算，避免了特征分解所需的大量計算時間和內存。作者采用一種新的回顧性得分檢驗并使用類似GRAMMAR-Gamma的方法對統計量進行校正。BOLT-LMM的改進部分借鑒了貝葉斯方法在動植育種基因組選擇中的應用，認為大部分標記效應較小但存在少部分大效應標記，其假設更貼合實際。作者使用高斯混合分布作為標記效應的先驗分布來擬合貝葉斯線性回歸，通過快速變分法得到近似的表型殘差。最后，基于表型殘差使用同樣的回顧性得分對每個標記進行檢驗并使用LD分數回歸（LD score regression，LDSC）[27]對統計量進行調整。使用23 294個人樣本的脂質、身高、體重指數和血壓等定量性狀，作者證明BOLT-LMM具有更高的統計功效，且運行效率相較FaST-LMM-Select、GEMMA和EMMAX等更高。

Jiang等[7]開發了基于MLM的GWAS工具fastGWA，采用高效的基于網格搜索的約束性最大似然（restricted maximum likelihood，REML）算法fastGWA-REML。在對GRM稀疏化的基礎上，方差組分估計時對矩陣使用Cholesky分解避免對其求逆。fastGWA使用了與GRAMMAR-Gamma類似的統計量檢驗關聯性。fastGWA方法比其他基于MLM的工具快幾個數量級，其內存使用量也極大降低。作者抽取了英國生物樣本庫（UK Biobank，UKB）中400 000個樣本，分別使用fastGWA與BOLT-LMM對體重指數進行GWAS分析。fastGWA的運行時間為BOLT-LMM的1.22%，內存使用量僅為BOLT-LMM的5%，使得Biobank級的GWAS運算成為可能。

上述方法更多聚焦于MLM，尤其是針對模型中作為隨機效應的GRM的計算優化，個別方法在構建GRM時做了一些簡化，但是總體思路接近，方法的統計功效依然與傳統的MLM類似。MLM較好地控制了假陽性，但是人們意識到，在一定程度上其存在標記效應和控制效應的混雜，造成了一定程度的假陰性結果。于是一些新的方法嘗試在模型中剝離隨機效應，從而解決這類混雜問題。

多位點混合模型（multi-locus mixed-model，MLMM）法[28]采用前向后向選擇組合的逐步選擇策略對常規的MLM進行改進。在前向選擇中，每一步都首先估計方差組分，將最顯著的標記作為固定效應的協變量逐步加入模型，用來壓縮模型中的隨機效應。過程中持續更新模型的方差組分，直到隨機效應解釋的變異接近于零或達到指定的循環數后結束前向選擇。類似地，在隨后的后向選擇中逐步把最不顯著的標記協變量從模型中剔除。為了避免近端污染（proximal contamination），作者建議模型中作為協變量的標記不參與GRM的構建。

MLMM法使用向前向后逐步回歸消除了一部分檢測標記效應的混雜問題，增強了GWAS的統計功效。在此基礎上，Liu等[14]提出交替使用固定效應和隨機效應（fixed and random model circulating probability unification，FarmCPU）來進一步解決模型中混雜問題的方法。對每一個標記，該方法把篩選后可能的關聯位點作為協變量加入固定效應模型來進行檢測。而候選關聯位點的篩選是在獨立的隨機效應模型中進行，更新候選關聯位點后重新進入固定效應模型對標記進行檢測，循環往復直到沒有新的候選關聯位點加入到固定效應模型中。由于模型中沒有同時出現固定和隨機效應，FarmCPU方法避免了不同效應的混雜，并同時控制了假陽性和假陰性。作者團隊針對FarmCPU的改進方法BLINK（Bayesian-information and linkage-disequilibrium iteratively nested keyway）法[29]更是不再使用隨機效應模型，而是使用連鎖不平衡和貝葉斯信息準則（Bayesian information criterion，BIC）來入選和篩選可能的關聯位點。作者使用模擬數據對BLINK和FarmCPU進行評估，發現BLINK的運行速度是FarmCPU的3—4倍，且BLINK更具有發掘額外遺傳位點的潛力。本文在單個處理器上（Intel(R) Xeon(R) Gold 6230 CPU @ 2.10GHz）比較了上述模型的單核運算時間。從Galbase數據庫[30]下載928個樣本的基因型數據，隨機抽取GGA 11上443 218個位點用于GWAS測試。表型數據使用GCTA[31]進行模擬，選擇100、200、400、600、800和928個樣本進行不同樣本量梯度的測試。經測試，EMMA、CMLM、SUPER運算時間依次減少，但遠高于其他模型（圖2）。FastGWA運算速度最快，其次為GEMMA、BLINK、FarmCPU，其運算時間在樣本量較大時明顯低于其他模型。BOLT-LMM、MLMM、FaST-LMM、EMMAX運算時間高于上述4個模型，可能由于樣本量和標記數目的限制，未能觀察到明顯差距（圖2）。

表1 GWAS中MLM的優化模型

圖2 混合模型算法概述與運算速度比較

2.3 GWAS中的廣義混合線性模型算法

隨著許多大型生物樣本庫和隊列的建立，基于MLM的GWAS在數量性狀中取得了巨大的成功。然而大部分人類遺傳病是病例對照（case-control）研究，屬于二分類的質量性狀，患病率較低的疾病還可能出現病例-對照比極不平衡的現象，這樣的數據不滿足MLM中殘差具有方差均等的假設。當存在群體分層，尤其病例-對照比不平衡時，基于MLM的GWAS可能無法有效地控制Ⅰ型錯誤率（假陽性），因此，研究者們對MLM做出擴展，將更適合此類情況的廣義混合線性模型（generalized linear mixed model，GLMM）應用于GWAS中，表2匯總了GWAS中GLMM的優化模型。

Chen等[32]提出廣義線性混合模型關聯檢驗（generalized linear mixed model association test，GMMAT），該方法基于殘差異方差的logistic混合模型。方法首先構建不含標記效應的零模型，模型擬合中作者使用懲罰擬似然法（penalized quasi-likelihood，PQL）和平均信息約束最大似然法（average information restricted maximum likelihood，AI-REML）進行參數估計。估計得到的模型參數被固定下來用于所有標記，并在此基礎上使用得分檢驗對每個標記效應進行估計。作者比較了GMMAT和SAS的PROC GLIMMIX過程，發現GMMAT在擬合具有一個方差分量的logistic混合模型時，其運算時間為SAS PROC GLIMMIX的1.5%。此外，在存在不平衡的病例-對照比的情況下，GMMAT更好地控制了假陽性。

Zhou等[33]提出了可擴展的精確廣義混合模型（scalable and accurate implementation of generalized mixed model，SAIGE）法，用于處理病例對照比失衡的大規模質量性狀數據。SAIGE法的步驟與GMMAT類似，但是過程中大量使用了計算優化。SAIGE法第一步使用AI-REML估計方差組分等參數，使用預處理共軛梯度法（preconditioned conjugate gradient，PCG）替代矩陣特征分解，節省了GRM相關的計算成本。第二步使用方差比來校準得分統計量的方差，并借助鞍點近似法（saddlepoint approximation，SPA）克服二分類性狀中病例對照比失衡的問題，獲得準確度高的值。作者從UKB中隨機抽取了冠狀動脈疾病樣本，分別使用BOLT-LMM、GMMAT和SAIGE進行GWAS分析。結果表明，SAIGE的假陽性率顯著優于前兩者。SAIGE的時間復雜度與BOLT-LMM相同，但由于logistic混合模型的迭代步驟多于MLM，導致其運算速度低于BOLT-LMM，但依然快于GMMAT。

同樣使用SPA來校正病例對照比失衡帶來的假陽性問題，Jiang等[34]開發的fastGWA-GLMM工具沿用了前序fastGWA法中基于網格搜索的算法估計方差組分，并使用GRM稀疏化提高計算效率。作者通過從UKB中隨機抽取樣本，將fastGWA-GLMM與SAIGE進行比較，在樣本量為400 000時，fastGWA-GLMM的運行速度為SAIGE的36.8倍，極大地節省了運算成本。

除病例-對照研究產生的二分類質量性狀外，臨床上常使用有序分類測量來衡量疾病的嚴重程度，如從1—9對疾病的感染程度打分，1為基本不感染，9為嚴重感染。把有序分類變量當作連續變量或者降級為二分類變量使用都不合適，基于此，Bi等[35]提出了比例優勢logistic混合模型（proportional odds logistic mixed model，POLMM）法。該方法使用和上述方法相似的統計架構，將logistic模型應用于有序分類表型，使用PQL與AI-REML對零模型進行擬合并進行參數估計，也通過SPA校準第二步得分檢驗的值。根據使用的GRM的不同，POLMM法提供DensePOLMM與FastPOLMM 2種實現，DensePOLMM使用稠密矩陣并通過PCG進行矩陣加速運算，FastPOLMM使用稀疏矩陣，在計算速度上有優勢但統計功效略低于DensePOLMM。作者使用BOLT- LMM、FastPOLMM-NoSPA和FastPOLMM對UKB中4種食品偏好的有序分類變量進行GWAS。當表型分布平衡時，BOLT-LMM獲得與FastPOLMM一致的結果；當表型分布不平衡時，FastPOLMM-NoSPA優于BOLT-LMM，而FastPOLMM較前兩者更好地控制了假陽性。

表2 GWAS中GLMM的優化模型

使用上述數據在單個處理器上（Intel(R) Xeon(R) CPU E3-1245 v6 @ 3.70GHz）比較不同樣本量下GLMM模型的單核運算時間。與MLM測試相似，fastGWA-GLMM具有最快的運算速度，FastPOLMM幾乎與其相同，SAIGE在樣本量較小時速度較快，在樣本量大于800時速度顯著降低，DensePOLMM在4種模型中使用了最多的計算時間（圖2）。

3 GWAS的應用

隨著MLM在GWAS中的引入，越來越多和動植物重要性狀關聯的基因和位點被發現，推動了遺傳學的發展。

3.1 GWAS在植物中的應用

自從Hansen等[36]使用BSA分析（bulked segregant analysis）首次將GWAS用于海甜菜后，GWAS方法已成功用于植物育種、逆境生長、抵御病蟲害等方面的候選基因鑒定。Atwell等[37]使用EMMA對擬南芥自交系107種表型進行了GWAS分析，鑒定到6種受單基因控制的表型，包括抗病響應、開花基因表達等。

水稻與人類生活息息相關，Zhao等[38]對不同國家的413種水稻品種開展GWAS分析，利用EMMA鑒定到34個與植物生理、生長發育以及植物形態發生相關的候選基因。Huang等[39]采用低深度測序數據進行基因分型，使用CMLM法將GWAS應用于秈稻的14種育種相關性狀，共篩選出37個與開花日期和增產性狀顯著相關的QTL，該分析促進了水稻優良農藝性狀的遺傳解析以及品種的選育。在全球950個水稻品種中，Huang等[40]使用CMLM分別對其中的地方粳稻群體、地方秈稻群體與全部品種進行GWAS分析，定位到32個與開花日期和農藝性狀顯著相關的QTL。產量是備受關注的農業性狀之一，研究人員在利用TASSEL開展的基于混合線性模型的14種高產水稻品種的GWAS分析中，找到1 152個和高產相關的顯著位點，同時檢測到8個高產相關滲入區段，揭示了高產相關基因由于受到人工馴化而發生了遺傳結構的改變[41]。

作為另一種重要的農作物，玉米中的GWAS研究也已廣泛開展。Wang等[42]通過TASSEL軟件在114種玉米自交系品種中定位到18個與抗黑穗病顯著相關的變異位點，同時揭示了玉米抗黑穗病的遺傳機理。Li等[43]利用CMLM在368種玉米自交系的103萬個變異位點中鑒定出74個與籽粒油分和脂肪酸形成相關的候選基因，結合eQTL（expression QTL）與共表達分析發現有三分之一的候選基因與油分的代謝途徑相關。同樣，研究人員使用CMLM進行GWAS分析，發現玉米的早花性狀與220個遺傳標記顯著關聯，其中大部分關聯位點與擬南芥同源[44]。綠色保持是植物延緩衰老的一種表現形式，Sekhon等[45]通過測定葉綠素熒光參數v/m指示植物延緩衰老的程度，使用GAPIT軟件進行關聯分析，鑒定到64個與其顯著相關的候選基因，其中14個基因已被證明與衰老過程相關。作者同時構建了延緩衰老與正常衰老基因的共表達網絡，對候選基因的生物學功能作出系統解釋。Chao等[46]通過TASSEL軟件對玉米內核鋅濃度進行GWAS分析，鑒定到鋅-煙酰胺轉運蛋白基因，并驗證了該基因過表達能夠使玉米粒中鋅濃度增加31.6%，實現玉米中鋅的生物強化。

3.2 GWAS在動物中的應用

GWAS也是動物遺傳分析的主要研究方法之一，通過定位影響家畜重要經濟性狀的變異位點和候選基因，幫助研究者更好地理解性狀的遺傳機制，在家畜育種和改良方面起到重要的推動作用。

在1 027頭杜洛克和二花臉雜交的F2代群體的GWAS分析中，Ren等[47]發掘出多個與復雜性狀顯著相關的候選基因：使用SAS的GLM模型發現豬耳朵大小受到PPARD進化保守區域變異的影響；使用R包GenABEL發現在糖原分解的級聯激活中發揮作用，其突變導致骨骼肌糖酵解潛能增加，進而影響豬肉產量[48]。Wang等[49]在82頭母豬中使用CMLM進行GWAS分析，鑒定到266個與豬出生重顯著相關的QTL。Guo等[50]分別使用MLM與貝葉斯混合模型，共定位到15個與仔豬數目與死亡率性狀有關的候選基因。Gozalo-Marcilla等[51]使用GEMMA對來自8個具有不同遺傳背景的275 590頭豬的背腩厚度性狀進行GWAS分析，共定位到264個位點與背腩厚度顯著關聯，鑒定了64個與脂肪代謝相關的候選基因。

雞的GWAS分析大多聚焦于生長、產蛋與抗病性狀。Gu等[52]使用烏骨雞與白洛克雞雜交的F2群體，利用PLINK軟件定位到基因組中3個與雞體重相關的區段。該課題組后續進行了多個種雞不同表型性狀的GWAS研究，包括羽毛形態、胡須與雞冠形態等方面[53-54]，同時使用祖先單倍型對雞的生長性狀進行精細定位，篩選出9個關鍵的候選區段供進一步研究[55]。在400只中國雞的生殖性狀GWAS分析中，Fan等[56]利用PLINK軟件鑒定到19個與蛋重、產蛋數等性狀相關的QTL，定位到、、等17個候選基因。Li等[57]使用GEMMA對839只母雞的多個蛋殼晶體結構相關性狀進行GWAS分析，鑒定到GGA1上55.6—69.1 Mb區間內部621個顯著信號，注釋到、、、、和等參與調節胞質鈣離子濃度生物過程的基因。Guo等[58]通過建立肉雞外翻-內翻畸形（valgus-varus deformity，VVD）病例組與對照組，基于233個樣本使用GEMMA進行GWAS分析，篩選到5個與VVD顯著相關的變異，并將GWAS與轉錄組分析整合，定位到重要的易感基因。

羊是中國重要的農業動物，其種類繁多，且具有豐富的遺傳資源[59]。Demars等[60]使用PLINK對2種羔羊開展GWAS后定位到與高產表型和排卵率顯著關聯，揭示了在羔羊卵巢發育功能中的關鍵作用，為探索生育障礙提供了重要依據。He等[61]將GWAS分析用于3個不同品種的中國本土綿羊，使用CMLM在雙角羊和四角羊中鑒定到和2個候選基因，幫助理解綿羊角發育的分子調控。

水產養殖的遺傳收益總體上高于陸地農業動物[62]，水產的抗病育種具有重要的經濟意義。大西洋鮭魚普遍受到細菌性腎病的侵擾，Holborn等[63]使用GenABEL對507只大西洋鮭魚對細菌性腎病下的抗性進行GWAS分析，發現其抗性為多基因性狀，并定位到2個與細菌性腎病抗性顯著相關的QTL。Peng等[64]使用PLINK對黃河鯉魚體重、體長和胴體重量等性狀進行QTL定位，鑒定出多個與神經發育，基礎代謝相關的基因，為黃河鯉魚生長性狀的選育提供了遺傳材料。黃花魚同樣是我國的傳統養殖魚類之一，黃花魚具有明顯的性別二態性，雌性黃花魚的生長性狀顯著優于雄性。LIN等[65]對905只黃花魚（463只雌性，442只雄性）開展GWAS分析，使用TASSEL在第21染色體處發現22個QTL與性別決定顯著相關，鑒定到包括調控精子發育、雌性激素代謝功能的14個候選基因。同時，對不同性別性腺指數的GWAS確定了第18染色體與雄性性腺發育相關的區段，鑒定到相關調節基因、和。

3.3 GWAS在微生物中的應用

微生物與植物的整個生長過程息息相關。植物病原微生物在定植后導致植物產生特定疾病，從而影響植物生長。Davila等[66]使用Fast-LMM在350株擬南芥中鑒定出干旱條件下與灰霉菌（）真菌病原體相關的轉錄因子，其與耐寒及壞死真菌抵抗相關。Zhang等[67]使用701種不同水稻種質和23種不同水稻白葉枯病菌（pv，）菌株進行跨物種GWAS，使用EMMAX篩選出47個毒力相關基因和318個水稻不完全抗性基因，并對毒力相關基因與抗病基因之間的基因互作加以闡述。Martins等[68]使用BLINK對豌豆的派倫霉菌（）抗性以及莖直徑等生長性狀進行GWAS，發現具有抗病性的等位基因導致了較低的株高，印證了植物抗病以犧牲自身生長為代價這一結論[69]。de Ronne等[70]同時使用ECMLM、FarmCPU和BLINK在357個大豆品種中發現了新的大豆疫霉菌（）抗性QTL，該QTL的LD區塊內包含與病原體抗性相關的乳膠蛋白編碼基因。

與病原微生物不同，根際微生物中存在一部分能夠促進植物營養吸收的共生體[71]。根際微生物正向促進了植物的營養吸收能力與范圍，同時調節植物的生長發育與抗性反應[72]。Bergelson等[73]使用EMMAX基于擬南芥的細菌豐度數據鑒定到與，其分別調控植物免疫與側根形成；通過真菌的豐度數據發現SNARE蛋白的靶點，其參與根毛蛋白質轉運。研究發現細菌與真菌豐富度的GWAS結果基本沒有重疊，表明細菌與真菌群落豐富度受到不同基因的影響。Deng等[74]使用GEMMA研究了高粱遺傳位點與根際微生物豐度的關系，并使用高粱遺傳信息成功預測根際微生物組成情況。

腸道益生菌能夠提高動物的飼料轉化率，增加動物體重，以及增產牛奶或雞蛋等農業產品，實現更多經濟價值[75-78]。Crespo-Piazuelo等[79]使用GEMMA對285只豬腸道微生物中18個屬的相對豐度進行GWAS分析，發現、、、、和與基因型存在顯著關聯，定位到包括免疫應答與生理調節相關的多個候選基因。Bergamaschi等[80]在1 028頭豬在斷奶期、生長中期、生長末期3個生長過程的糞便中提取微生物樣本，提取腸道微生物Alpha多樣性與分類操作單元（operational taxonomic unit，OTU）作為表型數據，使用EMMAX鑒定到候選基因。該基因在腸道組織中高度表達，與細胞增殖相關。

3.4 基于單倍型的GWAS應用

復雜性狀GWAS研究有助于了解復雜性狀的遺傳機制，定位到的候選基因為進一步的研究提供了指導方向。單倍型是染色體上共同遺傳的多個等位基因的組合，包含了等位基因間的連鎖信息。單位點GWAS每次檢驗一個SNP，單倍型GWAS（haplotype- based genome-wide association study，hGWAS）把整個單倍型區塊用于GWAS分析，檢測與性狀顯著關聯的單倍型區塊。由于單倍型可能包含了位點之間的互作信息，一些發表的基于單倍型的GWAS分析證明hGWAS在定位效果以及統計意義上均強于基于單位點的GWAS[81-83]。

在水稻中，Yano等[84]把單倍型作為固定效應進行GWAS分析，篩選出4個與水稻農藝性狀相關的候選基因。Ogawa等[85]使用日本的8個高產水稻品種構建了日本-多親高代雜交系（Japan-multi-parent advanced generation inter-cross，JAM），通過8個祖先群體進行hGWAS分析，鑒定到控制糯性胚乳和糯素長度性狀的QTL。與上述方法不同，Zhang等[86]開發了GLASCOW軟件將祖先單倍型加入GLMM作為第二個隨機效應組分，即GRM控制多基因效應而祖先單倍型控制群體分層，在3種單基因隱性疾病中獲得了EMMAX未檢測到的顯著結果。在瘦雞和肥雞兩個肉雞品種的GWAS分析中，Zhang等[87]定位到與腹部脂肪重量顯著關聯的132個單倍型區塊，篩選出7個可能在控制腹部脂肪含量中產生影響的候選基因。Howard等[88]在18 773個蘇格蘭家庭中發現了2種與重性抑郁障礙顯著關聯的單倍型，其中包括與雙相情感障礙相關的單倍型區域，該結果通過25 035個UKB中的個體加以驗證，為揭示重度抑郁癥的遺傳機制提供了思路。

4 展望

從MLM被引入GWAS起，研究人員持續對其進行優化。目前MLM在GWAS中的優化主要包括2種：其一為對隨機效應求解中GRM相關的計算優化；其二為對隨機效應中GRM構建的優化。方法的優化大幅控制了計算結果中的假陽性，并顯著提升了計算速度。對于計算優化部分，fastGWA作為目前最快的MLM算法實現，已經將GWAS中混合線性模型的時間復雜度降低至O（MN），使大規模數據的快速GWAS計算成為現實。伴隨著GWAS的廣泛使用，累積了越來越多的匯總數據（summary statistics），具有高統計功效和計算速度的META-GWAS-MLM算法有待開發以利用這些匯總數據，實現GWAS的聯合分析。另一方面，MLM控制了假陽性卻帶來了假陰性問題，降低了統計功效。近期發表的方法FarmCPU與BLINK嘗試剝離隨機效應，以解決MLM導致的假陰性問題。提高關聯分析的統計功效可以考慮多種方式，使用單倍型的hGWAS可能是提高功效的方法之一。單倍型包含了可能的標記互作信息，比單位點信息量更豐富。在GWAS中引入貝葉斯思想同樣能夠增加GWAS的統計功效，因此，MLM與先驗信息的結合也可能是未來控制假陰性的策略之一。協變量的優化同樣可能是混合模型GWAS的優化內容之一，更好地捕獲復雜數據全局和局部信息能更好地應對群體分層。對于復雜性狀而言，基因多效性和多基因效應是影響表型形成的關鍵因素，目前這方面的研究相對較少。開發和完善多位點混合模型GWAS的方法對于遺傳定位，以及研究復雜性狀的形成機制具有重要作用。

基于混合模型的GWAS已經廣泛應用于植物、動物和微生物的遺傳研究，為生物育種提供了理論基礎和新的思路。精心設計的GWAS試驗有助于更精確地定位到候選基因，比如使用重組自交系（recombinant inbred line，RIL）或深度雜交系（advanced intercross line，AIL）的實驗群體能夠更好地定位候選基因。除此之外，基因與環境互作影響動植物優質性狀形成的分子機制解析是未來GWAS研究的一個重要方向。目前已經廣泛認識到環境對于植物生長性狀具有重要的影響，適宜的環境可以顯著增加作物產量。通過挖掘并改良作物的環境互作基因，以及表型可塑性的相關基因，有助于了解作物的育種潛力，并幫助作物在適宜環境中增加產量，在惡劣環境中維持產量，在育種中有著重要作用。而這些在動物的環境適應性研究中同樣適用。

[1] BOTSTEIN D, RISCH N. Discovering genotypes underlying human phenotypes: past successes for mendelian disease, future approaches for complex disease. Nature Genetics, 2003, 33(3): 228-237.

[2] VISSCHER P M, BROWN M A, MCCARTHY M I, YANG J. Five years of GWAS discovery. The American Journal of Human Genetics, 2012, 90(1): 7-24.

[3] VISSCHER P M, WRAY N R, ZHANG Q, SKLAR P, MCCARTHY M I, BROWN M A, YANG J. 10 years of GWAS discovery: Biology, function, and translation. The American Journal of Human Genetics, 2017, 101(1): 5-22.

[4] YU J, PRESSOIR G, BRIGGS W H, BI I V, YAMASAKI M, DOEBLEY J F, MCMULLEN M D, GAUT B S, NIELSEN D M, HOLLAND J B, KRESOVICH S, BUCKLER E S. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Nature Genetics, 2006, 38(2): 203-208.

[5] LIPPERT C, LISTGARTEN J, LIU Y, KADIE C M, DAVIDSON R I, HECKERMAN D. FaST linear mixed models for genome-wide association studies. Nature Methods, 2011, 8(10): 833-835.

[6] LOH P R, TUCKER G, BULIK-SULLIVAN B K, VILHJáLMSSON B J, FINUCANE H K, SALEM R M, CHASMAN D I, RIDKER P M, NEALE B M, BERGER B, PATTERSON N, PRICE A L. Efficient Bayesian mixed-model analysis increases association power in large cohorts. Nature Genetics, 2015, 47(3): 284-290.

[7] JIANG L, ZHENG Z, QI T, KEMPER K E, WRAY N R, VISSCHER P M, YANG J. A resource-efficient tool for mixed model association analysis of large-scale data. Nature Genetics, 2019, 51(12): 1749-1755.

[8] 卜李那, 趙毅強. 全基因組關聯分析及其擴展方法的研究進展. 農業生物技術學報, 2019, 27(1): 150-158.

BU L N, ZHAO Y Q. Research progress of genome-wide association study and its extension methods. Journal of Agricultural Biotechnology, 2019, 27(1): 150-158. (in Chinese)

[9] CARDON L R, PALMER L J. Population stratification and spurious allelic association. The Lancet, 2003, 361(9357): 598-604.

[10] DEVLIN B, ROEDER K. Genomic control for association studies. Biometrics, 1999, 55(4): 997-1004.

[11] PRICE A L, PATTERSON N J, PLENGE R M, WEINBLATT M E, SHADICK N A, REICH D. Principal components analysis corrects for stratification in genome-wide association studies. Nature Genetics, 2006, 38(8): 904-909.

[12] SHAM P C, PURCELL S M. Statistical power and significance testing in large-scale genetic studies. Nature Reviews Genetics, 2014, 15(5): 335-346.

[13] GAO X, BECKER L C, BECKER D M, STARMER J D, PROVINCE M A. Avoiding the high Bonferroni penalty in genome-wide association studies. Genetic Epidemiology, 2010, 34(1): 100-105.

[14] LIU X L, HUANG M, FAN B, BUCKLER E S, ZHANG Z. Iterative usage of fixed and random effect models for powerful and efficient genome-wide association studies. Plos Genetics, 2016, 12(2): 1-24.

[15] ZHAO K Y, ARANZANA M J, KIM S, LISTER C, SHINDO C, TANG C, TOOMAJIAN C, ZHENG H G, DEAN C, MARJORAM P, NORDBORG M. Anexample of association mapping in structured samples. Plos Genetics, 2007, 3(1): 71-82.

[16] XIAO Y J, LIU H J, WU L J, WARBURTON M, YAN J B. Genome-wide association studies in maize: Praise and stargaze. Molecular Plant, 2017, 10(3): 359-374.

[17] 溫陽俊, 馮建英, 張瑾. 多位點關聯分析方法學的研究進展. 南京農業大學學報, 2022, 45(1): 1-10.

WEN Y J, FENG J Y, ZHANG J. Research progress of multi-locus genome-wide association study. Journal of Nanjing Agricultural University, 2022, 45(1): 1-10. (in Chinese)

[18] KANG H M, ZAITLEN N A, WADE C M, KIRBY A, HECKERMAN D, DALY M J, ESKIN E. Efficient control of population structure in model organism association mapping. Genetics, 2008, 178(3): 1709-1723.

[19] KANG H M, SUL J H, SERVICE S K, ZAITLEN N A, KONG S Y, FREIMER N B, SABATTI C, ESKIN E. Variance component model to account for sample structure in genome-wide association studies. Nature Genetics, 2010, 42(4): 348-354.

[20] ZHANG Z W, ERSOZ E, LAI C Q, TODHUNTER R J, TIWARI H K, GORE M A, BRADBURY P J, YU J M, ARNETT D K, ORDOVAS J M, BUCKLER E S. Mixed linear model approach adapted for genome-wide association studies. Nature Genetics, 2010, 42(4): 355-360.

[21] LI M, LIU X L, BRADBURY P, YU J M, ZHANG Y M, TODHUNTER R J, BUCKLER E S, ZHANG Z W. Enrichment of statistical power for genome-wide association studies. BMC Biology, 2014, 12(1): 1-10.

[22] AULCHENKO Y S, de KONING D J, HALEY C. Genomewide rapid association using mixed model and regression: a fast and simple method for genomewide pedigree-based quantitative trait loci association analysis. Genetics, 2007, 177(1): 577-585.

[23] SVISHCHEVA G R, AXENOVICH T I, BELONOGOVA N M, van DUIJN C M, AULCHENKO Y S. Rapid variance components-based method for whole-genome association analysis. Nature Genetics, 2012, 44(10): 1166-1170.

[24] ZHOU X, STEPHENS M. Genome-wide efficient mixed-model analysis for association studies. Nature Genetics, 2012, 44(7): 821-824.

[25] LISTGARTEN J, LIPPERT C, KADIE C M, DAVIDSON R I, ESKIN E, HECKERMAN D. Improved linear mixed models for genome-wide association studies. Nature Methods, 2012, 9(6): 525-526.

[26] WANG Q S, TIAN F, PAN Y C, BUCKLER E S, ZHANG Z W. A SUPER powerful method for genome wide association study. Plos One, 2014, 9(9): 1-9.

[27] BULIK-SULLIVAN B K, LOH P R, FINUCANE H K, RIPKE S, YANG J, PATTERSON N, DALY M J, PRICE A L, NEALE B M. LD Score regression distinguishes confounding from polygenicity in genome-wide association studies. Nature Genetics, 2015, 47(3): 291-295.

[28] SEGURA V, VILHJáLMSSON B J, PLATT A, KORTE A, SEREN ü, LONG Q, NORDBORG M. An efficient multi-locus mixed-model approach for genome-wide association studies in structured populations. Nature Genetics, 2012, 44(7): 825-830.

[29] HUANG M, LIU X L, ZHOU Y, SUMMERS R M, ZHANG Z W. BLINK: a package for the next level of genome-wide association studies with both individuals and markers in the millions. GigaScience, 2019, 8(2): 1-12.

[30] FU W W, WANG R, XU N Y, WANG J X, LI R, NANAEI H A, NIE Q H, ZHAO X, HAN J L, YANG N, JIANG Y. Galbase: A comprehensive repository for integrating chicken multi-omics data. BMC Genomics, 2022, 23(1): 1-11.

[31] YANG J, LEE S H, GODDARD M E, VISSCHER P M. GCTA: A tool for genome-wide complex trait analysis. American Journal of Human Genetics. 2011, 88(1): 76-82.

[32] CHEN H, WANG C L, CONOMOS M P, STILP A M, LI Z L, SOFER T, SZPIRO A A, CHEN W, BREHM J M, CELEDON J C, REDLINE S, PAPANICOLAOU G J, THORNTON T A, LAURIE C C, RICE K, LIN X H. Control for population structure and relatedness for binary traits in genetic association studies via logistic mixed models. The American Journal of Human Genetics, 2016, 98(4): 653-666.

[33] ZHOU W, NIELSEN J B, FRITSCHE L G, DEY R, GABRIELSEN M E, WOLFORD B N, LEFAIVE J, VANDEHAAR P, GAGLIANO S A, GIFFORD A, BASTARACHE L A, WEI W Q, DENNY J C, LIN M X, HVEEM K, KANG H M, ABECASIS G R, WILLER C J, LEE S. Efficiently controlling for case-control imbalance and sample relatedness in large-scale genetic association studies. Nature Genetics, 2018, 50(9): 1335-1341.

[34] JIANG L D, ZHENG Z L, FANG H L, YANG J. A generalized linear mixed model association tool for biobank-scale data. Nature Genetics, 2021, 53(11): 1616-1621.

[35] BI W J, ZHOU W, DEY R, MUKHERJEE B, SAMPSON J N, LEE S. Efficient mixed model approach for large-scale genome-wide association studies of ordinal categorical phenotypes. The American Journal of Human Genetics, 2021, 108(5): 825-839.

[36] HANSEN M, KRAFT T, GANESTAM S, S?LL T, NILSSON N O. Linkage disequilibrium mapping of the bolting gene in sea beet using AFLP markers. Genetical Research, 2001, 77(1): 61-66.

[37] ATWELL S, HUANG Y S, VILHJáLMSSON B J, WILLEMS G, HORTON M, LI Y, MENG D Z, PLATT A, TARONE A M, HU T T, JIANG R, MULIYATI N W, ZHANG X, AMER M A, BAXTER I, BRACHI B, CHORY J, DEAN C, DEBIEU M, de MEAUX J, ECKER J R, FAURE N, KNISKERN J M, JONES J D G, MICHAEL T, NEMRI A, ROUX F, SALT D E, TANG C L, TODESCO M, TRAW M B, WEIGEL D, MARJORAM P, BOREVITZ J O, BERGELSON J, NORDBORG M. Genome-wide association study of 107 phenotypes ininbred lines. Nature, 2010, 465(7298): 627-631.

[38] ZHAO K Y, TUNG C W, EIZENGA G C, WRIGHT M H, ALI M L, H PRICE A, NORTON G J, ISLAM M R, REYNOLDS A, MEZEY J, MCCLUNG A M, BUSTAMANTE C D, MCCOUCH S R. Genome-wide association mapping reveals a rich genetic architecture of complex traits in. Nature Communications, 2011, 2(1): 1-10.

[39] HUANG X H, WEI X H, SANG T, ZHAO Q, FENG Q, ZHAO Y, LI C Y, ZHU C R, LU T T, ZHANG Z W, LI M, FAN D L, GUO Y L, WANG A H, WANG L, DENG L W, LI W J, LU Y Q, WENG Q J, LIU K Y, HUANG T, ZHOU T Y, JING Y F, LI W, LIN Z, BUCKLER E S, QIAN Q, ZHANG Q F, LI J Y, HAN B. Genome-wide association studies of 14 agronomic traits in rice landraces. Nature Genetics, 2010, 42(11): 961-967.

[40] HUANG X H, ZHAO Y, WEI X H, LI C Y, WANG A H, ZHAO Q, LI W J, GUO Y L, DENG L W, ZHU C R, FAN D L, LU Y Q, WENG Q J, LIU K Y, ZHOU T F, JING Y F, SI L Z, DONG G J, HUANG T, LU T T, FENG Q, QIAN Q, LI J Y, HAN B. Genome-wide association study of flowering time and grain yield traits in a worldwide collection of rice germplasm. Nature Genetics, 2012, 44(1): 32-39.

[41] YONEMARU J I, MIZOBUCHI R, KATO H, YAMAMOTO T, YAMAMOTO E, MATSUBARA K, HIRABAYASHI H, TAKEUCHI Y, TSUNEMATSU H, ISHII T, OHTA H, MAEDA H, EBANA K, YANO M. Genomic regions involved in yield potential detected by genome-wide association analysis in Japanese high-yielding rice cultivars. BMC genomics, 2014, 15(1): 1-12.

[42] WANG M, YAN J B, ZHAO J R, SONG W, ZHANG X B, XIAO Y N, ZHENG Y L. Genome-wide association study (GWAS) of resistance to head smut in maize. Plant Science, 2012, 196(1): 125-131.

[43] LI H, PENG Z Y, YANG X H, WANG W D, FU J J, WANG J H, HAN Y J, CHAI Y C, GUO T T, YANG N, LIU J, WARBURTON M L, CHENG Y B, HAO X M, ZHANG P, ZHAO J Y, LIU Y J, WANG G Y, LI J S, YAN J B. Genome-wide association study dissects the genetic architecture of oil biosynthesis in maize kernels. Nature Genetics, 2013, 45(1): 43-50.

[44] LI Y X, LI C H, BRADBURY P J, LIU X L, LU F, ROMAY C M, GLAUBITZ J C, WU X, PENG B, SHI Y S, SONG Y C, ZHANG D F, BUCKLER E S, ZHANG Z W, LI Y, WANG T Y. Identification of genetic variants associated with maize flowering time using an extremely large multi-genetic background population. The Plant Journal, 2016, 86(5): 391-402.

[45] SEKHON R S, SASKI C, KUMAR R, FLINN B S, LUO F, BEISSINGER T M, ACKERMAN A J, BREITZMAN M W, BRIDGES W C, DE LEON N, KAEPPLER S M. Integrated genome-scale analysis identifies novel genes and networks underlying senescence in maize. The Plant Cell, 2019, 31(9): 1968-1989.

[46] CHAO Z F, CHEN Y Y, JI C, WANG Y L, HUANG X, ZHANG C Y, YANG J, SONG T, WU J C, GUO L X, LIU C B, HAN M L, WU Y R, YAN J B, CHAO D Y. A genome-wide association study identifies a transporter for zinc uploading to maize kernels. Embo Reports, 2023, 24(1): 1-19.

[47] REN J, DUAN Y Y, QIAO R M, YAO F, ZHANG Z Y, YANG B, GUO Y M, XIAO S J, WEI R X, OUYANG Z X, DING N S, AI H S, HUANG L S. A missense mutation in PPARD causes a major QTL effect on ear size in pigs. Plos Genetics, 2011, 7(5): 1-10.

[48] MA J W, YANG J, ZHOU L S, REN J, LIU X X, ZHANG H, YANG B, ZHANG Z Y, MA H B, XIE X H, XING Y Y, GUO Y M, HUANG L S. A splice mutation in the PHKG1 gene causes high glycogen content and low meat quality in pig skeletal muscle. Plos Genetics, 2014, 10(10): 1-13.

[49] WANG X M, LIU X L, DENG D D, YU M, LI X P. Genetic determinants of pig birth weight variability. BMC Genetics, 2016, 17(1): 41-48.

[50] GUO X Y, SU G S, CHRISTENSEN O F, JANSS L, LUND M S. Genome-wide association analyses using a Bayesian approach for litter size and piglet mortality in Danish Landrace and Yorkshire pigs. BMC Genomics, 2016, 17(1): 1-12.

[51] GOZALO-MARCILLA M, BUNTJER J, JOHNSSON M, BATISTA L, DIEZ F, WERNER C R, CHEN C Y, GORJANC G, MELLANBY R J, HICKEY J M, ROS-FREIXEDES R. Genetic architecture and major genes for backfat thickness in pig lines of diverse genetic backgrounds. Genetics, Selection, Evolution, 2021, 53(1): 1-14.

[52] GU X R, FENG C G, MA L, SONG C, WANG Y Q, DA Y, LI H F, CHEN K W, YE S H, GE C R, HU X X, LI N. Genome-wide association study of body weight in chicken F2resource population. Plos One, 2011, 6(7): 1-5.

[53] IMSLAND F, FENG C G, BOIJE H, BED'HOM B, FILLON V, DORSHORST B, RUBIN C J, LIU R R, GAO Y, GU X R, WANG Y Q, GOURICHON D, ZODY M C, ZECCHIN W, VIEAUD A, TIXIER-BOICHARD M, HU X X, HALLB??K F, LI N, ANDERSSON L. The Rose-comb mutation in chickens constitutes a structural rearrangement causing both altered comb morphology and defective sperm motility. Plos Genetics, 2012, 8(6): 1-12.

[54] GUO Y, GU X R, SHENG Z Y, WANG Y Q, LUO C L, LIU R R, QU H, SHU D M, WEN J, CROOIJMANS R P M A, CARLBORG ?, ZHAO Y Q, HU X X, LI N. A complex structural variation on chromosome 27 leads to the ectopic expression of hoxb8 and the muffs and beard phenotype in chickens. Plos Genetics, 2016, 12(6): 1-24.

[55] WANG Y Z, CAO X M, LUO C L, SHENG Z Y, ZHANG C Y, BIAN C, FENG C G, LI J X, GAO F, ZHAO Y Q, JIANG Z Q, QU H, SHU D M, CARLBORG ?, HU X X, LI N. Multiple ancestral haplotypes harboring regulatory mutations cumulatively contribute to a QTL affecting chicken growth traits. Communications Biology, 2020, 3(1): 1-13.

[56] FAN Q C, WU P F, DAI G J, ZHANG G X, ZHANG T, XUE Q, SHI H Q, WANG J Y. Identification of 19 loci for reproductive traits in a local Chinese chicken by genome-wide study. Genetics and Molecular Research, 2017, 16(1): 1-8.

[57] LI Q L, DUAN Z Y, SUN C J, ZHENG J X, XU G Y, YANG N. Genetic variations for the eggshell crystal structure revealed by genome-wide association study in chickens. BMC Genomics, 2021, 22(1): 1-12.

[58] GUO Y P, HUANG H T, ZHANG Z Z, MA Y C, LI J Z, TANG H H, MA H X, LI Z J, LI W T, LIU X J, KANG X T, HAN R L. Genome-wide association study identifies SNPs for growth performance and serum indicators inbroilers () using ddGBS sequencing. BMC Genomics, 2022, 23(1): 1-11.

[59] 張統雨, 朱才業, 杜立新, 趙福平. 羊重要性狀全基因組關聯分析研究進展. 遺傳, 2017, 39(06): 491-500.

ZHANG T Y, ZHU C Y, DU L X, ZHAO F P. Advances in genome-wide association studies for important traits in sheep and goats. Hereditas(Beijing), 2017, 39(6): 491-500. (in Chinese)

[60] DEMARS J, FABRE S, SARRY J, ROSSETTI R, GILBERT H, PERSANI L, TOSSER-KLOPP G, MULSANT P, NOWAK Z, DROBIK W, MARTYNIUK E, BODIN L. Genome-wide association studies identify two novel BMP15 mutations responsible for an atypical hyperprolificacy phenotype in sheep. Plos Genetics, 2013, 9(4): 1-13.

[61] HE X H, ZHOU Z K, PU Y B, CHEN X F, MA Y H, JIANG L. Mapping the four-horned locus and testing the polled locus in three Chinese sheep breeds. Animal Genetics, 2016, 47(5): 623-627.

[62] GJEDREM T. Genetic improvement for the development of efficient global aquaculture: A personal opinion review. Aquaculture, 2012, 344-349(1): 12-22.

[63] HOLBORN M K, ANG K P, ELLIOTT J A K, POWELL F, BOULDING E G. Genome wide association analysis for bacterial kidney disease resistance in a commercial North American Atlantic salmon () population using a 50K SNP panel. Aquaculture, 2018, 495(1): 465-471.

[64] PENG W Z, XU J, ZHANG Y, FENG J X, DONG C J, JIANG L K, FENG J Y, CHEN B H, GONG Y W, CHEN L, XU P. An ultra-high density linkage map and QTL mapping for sex and growth-related traits of common carp (). Scientific Reports, 2016, 6(1): 1-16.

[65] LIN H L, ZHOU Z X, ZHAO J, ZHOU T, BAI H Q, KE Q Z, PU F, ZHENG W Q, XU P. Genome-wide association study identifies genomic loci of sex determination and gonadosomatic index traits in large yellow croaker (). Marine Biotechnology, 2021, 23(1): 127-139.

[66] DAVILA OLIVAS N H, KRUIJER W, GORT G, WIJNEN C L, VAN LOON J J A, DICKE M. Genome-wide association analysis reveals distinct genetic architectures for single and combined stress responses in. New Phytologist, 2017, 213(2): 838-851.

[67] ZHANG F, HU Z Q, WU Z C, LU J L, SHI Y Y, XU J L, WANG X Y, WANG J P, ZHANG F, WANG M M, SHI X R, CUI Y R, VERA CRUZ C, ZHUO D L, HU D D, LI M, WANG W S, ZHAO X Q, ZHENG T Q, FU B Y, ALI J, ZHOU Y L, LI Z K. Reciprocal adaptation of rice andpv.: cross-species 2D GWAS reveals the underlying genetics. The Plant Cell, 2021, 33(8): 2538-2561.

[68] MARTINS L B, BALINT-KURTI P, REBERG-HORTON S C. Genome-wide association study for morphological traits and resistance to Peryonella pinodes in the USDA pea single plant plus collection. G3 Genes|Genomes|Genetics, 2022, 12(9): 1-8.

[69] KARASOV T L, CHAE E, HERMAN J J, BERGELSON J. Mechanisms to mitigate the trade-off between growth and defense. The Plant Cell, 2017, 29(4): 666-680.

[70] DE RONNE M, SANTHANAM P, CINGET B, LABBé C, LEBRETON A, YE H, VUONG T D, HU H F, VALLIYODAN B, EDWARDS D, NGUYEN H T, BELZILE F, BéLANGER R. Mapping of partial resistance to Phytophthora sojae in soybean PIs using whole-genome sequencing reveals a major QTL. The Plant Genome, 2022, 15(1): 1-16.

[71] LIU Q, CHENG L, NIAN H, JIN J, LIAN T X. Linking plant functional genes to rhizosphere microbes: A review. Plant Biotechnology Journal, 2023, 21(5): 902-917.

[72] BAI B, LIU W D, QIU X Y, ZHANG J, ZHANG J Y, BAI Y. The root microbiome: Community assembly and its contributions to plant fitness. Journal of Integrative Plant Biology, 2022, 64(2): 230-243.

[73] BERGELSON J, MITTELSTRASS J, HORTON M W. Characterizing both bacteria and fungi improves understanding of theroot microbiome. Scientific Reports, 2019, 9(1): 1-11.

[74] DENG S W, CADDELL D F, XU G, DAHLEN L, WASHINGTON L, YANG J L, COLEMAN-DERR D. Genome wide association study reveals plant loci controlling heritability of the rhizosphere microbiome. The ISME Journal, 2021, 15(11): 3181-3194.

[75] de FREITAS A S, de DAVID D B, TAKAGAKI B M, ROESCH L F W. Microbial patterns in rumen are associated with gain of weight in beef cattle. Antonie Van Leeuwenhoek, 2020, 113(9): 1299-1312.

[76] MALTECCA C, BERGAMASCHI M, TIEZZI F. The interaction between microbiome and pig efficiency: A review. Journal of Animal Breeding and Genetics, 2020, 137(1): 4-13.

[77] XUE M Y, SUN H Z, WU X H, LIU J X, GUAN L L. Multi-omics reveals that the rumen microbiome and its metabolome together with the host metabolome contribute to individualized dairy cow performance. Microbiome, 2020, 8(1): 1-19.

[78] ELOKIL A A, MAGDY M, MELAK S, ISHFAQ H, BHUIYAN A, CUI L, JAMIL M, ZHAO S, LI S. Faecal microbiome sequences in relation to the egg-laying performance of hens using amplicon-based metagenomic association analysis. Animal, 2020, 14(4): 706-715.

[79] CRESPO-PIAZUELO D, MIGURA-GARCIA L, ESTELLé J, CRIADO-MESAS L, REVILLA M, CASTELLó A, MU?OZ M, GARCíA-CASCO J M, FERNáNDEZ A I, BALLESTER M, FOLCH J M. Association between the pig genome and its gut microbiota composition. Scientific Reports, 2019, 9(1): 1-11.

[80] BERGAMASCHI M, MALTECCA C, SCHILLEBEECKX C, MCNULTY N P, SCHWAB C, SHULL C, FIX J, TIEZZI F. Heritability and genome-wide association of swine gut microbiome features with growth and fatness parameters. Scientific Reports, 2020, 10(1): 1-12.

[81] WANG Y T, SUNG P Y, LIN P L, YU Y W, CHUNG R H. A multi-SNP association test for complex diseases incorporating an optimal P-value threshold algorithm in nuclear families. BMC Genomics, 2015, 16(1): 1-10.

[82] WANG F, MEYER N J, WALLEY K R, RUSSELL J A, FENG R. Causal genetic inference using haplotypes as instrumental variables. Genetic Epidemiology, 2016, 40(1): 35-44.

[83] N'DIAYE A, HAILE J K, CORY A T, CLARKE F R, CLARKE J M, KNOX R E, POZNIAK C J. Single marker and haplotype-based association analysis of semolina and pasta colour in elite durum wheat breeding lines using a high-density consensus map. Plos One, 2017, 12(1): 1-24.

[84] YANO K, YAMAMOTO E, AYA K, TAKEUCHI H, LO P C, HU L, YAMASAKI M, YOSHIDA S, KITANO H, HIRANO K, MATSUOKA M. Genome-wide association study using whole-genome sequencing rapidly identifies new genes influencing agronomic traits in rice. Nature Genetics, 2016, 48(8): 927-934.

[85] OGAWA D, YAMAMOTO E, OHTANI T, KANNO N, TSUNEMATSU H, NONOUE Y, YANO M, YAMAMOTO T, YONEMARU J I. Haplotype-based allele mining in the Japan-MAGIC rice population. Scientific Reports, 2018, 8(1): 1-11.

[86] ZHANG Z, GUILLAUME F, SARTELET A, CHARLIER C, GEORGES M, FARNIR F, DRUET T. Ancestral haplotype-based association mapping with generalized linear mixed models accounting for stratification. Bioinformatics, 2012, 28(19): 2467-2473.

[87] ZHANG H, SHEN L Y, XU Z C, KRAMER L M, YU J Q, ZHANG X Y, NA W, YANG L L, CAO Z P, LUAN P, REECY J M, LI H. Haplotype-based genome-wide association studies for carcass and growth traits in chicken. Poultry Science, 2020, 99(5): 2349-2361.

[88] HOWARD D M, HALL L S, HAFFERTY J D, ZENG Y N, ADAMS M J, CLARKE T K, PORTEOUS D J, NAGY R, HAYWARD C, SMITH B H, MURRAY A D, RYAN N M, EVANS K L, HALEY C S, DEARY I J, THOMSON P A, MCINTOSH A M. Genome-wide haplotype-based association analysis of major depressive disorder in Generation Scotland and UK Biobank. Translational Psychiatry, 2017, 7(11): 1-9.

Principle, Optimization and Application of Mixed Models in Genome- Wide Association Study

College of Biological Sciences, China Agricultural University, Beijing 100193

Genome-wide association study (GWAS) is an effective method to locate genomic loci that are significantly associated with traits. With the accumulated phenotypic data, the continuous development of high-throughput genotyping technology, and the improved statistical methods, it promotes the wide application of GWAS in area of human disease and animal and plant genetics. False positives are one of the important concerns that impair the reliability of genome-wide association results. To control the false positives, in addition to correcting the-values, GWAS models have been continuously improved from the naive methods like ANOVA (for quantitative trait) or Chi-square test (for quality trait), to general linear model (GLM), which incorporates fixed-effect covariates, to the mixed linear model (MLM), which incorporates random effects. Fitting individual genetic effects into random effects defined by the genomic relationships matrix (GRM) is commonly adapted currently. Since the parameter estimation of MLM consumes a lot of computational resources, researchers have tried to optimize solving models and constructing GRM (which also improves computing efficiency), and the time complexity gradually decreased from O(MN3) to O(MN) for MLM-based methods, achieving a great leap in computational speed and statistical efficacy. For inflations caused by unbalanced case-control data, researchers further correct the generalized mixed linear model (GLMM). This paper comprehensively introduces the basic principles and development of GWAS, with specific emphasis on the model improvement and optimization details. We also list the applications of MLM in GWAS in agriculture, including progress on animals, plants and microbes, as well as the application of haplotype in GWAS. Finally, we give prospects on the future developments of GWAS from the viewpoints of further model optimization and experimental design.

genome-wide association study; complex traits; random effects; genomic relationships matrix; mixed linear model

2022-12-04；

2023-03-02

國家重點研發計劃（2022YFF1000204）

譚力治，E-mail：tanlizhi@cau.edu.cn。通信作者趙毅強，E-mail：yiqiangz@cau.edu.cn

10.3864/j.issn.0578-1752.2023.09.001

（責任編輯 李莉）