啊，全基因组关联分析（GWAS）的计算原理，了解一下？
前言
关于全基因组关联分析（GWAS）原理的资料，网上有很多。
这也是我写了这么多GWAS的软件教程，却从来没有写过GWAS计算原理的原因。
恰巧之前微博上某位小可爱提问能否写一下GWAS的计算原理。我一顺口就答应了。
后面一直很懒，不愿意动笔，但想着既然答应了，不写说不过去。
我写这段话的意思是，如果你有任何关于GWAS分析问题或者疑问，希望我能写一下的，可以跟我说。
如果我认为有价值，写出来对大家有帮助的话，会写的。
GWAS所涉及的公式：最小二乘法
首先，我们来一个知识点的回顾：最小二乘法。
看下图，熟不熟悉！
这可是我们中学时解了很多遍的算术题。
图片来源：http://kitsprout.blogspot.com/2015/11/algorithm-least-squares.html
公式可以写为：y = ax + b
y：我们研究的表型
x：基因型数据，这里指每一个SNP
a：SNP的系数
b：残差，可以是环境变量，或者除了SNP之外的影响表型的因素
来个例子给我们讲讲呗，公式怎么套进去
如图所示，假定有一个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
用公式表示就是：
b = cor(y,x)*Sd(y)/Sd(x)
a = (10+7+6+8+5+4+2+6)/8 - ((2+2+1+2+1+1+0+2)/8)*b
cor(y,x)表示x和y的相关系数
Sd(y),Sd(x)分别表示y和x的标准差
可以自己手算一下，也可以借助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
y:我们研究的表型
x:基因型数据，这里指每一个SNP
a：SNP的系数
z:年龄，性别等因素
β：年龄，性别等因素的系数
b：残差，除了我们纳入的SNP，性别年龄等因素外等其他可能影响表型的因素；
计算原理同上。