不同統(tǒng)計(jì)分析模型在大豆菌核病GWAS中的應(yīng)用

張 羽, 周婉瑩, FRAN?OIS Belzile

(1. 陜西理工大學(xué) 生物科學(xué)與工程學(xué)院, 漢中 723000；2. Département de Phytologie, Université Laval, Québec, QC G1V 0A6, Canada)

二代測(cè)序技術(shù)由于其準(zhǔn)確性高、簡(jiǎn)單、快捷、高通量等特點(diǎn)，產(chǎn)生的SNPs(Single nucleotide polymorphism)在全基因組關(guān)聯(lián)分析(Genome-wide association study，GWAS)中的應(yīng)用越來越普遍，尤其是人類疾病通過全基因組關(guān)聯(lián)研究已經(jīng)成為熱點(diǎn)[1]，近年來在動(dòng)植物QTLs(Quantitative trait locus)定位中也廣泛應(yīng)用[2-5]。而SNPs質(zhì)量、數(shù)量及測(cè)序樣本背景、樣本量是影響全基因組關(guān)聯(lián)研究的兩個(gè)重要因素。SNPs的數(shù)據(jù)質(zhì)量控制關(guān)系著將有哪些SNPs參與后續(xù)的關(guān)聯(lián)分析，數(shù)據(jù)質(zhì)量控制通常包括5個(gè)方面：某個(gè)樣本的所有SNPs基因型分型成功率(Genotype,%)、某個(gè)SNP在所有樣本中的分型成功率(call ratio,%)、孟德爾錯(cuò)誤率、HWE(Hardy-Weinberg equilibrium)檢驗(yàn)P值和MAF(Minor allele frequency)。個(gè)體SNP的缺失率是反映DNA樣本質(zhì)量的重要指標(biāo)，如果缺失多，則說明該個(gè)體的DNA樣品質(zhì)量差。常用0.01、0.02或0.05作為界值，剔除缺失率大于界值的個(gè)體。關(guān)于分型準(zhǔn)確性和分型率在測(cè)序公司內(nèi)部也是重要的質(zhì)量控制步驟，通常都能達(dá)到分析對(duì)數(shù)據(jù)質(zhì)量的要求。而HWE檢驗(yàn)P值和MAF取決于研究人員如何根據(jù)測(cè)序樣本的遺傳背景、群體結(jié)構(gòu)作出參數(shù)選擇，這兩個(gè)參數(shù)選擇對(duì)關(guān)聯(lián)分析影響很大[6]。在沒有進(jìn)化影響下，當(dāng)基因一代一代傳遞時(shí)，群體的基因頻率和基因型頻率將保持不變，群體滿足HWE。不滿足HWE的群體，說明可能存在近親交配、遺傳漂移、嚴(yán)重突變、群體分層等因素，這樣的群體代表性差，不能作進(jìn)一步分析，但由于疾病/病害的發(fā)生可能導(dǎo)致遺傳不平衡，特別是對(duì)于自交產(chǎn)生的簡(jiǎn)單群體，偏離HWE的位點(diǎn)很可能是疾病/病害易感位點(diǎn)[7]。本研究的供試材料為自交作物，不考慮HWE，用STRUCTURE(http://taylor0.biology.ucla.edu/structureHarvesteroybase.org/tools.php)分析也顯示群體結(jié)構(gòu)簡(jiǎn)單。在GWAS中，如果MAF很低時(shí)，一方面說明變異小，提供的與疾病/病害關(guān)聯(lián)信息少；另一方面，關(guān)聯(lián)性檢驗(yàn)的統(tǒng)計(jì)學(xué)效能很低。因此，GWAS中需剔除MAF較低的SNPs，目前大多數(shù)研究剔除MAF的界值常選為0.01～0.05[8-9]。

另一個(gè)影響關(guān)聯(lián)分析的重要因素為樣本選擇。從本質(zhì)上講，QTLs分析是一個(gè)統(tǒng)計(jì)意義座位，是以概率標(biāo)準(zhǔn)說明基因組的哪些區(qū)段或位點(diǎn)與哪些數(shù)量性狀緊密關(guān)聯(lián)，理論上樣本數(shù)越大、多態(tài)性位點(diǎn)越多的關(guān)聯(lián)分析可信度越高，但很多研究為了降低研究成本，首先進(jìn)行大樣本量的表型鑒定，然后從中選擇極端表型類型(0/1性狀)進(jìn)行測(cè)序分析和QTLs定位，即選擇基因分型技術(shù)，這樣可以大大降低測(cè)序費(fèi)用。其原理是雖然數(shù)量性狀在一個(gè)自然群體中是連續(xù)變異的，但如果淘汰大多數(shù)中間類型，則高值組和低值組兩種極端表型的個(gè)體就可以明確地區(qū)分開來，分成兩組來分析[10-11]，這種方法類似于人類復(fù)雜疾病分析中的case-control方法。對(duì)每個(gè)QTL而言，在高值表型組中應(yīng)存在較多的高值基因型，而低值組中應(yīng)存在較多的低值基因型。如果某個(gè)標(biāo)記與QTL有連鎖，那么該標(biāo)記與QTL之間就會(huì)發(fā)生一定程度的共分離，于是其基因型分離比例頻率分布會(huì)偏離孟德爾規(guī)律。用卡平方測(cè)驗(yàn)方法對(duì)其中一組檢驗(yàn)這種偏離，就能推斷該標(biāo)記是否與QTL連鎖。本研究用全部樣本和極端表型材料進(jìn)行關(guān)聯(lián)分析。

研究發(fā)現(xiàn)，生物在遺傳過程中通常是很多SNPs聯(lián)系在一起作為一個(gè)整體往下傳遞，是一種遺傳標(biāo)記的非隨機(jī)性組合，即連鎖不平衡(LD，Linkage Disequilibrium)，當(dāng)位于同一條染色體上的兩個(gè)位點(diǎn)/等位基因同時(shí)存在的概率大于群體中因隨機(jī)分布而同時(shí)出現(xiàn)的概率時(shí)，就稱這兩個(gè)位點(diǎn)處于LD狀態(tài)，也稱單體型(Haplotype)[12]。所以，如果在二代技術(shù)產(chǎn)生的數(shù)以萬計(jì)的SNPs基礎(chǔ)上，研究以單體型塊(Blocks)為單位，只要檢測(cè)幾個(gè)標(biāo)簽SNP(tagSNP)，就可以識(shí)別出相應(yīng)的單體型結(jié)構(gòu)，進(jìn)而確定是否與疾病/病害相關(guān)。本研究將數(shù)據(jù)質(zhì)量控制、不同統(tǒng)計(jì)分析模型、不同分析方法、不同樣本等方面進(jìn)行關(guān)聯(lián)研究比較，以期為植物全基因組關(guān)聯(lián)分析提供參考和大豆抗菌核病育種提供理論指導(dǎo)。

1 材料與方法

126個(gè)大豆品種(系)的測(cè)序數(shù)據(jù)和其對(duì)菌核病反應(yīng)的表型值來自于加拿大Laval大學(xué)Fran?ois實(shí)驗(yàn)室。

2 結(jié)果與分析

2.1 數(shù)據(jù)質(zhì)量控制關(guān)聯(lián)分析

用PLINKv1.07(http://www.softpedia.com/get/Science-CAD/PLINK. shtml)進(jìn)行測(cè)試。供試的126個(gè)大豆品種(系)的每個(gè)樣本的所有SNPs基因型分型成功率為100%；每個(gè)SNP在所有樣本中的分型成功率為100%；由于材料特點(diǎn)，孟德爾錯(cuò)誤率為0。因此這3個(gè)參數(shù)對(duì)分析結(jié)果沒有影響。如果不考慮群體情況，HWE檢驗(yàn)對(duì)關(guān)聯(lián)結(jié)果影響最大，但大豆為自交作物，不符合HWE的位點(diǎn)很可能是與表型關(guān)聯(lián)位點(diǎn)，因此本研究不考慮HWE[13-15]。由于自交作物群體樣本較小，遺傳效應(yīng)的存在與否依賴于MAF在全基因組的機(jī)會(huì)分布，為了降低假陽性率，我們重點(diǎn)對(duì)MAF進(jìn)行了測(cè)試，126品種(系)的MAF在全基因組的分部見圖1，大約46.21%的位點(diǎn)MAF小于0.1。MAF取值0.01時(shí)，剩余30 125個(gè)SNPs；MAF取值0.05時(shí)，剩余20 691個(gè)SNPs；MAF取值0.1時(shí)，剩余16 203個(gè)SNPs；取值0.25時(shí)，剩余9105個(gè)SNPs參與關(guān)聯(lián)分析(只算加性效應(yīng)/一般線性模型)。為了較好地比較幾種結(jié)果，我們列出了P值小于1E-05的關(guān)聯(lián)位點(diǎn)(表1)。MAF取0.01(圖2)、0.05(圖3)和0.1(圖4)關(guān)聯(lián)分析結(jié)果完全相同，強(qiáng)關(guān)聯(lián)依次在3-20-1-4-17號(hào)染色體上。MAF取0.25時(shí)(圖5)的強(qiáng)關(guān)聯(lián)依次在20-1-3-4-17染色體上。通過分析發(fā)現(xiàn)，雖然MAF取值0.25時(shí)的強(qiáng)關(guān)聯(lián)染色體次序有所變化，但它們?cè)谙嗤旧w上的關(guān)聯(lián)位點(diǎn)與MAF取值0.01、0.05和0.1是相同的。由于測(cè)序樣本量相對(duì)較小，而SNPs相對(duì)較多，為了降低假陽性率，我們以MAF取值0.01的數(shù)據(jù)為例進(jìn)行了Permutation Test測(cè)試(圖6)，測(cè)試后的強(qiáng)關(guān)聯(lián)染色體依次為1-3-20-4-17，與不進(jìn)行Permutation Test測(cè)試的關(guān)聯(lián)大小染色體排序3-20-1-4-17有區(qū)別，P值稍有增大，但在同一染色體上的關(guān)聯(lián)位點(diǎn)相同，具體見表1。

圖1 126個(gè)大豆品種(系)的次等位基因分布

在GWAS中，數(shù)據(jù)膨脹可能導(dǎo)致假陽性率升高，為了降低假陽性關(guān)聯(lián)，通常用Q-Q plot估計(jì)數(shù)量性狀觀測(cè)值與預(yù)測(cè)值之間的差異。我們比較了不同MAF取值和Permutation test下的Q-Q plot(圖7)，其中Permutation test膨脹最小，其次依次為MAF0.01、MAF0.05、MAF0.1，膨脹最大為MAF0.25。因此后續(xù)的顯著關(guān)聯(lián)位點(diǎn)/候選基因注釋應(yīng)參考Permutation test結(jié)果。

圖3 MAF取值0.05用PLINK分析的SNP-Trait關(guān)聯(lián)結(jié)果圖

圖4 MAF取值0.1用PLINK分析的SNP-Trait關(guān)聯(lián)結(jié)果圖

圖5 MAF取值0.25用PLINK分析的SNP-Trait關(guān)聯(lián)結(jié)果圖

2.2 Haplotype-Trait關(guān)聯(lián)分析結(jié)果

由于在單個(gè)SNP-Trait關(guān)聯(lián)研究中，MAF取值0.01、0.05和0.1的結(jié)果一致，在Haplotype-Trait關(guān)聯(lián)中，我們以MAF取值0.05為例分析。Haplotype用Block定義，較強(qiáng)關(guān)聯(lián)依次在17-1-3-4-20號(hào)染色體上(表2)，與單個(gè)SNP-Trait關(guān)聯(lián)大小染色體排序3-20-1-4-17有差異，但峰值SNP都包含在相應(yīng)的Haplotype中，即Haplotype的tagSNP與單個(gè)SNP-Trait關(guān)聯(lián)位點(diǎn)是相同的，但由于在某個(gè)Haplotype中測(cè)序的某些材料的個(gè)別SNPs位點(diǎn)是雜合的，導(dǎo)致Haplotype解釋的關(guān)聯(lián)P值稍有增大，使得Haplotype-Trait關(guān)聯(lián)與SNP-Trait關(guān)聯(lián)程度大小染色體排序有差異。

圖6 MAF取值0.01permutation 1 000 000PLINK分析的SNP-Trait關(guān)聯(lián)結(jié)果圖

圖7 Q-Q plot

表1 不同MAF取值和Permutation下SNP-Trait關(guān)聯(lián)結(jié)果

表2 基于Haplotype-Trait的關(guān)聯(lián)結(jié)果(P<1E-5)

2.3 極端類型關(guān)聯(lián)分析

126個(gè)大豆品種(系)對(duì)菌核病的表型值從28.6～192.4 mm，表型變異呈連續(xù)性，極端類型通常從連續(xù)變異表型的兩端各取15%～35%。本研究的極端表型從126個(gè)樣本中從兩頭各取10%(26個(gè)樣本)、20%(50個(gè)樣本)和30%(76個(gè)樣本)進(jìn)行分析。由于極端表型(0/1性狀)模擬質(zhì)量性狀，表型值用0、1和2表示，分別為：1代表unaffected；2代表affected；0代表unknown，中間類型無法判斷表型值，因此，用HAPLOVIEW4.2軟件(http://www.broadinstitute.org/haploview)分析極端類型關(guān)聯(lián)。

26個(gè)、50個(gè)和76個(gè)品種(系)的SNP-Trait的較強(qiáng)關(guān)聯(lián)比較見圖8。可以看出，用較少極端類型分析的關(guān)聯(lián)位點(diǎn)都包含在了用較多極端類型分析的關(guān)聯(lián)位點(diǎn)內(nèi)，樣本量越大，關(guān)聯(lián)出的位點(diǎn)越多，且P值越小。由于在HAPLOVIEW4.2里，把極端類型的基因型材料的表現(xiàn)值模擬為質(zhì)量性狀(非此即彼性狀)，因此關(guān)聯(lián)出的P值都偏低(表3)。

26個(gè)、50個(gè)和76個(gè)品種(系)的基于Haplotype-Trait的關(guān)聯(lián)分析見表4，Pvalue隨著分析的樣本量的增加而減小。

3 討論

通過全基因組關(guān)聯(lián)進(jìn)行QTLs定位的方法本質(zhì)上是一個(gè)統(tǒng)計(jì)意義座位，是以概率標(biāo)準(zhǔn)說明在基因組的哪些位點(diǎn)/區(qū)段可能存在影響哪些數(shù)量性狀的位點(diǎn)/候選基因，從理論上講樣本數(shù)越大、DNA多態(tài)性位點(diǎn)越多的關(guān)聯(lián)分析可信度越高。遺傳學(xué)中認(rèn)為變異大于1%的為多態(tài)性，如果MAF取值過小會(huì)把突變作為多態(tài)性，取值過大可能會(huì)漏掉現(xiàn)實(shí)存在的多態(tài)性位點(diǎn)，從研究分析結(jié)果看MAF取值0.01～0.05最為合適，驗(yàn)證了前人的結(jié)果。

圖8 不同極端類型基于SNP-Trai關(guān)聯(lián)圖

極端材料在QTLs定位中已經(jīng)有報(bào)道。從本研究的分析可以看出，全部樣本關(guān)聯(lián)出的位點(diǎn)包含了用極端材料關(guān)聯(lián)出的位點(diǎn)，用來進(jìn)行QTLs定位的材料數(shù)目越多，關(guān)聯(lián)位點(diǎn)也越多，因此，用極端材料進(jìn)行QTLs定位有可能把有的關(guān)聯(lián)位點(diǎn)漏掉，同時(shí)，在選取極端材料時(shí)，需要提前進(jìn)行大量的表型鑒定和選擇，才能保證選到真正的極端材料進(jìn)行全基因組關(guān)聯(lián)研究，而對(duì)于表型鑒定過程復(fù)雜的一些性狀，會(huì)花費(fèi)更多的人力物力。另外，如果研究的性狀屬于數(shù)量性狀遺傳，那么在一個(gè)大的群體中，性狀變異呈一個(gè)連續(xù)的變異，大豆為自交作物，理論上中間類型會(huì)越來越少，兩端極端類型越來越多且頻率接近，因此，需要選取較多的極端類型進(jìn)行分析。從本研究看，30%的極端類型關(guān)聯(lián)出的位點(diǎn)出現(xiàn)在全部材料關(guān)聯(lián)出的較強(qiáng)關(guān)聯(lián)位點(diǎn)里。

表3 不同極端類型基于SNP-Trait的關(guān)聯(lián)結(jié)果

表4 不同極端類型基于Haplotype-Trait的關(guān)聯(lián)結(jié)果

在GWAS中，Haplotype發(fā)揮了重要的作用。一方面Haplotype包含了更多的SNPs，另一方面降低了分析的自由度。在Haplotype應(yīng)用中，分型和頻率推斷是關(guān)鍵，雖然現(xiàn)有的一些軟件如TASSEL5.0和FLAPJACK等有可視性功能，能讓研究人員形象的看到一些Haplotype的結(jié)構(gòu)，進(jìn)行關(guān)聯(lián)位點(diǎn)的直觀驗(yàn)證，但距離較遠(yuǎn)的及復(fù)雜的Haplotype還需軟件分析，不同軟件有不同的Haplotype定義方法，Blocks的劃分需要研究人員根據(jù)材料考慮。

不同統(tǒng)計(jì)分析模型關(guān)聯(lián)結(jié)果有差異，所以建議在關(guān)聯(lián)分析時(shí)多用幾種模型分析，然后優(yōu)先確定共同發(fā)現(xiàn)的QTLs。在一個(gè)基因組區(qū)域沒有研究清楚之前，由于生物轉(zhuǎn)錄模板鏈的選擇及選擇性剪切等因素影響，在假定候選基因關(guān)聯(lián)出來后再根據(jù)已知的生物信息學(xué)知識(shí)推導(dǎo)可能的候選基因，然后實(shí)驗(yàn)加以驗(yàn)證。

4 結(jié)論

用2種方法(基于單個(gè)SNP-Trait和Haplotype-Trait)和2種表型(全部表型和極端表型)的關(guān)聯(lián)分析，共同定位的較強(qiáng)關(guān)聯(lián)位點(diǎn)(只顯示峰值SNP)有Chr1：5589567；Chr3：34387780；Chr4：6353873；Chr13：3784700；Chr17：5734897；Chr20：42091969。這些可能的顯著SNPs位點(diǎn)還需多群體、大樣本的重復(fù)驗(yàn)證，最終找到穩(wěn)定的SNPs標(biāo)記為大豆抗菌核病育種工作服務(wù)。