anx_gz <- "ANX_2026_daner.gz"
if (!file.exists(anx_gz)) {
download.file("https://ndownloader.figshare.com/files/62089324", destfile = anx_gz, mode = "wb")
}
dme <- fread("../GRCH37_conversion/gwas_GRCh37.csv",
select = c("CHROM_GRCh37","GENPOS_GRCh37","ALLELE0","ALLELE1","BETA","SE","P","N"))
setnames(dme, c("CHROM_GRCh37","GENPOS_GRCh37","ALLELE0","ALLELE1","BETA","SE", "P"),
c("chrom","pos","dme_a0","dme_a1","dme_beta","dme_se", 'dme_p'))
anx <- fread(anx_gz, select = c("CHR","SNP","BP","A1","A2","OR","SE","P","Nca","Nco"))
setnames(anx, c("CHR","BP","A1","A2","OR","SE", "P"),
c("chrom","pos","anx_a1","anx_a0","anx_or","anx_se", "anx_p"))
anx[, anx_beta := log(anx_or)]
anx <- anx[, .(chrom, pos, anx_a0, anx_a1, anx_beta, anx_se, anx_p, Nca, Nco)]
# The txt file gives this info and says they filtered for minimum 70% of cases
anx_total_n <- 122083 + 729602
anx <- anx[(Nca + Nco) >= 0.7 * anx_total_n]
# The top hit for DME on chromosome 6
chrom_i = 6
pos_i = 26239404
window = 5e5
lower_limit = pos_i - window
upper_limit = pos_i + window
d1 <- dme[chrom == chrom_i & pos > lower_limit & pos < upper_limit]
d2 <- anx[chrom == chrom_i & pos > lower_limit & pos < upper_limit]
# Only keep SNPs that are in both GWAS
m <- merge(d1, d2, by = c("chrom", "pos"))
# Leave out palindromic allelles (AT, TA, CG, GC)
is_palindromic <- (m$dme_a0 == "A" & m$dme_a1 == "T") | (m$dme_a0 == "T" & m$dme_a1 == "A") |
(m$dme_a0 == "C" & m$dme_a1 == "G") | (m$dme_a0 == "G" & m$dme_a1 == "C")
m <- m[!is_palindromic]
# Align both GWAS
same_direction <- m$dme_a1 == m$anx_a1 & m$dme_a0 == m$anx_a0
opposite_direction <- m$dme_a1 == m$anx_a0 & m$dme_a0 == m$anx_a1
m <- m[same_direction | opposite_direction]
m[dme_a1 == anx_a0, anx_beta := -anx_beta]
# Give our SNPs a name for coloc
snp <- paste(m$chrom, m$pos, sep = ":")
d1 <- list(
snp = snp,
beta = m$dme_beta,
varbeta = m$dme_se^2,
type = "cc"
)
d2 <- list(
snp = snp,
beta = m$anx_beta,
varbeta = m$anx_se^2,
type = "cc"
)
res <- coloc.abf(d1, d2)
res$summary