# GWAS

这两天开始折腾 GWAS，又翻了翻 GWAS 的原理，好多都不懂，记录一下自己的心路历程。

这次参考了 Zhou (2022) 番茄图泛基因组的 GWAS 流程 YaoZhou89/TGG: tomato graph pangenome，使用了 MLM模型 与 LOCO (Leava one choromosome out) 的方法。流程主要是 LDAK计算亲缘关系矩阵 + gcta进行GWAS分析 ，一上来直接给我干懵了，这里记录一下我这两天折腾的过程，对这次的分析流程的部分原理进行介绍。

# MLM

Yu 等在玉米中使用 MLM 的模型，向 GWAS 模型中纳入了 亲缘关系矩阵 与 群体结构矩阵 ，Yang 等在 gcta 中实现了这一方法。根据 Yang 等文章中的描述。

$$y=X\beta+S\alpha+Qv+Zu+e$$

其中$y$ (n,1) 为表型向量；

$X\beta$ 为除 SNP 与群体结构以外的固定效应，$\beta$ (p,1) 为除了 SNP 或群体结构效应之外的其他固定效应的向量；

$S\alpha$ 为 SNP 固定效应，$\alpha$ 是一个 SNP 固定效应的向量；

$Qv$ 为群体结构固定效应，$v$ 为群体结构向量 (eg. PCA)；

$Zu$ 代表亲缘关系效应，$u$ 为随机多基因背景效应的向量；$Var(u)=2KV_g$，K 为 (n,n) 的亲缘关系矩阵，定义个体之间遗传协方差；$V_g$ 为遗传方差；

$X,S,Q,Z$ 分别为$\beta,\alpha,v,u$ 的关联矩阵；

$e$ 是一个残差效应的向量；$Var(e)=RV_R$，R (n,n) 非对角线元素为 0，对角线元素是获得每个表型数据点的观测数的倒数；$V_R$ 为残差；

# BLUE

上述模型中，对$\beta$ 的估计使用 BLUE 的方法，即利用 GLS (最小二乘法) 找到$\beta$ 使得$Var(y)$ 最小。

详见 Livro-Applications-of-Linear-Models-in-Animal-Breeding. 第三章，拉格朗日算子求条件极值，https://math.stackexchange.com/questions/1104376/how-to-set-up-lagrangian-for-matrix-constraints

或者可使用 EM算法 迭代求得模型最优参数$\beta$ (文中的$Y = X\beta + Zu + \epsilon$ 的$X\beta$ 可展开为$X\beta+S\alpha+Qv$，可见 Livro-Applications-of-Linear-Models-in-Animal-Breeding , 第二章)

# P

https://www.mv.helsinki.fi/home/mjxpirin/GWAS_course/material/GWAS2.html

有了估计模型，则可以对每个变异位点进行差异检验，上述文中使用 t 测对 SNP 的显著性进行检验。

# QQ plot

QQ plot 反应了期望 P 值与观测 P 值之间的关系，如下图所示（QQ plot 图 —— 评价你的统计模型是否合理 - 基迪奥生物）

图 a 中，Pvalue 观察值和期望值相同，说明分析模型是合理的。但所有的 P value 观测值都没有明显超过期望值，说明分析结果没有找到（与性状）显著关联的位点，可能原因包括：性状由微效多基因控制，效应太弱；群体大小不够等，这里先不展开详述。

图 b 是我们最期望看到的结果类型。在散点图的左下角是显著性低的位点，即确定与性状不关联的位点，这些位点的 P value 观测值应该与期望值一致。而图中这些点的确位于对角线上，说明分析模型是合理的。而在图形的右上角则是显著性较高的位点，是潜在与性状相关的候选位点。这些点位于对角线的上方，即位点的 P value 观测值超过了期望值，说明这些位点的效应超过了随机效应，进而说明这些位点是与性状显著相关的。小结了一下：这个图形的左下角说明了模型的合理性，右上角则说明找了关联位点，所以这是最理想的结果。（备注：在有显著关联位点的情况下，结合曼哈顿图进行展示，会更加醒目）

图 c 是大部分点位于对角线的下方，则说明大部分位点的 P value 观察值小于期望值。主要原因包括两种情况：（1）模型不合理，P value 被过度校正，导致 P value 显著性过低；（2）群体中大量 SNP 位点间存在连锁不平衡，有效位点数（相互间不存在连锁不平衡的位点）明显低于实际位点数，所以 P value 的期望值被低估了（即期望值的 - log10（P value）被高估了），也会出现这种情况。（3）模型中过度拟合，没有亲缘关系的群体在拟合过程中添加了亲缘关系矩阵。（4）个体过少。

图 d，则相反：大部分点位于对角线的上方，则说明大部分位点的 P value 观察值超过期望值。按照统计学的逻辑推导，就是大部分位点与某个性状显著相关。这显然是不符合生物学逻辑的，那么这只有一种可能：分析模型不合理，数据的假阳性过大，P value 观测值的显著性被高估了。

# Reference

A unified mixed-model method for association mapping that accounts for multiple levels of relatedness | Nature Genetics

Advantages and pitfalls in the application of mixed-model association methods | Nature Genetics

https://www.nature.com/articles/ng.546#Sec6