GWAS找到顯著信號位點(diǎn)后叭爱，需要解釋顯著信號位點(diǎn)如何影響表型挟憔。
常見的一個(gè)解釋方法是共定位分析。
主流的共定位分析包括：
- 1）GWAS和eQTL共定位峦失；
- 2）GWAS和sQTL共定位悴灵；
- 3）GWAS和meQTL共定位；
- 4）GWAS和pQTL共定位腿短；
其中屏箍，GWAS和eQTL共定位應(yīng)用最為廣泛。
具體來說橘忱，當(dāng)檢測到GWAS信號和eQTL共定位時(shí)赴魁，我們會認(rèn)為GWAS信號上的位點(diǎn)可能通過改變基因表達(dá)的生物學(xué)過程影響表型。
共定位分析有四種設(shè)想：
第一種設(shè)想 H0: 表型1（GWAS）和 表型2 （以eQTL為例）與某個(gè)基因組區(qū)域的所有SNP位點(diǎn)無顯著相關(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ū)動凝颇；
第四種設(shè)想 H4: 表型1（GWAS）和 表型2 （以eQTL為例）與某個(gè)基因組區(qū)域的SNP位點(diǎn)顯著相關(guān)潘拱，且由同一個(gè)因果變異位點(diǎn)驅(qū)動；
基于以上四種設(shè)想拧略，我們希望第四種設(shè)想 H4 在統(tǒng)計(jì)學(xué)上概率更高芦岂，這樣就能解釋顯著信號位點(diǎn)如何影響表型；
所以共定位分析垫蛆，本質(zhì)上是在檢驗(yàn)第四種的后驗(yàn)概率禽最；
講完共定位分析的相關(guān)概念，接下來我們以“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 下載測試數(shù)據(jù)
測試數(shù)據(jù)請?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編號rs_id懦趋，P值pval_nominal等：
注意事項(xiàng)：如果表型是二分類變量（case和control），輸入文件二選一：
1）rs編號
rs_id疹味、P值pval_nominal仅叫、SNP的效應(yīng)值beta帜篇、效應(yīng)值方差varbeta；2）rs編號
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編號rs_id，基因gene_id遂跟，次等位基因頻率maf、P值pval_nominal等：
注意事項(xiàng)：如果表型是連續(xù)型變量婴渡，輸入文件三選一：
1）rs編號
rs_id幻锁、P值pval_nominal、表型的標(biāo)準(zhǔn)差sdY边臼；2）rs編號
rs_id哄尔、P值pval_nominal、效應(yīng)值beta,效應(yīng)值方差varbeta, 樣本量N,次等位基因頻率MAF柠并；3）rs編號
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)
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)會放松要求到0.75。
這里假定后驗(yàn)概率大于0.95為共定位位點(diǎn)：
library(dplyr)
need_result=result$results %>% filter(SNP.PP.H4 > 0.95)
結(jié)果如下：
從圖上可以看出缰雇，SNP.4811位點(diǎn)的后驗(yàn)概率為1入偷。怎么找到這個(gè)位點(diǎn)呢？可以通過對應(yīng)的P值（1.81e-42）找到這個(gè)位點(diǎn)的rs號械哟；
4 結(jié)果解讀
對于輸出結(jié)果疏之，我們只需要關(guān)注最后一列信息SNP.PP.H4，對應(yīng)推文前面提到的第四種設(shè)想 H4暇咆。
SNP.PP.H4表示的是GWAS顯著信號和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ì)算沃呢，所以請不要把全基因組的信息都丟進(jìn)去（偷懶該打）年栓，這么做一個(gè)是沒意義，另外一個(gè)特別耗時(shí)薄霜，給計(jì)算機(jī)增加工作負(fù)擔(dān)某抓；
2）雖然我們沒必要把基因組的全部信息丟進(jìn)去纸兔，但也不意味著只放一個(gè)變異位點(diǎn)信息就行。
3）因此否副，正確的做法是汉矿，先提取GWAS中顯著變異位點(diǎn)上下游區(qū)域（這個(gè)區(qū)域多大自己定，沒有金標(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)纽乱，請考慮其他工具（允許多個(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
致謝橙子牛奶糖（陳文燕），請用參考模版：We thank the blogger (orange_milk_sugar, Wenyan Chen) for XXX
感謝小可愛們多年來的陪伴昆箕， 我與你們一起成長~
