【孟德尔随机化】文献复现（五）

生信菜鸟团

发布于 2024-02-26 10:24:18

4160

发布于 2024-02-26 10:24:18

文章被收录于专栏：生信菜鸟团生信菜鸟团

因为自己复现的结果始终与文献不一致，所以尝试联系了作者并找到了作者提供的代码，这次我们尝试用官方代码重跑一遍~

作者提供的数据和代码在这里：

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