GWAS课程学习笔记

参考课程： 基迪奥GWAS课程：https://www.omicshare.com/class/Home/Index/singlev?id=3

自然群体利用了进化过程中的染色体重组，容易进行基因定位。

1、GWAS分析常用的软件

（1）TASSEL

• 植物类项目应用较多，可以矫正群体结构和系谱关系（Trait Analysis by aSSociation, Evolution and Linkage; Bradbury et al, 2007, Bioinformatics 23:2633-2635）
  （2）EMMA
• 动物类项目应用较多，可矫正系谱关系（Kang et al, 2008, Genetics 178:1709-1723）
  （3）Plink
• 使用较为简单（Purcell et al, American Journal of Human Genetics, 2007, 81）

2、表型的处理：线性表型性状

• 正态性判断：R语言的shapiro.test(x)检验
• 如果是仅个别样本异常，建议剔除。如极端值、离开均值大于4倍SD的
• 若整体偏离散（如基因表达量值），建议取log2后，重新检验正态性。

3、材料的选择

主要从两方面考虑，一是其LD衰减和重组情况如何、二是群体结构如何。

（1）群体的选择

• 野生品种、地方品种、培育品种
  不同群体关联分析的效果不同

  
  
• 标记的效应越弱，要检测到这个标记所需的样本数目就更大，因此要先考虑研究的性状是偏质量的、还是偏主效基因的，or前人报道的没有主效基因的
• 如果是前人报道的无主效基因的，就要考虑增加样本，或确实定位不到位点

（2）基因型是否完全覆盖

• GWAS分析的基础就是基因与标记之间的LD是否连锁，不同群体的LD衰减距离不同，可以用hyploview进行计算。
  
• 通常当两个位点间R2>0.8时，认为两位点处于完全连锁不平衡，但这种连锁状态会随区域增加而不断降低。
  如何根据LD衰减距离判断做GWAS所需的标记个数？？
• 如果群体的LD衰减距离是100k，那么分析时就要保证每100k至少要有一个marker，那么1M就需要10个、1G就是10w个、3G就需要30w个
• 核心种质的LD衰减非常快，因此要增加标记密度

（3）基因型判断群体结构的影响（随机背景标记）

群体结构（Q矩阵）和个体检潜在的系谱关系（K矩阵），可能会导致假阳性（如下图）：

• 群体结构和性状分布恰好一致，会使人误以为只要是量群体特有的基因就都是与性状关联的，即将区分群体的背景标记认为是与性状相关的；
• 解决办法（2种）：
  ① 将群体结构作为协变量，引入到方程式里，将群体间的影响校正掉，剩下的效应可能是标记的效应；

• 计算群体结构（Q矩阵）：用structure或PCA分析的结果，作为群体结构的协变量，将其引入模型 ----- 具体操作见杨晓红老师GWAS操作教程课件
• 计算个体遗传关系（系谱关系，K矩阵）：用SPAGeDi软件

② 将两个群体分开，分别单独做GWAS，来敲除遗传背景对群体结构的影响

4、GWAS分析的多阶段设计

（1）什么是多阶段设计？

• 在人类疾病的GWAS研究中，常用两阶段法分析，比较严谨。第一阶段一般用覆盖全基因组的位点，第二阶段则聚焦在少量的候选位点的测序数据进行GWAS分析。

• 单阶段：一个群体进行关联分析 → 完成不严谨，一般为动植物类的研究
• 两阶段：
  （1）阶段1：找候选关联位点
  小样本（几百）全基因组关联分析，得到候选位点；
  （2）阶段2：候选位点的验证
  已有群体大样本（成千上万）或新的独立群体，只对候选位点的关联分析。
  （2）多阶段设计的优点
• 有验证的步骤：可靠；
• 降低成本：第二阶段的检测位点数较少；
• 解决潜在的多重检验校正的位点
  高密度芯片or全基因组重测序，SNP数量可达1M，多重检验过于严格。
  如：1M SNP，Bonferroni校正的adjusted p value阈值 = 0.05/110-6=510-8（太严格）
• 可以采用的方法：第一阶段放松过滤阈值，在第二阶段进行验证。由于第二阶段位点数较少，多重检验校正不会如此严格。

5、关联分析所需的模型

（1）模型原理

• 固定效应1：环境效应，如不同年份、不同地点数据
• 固定效应2：位点效应
• 固定效应3：群体结构效应，群体分层导致的，需要纠正，样本所属的亚群分类信息用Q矩阵表示
• 随机效应：潜在的系谱关系，K矩阵
• 随机误差
  关联分析时并不是说所有位点都要考虑，要结合自己的情况，选择合适的

（2）模型的选择

① 动物

• PCA分析初步判断；
• 一般而言，动物类样本在系谱清晰，且没有明显群体结构的情况下，可直接剔除离群样本；剔除离群样本后，再将剩下的个体做PCA分析，如果不再存在群体分层，即可用一般线性模型做关联分析；
• 若存在群体分层，再考虑使用Q矩阵进行矫正。

② 植物

• PCA分析初步判断；
• 植物（尤其作物）因品系间杂交更普遍（如玉米），故群体结构和不同品系间的系谱关系更普遍；分析时，同时使用一般线性模型和不同的混合线性模型，然后比较结果的好坏。

（3）如何判断模型是否合适？——qq图

① 正常的qq图：前贴后起

• GWAS分析后，p-value的-log10从低到高排序，看其与期望p-value之间的差别）

• 假如标记与性状完全不相关，则标记的p-value应该是正态分布，因此会一直沿着直线走，并且实际情况下，绝大部分标记确实是跟性状不相关。到了后期，标记的显著性增高，可能开始与性状之间存在相关，因此其观测到的p值会显著高于期望p值。

② 异常情况：过度矫正

• 过度矫正的可能原因：
  a. 群体结构或kinship矫正过于严格，导致观测值＜期望值；
  b. 期望p-value的随机分布是基于位点之间互相独立的假设，高通量测序or高密度芯片会导致很多相邻位点间存在连锁or相关关系，这样的话观测到的p值就不是完全随机的，若位点间实际存在

（4）关联分析的模型选择

• 做任何性状的关联分析时，都需要用至少2个模型进行模拟，判断最佳模型

  
  

（5）不同分析方法的最适范围：

6、示例：GWAS分析的一般步骤

step 1：通过进化树和PCA分析，看群体分层情况

step 2：不同模型的比较 —— 找出最佳模型

step 3：分群体和全群体分析 —— 当存在明显的群体分层时

Step 4：对定位到的位点的解读：优先解读可解读的，再去挖掘其他的

step 5： 结合RNA-seq或群体遗传学等其他方法来验证这个位点附近的基因可能是与性状相关的