不同统计分析模型在大豆菌核病GWAS中的应用

张 羽, 周婉莹, 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参与后续的关联分析，数据质量控制通常包括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的群体，说明可能存在近亲交配、遗传漂移、严重突变、群体分层等因素，这样的群体代表性差，不能作进一步分析，但由于疾病/病害的发生可能导致遗传不平衡，特别是对于自交产生的简单群体，偏离HWE的位点很可能是疾病/病害易感位点[7]。本研究的供试材料为自交作物，不考虑HWE，用STRUCTURE(http://taylor0.biology.ucla.edu/structureHarvesteroybase.org/tools.php)分析也显示群体结构简单。在GWAS中，如果MAF很低时，一方面说明变异小，提供的与疾病/病害关联信息少；另一方面，关联性检验的统计学效能很低。因此，GWAS中需剔除MAF较低的SNPs，目前大多数研究剔除MAF的界值常选为0.01～0.05[8-9]。

另一个影响关联分析的重要因素为样本选择。从本质上讲，QTLs分析是一个统计意义座位，是以概率标准说明基因组的哪些区段或位点与哪些数量性状紧密关联，理论上样本数越大、多态性位点越多的关联分析可信度越高，但很多研究为了降低研究成本，首先进行大样本量的表型鉴定，然后从中选择极端表型类型(0/1性状)进行测序分析和QTLs定位，即选择基因分型技术，这样可以大大降低测序费用。其原理是虽然数量性状在一个自然群体中是连续变异的，但如果淘汰大多数中间类型，则高值组和低值组两种极端表型的个体就可以明确地区分开来，分成两组来分析[10-11]，这种方法类似于人类复杂疾病分析中的case-control方法。对每个QTL而言，在高值表型组中应存在较多的高值基因型，而低值组中应存在较多的低值基因型。如果某个标记与QTL有连锁，那么该标记与QTL之间就会发生一定程度的共分离，于是其基因型分离比例频率分布会偏离孟德尔规律。用卡平方测验方法对其中一组检验这种偏离，就能推断该标记是否与QTL连锁。本研究用全部样本和极端表型材料进行关联分析。

研究发现，生物在遗传过程中通常是很多SNPs联系在一起作为一个整体往下传递，是一种遗传标记的非随机性组合，即连锁不平衡(LD，Linkage Disequilibrium)，当位于同一条染色体上的两个位点/等位基因同时存在的概率大于群体中因随机分布而同时出现的概率时，就称这两个位点处于LD状态，也称单体型(Haplotype)[12]。所以，如果在二代技术产生的数以万计的SNPs基础上，研究以单体型块(Blocks)为单位，只要检测几个标签SNP(tagSNP)，就可以识别出相应的单体型结构，进而确定是否与疾病/病害相关。本研究将数据质量控制、不同统计分析模型、不同分析方法、不同样本等方面进行关联研究比较，以期为植物全基因组关联分析提供参考和大豆抗菌核病育种提供理论指导。

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染色体上。通过分析发现，虽然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有区别，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。因此后续的显著关联位点/候选基因注释应参考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，表型变异呈连续性，极端类型通常从连续变异表型的两端各取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定位的方法本质上是一个统计意义座位，是以概率标准说明在基因组的哪些位点/区段可能存在影响哪些数量性状的位点/候选基因，从理论上讲样本数越大、DNA多态性位点越多的关联分析可信度越高。遗传学中认为变异大于1%的为多态性，如果MAF取值过小会把突变作为多态性，取值过大可能会漏掉现实存在的多态性位点，从研究分析结果看MAF取值0.01～0.05最为合适，验证了前人的结果。

图8 不同极端类型基于SNP-Trai关联图

极端材料在QTLs定位中已经有报道。从本研究的分析可以看出，全部样本关联出的位点包含了用极端材料关联出的位点，用来进行QTLs定位的材料数目越多，关联位点也越多，因此，用极端材料进行QTLs定位有可能把有的关联位点漏掉，同时，在选取极端材料时，需要提前进行大量的表型鉴定和选择，才能保证选到真正的极端材料进行全基因组关联研究，而对于表型鉴定过程复杂的一些性状，会花费更多的人力物力。另外，如果研究的性状属于数量性状遗传，那么在一个大的群体中，性状变异呈一个连续的变异，大豆为自交作物，理论上中间类型会越来越少，两端极端类型越来越多且频率接近，因此，需要选取较多的极端类型进行分析。从本研究看，30%的极端类型关联出的位点出现在全部材料关联出的较强关联位点里。

表3 不同极端类型基于SNP-Trait的关联结果

表4 不同极端类型基于Haplotype-Trait的关联结果

在GWAS中，Haplotype发挥了重要的作用。一方面Haplotype包含了更多的SNPs，另一方面降低了分析的自由度。在Haplotype应用中，分型和频率推断是关键，虽然现有的一些软件如TASSEL5.0和FLAPJACK等有可视性功能，能让研究人员形象的看到一些Haplotype的结构，进行关联位点的直观验证，但距离较远的及复杂的Haplotype还需软件分析，不同软件有不同的Haplotype定义方法，Blocks的划分需要研究人员根据材料考虑。

不同统计分析模型关联结果有差异，所以建议在关联分析时多用几种模型分析，然后优先确定共同发现的QTLs。在一个基因组区域没有研究清楚之前，由于生物转录模板链的选择及选择性剪切等因素影响，在假定候选基因关联出来后再根据已知的生物信息学知识推导可能的候选基因，然后实验加以验证。

4 结论

用2种方法(基于单个SNP-Trait和Haplotype-Trait)和2种表型(全部表型和极端表型)的关联分析，共同定位的较强关联位点(只显示峰值SNP)有Chr1：5589567；Chr3：34387780；Chr4：6353873；Chr13：3784700；Chr17：5734897；Chr20：42091969。这些可能的显著SNPs位点还需多群体、大样本的重复验证，最终找到稳定的SNPs标记为大豆抗菌核病育种工作服务。