張 羽, 周婉瑩, FRAN?OIS Belzile

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

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

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

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

1 材料與方法

126個大豆品種(系)的測序數據和其對菌核病反應的表型值來自于加拿大Laval大學Fran?ois實驗室。

2 結果與分析

2.1 數據質量控制關聯分析

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

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

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

圖3 MAF取值0.05用PLINK分析的SNP-Trait關聯結果圖

圖4 MAF取值0.1用PLINK分析的SNP-Trait關聯結果圖

圖5 MAF取值0.25用PLINK分析的SNP-Trait關聯結果圖

2.2 Haplotype-Trait關聯分析結果

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

圖6 MAF取值0.01permutation 1 000 000PLINK分析的SNP-Trait關聯結果圖

圖7 Q-Q plot

表1 不同MAF取值和Permutation下SNP-Trait關聯結果

表2 基于Haplotype-Trait的關聯結果(P<1E-5)

2.3 極端類型關聯分析

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

26個、50個和76個品種(系)的SNP-Trait的較強關聯比較見圖8?？梢钥闯?，用較少極端類型分析的關聯位點都包含在了用較多極端類型分析的關聯位點內，樣本量越大，關聯出的位點越多，且P值越小。由于在HAPLOVIEW4.2里，把極端類型的基因型材料的表現值模擬為質量性狀(非此即彼性狀)，因此關聯出的P值都偏低(表3)。

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

3 討論

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

圖8 不同極端類型基于SNP-Trai關聯圖

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

表3 不同極端類型基于SNP-Trait的關聯結果

表4 不同極端類型基于Haplotype-Trait的關聯結果

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

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

4 結論

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