張敬言，宋禹昕，張姮妤，張瑩，楊潤清，5

（1.黑龍江八一農(nóng)墾大學(xué)生命科學(xué)技術(shù)學(xué)院，大慶163319；2.南京農(nóng)業(yè)大學(xué)無錫漁業(yè)學(xué)院；3.黑龍江八一農(nóng)墾大學(xué)理學(xué)院；4.黑龍江八一農(nóng)墾大學(xué)動(dòng)物科學(xué)技術(shù)學(xué)院；5.中國水產(chǎn)科學(xué)研究院生物技術(shù)研究中心）

在全基因組關(guān)聯(lián)研究（GWAS）中，群體分層導(dǎo)致了檢測統(tǒng)計(jì)量的膨脹增加了關(guān)聯(lián)檢測的假陽性率，從而降低了對數(shù)量性狀核苷酸（QTNs）檢測的統(tǒng)計(jì)效力，它通常用個(gè)體在亞群中的比例表示[1-4]?；诩蚁禂?shù)據(jù)的GWAS 分析所關(guān)注的家系內(nèi)的信息不受分層影響，因?yàn)橥粋€(gè)家系內(nèi)遺傳和非遺傳等位基因都具有共同的遺傳祖先[5-8]。在非基于家系的關(guān)聯(lián)分析研究中，普遍使用全基因組遺傳標(biāo)記計(jì)算出的主成分（PCs）分析方法矯正群體分層的影響[9-10]。主成分分析（PCA）的GWAS 雖然只包括前幾個(gè)主成分帶來的群體分層混淆效應(yīng)，它通過親緣關(guān)系矩陣模擬個(gè)體間的遺傳協(xié)方差，所以這種矯正方法可以被認(rèn)為是主成分在表型值上的回歸[11-17]。在實(shí)際應(yīng)用基于主成分分析的GWAS 過程中，即使因?yàn)槔霉潭ㄐ?yīng)模型的協(xié)變量數(shù)只包含了前幾個(gè)主成分，即使這幾個(gè)主成分不和樣本數(shù)在同一個(gè)數(shù)量級，也仍然保持了合理的統(tǒng)計(jì)效力[18]，群體分層效應(yīng)可以在使用相關(guān)主成分較少的主成分分析的GWAS 模型中得到很好的糾正[19-20]。

間斷性狀如復(fù)雜疾病，抗病和異常行為是基因定位的重要目標(biāo)性狀，盡管這類可遺傳的性狀表型值由二進(jìn)制或不連續(xù)的數(shù)字表示，但它們的遺傳并不能通過簡單的孟德爾遺傳解釋，基因定位時(shí)同樣應(yīng)考慮除檢測標(biāo)記外的群體分層，家系結(jié)構(gòu)和親緣相關(guān)性[21-22]。許多動(dòng)植物的這類性狀可以以數(shù)量性狀的形式定量統(tǒng)計(jì)和描述，例如玉米對赤霉素感染的抗性可以用接種節(jié)間感染面積占總面積的比例表示，而水稻對稻瘟病菌感染的抗性則通過病損植株的數(shù)目占總植株的比例表示。然而大部分二元響應(yīng)變量的性狀無法定量表示這類性狀，與正態(tài)分布的數(shù)量性狀用剩余誤差作為矯正的表型不同，由于間斷性狀的自變量和因變量之間的尺度不同，無法簡單估計(jì)疾病表型的剩余誤差，這就需要將廣義線性模型（GLM）中的logit 回歸用于這類性狀的全基因組關(guān)聯(lián)分析。即使校正了一些臨床固定效應(yīng)協(xié)變量，logit 回歸仍然和簡單線性回歸一樣會受到群體分層導(dǎo)致檢驗(yàn)統(tǒng)計(jì)量的膨脹的影響。PLINK1.9 是一個(gè)被廣泛使用分析復(fù)雜疾病數(shù)據(jù)軟件，雖然1.9 版本的PLINK 絕大多數(shù)情況都是分析數(shù)據(jù)高效而穩(wěn)定，但當(dāng)使用上百個(gè)主成分矯正分層時(shí)，計(jì)算速度大大下降；同時(shí)PLINK 在水平很多的分類協(xié)變量時(shí)，需要轉(zhuǎn)換為多個(gè)二進(jìn)制虛變量，這也無形的增加了需要同時(shí)分析的協(xié)變量數(shù)目。最后如果分析的協(xié)變量和表型值之間相關(guān)極大時(shí)，PLINK 的回歸模型中協(xié)變量的回歸系數(shù)會遠(yuǎn)大于1，這導(dǎo)致標(biāo)記的效應(yīng)的估計(jì)值降低。針對這些問題描述了一種Logit 回歸能夠同時(shí)修正大量協(xié)變量的高效算法，我們擴(kuò)展了PLINK 功能，讓它可以高速實(shí)現(xiàn)。

1 材料和方法

1.1 材料

試驗(yàn)數(shù)據(jù)群體研究的核心包括Wellcome Trust Case Control Consortium（WTCCC）的4 種疾病，有：冠心病、類風(fēng)濕性關(guān)節(jié)炎、一型糖尿病和二型糖尿病。每一種疾病的病例樣本都是從英國各地隨機(jī)選擇的人群。這些試驗(yàn)的控制組有兩個(gè)來源：1 500 份是1958 年英國出生隊(duì)列的代表性樣本，1 500 份是英國聯(lián)邦三大國家血液服務(wù)機(jī)構(gòu)招募的獻(xiàn)血者。其中總樣本包含3 950 人，其中2 415 名男性，1 535 名女性，1 594 名患病個(gè)體和2 356 名對照個(gè)體。經(jīng)過質(zhì)量控制，關(guān)聯(lián)分析了313 518 個(gè)SNP。

模擬的基因型數(shù)據(jù)集是來自康奈爾大學(xué)的2 468 株玉米群體，經(jīng)過質(zhì)量控制后有300 000 SNP位點(diǎn)，模擬由對系譜中的每個(gè)個(gè)體，在一定長度的染色體片段上，按照重組率生成中等密度和高密度的遺傳標(biāo)記基因型，同時(shí)，在給定多個(gè)QTL 或標(biāo)記的位置和對數(shù)量性狀的效應(yīng)值、剩余多基因、永久環(huán)境和誤差方差組分的條件下，隨機(jī)產(chǎn)生個(gè)體目標(biāo)性狀的表型值。由于試驗(yàn)主要研究在QTNs 的位置和遺傳效應(yīng)發(fā)生變化時(shí)，新方法的性能如何在每一輪的模擬過程中，模擬基因組遺傳力固定為0.8 不變，保持固定的樣本發(fā)病率比例為50%，QTN 的數(shù)目分別設(shè)定為20 個(gè)、200 個(gè)和2 000 個(gè)三個(gè)水平，不同數(shù)目的QTNs 被隨機(jī)地放置在模擬基因組的全部多態(tài)SNP上，而且從形狀和尺度參數(shù)分別為1.66 和0.4 的伽馬分布中，隨機(jī)抽取每個(gè)模擬QTN 的遺傳效應(yīng)。

1.2 方法

1.2.1 基因組廣義線性模型

根據(jù)廣義線性模型理論，連接函數(shù)建立起了疾病性狀表型向量y 和被檢驗(yàn)標(biāo)記效應(yīng)的關(guān)系：

這就是二岐性狀基因組logit 回歸模型。其中，μ是每個(gè)個(gè)體疾病性狀的平均值組成的向量，β 是固定效應(yīng)向量，包含被檢驗(yàn)標(biāo)記的遺傳效應(yīng)βSNP；X 是由β 所對應(yīng)的系數(shù)矩陣。

在μ=Xβ 處，對似然函數(shù)進(jìn)行二階泰勒展開并化簡為：

根據(jù)迭代重加權(quán)最小二乘法估計(jì)回歸系數(shù)的估計(jì)值向量β：

和正態(tài)隨機(jī)變量的誤差方差σe2：

其中，dfe是誤差自由度，W=μ（1-μ）是權(quán)重，y*=是由表型值y 和回歸系數(shù)組成的新的因變量。

采用推斷卡方統(tǒng)計(jì)量，推斷SNP 與疾病性狀關(guān)聯(lián)：

此處，xSNPβ?SNP是被檢驗(yàn)標(biāo)記的回歸項(xiàng)，這個(gè)統(tǒng)計(jì)量服從自由度為1 的卡方分布。

1.2.2 考慮多個(gè)協(xié)變量快速基因組廣義線性模型

考慮群體分層時(shí)，通常以親緣關(guān)系矩陣的前i 個(gè)特征向量U1∶i作為logit 模型的協(xié)變量去矯正由群體分層導(dǎo)致的分層效應(yīng)，基因組廣義線性模型變?yōu)椋?/p>

由于特征向量間相互獨(dú)立，每一個(gè)協(xié)變量的回歸系數(shù)可以分別估計(jì)為不像正態(tài)分布的數(shù)量性狀那樣用剩余誤差作為矯正的表型，因?yàn)樽宰兞亢鸵蜃兞恐g的尺度不同無法估計(jì)疾病表型的剩余誤差。因此，我們把作為已知的協(xié)變量，只需要用去代替y*后，再利用廣義線性模型估計(jì)每個(gè)SNP 的遺傳效應(yīng)和標(biāo)準(zhǔn)誤。

1.2.3 復(fù)雜群體分層的快速矯正步驟

由于每個(gè)主成分是之間是相互獨(dú)立的，我們可以使用分別求解出每個(gè)加權(quán)主成分對應(yīng)的分層效應(yīng)值計(jì)算作為一個(gè)已知的協(xié)變量。我們提供兩種擴(kuò)展方式：

（1）基于PLINK，利用PLINK 軟件自身的協(xié)變量回歸算法，這時(shí)幾百個(gè)協(xié)變量回歸模型等價(jià)轉(zhuǎn)化成單個(gè)協(xié)變量回歸模型。

1.2.4 主成分選擇標(biāo)準(zhǔn)

為了綜合評價(jià)幾種方法的優(yōu)劣性，以及比較優(yōu)化基因組控制的廣義線性混合模型主成分回歸法自身的穩(wěn)健性，對于統(tǒng)計(jì)特性評價(jià)采用如下兩個(gè)標(biāo)準(zhǔn)驗(yàn)證：

（1）基因組控制效果：通過全基因組統(tǒng)計(jì)量的平均卡方值控制值（gc 值）接近1 的程度，并結(jié)合QQ 圖和理論線的偏離程度判斷基因組數(shù)據(jù)中存在的群體分層的校正效果。

（2）QTL 檢測個(gè)數(shù)：同一群體，不存在明顯的假陽性的條件下，以Bonferroni 閾值為標(biāo)準(zhǔn)比較QTN的檢測個(gè)數(shù)。

2 結(jié)果與分析

2.1 模擬結(jié)果

模擬方法包括4 種：直接運(yùn)行PLINK 的GLM 模型算法，PLINK 自帶的矯正PC 協(xié)變量的GLM 模型算法（PLINK-PCs），基于PLINK 矯正PC 協(xié)變量的GLM 模型算法（Bs-PLINK-PCs）和擴(kuò)展PLINK 矯正PC 協(xié)變量的GLM 模型算法（Ex-PLINK-PCs）。從圖1 上半部分的Q-Q 圖可以看出直接運(yùn)行PLINK 的GLM 模型算法的黃色線明顯從起始位置就開始上揚(yáng)，并且平均基因組控制值分別等于1.94，1.71 和1.69，這兩點(diǎn)表明群體內(nèi)存在明顯的群體分層需要矯正。而后三種方法的QQ 圖幾乎完全一致的重合且前端都很好的貼合理論線，基因組控制值也都在1 附近，這說明這幾種方法均沒有出現(xiàn)明顯的假陽性或者假陰性。在檢測效力上除了在模擬20 個(gè)QTNs 時(shí)擴(kuò)展PLINK 矯正PC 協(xié)變量的GLM 模型算法（Ex-PLINK-PCs）的藍(lán)色會有極其微小的提高外其他情況三種方法均完全重合，這說明了只采用主成分作為協(xié)變量時(shí)，我們提出的兩種方案均可以替代PLINK自帶的矯正PC 協(xié)變量的GLM 模型算法，且在只有主成分作為協(xié)變量時(shí)，我們提出的兩種方案結(jié)果也幾乎沒有差別。

圖1 使用玉米基因組模擬100 輪后2 種方法表現(xiàn)的Q-Q 圖和QTNs 檢測效力圖，其中黃色代表直接運(yùn)行PLINK 的GLM 模型算法，綠色代表基于PLINK 矯正PC 協(xié)變量的GLM 模型算法（Bs-PLINK-PCs）。Fig.1 Q-Q and Power plots of two methods after simulated 100 times using maize dataset.Yellow represents for PLINK GLM model algorithm，green represents corrected PC covariate GLM model algorithm by our method based on PLINK（Bs-PLINK-PCs）.

2.2 實(shí)際資料結(jié)果

表1 為五種方法在4 組數(shù)據(jù)的基因組控制值的比較，從表1 中可看出：無論是哪個(gè)群體，在不放任何協(xié)變量的廣義線性回歸模型的基因組控制值均大于1，這表明各個(gè)群體中都有一定的群體分層現(xiàn)象存在，這一現(xiàn)象導(dǎo)致了統(tǒng)計(jì)量的膨脹。同時(shí)隨著引入的主成分協(xié)變量增多，基因組控制值逐漸降低，群體分層的影響也逐漸變小，但是無論是將1 個(gè)，5 個(gè)還是10 個(gè)主成分作為協(xié)變量都無法完全的矯正統(tǒng)計(jì)量的膨脹，然而當(dāng)主成分協(xié)變量的數(shù)目增多到200 個(gè)時(shí)，基因組控制值達(dá)到1 左右，良好的矯正了群體分層效應(yīng)。由于前4 種方法的檢驗(yàn)統(tǒng)計(jì)量造成的假陽性，檢測到的潛在QTN 不具有可比性，所以只比較200 個(gè)主成分作為協(xié)變量的廣義線性模型回歸QTN 檢測能力，見表2。

表1 五種方法的基因組在4 組數(shù)據(jù)的基因組控制值Table 1 Genome control values of five methods using four datasets

表2 基于PLINK 矯正PC 協(xié)變量的GLM 模型算法在4 組數(shù)據(jù)定位p 值超過閾值潛在QTN 的個(gè)數(shù)Table 2 Number of potential QTNs which p-value above the Bonferroni threshold for corrected PC covariate GLM model algorithm by our method based on PLINK

2.2.1 冠狀動(dòng)脈粥樣硬化性心臟病疾病群體（CAD）結(jié)果

在冠狀動(dòng)脈粥樣硬化性心臟病疾病群體中，采用200 個(gè)PC 作為協(xié)變量明顯的矯正了檢驗(yàn)統(tǒng)計(jì)量的膨脹，這令我們的結(jié)果的基因組控制值十分接近1。

2.2.2 類風(fēng)濕關(guān)節(jié)炎（RA）結(jié)果

類風(fēng)濕關(guān)節(jié)炎群體中，檢驗(yàn)統(tǒng)計(jì)量膨脹不如冠狀動(dòng)脈粥樣硬化性心臟病疾病群體明顯，但是采用10 個(gè)主成分作為協(xié)變量仍然無法完全矯正群體分層，采用200 個(gè)主成分作為協(xié)變量雖然矯正了檢驗(yàn)統(tǒng)計(jì)量的膨脹，但基因組控制值為0.999 7，略低于1有假陰性趨勢，可能進(jìn)一步需要減少主成分的個(gè)數(shù)，可能影響潛在QTN 總數(shù)。

2.2.3 一型糖尿?。═1D）結(jié)果

一型糖尿病群體結(jié)果和類風(fēng)濕性關(guān)節(jié)炎群體類似，統(tǒng)計(jì)量都是只存在略微的膨脹，這令200 個(gè)PC作為協(xié)變量的基因組控制值收縮到0.997 5。

2.2.4 二型糖尿?。═2D）結(jié)果

二型糖尿病群體統(tǒng)計(jì)量膨脹適中，采用200 個(gè)主成分在4 個(gè)群體中的基因組控制值為0.999 9，最接近1。

圖2 使用廣義線性模型法分析冠狀動(dòng)脈粥樣硬化性心臟病疾病群體時(shí)的統(tǒng)計(jì)量膨脹在Q-Q 圖上的體現(xiàn)，左圖為右圖淺紅色區(qū)域的放大Fig.2 The representation of false positives on Q-Q plot using generalized linear model in the CAD disease population and the left image is an enlargement of the light red area of the right

圖3 冠狀動(dòng)脈粥樣硬化性心臟病疾病使用200 個(gè)主成分作為協(xié)變量的GLM 模型法的曼哈頓圖與Q-Q 圖Fig.3 Manhattan plot and Q-Q plot of CAD using generalized linear model with 200 principal components as covariates

圖4 使用廣義線性模型法分析類風(fēng)濕關(guān)節(jié)炎疾病群體時(shí)的統(tǒng)計(jì)量膨脹在Q-Q 圖上的體現(xiàn)，左圖為右圖淺紅色區(qū)域的放大Fig.4 The representation of false positives on Q-Q plot using generalized linear model in the RA disease population and the left image is an enlargement of the light red area of the right

圖5 類風(fēng)濕關(guān)節(jié)炎疾病使用200 個(gè)主成分作為協(xié)變量的GLM 模型Fig.5 Manhattan plot and Q-Q plot of RA disease by using generalized linear model with 200 principal components as covariates

圖6 使用廣義線性模型法分析一型糖尿病疾病群體時(shí)的統(tǒng)計(jì)量膨脹在Q-Q 圖上的體現(xiàn)，左圖為右圖淺紅色區(qū)域的放大Fig.6 The representation of false positives on Q-Q plot using generalized linear model in the T1d disease population and the left image is an enlargement of the light red area of the right

圖7 一型糖尿病疾病使用200 個(gè)主成分作為協(xié)變量的GLM 模型的曼哈頓圖與Q-Q 圖Fig.7 Manhattan and Q-Q plot of T1d disease using generalized linear model with 200 principal components as covariates

圖8 使用廣義線性模型法分析二型糖尿病疾病群體時(shí)的統(tǒng)計(jì)量膨脹在Q-Q 圖上的體現(xiàn)，左圖為右圖淺紅色區(qū)域的放大Fig.8 The representation of false positives on Q-Q plot using generalized linear model in the T2d disease population and the left image is an enlargement of the light red area of the right

圖9 二型糖尿病疾病使用200 個(gè)主成分作為協(xié)變量的GLM 模型的曼哈頓圖與Q-Q 圖Fig.9 Manhattan and Q-Q plot of T2D disease using generalized linear model with 200 principal components as covariates

2.2.5 計(jì)算時(shí)間和內(nèi)存消耗比較

這里借助PLINK 軟件分析比較了使用直接使用PLINK 軟件進(jìn)行廣義線性模型分析，在已知主成分的情況下，使用PLINK 軟件自帶的功能進(jìn)行同時(shí)考慮100 個(gè)和200 個(gè)主成分作為協(xié)變量的廣義線性模型分析，經(jīng)過我們拓展的考慮200 個(gè)主成分作為協(xié)變量的廣義線性模型分析在同等情況下：Intel Core i5-8350U 1.70GHz，16GBRAM 的個(gè)人筆記本電腦上比較使計(jì)算所需時(shí)間和軟件運(yùn)行峰值所占用的內(nèi)存。這里其中總樣本包含3 950 人，關(guān)聯(lián)分析了313 518 個(gè)SNP。由表3 可知，無論何種情況，使用廣義線性模型所需的內(nèi)存均不大，并且使用PLINK 自帶的考慮主成分作為協(xié)變量的廣義線性模型分析時(shí)，隨著同時(shí)考慮的主成分增加，雖然內(nèi)存占用沒有增多，但方法所需時(shí)間會隨之大大增加。新方法即使同時(shí)考慮200 個(gè)協(xié)變量，速度也和不考慮任何協(xié)變量的PLINK 廣義線性模型速度和內(nèi)存差別不大，比PLINK 自帶功能速度快200 多倍。因?yàn)榈讓映绦蜻壿嫴煌?，R 語言的程序慢于C 語言編寫的PLINK 軟件，所以表中沒有對比這種方案和PLINK 軟件的速度。

表3 不同方法在關(guān)聯(lián)分析3 950 人的313 518 個(gè)SNP 所需的計(jì)算時(shí)間和峰值內(nèi)存Table 3 Calculation time and peak memory occupied（RAM）using different methods to analysis 3 950 people with 313 518 SNPs

3 討論

通過以上的研究，展示了快速的廣義線性模型的主成分回歸法可以很好的矯正關(guān)聯(lián)分析中的群體分層，結(jié)果顯示無論是評價(jià)指標(biāo)中的Q-Q 圖還是基因組控制值，都顯示只要考慮作為協(xié)變量的主成分?jǐn)?shù)目恰當(dāng)，就可以很好地控制全基因組關(guān)聯(lián)分析中的假陽性。提出的兩種方法檢測到的QTNs 都幾乎一致等于PLINK 自帶方法，并且方法在實(shí)際測試中分析200 個(gè)PCs 作為協(xié)變量時(shí)，速度最多比PLINK 快200 倍以上。人類的實(shí)際數(shù)據(jù)分析結(jié)果中有兩組數(shù)據(jù)的基因組控制值都略低于1，導(dǎo)致這種情況的原因可能是一致的選擇200 個(gè)主成分作為協(xié)變量，但是不同群體的所需矯正的主成分協(xié)變量數(shù)目與群體本身有關(guān)，所以200 個(gè)主成分并非是所有群體的最優(yōu)的主成分?jǐn)?shù)目。在后續(xù)的研究中可以以研究不同群體矯正群體分層的最優(yōu)主成分?jǐn)?shù)目為目標(biāo)再拓展該方法的實(shí)用性。方法雖然借助PLINK 軟件，但是保持PLINK 計(jì)算內(nèi)存消耗不大的優(yōu)勢下，又比其自帶多協(xié)變量回歸速度上也極大提高。同時(shí)主成分是實(shí)現(xiàn)親緣關(guān)系矩陣普分解得到的特征向量，實(shí)現(xiàn)親緣關(guān)系矩陣的計(jì)算方法有很多，在后續(xù)研究中親緣關(guān)系矩陣計(jì)算上也可以尋求簡化算法。同時(shí)從模擬中發(fā)現(xiàn)，提出的兩種方案在僅僅包含主成分作為協(xié)變量時(shí)差別不大，這是由于主成分和表型值間的關(guān)系不大，但是假如協(xié)變量和表型值間的回歸系數(shù)遠(yuǎn)大于1時(shí)會極大的影響QTNs 的檢測效力。提出的第二種方案現(xiàn)在慢于PLINK 軟件，但是并非算法不簡捷，相反的少求解一個(gè)回歸系數(shù)會令方程更簡單，所以這種方法如果擴(kuò)展到PLINK 的快速回歸計(jì)算中，不但會比用R 語言編寫的程序計(jì)算效率更高，理論上會比PLINK 軟件本身的速度進(jìn)一步提高。在日后的研究里，可以嘗試適當(dāng)擴(kuò)大整體樣本數(shù)，這可以增加結(jié)果的精確程度。同時(shí)，也可以嘗試對方法進(jìn)行改良和優(yōu)化，期望新方法可以在不同群體中都選擇最合適的主成分?jǐn)?shù)目。

4 結(jié)論

1.9 版本的PLINK 也是一個(gè)通常用于研究疾病對照試驗(yàn)基因定位的工具，在性能和兼容性方面有了顯著的改進(jìn)，雖然它可以糾正數(shù)百個(gè)協(xié)變量，也可以糾正來自多級分類協(xié)變量的虛變量，但隨著協(xié)變量數(shù)目增加到一定數(shù)目，掃描的速度也會大大降低。在此基礎(chǔ)上拓展了PLINK 軟件，讓它同時(shí)快速的修正大量協(xié)變量去矯正關(guān)聯(lián)分析中存在的群體分層，極大地提高了間斷性狀廣義線性模型全基因組關(guān)聯(lián)分析的運(yùn)算效率，操作簡單方便，易于理解。