极智课堂 01期丨关于GWAS那些事儿~

谈到GWAS你是否仍然对一些问题懵懵懂懂？是否对某些概念有只知其名不解其来历的困惑？本期我们整理了GWAS相关的一些问题，内含表型分布正态性检验小tips、单体型（haplotype）概念解析，一起来看看吧~~

Q1. GWAS分析的优势？

GWAS是自然群体性状定位常用的分析方法，它的优势是不用专门构建遗传群体，直接采集自然的群体材料即可，并且可以同时满足多性状定位。

Q2. 表型采集建议？

表型性状一般可以分为数量性状、质量性状、分级性状几类，不同类型的表型数据在采集时有不同的注意事项。

（1） 数量性状：

一般由多基因控制，能够测量得到具体数值，数据分布是连续的，受环境影响大。

① 尽量保证样本材料在同一环境下培育或养殖；

② 整体表型结果应尽量符合正态分布；

③ 如有多年多点数据，可估算BLUE/BLUP值再进行关联分析；

（2） 质量性状

单基因控制，无法用具体数值衡量，可转换成0、1表示。建议采集的两种性状的样本数应当尽量相近。

（3） 分级性状

表型分布类似质量分布，但实际受多基因控制，如抗性性状。

① 可将分级数量化，如果根据抗性强弱赋值0、1、2、3、4等，等级顺序必须固定且有意义，不同的顺序会获得不同的结果，随意的顺序会导致结果难以解释；

② 采用一定的测定方法将其数量化，如使用存活率、死亡率、病斑大小等表示抗性的强弱。

注：群体中某样本在某位点缺失表型，用NA表示，建议单个群体中，缺失表型的个数尽量不要多于群体的1/3。

Q3. 采集统计完的表型就可以直接用于分析吗？

在GWAS分析前，一般需要先对表型数据进行正态性检验，极端异常值去除，并对多年多点表型值进行处理。

GWAS分析属于线性模型，要求数据尽量符合正态分布。首先绘制频率分布图，观测数据分布情况，如果得到的表型数据不符合正态分布，将对表型数据进行正态化处理。

这些工作一般需要生物信息学技术完成处理，但是作为第一手表型数据的采集人，是否可以自己判断表型符合分析需求呢？小编探索发现，利用Office excel的数据分析功能，只需3步就可以实现表型数据的正态性分布初步检验！这里为大家详细介绍一下：

1）以Trait 1为例，将统计好的表型数据记录为一列，先利用【数据】-【筛选】功能查看表型数据的分布范围。也可以利用排序功能完成这步，我这组数据是0.15-18.68；

2）分段。根据自己的数据分布范围确定横坐标范围，一般可以按20个小段。我这里可以按0-20整体范围，1以间隔分为20个区间。在旁边的单元格填写分段范围。

3）绘制频率分布直方图。在【数据】工具栏点击【数据分析】，跳出一个数据分析的选择框，选择“直方图”。输入区域选择表型数据所在的所有单元格，接收区域选择设置好的分段范围所在的单元格，输出区域是生成结果的位置，可以随意选择。勾选上图表输出，点击确定，这组表型数据的频率分布直方图就绘制好啦，是不是so easy？

数据工具栏中找不到【数据分析】的小伙伴不要着急，在【文件】-【选项】中找到【加载项】，点击【转到】，跳转出加载宏的小窗口，勾选分析工具库之后就会出现啦。

当然，有条件的话还是建议大家选择更加专业的分析工具，例如SPSS，SPSS的优点是可以针对数据生成正态分布的模拟曲线，更加清晰地看出表型数据分布与理想的正态分布之前的差异。操作也非常简单，同样以Trait 1的表型数据为例，为大家介绍一下。

1）把表型数据粘贴到SPSS表格中，注意粘贴前需要把数据保存为数值格式；

2）【分析】-【描述统计】-【频率】，跳转出频率小窗口，点击中间的按钮，把数据加到变量栏。

3）分别点击【统计量】和【图表】，在【统计量】窗口中，勾选分布的偏度和峰度，在【图表】窗口中，勾选直方图以及在直方图上显示正态曲线，点击确定即在新窗口中生成结果。

对于不符合正态分布的数据，可以通过去除极端异常值、多年多点表型值处理使数据正态化，处理完再重新检验确定。

Q4.GWAS分析模型/软件怎么选择？

GWAS分析中的使用的模型主要为一般线性模型（GLM，General Linear Model）和混合线性模型（MLM，Mixed Linear Model）两种：

相比较而言MLM在GLM的基础上考虑了亲缘关系的影响，可以消除群体内个体间亲缘关系造成的假阳性过高的情况，比之更加全面。

随着测序技术的发展，测序分析成本降低，用于GWAS分析的样本量和标记量越来越多，用原始的MLM模型分析所需的时间非常长。人们尝试了不同的优化方法，发展出以混合线性模型为核心的一系列分析模型，如EMMAX、FaST-LMM、GEMMA、FarmCPU等，这些模型被打包开发成相应的分析软件。

其中EMMAX和GEMMA由于计算速度快、检测能力也较好，被广泛应用于动植物的GWAS分析。

Q5.GWAS分析结果如何解读?

GWAS分析完成后得到的主要结果是，所有被分析的变异位点与目标性状的关联性P 值、效应大小及其方向的数据。这些数据通常使用Manhattan图和QQ图展示。

Manhattan图和QQ图具体该怎么看怎么理解呢？一起往下看。

曼哈顿图是GWAS分析得到的最直观的结果，其实质是点图，比较多的点堆叠在一起的效果。横坐标是染色体，纵坐标是-log10(P)，变异位点与目标性状之间相关性越高，原始P值越小，-log10转换后的值越大，在曼哈顿图中就越高。曼哈顿图中超过阈值线的点表示相应变异位点与目标性状显著相关。在动植物中，影响表型性状的位点一般处于连锁不平衡状态，也就是说通过GWAS分析筛选到的显著位点一般成簇存在，如下方曼哈顿图中的红框所示，我们认为这些位点具有更高的可靠性。当然，筛选到的非成簇的显著SNP位点，如果在多年、多点数据中重复出现，也被认为是更可靠的。针对这些显著相关的位点，可进一步关注其注释信息，结合其位置、基因型、所在基因、所在基因功能等分析其生物学意义。

QQ图（quantile-quantile plot），也叫做分位图，是判断GWAS分析结果假阳性、假阴性的重要指标。

Q6. GWAS分析结果中的假阳性、假阴性怎么理解呢？

一般认为关联P 值达到显著性，就说明相应的变异位点对表型有显著影响。QQ图的实质是比较观测P 值（Y轴）和期望P 值（X轴）的一致性。观测P 值可以理解为与随机遗传漂变有关。随机漂变，是随机出现在染色体上突变，它符合均匀分布。对于GWAS分析中的变异位点来说，其观测P 值应该为随机漂变与基因效应的加和。

我们来结合实际分析中几种QQ图理解一下：

1）QQ图中，点都分布在从原点发出的45°线上，观测P 值与期望P 值几乎完全一致，说明几乎没有基因效应发挥作用，GWAS分析也就没有找到与目标性状显著关联的位点。

2）P 值低的位点，即与性状无显著关联的位点，观测值与期望值一致。P 值高的位点，观测值超过期望值，说明这些位点的基因效应明显使其超过了随机效应，进而说明这些位点是与性状显著相关的，也就是我们的目标位点。

3）如果大部分位点P 值观测值低于期望值，也就是可能存在假阴性的情况。说明可能存在模型设置不合适，P 值被过度校正，例如群体分化不明显时。

4）如果很多位点P 值观测值都高于期望值，即很多位点都与某一个性状显著相关，这显然是不合理的，说明假阳性过高，模型选择不合理，或者数据不符合GWAS分析要求。

Q7. 单体型（haplotype）构建及解读

由于动植物中控制性状的位点具有连锁不平衡性，可以通过单倍型分析，确定与目标表型相关的控制区域范围。

单倍型：指在同一染色体上一定区域内若干个决定同一性状的且紧密连锁的位点，可以是两个基因座或整条染色体。单倍型确定后，绘制单倍型图可用于关联分析。

单倍型块(Haplotype Block)：指染色体上存在着的连续的、稳定的、几乎不被重组所打断的单倍型区域。以单倍型块为基础，进行目标性状与基因的关联分析是目前最经济的连锁不平衡(LD) 分析方法。

如何理解LD block和haplotype呢？

以上图为例，最上方蓝色粗线表示一条染色体上的一段区域，上面一道道白色的竖线表示用于分析的变异位点，以SNP为例。计算两两SNP之间的LD值，以0-100（百分数）表示，数值越高颜色越深，形成下方的热图。

从上图中每两个SNP之间的LD值都可以清楚找到，例如SNP8和SNP11的数值为98。因此，SNP8和SNP11的连锁相关性为0.98，彼此处于相当高的连锁不平衡状态。这些位于同一条染色体上处于连锁不平衡状态的SNP所在的一段连续的区域可以划分为一个LD block，也称为haplotype block。Block的划分可以用软件计算出来，划分时不同的软件阈值设置会有不同的block划分结果。例如上图即为以0.8为阈值划分的block，所以block3没有包含SNP6和SNP19等。

通过单倍型和Block分析，可以确定显著位点的关联范围，发现更多可能与目标表型相关的关键变异位点或基因。

Q8. 根据GWAS结果如何展开后续分析？

根据单倍型分析找到可靠性高的SNP，结合注释信息寻找关键SNP，对关键SNP细化分析，例如表型解释率PVE、基因型表型相关性（box图）、优势等位基因型判定等。

我们可以从发表的文章中看到类似的思路，例如2020年发表在Plant Biotechnology Journal上的“Identification of loci controlling adaptation in Chinese soyabean landraces via a combination of conventional and bioclimatic GWAS”。利用两个地点的开花时间数据做GWAS分析找到了两个强关联位点Chr12：5470311和Chr12：5914898。对这两个位点进行block分析，找到了开花时间的候选基因，而且这个候选基因与拟南芥中已知的一个开花基因同源。进一步分析Chr12:5470311和Chr12:5914898 SNP位点的特定基因型与3个亚群开花时间的一致性（图C，D）、Fst、等位基因频率分析鉴定（图E），研究证明GWAS分析筛选到的两个显著SNP位点均参与了开花时间的调控。