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

生信菜鸟团

发布于 2024-02-28 12:50:36

3680

发布于 2024-02-28 12:50:36

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

话不多说，今天咱们的目标是这张图里的male部分，也是文章最重要的结论来源~

Instrument selection

根据工具变量的条件过滤以后，得到的工具变量——

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_dsbp<-convert_outcome_to_exposure(exp_dsbp)

Mendelian randomization

对应的代码如下：

matrix<-ld_matrix(exp_dsbp$SNP)
a<-unlist(stringr::str_split(colnames(matrix),"_"))
colnames(matrix)<-a [seq(from=1,to=length(a), by =3)]
a<-unlist(stringr::str_split(rownames(matrix),"_"))
rownames(matrix)<-a [seq(from=1,to=length(a), by =3)]

mr<-function(exp,out){
  #outcome<-input[input$SNP %in% exp$SNP,]
  dat<-harmonise_data(exp,out)
  #ld<-matrix[dat$SNP,dat$SNP]
  IVW<-IVWcorrel(betaYG=dat$beta.outcome,  sebetaYG=dat$se.outcome,  betaXG=dat$beta.exposure,  sebetaXG=dat$se.exposure,  rho=matrix[dat$SNP,dat$SNP])
  return(IVW)
}

主要的结局

Number of children

kids_men<-extract_outcome_data(exp_dsbp$SNP ,'ukb-b-2227', proxies = F)  

mr(exp_dsbp,kids_men)

beta_IVWcorrel se_IVWcorrel.random p_IVWcorrel.random n_snp   F_stat
[1,]    -0.03797464          0.00809483       2.715781e-06     5 25.58036

# 置信区间：95%CI = (ES-1.96SE, ES+1.96SE)lCI <- 0.03797464*5.5*1.32269-1.96*0.00809483*5.5*1.32269;lCI # 0.1608368
hCI <- 0.03797464*5.5*1.32269+1.96*0.00809483*5.5*1.32269;hCI # 0.3916786

大家注意到这里为什么要乘5.5再乘1.32269了吗？【也是问了作者大大才注意到这个细节的，( ╯□╰ )】

Data-Field 2405 (ox.ac.uk)

Number of sexual partners

##读取结局数据
dsbp_dat <- read_outcome_data(filename = "../no_sex_men_imputed.txt.gz",
                             snps =exp_dsbp$SNP,
                             sep = "\t", ##很关键！否则容易报错
                             phenotype_col = "Outcome", ##要有这一列
                             snp_col = "SNP",
                             beta_col = "BETA",
                             se_col = "SE",
                             eaf_col = "A1FREQ",
                             effect_allele_col = "ALLELE1",
                             other_allele_col = "ALLELE0",
                             pval_col = "P_LINREG",
                             chr_col = "CHR",
                             pos_col = "BP")
mr(exp_dsbp,dsbp_dat)


# 置信区间：95%CI = (ES-1.96SE, ES+1.96SE)
lCI <- 0.2082813*5.5-1.96*0.9169092*5.5;lCI #-8.738734
hCI <- 0.2082813*5.5+1.96*0.9169092*5.5;hCI #11.02983

Probability of never having had sexual intercourse

dsbp_dat <- read_outcome_data(filename = "../vergin_men_imputed.txt.gz",
                              snps =exp_dsbp$SNP,
                              sep = "\t", ##很关键！否则容易报错
                              phenotype_col = "Outcome", ##要有这一列
                              snp_col = "SNP",
                              beta_col = "BETA",
                              se_col = "SE",
                              eaf_col = "A1FREQ",
                              effect_allele_col = "ALLELE1",
                              other_allele_col = "ALLELE0",
                              pval_col = "P_LINREG",
                              chr_col = "CHR",
                              pos_col = "BP")
mr(exp_dsbp,dsbp_dat)

# 置信区间：95%CI = (ES-1.96SE, ES+1.96SE)
lCI <- -3.969715e-06*5.5-1.96*0.001066891*5.5;lCI #-0.01152292
hCI <- -3.969715e-06*5.5+1.96*0.001066891*5.5;hCI # 0.01147925

Subjective wellbeing

##读取结局数据
dsbp_dat <- read_outcome_data(filename = "../wellbeing_men_imputed.txt.gz",
                              snps =exp_dsbp$SNP,
                              sep = "\t", ##很关键！否则容易报错
                              phenotype_col = "Outcome", ##要有这一列
                              snp_col = "SNP",
                              beta_col = "BETA",
                              se_col = "SE",
                              eaf_col = "A1FREQ",
                              effect_allele_col = "ALLELE1",
                              other_allele_col = "ALLELE0",
                              pval_col = "P_LINREG",
                              chr_col = "CHR",
                              pos_col = "BP")
mr(exp_dsbp,dsbp_dat)

同上，有个问题是这里主观幸福感以标准偏差单位衡量。但SD信息查询不到……

最后统一校正p值

p.adjust(c(2.71578080673037e-06,0.820302536933148,0.997031217792966,0.692947161539899),method="fdr") #dbp

总体上来说，目前文章结果基本是能够复现的，而且通过学习这样的文章，对我们尝试不一样的MR分析也会有帮助，希望大家都能有所收获~

碎碎念——

这篇文章【A drug target for erectile dysfunction to help improve fertility, sexual activity, and wellbeing: mendelian randomisation study】的复现基本就告一段落啦~感谢大家一路以来的支持，特别要感谢这篇文章的作者Benjamin博士不厌其烦地答疑解惑和小伙伴冯雨瑶的交流讨论。

好的文章确实是值得反复琢磨的，在反复咀嚼的过程中不断提高，对我来说也是弥足珍贵的体验。

下一期我会找找看有没有孟德尔和单细胞的跨界合作文章，为大家提供一些不一样的东西……

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2024-02-28，如有侵权请联系 cloudcommunity@tencent.com 删除

变量