MR-Analysis 2-26-2026

Mendelian Randomization

Import Exposure

library(tidyverse)
library(TwoSampleMR)
library(LDlinkR)      # LD and proxy snps
library(ieugwasr)
library(fixtuRes)

exposure_path_A2 = "gwas_resources/A2EUR_GRCH38.tsv.gz"
exposure_c1 <- read_tsv(exposure_path_A2,comment = "##", 
                        col_select = c(1,2,3,4,5,6,7,8,9,14,15))
colnames(exposure_c1)<- c("CHR","POS","effect_allele.exposure","other_allele.exposure","SNP","N","beta.exposure","se.exposure","pval","eaf.exposure","rsid")
exposure_c1$exposure <- rep(c("A2_severe_respiratory_COVID19"),times=nrow(exposure_c1));
clump_id<- "a2grch"
#calculate z value
exposure_c1$z.exposure <- exposure_c1$beta.exposure/exposure_c1$se.exposure;
exposure_c1$id.exposure <- rep(c(clump_id),times=nrow(exposure_c1));


exposure_dat1 <- filter(exposure_c1, pval < 5e-8)
exposure_clump1 <-ld_clump(
    #dplyr::tibble(rsid=exposure_c1$rsid, pval=exposure_c1$pval, id=exposure_c1$trait_id),
    dat = exposure_dat1,
    plink_bin = genetics.binaRies::get_plink_binary(),
    bfile = "1kgv3/EUR"
    
)

write_csv(exposure_clump1, "data/exposure_A2GRCH38.csv", append = FALSE)

exposure_path_B2 = "gwas_resources/B2EUR_GRCH38.tsv.gz"
exposure_c2 <- read_tsv(exposure_path_B2,comment = "##", 
                        col_select = c(1,2,3,4,5,6,7,8,9,14,15))
colnames(exposure_c2)<- c("CHR","POS","effect_allele.exposure","other_allele.exposure","SNP","N","beta.exposure","se.exposure","pval","eaf.exposure","rsid")
exposure_c2$exposure <- rep(c("B2_hospitalized_COVID19"),times=nrow(exposure_c2));
clump_id<- "b2grch"
#calculate z value
exposure_c2$z.exposure <- exposure_c2$beta.exposure/exposure_c2$se.exposure;
exposure_c2$id.exposure <- rep(c(clump_id),times=nrow(exposure_c2));


exposure_dat2 <- filter(exposure_c2, pval < 5e-8)
exposure_clump2 <-ld_clump(
    #dplyr::tibble(rsid=exposure_c1$rsid, pval=exposure_c1$pval, id=exposure_c1$trait_id),
    dat = exposure_dat2,
    plink_bin = genetics.binaRies::get_plink_binary(),
    bfile = "1kgv3/EUR"
    
)

write_csv(exposure_clump2, "data/exposure_B2GRCH38.csv", append = FALSE)

exposure_path_C2 = "gwas_resources/C2EUR_GRCH38.tsv.gz"
exposure_c3 <- read_tsv(exposure_path_C2,comment = "##", 
                        col_select = c(1,2,3,4,5,6,7,8,9,14,15))
colnames(exposure_c3)<- c("CHR","POS","effect_allele.exposure","other_allele.exposure","SNP","N","beta.exposure","se.exposure","pval","eaf.exposure","rsid")
exposure_c3$exposure <- rep(c("C2_COVID19"),times=nrow(exposure_c3));
clump_id<- "c2grch"
#calculate z value
exposure_c3$z.exposure <- exposure_c3$beta.exposure/exposure_c3$se.exposure;
exposure_c3$id.exposure <- rep(c(clump_id),times=nrow(exposure_c3));


exposure_dat3 <- filter(exposure_c3, pval < 5e-8)
exposure_clump3 <-ld_clump(
    #dplyr::tibble(rsid=exposure_c1$rsid, pval=exposure_c1$pval, id=exposure_c1$trait_id),
    dat = exposure_dat3,
    plink_bin = genetics.binaRies::get_plink_binary(),
    bfile = "1kgv3/EUR"
    
)

write_csv(exposure_clump3, "data/exposure_C2GRCH38.csv", append = FALSE)

Import Outcome

outcome1 <- read_tsv('gwas_resources/AKIJiangL2021h.tsv.gz', col_select = c(2,3,4,5,6,7,11,12,18, 21,22))
colnames(outcome1)<- c("rsid","CHR","POS","other_allele.outcome","effect_allele.outcome","beta.outcome","eaf.outcome","code","N","se.outcome","pval.outcome")
outcome1$outcome <- rep(c("PheCode 585.1"),times=nrow(outcome1));
clump_id <- "jiangl2021";
outcome1$id.outcome <- rep(c(clump_id),times=nrow(outcome1));

exposure_snps <- read_csv("data/exposure_A2GRCH38.csv") 
#filter outcome on SNPs in the exposure
outcome_clumped <- filter(outcome1, rsid %in% exposure_snps$rsid)
#harmonise data
mr_dat <- harmonise_data(dplyr::tibble(SNP=exposure_snps$rsid, beta.exposure=exposure_snps$beta.exposure, se.exposure=exposure_snps$se.exposure, effect_allele.exposure = exposure_snps$effect_allele.exposure, other_allele.exposure = exposure_snps$other_allele.exposure,eaf.exposure = exposure_snps$eaf.exposure, id.exposure = exposure_snps$id.exposure, exposure = exposure_snps$exposure), dplyr::tibble(SNP=outcome_clumped$rsid, beta.outcome=outcome_clumped$beta.outcome, se.outcome=outcome_clumped$se.outcome, effect_allele.outcome = outcome_clumped$effect_allele.outcome, other_allele.outcome = outcome_clumped$other_allele.outcome,eaf.outcome = outcome_clumped$eaf.outcome , outcome = outcome_clumped$outcome, id.outcome = outcome_clumped$id.outcome), action = 2)
write_csv(mr_dat, "harmonized/harmonized_A2_JiangL")

outcome2 <- read_tsv('gwas_resources/DouvilleNJModel1h.tsv.gz', col_select = c(1,2,3,4,5,6,7,8,9,13))
colnames(outcome2)<- c("CHR","POS","effect_allele.outcome","other_allele.outcome","beta.outcome","se.outcome","eaf.outcome","pval.outcome","rsid","N")
outcome2$outcome <- rep(c("Sepsis AKI"),times=nrow(outcome2));
clump_id <- "DouvilleNJ2025";
outcome2$id.outcome <- rep(c(clump_id),times=nrow(outcome2));


exposure_snps <- read_csv("data/exposure_A2GRCH38.csv") 
#filter outcome on SNPs in the exposure
outcome_clumped <- filter(outcome2, rsid %in% exposure_snps$rsid)
#harmonise data
mr_dat <- harmonise_data(dplyr::tibble(SNP=exposure_snps$rsid, beta.exposure=exposure_snps$beta.exposure, se.exposure=exposure_snps$se.exposure, effect_allele.exposure = exposure_snps$effect_allele.exposure, other_allele.exposure = exposure_snps$other_allele.exposure,eaf.exposure = exposure_snps$eaf.exposure, id.exposure = exposure_snps$id.exposure, exposure = exposure_snps$exposure), dplyr::tibble(SNP=outcome_clumped$rsid, beta.outcome=outcome_clumped$beta.outcome, se.outcome=outcome_clumped$se.outcome, effect_allele.outcome = outcome_clumped$effect_allele.outcome, other_allele.outcome = outcome_clumped$other_allele.outcome,eaf.outcome = outcome_clumped$eaf.outcome , outcome = outcome_clumped$outcome, id.outcome = outcome_clumped$id.outcome), action = 2)
write_csv(mr_dat, "harmonized/harmonized_A2_Douville")

outcome3 <- read_tsv('gwas_resources/AKIVermaAh.tsv.gz', col_select = c(1,2,3,4,5,6,7,8,9,13))
colnames(outcome3)<- c("CHR","POS","effect_allele.outcome","other_allele.outcome","oddsratio.outcome","se.outcome","eaf.outcome","pval.outcome","rsid","N")
outcome3$outcome <- rep(c("PheCode 585.1 (Verma)"),times=nrow(outcome3));
clump_id <- "vermaA2024";
outcome3$id.outcome <- rep(c(clump_id),times=nrow(outcome3));
outcome3$beta.outcome <- log(outcome3$oddsratio.outcome)


exposure_snps <- read_csv("data/exposure_A2GRCH38.csv") 
#filter outcome on SNPs in the exposure
outcome_clumped <- filter(outcome3, rsid %in% exposure_snps$rsid)
#harmonise data
mr_dat <- harmonise_data(dplyr::tibble(SNP=exposure_snps$rsid, beta.exposure=exposure_snps$beta.exposure, se.exposure=exposure_snps$se.exposure, effect_allele.exposure = exposure_snps$effect_allele.exposure, other_allele.exposure = exposure_snps$other_allele.exposure,eaf.exposure = exposure_snps$eaf.exposure, id.exposure = exposure_snps$id.exposure, exposure = exposure_snps$exposure), dplyr::tibble(SNP=outcome_clumped$rsid, beta.outcome=outcome_clumped$beta.outcome, se.outcome=outcome_clumped$se.outcome, effect_allele.outcome = outcome_clumped$effect_allele.outcome, other_allele.outcome = outcome_clumped$other_allele.outcome,eaf.outcome = outcome_clumped$eaf.outcome , outcome = outcome_clumped$outcome, id.outcome = outcome_clumped$id.outcome), action = 2)
write_csv(mr_dat, "harmonized/harmonized_A2_VermaA")

MR Analysis

library(tidyverse)    # Data wrangling 
library(knitr)
library(TwoSampleMR)  # MR 
library(gt)
library(LDlinkR)      # LD and proxy snps
library(RadialMR)     # Radial MR sensetivity analysis 
library(phenoscanner)

mr_dat <- read_csv("harmonized/harmonized_A2_JiangL")

Rows: 28 Columns: 22── Column specification ──────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (9): SNP, effect_allele.exposure, other_allele.exposure, effect_allele.outcome, other_allele.outcome, ...
dbl (8): beta.exposure, beta.outcome, eaf.exposure, eaf.outcome, se.outcome, se.exposure, action, SNP_index
lgl (5): remove, palindromic, ambiguous, mr_keep, samplesize.outcome
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

mr_pleiotrophy_result <- mr_pleiotropy_test(mr_dat)
View(mr_pleiotrophy_result)
mr_res <- mr(mr_dat, method_list = c(
  "mr_two_sample_ml",  "mr_egger_regression", 
  "mr_simple_median", "mr_weighted_median", 
  "mr_ivw_fe", "mr_ivw_mre",
  "mr_simple_mode", "mr_weighted_mode"
  ))

Analysing 'a2grch' on 'jiangl2021'

res_single <- mr_singlesnp(mr_dat, all_method = c("mr_ivw_fe", "mr_egger_regression", "mr_weighted_median", "mr_weighted_mode")) %>% as_tibble()
View(res_single)
odd <- generate_odds_ratios(mr_res)
print("Finished Analyzing")

[1] "Finished Analyzing"

write.csv(odd,"results/A2/jiangla2_oddsratios.csv")
write.csv(res_single,"results/A2/jiangla2_res_single.csv")
scatter_p <- mr_scatter_plot(mr_res, mr_dat) 
scatter_out_p <- scatter_p[[1]] + theme_bw() + 
  guides(color=guide_legend(ncol =1)) + 
  theme(
    text = element_text(size = 8), 
  )

scatter_out_p
ggsave("results/A2/scatterplotA2JiangL.png", plot = scatter_out_p, units = 'in', height = 4, width = 9)

forrest_p <- mr_forest_plot(res_single)
forrest_p[[1]]
ggsave("results/A2/forrestplotA2JiangL.png", plot =forrest_p[[1]] + theme(text = element_text(size = 8)), units = 'in', height = 9, width = 9)

res_loo <- mr_leaveoneout(mr_dat, method = mr_ivw_fe) %>% as_tibble()
loo_p <- mr_leaveoneout_plot(res_loo)
loo_p[[1]]
ggsave("results/A2/looA2JiangL.png",plot = loo_p[[1]], units = 'in', height = 9, width = 9)

funnel_p <- mr_funnel_plot(res_single)
funnel_out_p <- funnel_p[[1]] + theme_bw() + 
  guides(color=guide_legend(ncol =1)) + 
  theme(
    text = element_text(size = 8), 
  )

funnel_out_p
ggsave("results/A2/funnelplotA2JiangL.png", plot = funnel_out_p, units = 'in', height = 4, width = 9)

mr_dat <- read_csv("harmonized/harmonized_A2_Douville")

Rows: 28 Columns: 22── Column specification ──────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (9): SNP, effect_allele.exposure, other_allele.exposure, effect_allele.outcome, other_allele.outcome, ...
dbl (8): beta.exposure, beta.outcome, eaf.exposure, eaf.outcome, se.outcome, se.exposure, action, SNP_index
lgl (5): remove, palindromic, ambiguous, mr_keep, samplesize.outcome
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

mr_res <- mr(mr_dat, method_list = c(
  "mr_two_sample_ml",  "mr_egger_regression", 
  "mr_simple_median", "mr_weighted_median", 
  "mr_ivw_fe", "mr_ivw_mre",
  "mr_simple_mode", "mr_weighted_mode"
  ))

Analysing 'a2grch' on 'DouvilleNJ2025'

res_single <- mr_singlesnp(mr_dat, all_method = c("mr_ivw_fe", "mr_egger_regression", "mr_weighted_median", "mr_weighted_mode")) %>% as_tibble()
View(res_single)
odd <- generate_odds_ratios(mr_res)
print("Finished Analyzing")

[1] "Finished Analyzing"

write.csv(odd,"results/A2/douvillea2_oddsratios.csv")
write.csv(res_single,"results/A2/douvillea2_res_single.csv")

scatter_p <- mr_scatter_plot(mr_res, mr_dat) 
scatter_out_p <- scatter_p[[1]] + theme_bw() + 
  guides(color=guide_legend(ncol =1)) + 
  theme(
    text = element_text(size = 8), 
  )

scatter_out_p
ggsave("results/A2/scatterplotA2Douville.png", plot = scatter_out_p, units = 'in', height = 4, width = 9)

forrest_p <- mr_forest_plot(res_single)
forrest_p[[1]]
ggsave("results/A2/forrestplotA2Douville.png", plot =forrest_p[[1]] + theme(text = element_text(size = 8)), units = 'in', height = 9, width = 9)

res_loo <- mr_leaveoneout(mr_dat, method = mr_ivw_fe) %>% as_tibble()
loo_p <- mr_leaveoneout_plot(res_loo)
loo_p[[1]]
ggsave("results/A2/looA2Douville.png",plot = loo_p[[1]], units = 'in', height = 9, width = 9)

ggman(filter(tc_ss, p_value < 0.05 & p_value > 1e-100), snp = "SNPID", bp = "base_pair_location", chrom = "chromosome", pvalue = "p_value", relative.positions = TRUE) + 
  theme_classic()

Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as of ggplot2 3.3.4.