因为自己复现的结果始终与文献不一致,所以尝试联系了作者并找到了作者提供的代码,这次我们尝试用官方代码重跑一遍~
作者提供的数据和代码在这里:
https:// doi.org/10.17605/OSF.IO/MUERZ.
今天的内容建议参考【孟德尔随机化】文献复现(三),比较一下代码的异同以便进一步思考~
我们先来看table 1部分:
chr<-4
window <- 0 #using a more stringent window for the anlysis.
gene_start<-120415550#119494395#
gene_end<-120550146#119628991#
chrpos <- paste0(chr, ":", gene_start - window, "-", gene_end + window)
获取PDE5的eqtl数据:
eqtl<-extract_instruments('eqtl-a-ENSG00000138735' ,clump = F)
eqtl<-eqtl$SNP[eqtl$pval.exposure< 5*10^-8 & eqtl$pos.exposure>gene_start-window & eqtl$pos.exposure<gene_end+window & eqtl$chr.exposure==chr]
misssense<-c("rs3733526","rs139979143", "rs114886951","rs149385790") #missesnse varailbes
获取收缩压数据:
exp_dsbp<- extract_outcome_data(c(eqtl,misssense),
c("ieu-b-39"),
proxies = F,
rsq = 0.8,
align_alleles = 1,
palindromes = 1,
maf_threshold = 0.3,
access_token = ieugwasr::check_access_token(),
splitsize = 10000,
proxy_splitsize = 500)
在线clump:
exp_dsbp<-convert_outcome_to_exposure(exp_dsbp)
exp_dsbp<-clump_data(exp_dsbp,clump_r2 = 0.35, clump_kb = 10000)
可以发现这里的snp其实和文献中的并不完全重合,反而和之前俺用本地clump复现的结果一致,除了beta值的正负恰好相反(……)
向作者请教以后,作者给的答复如下
所以还是可能与clump的方法有关,本地和在线竟然差这么大!
但为了尽量重复原文结果,我们还是用文献中的5个snp接着分析——
exp_dsbp<- extract_outcome_data(c("rs80223330", "rs12646525", "rs17355550", "rs10050092", "rs66887589"), c("ieu-b-39"), proxies = F, rsq = 0.8, align_alleles = 1, palindromes = 1, maf_threshold = 0.3, access_token = ieugwasr::check_access_token(),splitsize = 10000, proxy_splitsize = 500)
exp<-convert_outcome_to_exposure(exp_dsbp)
然后来看看困扰小编很久的S table1是如何实现的:
ED<- extract_outcome_data(exp$SNP, c("ebi-a-GCST006956","finn-b-I9_HYPTENSPUL","finn-b-ERECTILE_DYSFUNCTION"), proxies = F, rsq = 0.8, align_alleles = 1, palindromes = 1, maf_threshold = 0.3, access_token = ieugwasr::check_access_token(),splitsize = 10000, proxy_splitsize = 500)
接着还是平平无奇的harmonise
dat <- harmonise_data(exp, ED, action = 2)
重头戏在这里:
for (i in dat$SNP[dat$outcome.deprecated=="Erectile dysfunction || || "&dat$id.exposure=="ieu-b-39"&dat$id.outcome=="finn-b-ERECTILE_DYSFUNCTION"]){
mED<-meta::metagen(TE=dat$beta.outcome[dat$SNP==i&dat$outcome.deprecated=="Erectile dysfunction || || "&dat$id.exposure=="ieu-b-39"],seTE=dat$se.outcome[dat$SNP==i&dat$outcome.deprecated=="Erectile dysfunction || || "&dat$id.exposure=="ieu-b-39"])
dat$beta.outcome[dat$SNP==i&dat$outcome.deprecated=="Erectile dysfunction || || "&dat$id.exposure=="ieu-b-39"&dat$id.outcome=="ebi-a-GCST006956"]<-mED$TE.fixed
dat$beta.outcome[dat$SNP==i&dat$outcome.deprecated=="Erectile dysfunction || || "&dat$id.exposure=="ieu-b-39"&dat$id.outcome=="finn-b-ERECTILE_DYSFUNCTION"]<-NA
dat$se.outcome[dat$SNP==i&dat$outcome.deprecated=="Erectile dysfunction || || "&dat$id.exposure=="ieu-b-39"&dat$id.outcome=="ebi-a-GCST006956"]<-mED$seTE.fixed
dat$se.outcome[dat$SNP==i&dat$outcome.deprecated=="Erectile dysfunction || || "&dat$id.exposure=="ieu-b-39"&dat$id.outcome=="finn-b-ERECTILE_DYSFUNCTION"]<-NA
}
dat<-dat[!is.na(dat$beta.outcome),]
for (i in c("ebi-a-GCST006956","finn-b-I9_HYPTENSPUL" )){#"finn-b-ERECTILE_DYSFUNCTION"
for (j in c("ieu-b-39")){#,"ieu-b-38"
print(j)
print(i)
IVW<-IVWcorrel(
betaYG=dat$beta.outcome[dat$id.exposure==j&dat$id.outcome==i],
sebetaYG=dat$se.outcome[dat$id.exposure==j&dat$id.outcome==i],
betaXG=dat$beta.exposure[dat$id.exposure==j&dat$id.outcome==i],
sebetaXG=dat$se.exposure[dat$id.exposure==j&dat$id.outcome==i],
rho=matrix[dat$SNP[dat$id.exposure==j&dat$id.outcome==i],dat$SNP[dat$id.exposure==j&dat$id.outcome==i]]
)
print(IVW)
}}
[1] "ieu-b-39"
[1] "ebi-a-GCST006956"
beta_IVWcorrel se_IVWcorrel.random p_IVWcorrel.random n_snp F_stat
[1,] 0.1276879 0.05520846 0.02073184 5 25.58036
[1] "ieu-b-39"
[1] "finn-b-I9_HYPTENSPUL"
beta_IVWcorrel se_IVWcorrel.random p_IVWcorrel.random n_snp F_stat
[1,] 3.293556 0.5730375 9.055132e-09 5 25.58036
努力破禁锢,下周继续~