話不多說(shuō)，直接上代碼

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ì)性檢測(cè)結(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

大家來(lái)點(diǎn)贊啊，點(diǎn)贊超過100，分享day-4代碼，多個(gè)工具變量F值的計(jì)算
另外，還有一些reviewer關(guān)心統(tǒng)計(jì)功效，power計(jì)算的R代碼on the way