共定位分析是一種用貝葉斯檢驗來分析兩種可能相關的表型是否在給定區域內共享共同的遺傳因果變異的方法，本文中，我們使用 coloc 包來進行共定位分析。
數據的准備
MAF的計算
我們的共定位的目的是為了找出caQTL、eQTL與meQTL的共定位，主角是我們的caQTL數據。coloc需要SNP的次等位基因頻率（MAF）以及樣本量信息，這兩部分信息需要從基因型文件中提取，因此，我們的第一步就是從基因型數據中提取caQTL中的caSNP的次等位基因頻率和樣本量信息，由於基因型矩陣分別使用0、1、2對AA、Aa和aa進行了編碼，因此對於較小頻率的突變等位基因a，其等位基因頻率為：
\[AF_{a} = \frac{0N_{AA} + 1N_{Aa} + 2N_{aa}}{2N_{AA} + 2N_{Aa} + 2N_{aa}}，\ 其中N_{AA}、N_{Aa}、N_{aa}分別代表三種基因型的樣本數 \]
上面的這個公式實質上等價於基因型矩陣每一行的均值的一半，因此我們可以使用 mean(row) / 2 來計算每一個SNP的alt等位基因頻率，然后再判斷此頻率是否在5%到95%之間。
library("readr")
library("tidyverse")
# 導入原始的caQTL文件
caqtl_raw = read_tsv("caqtl_raw.tsv") %>% distinct(snps, .keep_all = TRUE)
# 導入eQTL文件
eqtl = read_tsv("eqtl.tsv", col_types = cols())
eqtl_unique = eqtl %>% arrange(`p-value`) %>% distinct(snp, .keep_all = TRUE)
# 計算等位基因頻率的函數
calculate_AF = function(row) {
x = as.numeric(row)
AF = round(mean(x[!is.na(x)]) / 2, 3)
return(min(AF, 1-AF))
}
# 導入基因型文件並計算MAF
genotype = read_tsv("BRCA_genotype")
genotype$maf = apply(genotype, 1, calculate_AF)
# 篩選caQTL和eQTL中出現的SNP
genotype_caqtl = genotype %>% filter(ID %in% caqtl_raw$snps)
genotype_caqtl = genotype_caqtl[match(caqtl_raw$snps, genotype_caqtl$ID),]
genotype_caqtl$maf = apply(genotype_caqtl, 1, calculate_AF)
genotype_eqtl = genotype %>%
mutate(ID = gsub(":.*", "", ID)) %>%
filter(ID %in% intersect(ID, eqtl_unique$snp)) %>%
distinct(ID, .keep_all = TRUE)
eqtl_unique = eqtl_unique %>% filter(snp %in% intersect(snp, genotype_eqtl$ID))
genotype_eqtl = genotype_eqtl[match(eqtl_unique$snp, genotype_eqtl$ID),]
genotype_eqtl$maf = apply(genotype_eqtl, 1, calculate_AF)
處理兩種QTL的數據為coloc需要的格式
coloc不接受數據表格式，因此需要將數據包格式的數據轉換為列表的格式，對於定量數據，需要以列表的形式指定額外的屬性：
type: "quant"：數據的類型MAF：SNP的次等位基因頻率N：數據的樣本量
數據的處理方式如下：
# 處理caQTL數據為coloc需要的格式
caqtl = read_tsv("caQTL.tsv", col_types = cols())
caqtl_clean = caqtl %>%
select(caSNP_ID,caSNP_chr,caSNP_pos,std_error,beta,pvalue) %>%
mutate(varbeta = std_error^2) %>%
select(-std_error) %>%
rename("snp" = caSNP_ID, "p" = pvalue) %>%
unite(col = "position", caSNP_chr,caSNP_pos,sep = ":") %>%
arrange(p) %>%
distinct(snp, .keep_all = TRUE) %>%
as.list()
# 處理eQTL數據為coloc需要的格式
eqtl_clean = eqtl_unique %>%
rename("tstat" = "t-stat", "p" = "p-value") %>%
mutate(varbeta = (beta / tstat) ^ 2) %>%
select(snp, chr, position, beta, varbeta, p) %>%
unite(col = "position", chr, position, sep = ":") %>%
arrange(p) %>%
distinct(snp, .keep_all = TRUE) %>%
as.list()
# 轉換為列表格式並加入上一部分的結果
caqtl_clean$type = "quant"
caqtl_clean$MAF = genotype_caqtl$maf
caqtl_clean$N = 74
eqtl_clean$type = "quant"
eqtl_clean$MAF = genotype_eqtl$maf
eqtl_clean$N = 1093
# 檢查數據表是否有問題
check_dataset(caqtl_clean) # NULL
check_dataset(eqtl_clean) # NULL
coloc分析以及結果解讀
如果數據的檢查函數返回 NULL 說明數據沒有問題，接着就可以進行共定位分析了，使用的函數是 coloc.abf
coloc_results = coloc.abf(eqtl_clean, caqtl_clean)
根據共定位算法，結果中的后驗概率類型包含以下5種，分別是：
H0：兩種性狀在區域內都沒有遺傳學關聯H1：只有性狀1在區域內具有遺傳學關聯H2：只有性狀2在區域內具有遺傳學關聯H3：兩種性狀在區域內具有不同的遺傳學關聯H4：兩種性狀在區域內具有同一的遺傳學關聯（共定位證據）
我們選取后驗概率比較高的點進行查看
coloc_results$results %>% filter(SNP.PP.H4 > 0.75) %>% select(snp)
發現只有 rs2031947 一個結果，其余的SNP都不顯著，我們在caQTL和eQTL數據庫中進行尋找，結果如下：
可以看到caQTL和eQTL的趨勢比較一致，而與meQTL相反，而且這個SNP對應的基因 GSTM1 曾被預測與多種癌症相關。