-
Notifications
You must be signed in to change notification settings - Fork 0
/
coloc_analysis.R
36 lines (29 loc) · 919 Bytes
/
coloc_analysis.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
library("coloc")
library(dplyr)
# Import Phenotype 1 (GWAS) data:
gwas <- read.table(file="/~path/GWAS.txt", header=T);
# Import Phenotype 2 (eQTL) data:
eqtl <- read.table(file="/~path/eQTL.txt", header=T);
# Merging GWAS and eQTL data
merged_data <- merge(gwas, eqtl, by="rs_id")
# Colocalization analysis
result <- coloc.abf(
dataset1 = list(
snp = merged_data$rs_id,
pvalues = merged_data$pval_nominal.x,
beta = log(merged_data$OR),
type = "cc",
s = 0.33,
N = 50000
),
dataset2 = list(
snp = merged_data$rs_id,
pvalues = merged_data$pval_nominal.y,
type = "quant",
N = 10000
),
MAF = merged_data$MAF
)
# Screen for colocalized sites
library(dplyr)
need_result=result$results %>% filter(SNP.PP.H4 > 0.8)