话不多说,今天咱们的目标是这张图里的male部分,也是文章最重要的结论来源~
根据工具变量的条件过滤以后,得到的工具变量——
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)
对应的代码如下:
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)
}
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)
##读取结局数据
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
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
##读取结局数据
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博士不厌其烦地答疑解惑和小伙伴冯雨瑶的交流讨论。
好的文章确实是值得反复琢磨的,在反复咀嚼的过程中不断提高,对我来说也是弥足珍贵的体验。
下一期我会找找看有没有孟德尔和单细胞的跨界合作文章,为大家提供一些不一样的东西……