趙敬麗，李淑玲，高　進(jìn)，楊潤清,3*

（1.南京農(nóng)業(yè)大學(xué)無錫漁業(yè)學(xué)院，江蘇　無錫　214081;2.東北農(nóng)業(yè)大學(xué)生命科學(xué)學(xué)院，哈爾濱　150030;3.中國水產(chǎn)科學(xué)研究院生物技術(shù)研究中心，北京　100141）

在遺傳學(xué)領(lǐng)域，全基因組關(guān)聯(lián)分析（Genomewide association study，GWAS）借助單核苷酸多態(tài)性（Single nucleotide ploymorphism，SNP）分子遺傳標(biāo)記，在全基因組范圍作表型與基因型間相關(guān)性分析，發(fā)現(xiàn)影響復(fù)雜性狀基因位點(diǎn)（Quantitative trait nucleotide，QTN）或主效基因。近年來，GWAS在人類復(fù)雜疾病性狀和動物重要經(jīng)濟(jì)性狀中應(yīng)用廣泛。受群體分層和個體間親緣關(guān)系等混雜因素影響，表型與基因型間簡單回歸分析難以適應(yīng)復(fù)雜GWAS。因此，將線性混合模型（Linear mixed model,LMM）引入GWAS[1]。通過剖析由群體分層和個體間親緣關(guān)系所致的混雜偏差，LMM從大量數(shù)據(jù)中分離真實(shí)信號，有效控制QTN定位的Ⅰ型錯誤率（TypeⅠerror）和Ⅱ型錯誤率（TypeⅡ error）。使用全基因組標(biāo)記[2]計(jì)算的實(shí)現(xiàn)親緣關(guān)系矩陣（Realised relationship matrix，RRM）[1]估計(jì)剩余微效多基因方差，基于約束最大似然法（Restricted maximum likelihood,REML）[3]的LMM求解復(fù)雜。因此GRAMMAR[4]、 EMMA[5]、 EMMAX[6]、 CMLM[7]、 P3D[7]、FaST-LMM[8]、 GEMMA[9]、 GRAMMAR-Gamma[10]和BOLT-LMM[11]等簡化算法相繼出現(xiàn)。

大多數(shù)改進(jìn)算法均通過降階LMM實(shí)現(xiàn)簡化計(jì)算。例如：CMLM法根據(jù)實(shí)現(xiàn)親緣關(guān)系矩陣將個體聚類成組，將動物模型轉(zhuǎn)化成公畜模型。EMMAX法利用由基因組最佳線性無偏預(yù)測（Genomic best linear unbiased prediction,GBLUP）[12]計(jì)算的全基因組遺傳方差估計(jì)剩余微效多基因方差，將線性混合模型轉(zhuǎn)化成當(dāng)前檢驗(yàn)標(biāo)記效應(yīng)的剩余誤差方差不齊的線性回歸模型。GRAMMAR法則將GBLUP估計(jì)剩余項(xiàng)作為新表型值并對每個SNP標(biāo)記作回歸分析，此后將剩余微效多基因方差固定為基因組遺傳方差，GRAMMAR-Gamma和BOLT-LMM法為對該算法改進(jìn)。

與REML相比，EMMA法通過對表型值和標(biāo)記指示變量作譜分解，避免每輪迭代中似然函數(shù)計(jì)算涉及冗余矩陣運(yùn)算，成倍量級地提高線性混合模型求解速度。與譜分解每個檢驗(yàn)標(biāo)記EMMA法不同，F(xiàn)aST-LMM法僅需一次譜分解便可檢驗(yàn)所有標(biāo)記。GEMMA算法在譜分解基礎(chǔ)上對目標(biāo)似然函數(shù)作二階求導(dǎo)，尋求全局最優(yōu)。所有算法均通過預(yù)先估計(jì)基因組遺傳力固定剩余微效多基因方差，提高全基因組高通量標(biāo)記分析計(jì)算效率。雖然與完整（非簡化）LMM統(tǒng)計(jì)效力相同，但過高估計(jì)剩余微效多基因效應(yīng)，降低表型擬合優(yōu)度。本研究用剩余多基因遺傳力代替方差比值，將多基因遺傳力求解限制在（0,1）區(qū)間。在FaST-LMM框架內(nèi)，引入R語言RcppArmadillo程序包[13]中極速線性模型擬合函數(shù)（fastLmPure函數(shù)）快速估計(jì)SNP效應(yīng)和完整LMM最大似然值。可通過http://th.archive.ubuntu.com/cran/web/packages/RcppArmadillo/RcppArmadillo.pdf查找fastLmPure函數(shù)具體用法。

1　理論方法

1.1　基因組混合模型

在校正群體分層等非遺傳固定效應(yīng)后，構(gòu)建如下含有當(dāng)前檢驗(yàn)SNP加性遺傳效應(yīng)一般線性模型：

y為n個個體校正表型值向量，μ為群體均值，β為被檢驗(yàn)標(biāo)記加性遺傳效應(yīng)，u為隨機(jī)剩余多基因效應(yīng)，且u～N(0,)，其中K為由遺傳標(biāo)記構(gòu)建RRM，為未知剩余微效多基因方差，ε為剩余誤差效應(yīng)向量，且 ε～N(0,)，其中I為單位矩陣，為誤差方差；1為n維列向量，X和Z分別為β和u指示變量矩陣。

在隨機(jī)效應(yīng)分布假設(shè)下，校正表型值協(xié)方差如下:

對K作譜分解即K=USUT，其中S為包含K按降序排列的特征根對角矩陣，U為特征根對應(yīng)特征向量矩陣。根據(jù)UUT=I，協(xié)方差矩陣如下形式：

對表型值y和指示變量矩陣X作譜分解轉(zhuǎn)換可得y?=UTy和 X?=UT[1 X]，則（1）式將變?yōu)槿缦翷MM:

因此，通過加權(quán)最小二乘法獲得β和σ2ε最大似然估計(jì)值，表示如下：

對數(shù)似然函數(shù)進(jìn)一步簡化：

此式關(guān)于h2函數(shù)。因此，通過對h2在區(qū)間（0,1）內(nèi)一維掃描，優(yōu)化上述與h2有關(guān)對數(shù)似然函數(shù)，找到h2最大似然估計(jì)值。同時，利用與優(yōu)化遺傳力h2相對應(yīng)β和，統(tǒng)計(jì)推斷當(dāng)前檢驗(yàn)SNP遺傳效應(yīng)。

1.2　極速回歸掃描法

通過搜索剩余多基因遺傳力最優(yōu)估計(jì)值，單個SNP的LMM可采用重復(fù)加權(quán)最小二乘方法實(shí)現(xiàn)QTN統(tǒng)計(jì)推斷。對于高通量SNPs而言，逐個SNP的GWMMAS變成優(yōu)化每個SNP對應(yīng)的剩余多基因遺傳力全基因組回歸掃描問題。為提高GWMMAS計(jì)算效率，在求解上述模型（5）時引入極速線性模型擬合函數(shù)（fastLmPure函數(shù)），加快對式（6）中當(dāng)前檢驗(yàn)SNP效應(yīng)和LMM極大似然值重復(fù)加權(quán)最小二乘估計(jì)。與常規(guī)線性模型擬合函數(shù)（lm函數(shù)）相比，fastLmPure函數(shù)運(yùn)行速度提高幾十倍。fastLmPure函數(shù)僅輸出當(dāng)前檢驗(yàn)SNP遺傳效應(yīng)和標(biāo)準(zhǔn)差，，2logL，t值和p值等重要統(tǒng)計(jì)量需要運(yùn)行fastLmPure函數(shù)后間接計(jì)算。

理論上，當(dāng)前檢驗(yàn)SNP剩余多基因遺傳力等于性狀基因組遺傳力和該SNP遺傳力（由SNP效應(yīng)解釋的表型方差比例）之間差異。盡管剩余多基因遺傳力在高通量SNPs間不同，但仍接近性狀基因組遺傳力，因?yàn)槌齉TN外的大多數(shù)SNPs對數(shù)量性狀無作用。因此，性狀基因組遺傳力必須在無效模型（不含SNP效應(yīng)LMM）下預(yù)先估計(jì)。在GWMMAS中，從基因組遺傳力處出發(fā)向下一步或幾步掃描，便可快速找到最優(yōu)剩余多基因遺傳力。

每個SNP剩余多基因遺傳力被固定為基因組遺傳力后，前文中快速回歸掃描被簡化為EMMAX算法[6]，其全基因組掃描速度在不優(yōu)化剩余多基因遺傳力情況下，通過fastLmPure函數(shù)達(dá)最高值。由于基因組遺傳力中已涵蓋當(dāng)前檢驗(yàn)SNP變異，因此EMMAX法SNP檢測效力略降。該方法估計(jì)P值可作為每個SNP快速回歸掃描參考。本研究僅選擇大效應(yīng)或較低顯著水平（0.05或0.01）SNPs以進(jìn)一步提高GWMMAS計(jì)算效率，優(yōu)化其剩余多基因遺傳力。至此，計(jì)算時間復(fù)雜度變?yōu)镺（imn），其中i為全基因組回歸掃描迭代次數(shù)（1＜i≤2）。

1.3　程序設(shè)計(jì)

在R環(huán)境中，按照上述基因組混合模型求解步驟，使用fastLmPure函數(shù)執(zhí)行GWMMAS，執(zhí)行式（5）線性模型求解。實(shí)現(xiàn)親緣關(guān)系矩陣K被計(jì)算為K=scale(X)T?scale(X)，其中m為標(biāo)記個數(shù)。求解K特征根S和特征向量，將y和X分別譜變換為y?和X?。給定一個剩余多基因遺傳力h2，可產(chǎn)生fastLmPure函數(shù)因變量y*=和自變量X*=。設(shè)計(jì)極速回歸求解LMM子程序，即fastLmPure函數(shù)執(zhí)行GWMMAS子程序如下：

lmm＜-function(ystar,xstar,w){

fit＜-fastLmPure(y=ystar,X=xstar)

effects＜-fit$coefficients

residual＜-ystar-xstar*effects

ve＜-var(residual)

std＜-fit$stderr

t＜-effects[2]/std[2]

p＜-2*pnorm(abs(t),lower.tail=FALSE)

logL＜-log(det(w))+nobs*log(ve)

}

對每一個SNP，可以從基因組遺傳力估計(jì)值出發(fā)，向下快速搜尋SNP剩余多基因遺傳力最大似然估計(jì)值?；蚪M遺傳力也適用于EMMAX和GRAMMAR算法?；跓o效模型y=1μ+Zu+ε，基因組遺傳力快速估計(jì)程序由fastLmPure函數(shù)編碼。如果目標(biāo)數(shù)量性狀遺傳力來自系譜，基因組遺傳力可通過在給定遺傳力周圍幾步掃描快速獲得?；蚪M估計(jì)育種值可預(yù)測為GV-1(y-1μ)，同樣適用于GRAMMAR算法剩余項(xiàng)。

2　模擬驗(yàn)證

為評價(jià)極速回歸掃描法性能，利用玉米基因組數(shù)據(jù)集模擬驗(yàn)證，數(shù)據(jù)集可從URL:http://www.panzea.org/#！genotypes/cctl下載。該數(shù)據(jù)集提供2 648個體的681 258個多態(tài)SNPs標(biāo)記分型結(jié)果。從這些標(biāo)記中隨機(jī)抽取300 000 SNPs組成標(biāo)記基因型數(shù)據(jù)集。將100和1 000 QTNs隨機(jī)放置在除6和8號以外染色體SNP上。假設(shè)被模擬QTNs可解釋60%表型變異，從shape=1.66和scale=0.4伽馬分布中抽取其遺傳效應(yīng)。同時給定群體均值為0，誤差方差為5，再根據(jù)模擬QTL基因型隨機(jī)產(chǎn)生目標(biāo)性狀個體表型值。所有模擬均在一個CentOS 6.5操作系統(tǒng)展開，該系統(tǒng)運(yùn)行于2.60 GHz的Intel Xeon E5-2660 Opteron（tm）處理器，512 GB內(nèi)存和20 TB硬盤的服務(wù)器。去除數(shù)據(jù)輸入、RRM計(jì)算與表型值和標(biāo)記指示變量譜分解所需時間，極速回歸掃描法作回歸掃描所需時間為1.986 min，明顯低于R語言中l(wèi)m函數(shù)線性模型求解時間（21.289 min）。

2.1　群體分層判別

系譜群體GWAS中群體分層造成目標(biāo)性狀與無關(guān)基因間假關(guān)聯(lián),GWAS假陽性率較高。所謂群體分層（又稱為群體結(jié)構(gòu)）是指一個群體中亞群之間在等位基因頻率上存在系統(tǒng)性差異。假陽性率是指被錯誤鑒定為真QTN的假Q(mào)TN數(shù)與實(shí)際假Q(mào)TN數(shù)目比率。膨脹系數(shù)（基因組控制λGC）可判別群體分層是否廣泛適用。λGC是所有SNP檢驗(yàn)統(tǒng)計(jì)量中位數(shù)（或均數(shù)）與理論分布中位數(shù)比值，實(shí)際研究中被定義為所有SNP卡方統(tǒng)計(jì)量均值[14]。當(dāng)λGC≈1時，表明群體分層不存在；當(dāng)λGC＞1時表明存在群體分層或其他混雜變量如家系結(jié)構(gòu)或隱含親緣關(guān)系。Q-Q（Quantile-Quantile）圖是將檢驗(yàn)統(tǒng)計(jì)量可視化的一種標(biāo)準(zhǔn)圖形技術(shù)，用于判斷兩個數(shù)據(jù)集是否來自具有共同分布種群，反映分析模型的統(tǒng)計(jì)性質(zhì)（見圖1）。

圖1　100和1 000 QTN設(shè)置下Q-Q圖Fig.1　Q-Q plots under the 100 and 1 000 QTN settings

由圖1a和b可知，黃色和綠色線嚴(yán)重偏離理論線，表明群體分層存在；對于藍(lán)色和紅色線而言，大多數(shù)標(biāo)記觀測P值可較好擬合理論閾值線，少數(shù)標(biāo)記嚴(yán)重偏離理論線，表明無群體分層影響，且有QTN存在。

2.2　模型統(tǒng)計(jì)性質(zhì)比較

Liu等介紹3種混雜水平設(shè)置下無效分布[15]。本研究采取第3種方法，將遺傳標(biāo)記劃分為QTN和非QTN區(qū)標(biāo)記，位于非QTN區(qū)標(biāo)記推導(dǎo)無效分布。在模擬中，將6和8號染色體設(shè)置為非QTN區(qū)域，該區(qū)域包括50 000個SNPs。借助非QTN區(qū)域Q-Q圖評估是否出現(xiàn)假陰性或假陽性。圖1分別展示QTN和非QTN區(qū)域Q-Q關(guān)系。不論QTN區(qū)還是非QTN區(qū)，極速回歸掃描法和EMMAX法優(yōu)于PLINK和GRAMMAR法，表明極速回歸掃描法和EMMAX法可有效擬合群體分層效應(yīng)。由非QTN區(qū)域Q-Q圖可見，PLINK法假陽性率明顯偏高，而GRAMMAR法假陰性率較高。這些錯誤率隨模擬QTNs數(shù)目增加而升高。

2.3　統(tǒng)計(jì)檢測效力比較

QTN統(tǒng)計(jì)檢測效力被定義為檢測到QTN數(shù)占標(biāo)記總數(shù)比例。圖2描述QTN統(tǒng)計(jì)檢測效力與Ⅰ型錯誤率關(guān)系。PLINK法檢測效力與極速回歸掃描法相近，但其假陽性率較高；GRAMMAR法檢測效力低于極速回歸掃描法。盡管EMMAX法與極速回歸掃描法Q-Q圖幾乎重疊，但圖2 c和d展示極速回歸掃描法檢測到比EMMAX法更多QTNs。本研究未與GRAMMAR-Gamma和BOLT-LMM算法比較，因兩算法雖可獲得與極速回歸掃描法相近的統(tǒng)計(jì)檢測效力，但計(jì)算效率較GRAMMAR算法低。

圖2　100和1 000 QTNs設(shè)置下不同一類錯誤率統(tǒng)計(jì)檢測效力Fig.2　Powers against different levels of Type I error under the 100 and 1 000 QTN settings

3　牙鲆生長性狀GWMMAS

3.1　GWAS在魚類經(jīng)濟(jì)性狀中應(yīng)用

牙鲆（Paralichthys olivaceus）是我國北方重要海水養(yǎng)殖魚類，其生長和耐溫等多數(shù)重要經(jīng)濟(jì)性狀，均屬于數(shù)量性狀。探索與牙鲆重要經(jīng)濟(jì)性狀關(guān)聯(lián)基因或QTL，可了解這些性狀遺傳機(jī)制，提高標(biāo)記輔助育種效率。近年來，GWAS在魚類等低等脊椎動物經(jīng)濟(jì)性狀研究中逐漸增多，主要集中于抗病性、抗逆性、發(fā)育、性成熟和生長特性等方面[16]。盡管相關(guān)研究尚處于初步階段[17]，GWAS已應(yīng)用于大西洋鮭、虹鱒等少數(shù)魚類生長研究，如Gutierrez等研究大西洋鮭生長相關(guān)GWAS，篩選到1個位于Ssa13連鎖群SNP位點(diǎn)[18]；在虹鱒中，GWAS檢測到1個位于omy5位點(diǎn)，可分別解釋10、13月齡體重1.4%和1%變異[19]。此外，還有F1代雜交叉尾鮰高溫耐受性GWAS，通過EMMAX方法共檢測到3個顯著性SNP位點(diǎn)[20]。其中，位于14號連鎖群位點(diǎn)可解釋釋表型變異12.1%，另外兩個位于16號連鎖群SNPs分別可解釋表型變異11.3%和11.5%?，F(xiàn)有魚類GWAS研究主要針對單一性狀作主效QTN篩選，有關(guān)牙鲆鲆經(jīng)濟(jì)性狀GWAS研究未見報(bào)道。牙鲆基因組現(xiàn)已由Shao等測序并組裝完成[21]，利用基因組信息在全基因組范圍內(nèi)高精度定位牙鲆重要經(jīng)濟(jì)性狀基因并估計(jì)基因組育種值，可為牙鲆經(jīng)濟(jì)性狀分子遺傳解析提供有利條件。

3.2　牙鲆雙單倍體群體建立

圖3　牙鲆2個形態(tài)性狀Q-QFig.3　Q-Q plots of two morphological traits in Japanese flounder

本研究以威海圣航水產(chǎn)科技有限公司建立的牙鲆育種群作為驗(yàn)證極速回歸掃描法性能材料。自2014年開始利用不同來源牙鲆原種群，在兩代閉鎖繁育基礎(chǔ)上，通過靜水壓方法作雌核發(fā)育牙鲆誘導(dǎo)，將紫外線照射滅活牙鲆精子與正常卵受精，受精后靜水壓處理，誘導(dǎo)染色體加倍，培育雙單倍體（DH）牙鲆。按出生時間分為春季和秋季兩批，共162個體。4月齡時，對162尾牙鲆注射電子標(biāo)記，提取鰭條組織DNA。隨后將標(biāo)記后個體混養(yǎng)在6 m×6 m×1 m水泥池，自然光周期下，利用循環(huán)水系統(tǒng)養(yǎng)殖，水溫5～24℃。養(yǎng)殖期間，每天投喂兩次商業(yè)化餌料至飽食。自注射電子標(biāo)記起，電子秤定期稱量每個個體體重；在參考標(biāo)尺下，定期用數(shù)碼相機(jī)從一定高度向下垂直拍攝個體形態(tài)性狀。直到36月齡，測量2～9次。每次測量前，利用0.05%MS-222將待測個體麻醉，避免處理壓力。

根據(jù)拍攝圖像，利用ImageJ軟件獲得不同日齡全長、體長、頭長、體高和體厚等969個表型記錄。采用簡化基因組2b-RAD高通量標(biāo)記分型技術(shù)，獲得17 297個多態(tài)SNP標(biāo)記，參考牙鲆基因組建立多態(tài)標(biāo)記物理圖譜。

3.3　結(jié)果分析

分別采用PLINK、EMMAX、GRAMMAR和極速回歸掃描4種算法對牙鲆體長（BL）和體高（BH）作GWMMAS，4種算法統(tǒng)計(jì)性質(zhì)比較見圖3。與模擬結(jié)果類似，PLINK法假陽性率較高，EMMAX法產(chǎn)生一定程度假陰性，GRAMMAR法假陰性程度更高，極速回歸掃描法表現(xiàn)最佳。

在GWAS中，通常采用曼哈頓圖直觀反映所有標(biāo)記顯著性情況。圖4～5分別為4種GWAS方法對牙鲆BL和BH基因定位結(jié)果。對于BL而言，EMMAX與GRAMMAR法并未檢測到任何與該性狀相關(guān)聯(lián)基因位點(diǎn)，而PLINK法雖檢測到更多高顯著水準(zhǔn)標(biāo)記，但并未發(fā)現(xiàn)與性狀顯著關(guān)聯(lián)位點(diǎn)，極速回歸掃描法在2和20號染色體上鑒定到可能調(diào)控BL基因位點(diǎn)；對于BH，雖然PLINK法與極速回歸掃描法均檢測到顯著位點(diǎn)，但由圖5可知，極速回歸掃描法檢測標(biāo)記的顯著性明顯高于PLINK法。該數(shù)據(jù)充分論證極速回歸掃描法在基因定位方面優(yōu)勢，且模型可有效控制QTN定位假陽性率。

圖4　牙鲆體長4種比較方法曼哈頓圖Fig.4　Manhattan plots of body length in Japanese flounder for the four compared algorithms

圖5　牙鲆體高4種比較方法曼哈頓圖Fig.5　Manhattan plots of body height in Japanese flounder for the four compared algorithms

4　結(jié)論

在GWMMAS中，無論采用REML法還是譜變換法求解每個標(biāo)記對應(yīng)剩余多基因方差，隨標(biāo)記數(shù)目增多，計(jì)算時間增加。盡管線性混合模型簡化求解方法提高計(jì)算效率，但均降低剩余多基因方差估計(jì)精度及QTN檢測效率。本研究提出極速回歸掃描法以數(shù)量性狀基因組遺傳力為出發(fā)點(diǎn)，探究每個SNP最優(yōu)剩余多基因遺傳力，將GWMMAS轉(zhuǎn)換成高效回歸掃描。通過使用R語言Rcp-pArmadillo程序包中極速線性模型擬合函數(shù)（fastLmPure函數(shù)），縮短計(jì)算時間。