GWAS找到顯著信號(hào)位點(diǎn)后，需要解釋顯著信號(hào)位點(diǎn)如何影響表型。
常見的一個(gè)解釋方法是共定位分析。

主流的共定位分析包括：

• 1）GWAS和eQTL共定位；
• 2）GWAS和sQTL共定位；
• 3）GWAS和meQTL共定位；
• 4）GWAS和pQTL共定位；

其中，GWAS和eQTL共定位應(yīng)用最為廣泛。

具體來(lái)說(shuō)，當(dāng)檢測(cè)到GWAS信號(hào)和eQTL共定位時(shí)，我們會(huì)認(rèn)為GWAS信號(hào)上的位點(diǎn)可能通過(guò)改變基因表達(dá)的生物學(xué)過(guò)程影響表型。

共定位分析有四種設(shè)想：

  第一種設(shè)想 H0: 表型1（GWAS）和 表型2 （以eQTL為例）與某個(gè)基因組區(qū)域的所有SNP位點(diǎn)無(wú)顯著相關(guān)；  

  第二種設(shè)想 H1/H2: 表型1（GWAS）或表型2（以eQTL為例）與某個(gè)基因組區(qū)域的SNP位點(diǎn)顯著相關(guān)；

  第三種設(shè)想 H3: 表型1（GWAS）和 表型2 （以eQTL為例）與某個(gè)基因組區(qū)域的SNP位點(diǎn)顯著相關(guān)，但由不同的因果變異位點(diǎn)驅(qū)動(dòng)；

   第四種設(shè)想 H4: 表型1（GWAS）和 表型2 （以eQTL為例）與某個(gè)基因組區(qū)域的SNP位點(diǎn)顯著相關(guān)，且由同一個(gè)因果變異位點(diǎn)驅(qū)動(dòng)；

基于以上四種設(shè)想，我們希望第四種設(shè)想 H4 在統(tǒng)計(jì)學(xué)上概率更高，這樣就能解釋顯著信號(hào)位點(diǎn)如何影響表型；

所以共定位分析，本質(zhì)上是在檢驗(yàn)第四種的后驗(yàn)概率；

講完共定位分析的相關(guān)概念，接下來(lái)我們以“GWAS和eQTL共定位”為例講一下如何使用coloc進(jìn)行共定位分析。

1 下載、安裝coloc

if(!require("remotes"))
  install.packages("remotes")
  install.packages("dplyr")
library(remotes)
install_github("chr1swallace/coloc",build_vignettes=TRUE)
library("coloc")
library(dplyr)

2 下載測(cè)試數(shù)據(jù)

測(cè)試數(shù)據(jù)請(qǐng)?jiān)凇癰io生物信息”回復(fù)"coloc"；
或者用自己的數(shù)據(jù)分析（如果有的話）；

3 分析

3.1 導(dǎo)入表型1（GWAS）數(shù)據(jù)：

gwas <- read.table(file="E:/path_to_GWAS/GWAS.txt", header=T);

GWAS數(shù)據(jù)包括：rs編號(hào)rs_id，P值pval_nominal等：

image

注意事項(xiàng)：如果表型是二分類變量（case和control），輸入文件二選一：

1）rs編號(hào)rs_id、P值pval_nominal、SNP的效應(yīng)值beta、效應(yīng)值方差varbeta；

2）rs編號(hào)rs_id、P值pval_nominal、case在所有樣本中的比例s

3.2 導(dǎo)入表型2（eQTL）數(shù)據(jù)：

eqtl <- read.table(file="E:/path_to_eqtl/eQTL.txt", header=T);

eQTL數(shù)據(jù)包括：rs編號(hào)rs_id，基因gene_id，次等位基因頻率maf、P值pval_nominal等：

image

注意事項(xiàng)：如果表型是連續(xù)型變量，輸入文件三選一：

1）rs編號(hào)rs_id、P值pval_nominal、表型的標(biāo)準(zhǔn)差sdY；

2）rs編號(hào)rs_id、P值pval_nominal、效應(yīng)值beta,效應(yīng)值方差 varbeta, 樣本量N,次等位基因頻率 MAF；

3）rs編號(hào)rs_id、P值pval_nominal、次等位基因頻率 MAF；

3.3 合并GWAS和eQTL數(shù)據(jù)：

input <- merge(eqtl, gwas, by="rs_id", all=FALSE, suffixes=c("_eqtl","_gwas"))
head(input)

image

3.4 共定位分析

result <- coloc.abf(dataset1=list(pvalues=input$pval_nominal_gwas, type="cc", s=0.33, N=50000), dataset2=list(pvalues=input$pval_nominal_eqtl, type="quant", N=10000), MAF=input$maf)

dataset1的type="cc"指的是GWAS的表型是二分類（case和control）；

dataset2的type="quant"指的是eQTL的表型（基因表達(dá)量）是連續(xù)型

N指樣本量；

3.5 篩選共定位的位點(diǎn)

通常情況下，很多文獻(xiàn)認(rèn)為PPA > 0.95的位點(diǎn)是共定位位點(diǎn)，也有一些文獻(xiàn)會(huì)放松要求到0.75。

這里假定后驗(yàn)概率大于0.95為共定位位點(diǎn)：

library(dplyr)
need_result=result$results %>% filter(SNP.PP.H4 > 0.95)

結(jié)果如下：

image

從圖上可以看出，SNP.4811位點(diǎn)的后驗(yàn)概率為1。怎么找到這個(gè)位點(diǎn)呢？可以通過(guò)對(duì)應(yīng)的P值（1.81e-42）找到這個(gè)位點(diǎn)的rs號(hào)；

4 結(jié)果解讀

對(duì)于輸出結(jié)果，我們只需要關(guān)注最后一列信息SNP.PP.H4，對(duì)應(yīng)推文前面提到的第四種設(shè)想 H4。

SNP.PP.H4表示的是GWAS顯著信號(hào)和eQTL位點(diǎn)為同一個(gè)位點(diǎn)的后驗(yàn)概率，范圍在0-1之間，0表示概率為0%，1表示概率為100%。后驗(yàn)概率越高越好。

5 注意事項(xiàng)

1）由于共定位分析是基于某個(gè)基因組區(qū)域進(jìn)行計(jì)算，所以請(qǐng)不要把全基因組的信息都丟進(jìn)去（偷懶該打），這么做一個(gè)是沒(méi)意義，另外一個(gè)特別耗時(shí)，給計(jì)算機(jī)增加工作負(fù)擔(dān)；

2）雖然我們沒(méi)必要把基因組的全部信息丟進(jìn)去，但也不意味著只放一個(gè)變異位點(diǎn)信息就行。

3）因此，正確的做法是，先提取GWAS中顯著變異位點(diǎn)上下游區(qū)域（這個(gè)區(qū)域多大自己定，沒(méi)有金標(biāo)準(zhǔn)）的GWAS summary數(shù)據(jù)，理想情況下，提取后顯著變異位點(diǎn)所在基因組區(qū)域的SNP數(shù)量在1,000-10,000之間。

4）該方法的設(shè)想是只有一個(gè)causal 位點(diǎn)，所以如果表型1（GWAS）和 表型2 （以eQTL為例）在某個(gè)區(qū)域有多個(gè)顯著位點(diǎn)的話，用該方法是定位不出結(jié)果的，這是該方法的局限，所以如果某個(gè)基因組區(qū)域存在多個(gè)顯著位點(diǎn)，請(qǐng)考慮其他工具（允許多個(gè)causal 位點(diǎn)共定位的工具）

特別鳴謝：
https://chr1swallace.github.io/coloc/FAQ.html
https://hanruizhang.github.io/GWAS-eQTL-Colocalization/
https://rpubs.com/Charleen_D_Adams/743547

致謝橙子牛奶糖（陳文燕），請(qǐng)用參考模版：We thank the blogger (orange_milk_sugar, Wenyan Chen) for XXX

感謝小可愛們多年來(lái)的陪伴， 我與你們一起成長(zhǎng)~