話不多說师痕，直接上代碼
library(TwoSampleMR)
library(data.table)
#step 1. read exposure data
exposure_dat <- read_exposure_data('met-a-746_MR_format.txt', clump = F, sep = "\t",phenotype_col = "Phenotype",snp = "SNP", beta_col = "beta", se_col = "se", effect_allele_col = "effect_allele", other_allele_col = "other_allele", pval_col = "pval", samplesize_col = "samplesize", eaf_col = "eaf")
#step 2. exposure data clump, parameters are from the article（default parameter is clump_r2 = 0.001,clump_kb = 10000）
exposure_dat_clump <- clump_data(exposure_dat, clump_r2 = 0.1, pop = "EUR",clump_kb = 500)
#step 3. read outcome data
outcome_data <- fread('Alzheimer_GWAS_summary_example.txt')
outcome_dat <- format_data( dat=outcome_data, type = "outcome", snps = exposure_dat_clump$SNP, header = TRUE, phenotype_col = "Phenotype", snp_col ="SNP",beta_col ="beta",se_col ="se",effect_allele_col ="effect_allele",other_allele_col ="other_allele",pval_col ="pval",samplesize_col = "samplesize", eaf_col = "eaf")
#step 4. harmonise
dat <- harmonise_data(exposure_dat_clump, outcome_dat)
#step 5. caculate F-stat for each SNP
dat$EAF2 <- (1 - dat$eaf.exposure)
dat$MAF <- pmin(dat$eaf.exposure, dat$EAF2)
PVEfx <- function(BETA, MAF, SE, N){
pve <- (2*(BETA^2)*MAF*(1 - MAF))/((2*(BETA^2)*MAF*(1 - MAF)) + ((SE^2)*2*N*MAF*(1 - MAF)))
return(pve)
}
dat$PVE <- mapply(PVEfx, dat$beta.exposure, dat$MAF, dat$se.exposure, N = dat$samplesize.exposure)
dat$FSTAT <- ((dat$samplesize.exposure - 1 - 1)/1)*(dat$PVE/(1 - dat$PVE))
#step 6. heterogeneity test, heterogeneity (Inverse variance weighted) Q-pval = 0.3858222 > 0.05, then choose Inverse variance weighted (fixed effects) method
mr_results_het <- mr_heterogeneity(dat)
#step 7. MR analysis using Inverse variance weighted (fixed effects) method
res <- mr(dat, method_list = c("mr_ivw_fe"))
#step 8. Add OR and CI information
res <- generate_odds_ratios(res)
部分結(jié)果展示
step 5.
所有SNP的F-stat值都大于10试幽，因此都納入分析中
step 5 結(jié)果展示
step 6. 異質(zhì)性檢測結(jié)果
p = 0.3858222 > 0.05赴涵，因此選擇Inverse variance weighted (fixed effects) method匿值；如果 p < 0.05灿巧，選擇Inverse variance weighted (multiplicative random effects) method
step 6 結(jié)果展示
step 7. Inverse variance weighted (fixed effects) 結(jié)果
step 7 和 step 8結(jié)果展示
和原文章的對(duì)比蹦肴，原文章附件中Table S6包含了代謝物和 Alzheimer的所有MR結(jié)果僚碎，最后一行就是test的代謝物，Q estimate 均為5.25阴幌， P value for Q estimate 均為0.39勺阐，均選擇 Fixed-effect model 進(jìn)行分析，OR 值均為0.69矛双，95% CI均為0.57-0.84渊抽，P value 均為1.98×10-4。
復(fù)現(xiàn)完成议忽，喜大普奔......然后循環(huán)分析所有代謝物
原文章Table S6
大家來點(diǎn)贊啊懒闷，點(diǎn)贊超過100，分享day-4代碼，多個(gè)工具變量F值的計(jì)算
另外愤估，還有一些reviewer關(guān)心統(tǒng)計(jì)功效帮辟，power計(jì)算的R代碼on the way
