全基因組關聯分析（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，性別年齡等因素外等其餘可能影響表型的因素；

計算原理同上。