全基因组关联分析（GWAS）的计算原理

Posted 2023-04-27

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了全基因组关联分析（GWAS）的计算原理相关的知识，希望对你有一定的参考价值。

参考技术A

关于全基因组关联分析（GWAS）原理的资料，网上有很多。

这也是我写了这么多GWAS的软件教程，却从来没有写过GWAS计算原理的原因。

恰巧之前微博上某位小可爱提问能否写一下GWAS的计算原理。我一顺口就答应了。

后面一直很懒，不愿意动笔，但想着既然答应了，不写说不过去。

我写这段话的意思是，如果你有任何关于GWAS分析问题或者疑问，希望我能写一下的，可以跟我说。

如果我认为有价值，写出来对大家有帮助的话，会写的。

首先，我们来一个知识点的回顾：最小二乘法。

看下图，熟不熟悉！

这可是我们中学时解了很多遍的算术题。

公式可以写为： y = ax + b

如图所示，假定有一个SNP,叫 rs123: T>C

我们定义C为风险位点， 以加性模型为例 ，一个C=1,T=0

那么CC=2,CT=1,TT=0

根据上面的公式：

SNP对应的值x分别为：2,2,1,2,1,1,0,2
对应的表型y分别为10,7,6,8,5,4,2,6

回顾我们前面提到的公式：y = ax + b

现在我们有：

10= 2a+b

7= 2a+b

6= 1a+b

8= 2a+b

5= 1a+b

4= 1a+b

2= 0+b

6= 2a+b

转化一下，就是：

2a+b - 10 = 0

2a+b - 7 = 0

1a+b - 6 = 0

2a+b - 8 =0

1a+b - 5 = 0

1a+b - 4 = 0

0+b -2 = 0

2a+b -6 = 0

我们的任务就是，找到合适的a，b使得

（2a+b - 10)^2 + (2a+b - 7)^2 + (1a+b - 6)^2 + (2a+b - 8)^2 + (1a+b - 5)^2 + (1a+b - 4)^2 + (0+b -2)^2 + (2a+b -6)^2 最小。

a，b的求值涉及到最小二乘法的推导，感兴趣的看这篇文章： https://zhuanlan.zhihu.com/p/53556591

用公式表示就是：

可以自己手算一下，也可以借助R语言：

x=c(2,2,1,2,1,1,0,2)

y=c(10,7,6,8,5,4,2,6)

Ex=mean(x);Ex

Ey=mean(y);Ey

Sx=sd(x);Sx

Sy=sd(y);Sy

corn=cor(y,x) ; corn

b=corn*Sy/Sx ; b

a=Ey-b*Ex ; a

最后拟合的结果是：a的最优化为 2.8387， b的最优化为 2.0968 ，公式 y = 2.8387* x + 2.0968

R语言的lm函数也可以计算a和b，完全不需要我们自己一个个手动推导。lm函数结果的Intercept即为b值，x所在行对应的Estimate值即为a值

回归到我们的全基因组关联分析，这里的a即为beta（OR）值

所以搞明白全基因组关联分析的值是怎么来的了吗，这个就是它的计算原理

上面列出来的公式只是简单的计算基因型与表型之间的相关性。

实际上，我们在计算的时候，会加入其他的变量，比如性别，年龄，品系等。

这些因素也是影响表型的变量。

因此，当考虑其他变量存在时，计算公式会稍微改变一下：y = ax + zβ + b

计算原理同上。

全基因组关联分析（GWAS）:为何我的QQ图那么飘

参考技术A

前段时间有位小可爱问我，为什么她的QQ图特别飘，如果你不理解怎样算飘，请看下图：

理想的QQ图应该是这样的：

我当时的第一反应是：1）群体分层造成的；2）表型分布有问题。因此让她检查一下数据的群体分层情况，如果没有问题就看一下表型分布。

这段时间有空了，我觉得有必要梳理一下这个飘逸的QQ图，到底是怎么回事儿以及如何确定这么飘逸的QQ图有没有问题。

产生飘逸的qq图的原因有很多，比如我们喜闻乐见的： 基因多效性（polygenicity） 。也有可能是 混淆偏倚 ，比如群体分层，假如样本中混合了欧洲、非洲、亚洲等各个地方的群体，本身各个群体的SNP频率差异就很大，如果不加以群体分层控制，关联分析的时候就会产生很多偏离预期值的SNP位点。

很久以前，出现飘逸的QQ图的话，我们可以搭配膨胀系数（膨胀系数的计算）一起看，膨胀系数如果接近1（比如1.01，1.2这种不算接近1），那么也还算过得去。

但，膨胀系数接近1这种是比较理想的情况。实际情况是，很多人的 QQ图不仅飘逸，膨胀系数还老高 。

这就尴尬了，连膨胀系数都无法确定这个QQ图飘的正不正常了。

所以呢，接下来还有什么方法确定我们的基因组数据有没有问题呢？

接下来我要介绍一款神人工具 LDSC (LD SCore) ,全称是LD Score regression

这款工具就是帮我们实现 如何区分飘逸的QQ图正不正常 。

具体来说，就是通过LDSC工具计算基因组数据的LD回归截距，如果是基因多效性，那么截距会接近1（比如1.004），如果是群体分层等混淆因素引起的，那么LD回归截距就会远离1（比如1.30）。

除了看截距数值，我们还可以通过画LD Score的图来确定数据是否有问题。

如果画出来的LDscore图是下面这种形式，说明GWAS结果是没有问题的，QQ图飘逸就让它飘逸吧。

如果画出来的LDscore图是下面这种形式，说明QQ图是有问题的。

总之，看截距。

截距很接近1，就不用管QQ图好不好看了。