04_mvmr_run_analysis.Rmd

---
title: "MVMR analysis of BMIs and mediators to BC"
author: "Marina Vabistsevits"
date: "29/05/2020"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(readr)
library(vroom)

library(tidyr)
library(purrr)
library(tibble)
library(dplyr)

library(TwoSampleMR)
library(MVMR)
```

```{r message=F}
# set path for pre-calculated data, outside the code repo
# `local` / `remote` (reading data from RDSF)
currently_working_env = "local"
source("set_paths.R")
set_paths(currently_working_env)

# metafile
data_lookup<-read_csv(paste0("metadata/data_lookup.csv")) 
phenocor_values<-read_tsv(paste0("metadata/pheno_correlations.tsv")) 

# load functions
source("functions_mvmr.R")
source("functions.R")

# breast cancer dataset
# 1126 full, 1127 ER+, 1128 ER-
breast_cancer_id <- "ieu-a-1126" 
bc_data <- stringr::str_split(breast_cancer_id, "-")[[1]][3]
```

# Load all BMI files
```{r, cache=TRUE, message=F}
# Load BMI exposures
early_bmi_file <- data_lookup %>% filter(trait == "Childhood BMI") %>% pull(tophits_data)
early_bmi_exp <- read_tsv(paste0(data_path_tophits, early_bmi_file))
dim(early_bmi_exp) ## 115 in new

adult_bmi_file <- data_lookup %>% filter(trait == "Adult BMI") %>% pull(tophits_data)
adult_bmi_exp <- read_tsv(paste0(data_path_tophits, adult_bmi_file))
dim(adult_bmi_exp) # 173 in new

height_bmi_file <- data_lookup %>% filter(trait == "Childhood height") %>% pull(tophits_data)
height_bmi_exp <- read_tsv(paste0(data_path_tophits, height_bmi_file))
dim(height_bmi_exp) # 495 in new

# Load BMI outcomes
early_bmi_gwas_file <- data_lookup %>% filter(trait == "Childhood BMI") %>% pull(full_data)
early_bmi_gwas <- vroom(paste0(data_path_gwas, early_bmi_gwas_file))
dim(early_bmi_gwas) 

adult_bmi_gwas_file <- data_lookup %>% filter(trait == "Adult BMI") %>% pull(full_data)
adult_bmi_gwas <- vroom(paste0(data_path_gwas, adult_bmi_gwas_file))
dim(adult_bmi_gwas) 

height_bmi_gwas_file <- data_lookup %>% filter(trait == "Childhood height") %>% pull(full_data)
height_bmi_gwas <- vroom(paste0(data_path_gwas, height_bmi_gwas_file))
dim(height_bmi_gwas) 

```


# Run MVMR 0:  BMIs -> BC

```{r}
# put all exposure and full gwas dat into lists
exposure_list <- list(early_bmi_exp, adult_bmi_exp)#, height_bmi_exp)
full_gwas_list<- list(early_bmi_gwas, adult_bmi_gwas)#, height_bmi_gwas)

# create exposure_dat format
exposure_dat <- get_mv_exposures(exposure_list, full_gwas_list, clump_exposures = T)

#Next, also extract those SNPs from the outcome.
outcome_dat <- extract_outcome_data(exposure_dat$SNP, breast_cancer_id)

#Once the data has been obtained, harmonise so that all are on the same reference allele.
mvdat <- mv_harmonise_data(exposure_dat, outcome_dat)

#Finally, perform the multivariable MR analysis
res_bmis <- mv_multiple(mvdat)

mv_res_bmis<- res_bmis$result %>%
              split_outcome() %>%
              separate(outcome, "outcome", sep="[(]") %>% 
              generate_odds_ratios() %>% 
              select(-id.exposure, -id.outcome)

write_tsv(mv_res_bmis, paste0(results_path, "bmi/merged/mvmr_BMIs-BCAC_", bc_data,".tsv"))
#write_tsv(mv_res_bmis, paste0(results_path, "bmi/merged/mvmr_height_BMIs-BCAC_", bc_data,".tsv"))


# create MVMR package input
mvmr_input <- make_mvmr_input(exposure_dat, outcome.id.mrbase=breast_cancer_id)

# format data to be in MVMR package-compatiable df
mvmr_out <- format_mvmr(BXGs = mvmr_input$XGs %>% select(contains("beta")),  # exposure betas
                        BYG = mvmr_input$YG$beta.outcome,                        # outcome beta
                        seBXGs = mvmr_input$XGs %>% select(contains("se")),      # exposure SEs
                        seBYG = mvmr_input$YG$se.outcome,                        # outcome SEs
                        RSID = mvmr_input$XGs$SNP)                               # SNPs
#  estimate causal effects using method in MVMR package
mvmr_res <-ivw_mvmr(r_input=mvmr_out) %>% 
           tidy_mvmr_output() %>% 
           mutate(exposure = mvmr_input$exposures,
                  outcome = breast_cancer_id)

write_tsv(mvmr_res, paste0(results_path, "bmi/merged/mvmr_BMIs-BCAC_", bc_data, "_using_MVMR",".tsv"))
#write_tsv(mvmr_res, paste0(results_path, "bmi/merged/mvmr_height_BMIs-BCAC_", bc_data, "_using_MVMR",".tsv"))

## Sensitivity tests
# find phenotypic correlation value in table and estimate gencov
print(paste0("Using phenocor for calculating Fst"))
pheno_mat <- filter(phenocor_values, mediator %in% c('Childhood BMI', 'Adult BMI')) %>%  select(early_bmi, adult_bmi) %>%  as.matrix()
colnames(pheno_mat) = rownames(pheno_mat) = c("Childhood BMI", 'Adult BMI')
print(pheno_mat)
#extract SE matrix
se_matrix <- mvmr_out %>% as_tibble() %>% select(contains("sebetaX")) %>% as.data.frame()
#estimate gencov
gen_cov <- phenocov_mvmr(Pcov = as.matrix(pheno_mat), seBXGs = se_matrix)
  #Test for weak instruments
sres <- strength_mvmr(r_input=mvmr_out, gencov=gen_cov)
colnames(sres) = paste(c("Childhood BMI", 'Adult BMI'), "(Fst)")
print(sres)
#Test for horizontal pleiotropy
pres <- pleiotropy_mvmr(r_input=mvmr_out, gencov=gen_cov)
mvmr_sens_df <- sres
mvmr_sens_df$Qstat <- pres$Qstat
mvmr_sens_df$Qpval <- pres$Qpval
write_tsv(mvmr_sens_df, paste0(results_path, "bmi/merged/mvmr_sens_earlyBMI-adultBMI-to-BCAC_", bc_data,"_", Sys.Date(),".tsv"))

```



# Run MVMR for each mediator separately in a loop
# can do by mediators category 

```{r}
# specify group to process if data is in textfiles
current_trait_category <- "hormones"
mediators <- data_lookup %>% filter(trait_category == current_trait_category) %>% filter(!is.na(full_data))  %>%  pull(full_data)

current_trait_category <- "hormones_splitsample"
mediators <- data_lookup %>% filter(trait_category == current_trait_category) %>% filter(!is.na(full_data))  %>%  pull(full_data)

# specify group to process  if the data is in MRBase
current_trait_category <- "reproductive_traits"
mediators <- data_lookup %>% filter(trait_category == current_trait_category) %>% filter(!is.na(mrbase.id)) %>%  pull(mrbase.id)

current_trait_category <- "glycemic_traits"
# run 2 loops: mrbase and text data
mediators <- data_lookup %>% filter(trait_category == current_trait_category) %>% filter(!is.na(mrbase.id)) %>%  pull(mrbase.id)
mediators <- data_lookup %>% filter(trait_category == current_trait_category) %>% filter(!is.na(full_data)) %>%  pull(full_data)

# specify group to process if data is in textfiles
current_trait_category <- "physical_traits"
mediators <- data_lookup %>% filter(trait_category == current_trait_category) %>% filter(!is.na(full_data))  %>%  pull(full_data)
```

# Massive for loop to:
## Load full GWAs and instruments for a mediator
## Run MVMR 1: Childhood BMI and Adult BMI as exposures, Mediator as outcome
## Run MVMR 2: Childhood BMI and Mediator as exposures, Breast cancer as outcome
## Run MVMR 3: Childhood BMI, Adult BMI, and Mediator as exposures, Breast cancer as outcome
## Run MVMR 4: Childhood BMI, Childhood height, and Mediator as exposures, Breast cancer as outcome

```{r message= F}
results_path_sub <-   paste0(results_path, current_trait_category, "/")

run_analysis_1 = F
run_analysis_2 = T
run_analysis_3 = F
run_analysis_4 = F

for (i in 1:length(mediators)){
  
  if( mediators[i] %in% data_lookup$full_data ){
    # it's a text file
    format <-"textfile"
    mediator_name <- data_lookup %>% filter(full_data == mediators[i]) %>% pull(trait)
    mediator_file_name <- data_lookup %>% filter(full_data == mediators[i]) %>% pull(trait_file_name)

    
    # load mediator instruments
    tophits_file <- data_lookup %>% filter(full_data == mediators[i]) %>% pull(tophits_data)
    exposure_mediator <- vroom(paste0(data_path_tophits, tophits_file), 
                               col_types = cols(effect_allele.exposure = col_character())) # to avoid T being read as TRUE

  
    # load full GWAS data (outcomes) and subset to exposure SNPs
    outcome_mediator <- vroom(paste0(data_path_gwas, mediators[i]))


  } else if ( mediators[i] %in% data_lookup$mrbase.id ){
    # it's in mrbase
    format <- "mrbase"
    mediator_name <- data_lookup %>% filter(mrbase.id == mediators[i]) %>% pull(trait)
    mediator_file_name <- data_lookup %>% filter(mrbase.id == mediators[i]) %>% pull(trait_file_name)

    # load mediator instruments
    exposure_mediator <- extract_instruments(mediators[i])
    if (is.null(exposure_mediator)){ exposure_mediator <- extract_instruments(mediators[i], p1 = 10e-07)} # if no SNPs returned, try lower pval     

    exposure_mediator <- clump_data(exposure_mediator)
    exposure_mediator$exposure <- mediator_name
  }
  print(paste0("Currently processing ", mediator_name, " from " , format ))
  
  # make sure we have a place to write
  mediator_dir <- paste0(results_path_sub, mediator_file_name, "/")
  mediator_dir_w_backslash <- gsub(" ", "\\ ", mediator_dir, fixed=T) # create path vector escaping spaces, otherwise sytem call cant process it
  if(!dir.exists(mediator_dir)){ system(paste("mkdir -p", mediator_dir_w_backslash))}
  
  #
  #
  #
  
  
  # Analysis 1. Multivariable MR: Childhood BMI and Adult BMI as exposures, Mediator as outcome
  if (run_analysis_1) {
    print("=========== Running analysis 1 =============")
  
    # put all exposure and full gwas dat into lists
    exposure_list <- list(early_bmi_exp, adult_bmi_exp)
    full_gwas_list<- list(early_bmi_gwas, adult_bmi_gwas)
  
    # create exposure_dat format
    exposure_dat <- get_mv_exposures(exposure_list, full_gwas_list, clump_exposures = T)
  
    #Next, also extract those SNPs from the outcome.
    if (format == "mrbase"){
      outcome_dat <- extract_outcome_data(snps = exposure_dat$SNP, 
                                          outcomes = mediators[i])
    } else if (format == "textfile") {
      outcome_dat <- outcome_mediator %>% filter(SNP %in% exposure_dat$SNP)
    }
    #Once the data has been obtained, harmonise so that all are on the same reference allele.
    mvdat <- mv_harmonise_data(exposure_dat, outcome_dat)
    
    #Finally, perform the multivariable MR analysis
    res_bmis <- mv_multiple(mvdat)
    
    mv_res_bmis<- res_bmis$result %>%
                  split_outcome() %>%
                  generate_odds_ratios() %>% 
                  select(-id.exposure, -id.outcome)
    mv_res_bmis$outcome.full<-mediator_name
    
    write_tsv(mv_res_bmis, paste0(mediator_dir, "mvmr_BMIs-", mediator_file_name, "_using_2SMR_", Sys.Date(),".tsv"))
    
    rm(exposure_list, full_gwas_list, exposure_dat, outcome_dat)
  }  
  #
  #
  #
  
  
  # Analysis 2. Multivariable MR: Childhood BMI and Mediator as exposures, Breast cancer as outcome
  if (run_analysis_2) {
  print("=========== Running analysis 2 =============")
    
    # put all exposure and full gwas dat into lists
    exposure_list <- list(early_bmi_exp, exposure_mediator)
    
    
    if (format == "mrbase"){
      outcome_mediator <- extract_outcome_data(snps = exposure_list %>% 
                                              purrr::reduce(bind_rows) %>% pull(SNP), 
                                       outcomes = mediators[i])
      outcome_mediator$outcome <- mediator_name
    }
    
    full_gwas_list <- list(early_bmi_gwas, outcome_mediator)
    
    # create exposure_dat format
    exposure_dat <- get_mv_exposures(exposure_list, full_gwas_list, clump_exposures = T) 
  
    #Next, also extract those SNPs from the outcome.
    outcome_dat <- extract_outcome_data(exposure_dat$SNP, breast_cancer_id)
    
    #Once the data has been obtained, harmonise so that all are on the same reference allele.
    mvdat <- mv_harmonise_data(exposure_dat, outcome_dat)
    
    #Finally, perform the multivariable MR analysis
    res <- mv_multiple(mvdat)
    
    mv_res<- res$result %>%
                  split_outcome() %>% 
                  separate(outcome, "outcome", sep="[(]") %>% 
                  generate_odds_ratios() %>% 
                  select(-id.exposure, -id.outcome)
    write_tsv(mv_res, paste0(mediator_dir, "mvmr_earlyBMI-", mediator_file_name ,"-to-BCAC_", bc_data,"_using_2SMR_", Sys.Date(),".tsv"))
  
    ## sensitivity analysis
    
    # create MVMR package input
    mvmr_input <- make_mvmr_input(exposure_dat, outcome.id.mrbase=breast_cancer_id)
    
    # format data to be in MVMR package-compatiable df
    mvmr_out <- format_mvmr(BXGs = mvmr_input$XGs %>% select(contains("beta")),  # exposure betas
                        BYG = mvmr_input$YG$beta.outcome,                        # outcome beta
                        seBXGs = mvmr_input$XGs %>% select(contains("se")),      # exposure SEs
                        seBYG = mvmr_input$YG$se.outcome,                        # outcome SEs
                        RSID = mvmr_input$XGs$SNP)                               # SNPs

  #  estimate causal effects using method in MVMR package
    mvmr_res <-ivw_mvmr(r_input=mvmr_out) %>% 
               tidy_mvmr_output() %>% 
               mutate(exposure = mvmr_input$exposures,
                      outcome = breast_cancer_id)
    
    #write_tsv(mvmr_res, paste0(mediator_dir, "mvmr_earlyBMI-", mediator_file_name ,"-to-BCAC_", bc_data,"_using_MVMR_", Sys.Date(),".tsv"))
    
  
    
    # find phenotypic correlation value in table and estimate gencov
    if (mediator_name %in% phenocor_values$mediator){
      print(paste0("Using phenocor for calculating Fst"))
      phenocor_values <- filter(phenocor_values, mediator == mediator_name)
      pheno_mat <- matrix(c(1, phenocor_values$early_bmi,
                        phenocor_values$early_bmi, 1), nrow=2, ncol=2)
      colnames(pheno_mat) = rownames(pheno_mat) = c("Childhood BMI", mediator_name)
      print(pheno_mat)
      #extract SE matrixm ### NB USE MVMR_OUT now
      se_matrix <- mvmr_out %>% as_tibble() %>% select(contains("sebetaX")) %>% as.data.frame()
  
      #estimate gencov
      gen_cov <- phenocov_mvmr(Pcov = as.matrix(pheno_mat), seBXGs = se_matrix)
    } else{
      gen_cov <- 0
    }
    
      #Test for weak instruments
    sres <- strength_mvmr(r_input=mvmr_out, gencov=gen_cov)
    colnames(sres) = paste(c("Childhood BMI", mediator_name), "(Fst)")
    print(sres)
    
    #Test for horizontal pleiotropy
    pres <- pleiotropy_mvmr(r_input=mvmr_out, gencov=gen_cov)
    
    mvmr_sens_df <- sres
    mvmr_sens_df$Qstat <- pres$Qstat
    mvmr_sens_df$Qpval <- pres$Qpval
    
    write_tsv(mvmr_sens_df, paste0(mediator_dir, "mvmr_sens_earlyBMI-", mediator_file_name ,"-to-BCAC_", bc_data,"_", Sys.Date(),".tsv"))
  
    #if (any(as.vector(sres) < 10)){
    #  
    #  print(paste0("Fst is < 10 for one of the exposures; calculating Qhet"))
    #  mvmr_res_Qhet<-qhet_mvmr_tmp(r_input=mvmr_out, pcor=pheno_mat, CI=T, iterations=5) %>% 
    #                  mutate(exposure = unique(exposure_dat$exposure),
    #                  outcome = breast_cancer_id)
    #  write_tsv(mvmr_res_Qhet, paste0(mediator_dir, "mvmr_earlyBMI-", mediator_file_name ,"-to-BCAC_", bc_data,"_using_MVMR_Qhet_", Sys.Date(),".tsv"))
    #}
     
    rm(exposure_list, full_gwas_list, exposure_dat, outcome_dat)
  }
  #
  #
  #
  
 # Analysis 3. Multivariable MR: Childhood BMI, Adult BMI, and Mediator as exposures, Breast cancer as outcome

  if (run_analysis_3) {
  print("=========== Running analysis 3 =============")
  
    exposure_list <- list(early_bmi_exp, adult_bmi_exp, exposure_mediator)
    if (format == "mrbase"){
      outcome_mediator <- extract_outcome_data(snps = exposure_list %>% purrr::reduce(bind_rows)  %>% pull(SNP), 
                                       outcomes = mediators[i])
      outcome_mediator$outcome <- mediator_name
    }
    
    full_gwas_list<- list(early_bmi_gwas, adult_bmi_gwas,  outcome_mediator)
    
    # create exposure_dat format
    exposure_dat <- get_mv_exposures(exposure_list, full_gwas_list, clump_exposures = T)
  
    #Next, also extract those SNPs from the outcome.
    outcome_dat <- extract_outcome_data(exposure_dat$SNP, breast_cancer_id)
    
    #Once the data has been obtained, harmonise so that all are on the same reference allele.
    mvdat <- mv_harmonise_data(exposure_dat, outcome_dat)
    
    #Finally, perform the multivariable MR analysis
    res <- mv_multiple(mvdat)
    
    mv_3exp<- res$result %>%
                  split_outcome() %>% 
                  separate(outcome, "outcome", sep="[(]") %>% 
                  generate_odds_ratios() %>% 
                  select(-id.exposure, -id.outcome)
    
    write_tsv(mv_3exp, paste0(mediator_dir, "mvmr_adultBMI-earlyBMI-", mediator_file_name,"-to-BCAC_", bc_data,"_using_2SMR_", Sys.Date(),".tsv"))
    
    
     ## sensitivity analysis
    
    # create MVMR package input
    mvmr_input <- make_mvmr_input(exposure_dat, outcome.id.mrbase=breast_cancer_id)
    
    # format data to be in MVMR package-compatiable df
    mvmr_out <- format_mvmr(BXGs = mvmr_input$XGs %>% select(contains("beta")),  # exposure betas
                            BYG = mvmr_input$YG$beta.outcome,                        # outcome beta
                            seBXGs = mvmr_input$XGs %>% select(contains("se")),      # exposure SEs
                            seBYG = mvmr_input$YG$se.outcome,                        # outcome SEs
                            RSID = mvmr_input$XGs$SNP)                               # SNPs

    # find phenotypic correlation value in table and estimate gencov
    if (mediator_name %in% phenocor_values$mediator){
      print(paste0("Using phenocor for calculating Fst"))
      
      
      y <- phenocor_values %>% filter(mediator %in% c(mediator_name, "Childhood BMI", "Adult BMI")) %>% 
          select(mediator, adult_bmi, early_bmi) %>% # order matters
          mutate(mediator =ifelse(mediator == mediator_name, 'mediator', mediator)) %>% 
          arrange(mediator) %>%  # always bottom row
          column_to_rownames('mediator') 
        
      pheno_mat <- y %>% mutate(mediator = c(as.vector(as.matrix(y['mediator',])), 1)) %>% 
                         rename("Childhood BMI" = "early_bmi", "Adult BMI" = "adult_bmi")

      print(pheno_mat)
      
      #extract SE matrixm ### NB USE MVMR_OUT now
      se_matrix <- mvmr_out %>% as_tibble() %>% select(contains("sebetaX")) %>% as.data.frame()
  
      #estimate gencov
      gen_cov <- phenocov_mvmr(Pcov = as.matrix(pheno_mat), seBXGs = se_matrix)
    } else{
      gen_cov <- 0
    }
    
      #Test for weak instruments
    sres <- strength_mvmr(r_input=mvmr_out, gencov=gen_cov)
    colnames(sres) = paste(mvmr_input$exposures, "(Fst)")
    print(sres)
    
    #Test for horizontal pleiotropy
    pres <- pleiotropy_mvmr(r_input=mvmr_out, gencov=gen_cov)
    
    mvmr_sens_df <- sres
    mvmr_sens_df$Qstat <- pres$Qstat
    mvmr_sens_df$Qpval <- pres$Qpval
    
    write_tsv(mvmr_sens_df, paste0(mediator_dir, "mvmr_sens_adultBMI-earlyBMI-", mediator_file_name ,"-to-BCAC_", bc_data,"_", Sys.Date(),".tsv"))
  
    
    rm(exposure_list, outcome_mediator, full_gwas_list, exposure_dat, outcome_dat)
  }
  
  
  if (run_analysis_4) {
  print("=========== Running analysis 4 =============")
  
    exposure_list <- list(early_bmi_exp, height_bmi_exp, exposure_mediator)
    if (format == "mrbase"){
      outcome_mediator <- extract_outcome_data(snps = exposure_list %>% purrr::reduce(bind_rows)  %>% pull(SNP), 
                                       outcomes = mediators[i])
      outcome_mediator$outcome <- mediator_name
    }
    
    full_gwas_list<- list(early_bmi_gwas, height_bmi_gwas,  outcome_mediator)
    
    # create exposure_dat format
    exposure_dat <- get_mv_exposures(exposure_list, full_gwas_list, clump_exposures = T)
  
    #Next, also extract those SNPs from the outcome.
    outcome_dat <- extract_outcome_data(exposure_dat$SNP, breast_cancer_id)
    
    #Once the data has been obtained, harmonise so that all are on the same reference allele.
    mvdat <- mv_harmonise_data(exposure_dat, outcome_dat)
    
    #Finally, perform the multivariable MR analysis
    res <- mv_multiple(mvdat)
    
    mv_3exp<- res$result %>%
                  split_outcome() %>% 
                  separate(outcome, "outcome", sep="[(]") %>% 
                  generate_odds_ratios() %>% 
                  select(-id.exposure, -id.outcome)
    
    write_tsv(mv_3exp, paste0(mediator_dir, "mvmr_height-earlyBMI-", mediator_file_name,"-to-BCAC_", bc_data,"_using_2SMR_", Sys.Date(),".tsv"))
    
    
     ## sensitivity analysis
    
    # create MVMR package input
    mvmr_input <- make_mvmr_input(exposure_dat, outcome.id.mrbase=breast_cancer_id)
    
    # format data to be in MVMR package-compatiable df
    mvmr_out <- format_mvmr(BXGs = mvmr_input$XGs %>% select(contains("beta")),  # exposure betas
                            BYG = mvmr_input$YG$beta.outcome,                        # outcome beta
                            seBXGs = mvmr_input$XGs %>% select(contains("se")),      # exposure SEs
                            seBYG = mvmr_input$YG$se.outcome,                        # outcome SEs
                            RSID = mvmr_input$XGs$SNP)                               # SNPs

    # find phenotypic correlation value in table and estimate gencov
    if (mediator_name %in% phenocor_values$mediator){
      print(paste0("Using phenocor for calculating Fst"))
      
      
      y <- phenocor_values %>% filter(mediator %in% c(mediator_name, "Childhood height", "Childhood BMI")) %>% 
          select(mediator, early_bmi, child_height) %>%
          mutate(mediator =ifelse(mediator == mediator_name, 'mediator', mediator)) %>% 
          arrange(mediator) %>%  # always bottom row
          column_to_rownames('mediator')
      pheno_mat <- y %>% mutate(mediator = c(as.vector(as.matrix(y['mediator',])), 1)) %>% 
                         rename("Childhood BMI" = "early_bmi", "Childhood height" = "child_height")

      print(pheno_mat)
      
      #extract SE matrixm ### NB USE MVMR_OUT now
      se_matrix <- mvmr_out %>% as_tibble() %>% select(contains("sebetaX")) %>% as.data.frame()
  
      #estimate gencov
      gen_cov <- phenocov_mvmr(Pcov = as.matrix(pheno_mat), seBXGs = se_matrix)
    } else{
      gen_cov <- 0
    }
    
      #Test for weak instruments
    sres <- strength_mvmr(r_input=mvmr_out, gencov=gen_cov)
    colnames(sres) = paste(mvmr_input$exposures, "(Fst)")
    print(sres)
    
    #Test for horizontal pleiotropy
    pres <- pleiotropy_mvmr(r_input=mvmr_out, gencov=gen_cov)
    
    mvmr_sens_df <- sres
    mvmr_sens_df$Qstat <- pres$Qstat
    mvmr_sens_df$Qpval <- pres$Qpval
    
    write_tsv(mvmr_sens_df, paste0(mediator_dir, "mvmr_sens_earlyBMI-height-", mediator_file_name ,"-to-BCAC_", bc_data,"_", Sys.Date(),".tsv"))
  
    
    rm(exposure_list, outcome_mediator, full_gwas_list, exposure_dat, outcome_dat)
  }
  
  print("Finished.")
}
```


# Merge the MVMR results into one table for each trait category

```{r message=FALSE}
# select trait category
current_trait_category <- "hormones"
traits_subfolders <- list.files(path = paste0(results_path, current_trait_category), full.names = T)
traits_subfolders <- traits_subfolders[!grepl("merged", traits_subfolders)] # drop merged folder from list

# use this BC type
bc_data <- "1126"

# use this MVMR type
mvmr_method<- "using_2SMR" #"using_MVMR" # "using_2SMR"

mvmr_types <- c("mvmr_BMIs-", "mvmr_earlyBMI-", "mvmr_adultBMI-earlyBMI-", "mvmr_height-earlyBMI-")
mvmr_types <- c("mvmr_height-earlyBMI-")



# for each mvmr type, load all indiv trait files, merge, save
for (mvmr_type in mvmr_types) {

  print(paste0("Current trait category: ", current_trait_category))
  
  trait_files <- c()
  for ( folder_path in traits_subfolders){
    file_path <- list.files(path = folder_path, pattern = paste0(mvmr_method, "_2"), full.names = T)
    trait_files <- c(trait_files, file_path)
   
    # for mvmr with breast cancer, subset to the correct type
    if (mvmr_type != "mvmr_BMIs-"){
      trait_files <- trait_files[grepl(bc_data, trait_files)]}
    
    # merged files created today OR specify date 
    trait_files <- trait_files[grepl(Sys.Date(), trait_files)] ####### PAY ATTENTION HERE
    #trait_files <- trait_files[grepl("2020-07-29", trait_files)] ####### PAY ATTENTION HERE

    
  }
  print(paste0("Reading data from analysis: ", mvmr_type))
  
  # read all individual mediators
  l <- lapply(trait_files, read_tsv)
  df <- l %>% purrr::reduce(bind_rows) 
  
  if (mvmr_type == "mvmr_BMIs-") {
    df <- add_trait_type_out(df, current_trait_category)
    
    file_prefix <- paste0(mvmr_type,"to-", current_trait_category)
    
  } else if (mvmr_type == "mvmr_earlyBMI-") {
     df <- add_trait_type_exp(df, current_trait_category)
     trait_type_vec<- select(df, trait_type) %>% drop_na() %>% slice(rep(1:n(), each = 2))
     df$trait_type <- trait_type_vec$trait_type
     
     file_prefix <- paste0(mvmr_type, current_trait_category, "-to-BCAC_", bc_data)

  } else if (mvmr_type == "mvmr_adultBMI-earlyBMI-") {
     df <- add_trait_type_exp(df, current_trait_category)
     trait_type_vec<- select(df, trait_type) %>% drop_na() %>% slice(rep(1:n(), each = 3)) 
     df$trait_type <- trait_type_vec$trait_type
     
     file_prefix <- paste0(mvmr_type, current_trait_category, "-to-BCAC_", bc_data)
     
  } else if (mvmr_type == "mvmr_height-earlyBMI-") {
     df <- add_trait_type_exp(df, current_trait_category)
     trait_type_vec<- select(df, trait_type) %>% drop_na() %>% slice(rep(1:n(), each = 3)) 
     df$trait_type <- trait_type_vec$trait_type
     
     file_prefix <- paste0(mvmr_type, current_trait_category, "-to-BCAC_", bc_data)
  }
     
  write_tsv(df, paste0(results_path, current_trait_category, "/merged/merged_", file_prefix, "_", mvmr_method, ".tsv"))
  print(paste0(" -> Merged and saved."))

}
```


# Merge sensitivity analysis into one df

```{r message=FALSE}
# select trait category
current_trait_category <- "reproductive_traits"
traits_subfolders <- list.files(path = paste0(results_path, current_trait_category), full.names = T)
traits_subfolders <- traits_subfolders[!grepl("merged", traits_subfolders)] # drop merged folder from list

# use this BC type
bc_data <- "1126"

mvmr_types <- c("mvmr_sens_earlyBMI-")

# for each mvmr type, load all indiv trait files, merge, save
for (mvmr_type in mvmr_types) {

  print(paste0("Current trait category: ", current_trait_category))

  trait_files_sens <- c()
  for ( folder_path in traits_subfolders){
    file_path <- list.files(path = folder_path, pattern = paste0(mvmr_type, "*"), full.names = T)
    trait_files_sens <- c(trait_files_sens, file_path)
    
    # for mvmr with breast cancer, subset to the correct type
    trait_files_sens <- trait_files_sens[grepl(bc_data, trait_files_sens)]
    
    # merged files created today OR specify date 
    #trait_files_sens <- trait_files_sens[grepl(Sys.Date(), trait_files_sens)] ####### PAY ATTENTION HERE
    trait_files_sens <- trait_files_sens[grepl("2020-08-18", trait_files_sens)] ####### PAY ATTENTION HERE
    
  }
  print(paste0("Reading data from analysis: ", mvmr_type))
  
  read_custom <- function(file_list){
    x<-read_tsv(file_list)
    x$mediator <-colnames(x)[2]
    colnames(x)[1:2]<-c("Fst_BMI", "Fst_mediator")
    return(x)
  }
  
  # read all individual mediators
  l <- lapply(trait_files_sens, read_custom)
  df <- l %>% purrr::reduce(bind_rows) 
  
  file_prefix <- paste0(mvmr_type, current_trait_category, "-to-BCAC_", bc_data)
 
  write_tsv(df, paste0(results_path, current_trait_category, "/merged/merged_sens", file_prefix,"new.tsv"))
  print(paste0(" -> Merged and saved.")) 
  
  
  # check of any got Fst < 10, then collect Qhet results fro those
  if (any(df$Fst_BMI) <10 | any(df$Fst_mediator) <10 ){
   
    qhet_files <- c()
    for ( folder_path in traits_subfolders){
    file_path <- list.files(path = folder_path, pattern = paste0("Qhet"), full.names = T)
    qhet_files <- c(qhet_files, file_path)
    
    # for mvmr with breast cancer, subset to the correct type
    qhet_files <- qhet_files[grepl(bc_data, qhet_files)]
    
    # merged files created today OR specify date 
    #qhet_files <- qhet_files[grepl(Sys.Date(), qhet_files)] ####### PAY ATTENTION HERE
    qhet_files <- qhet_files[grepl("2020-07-29", qhet_files)] ####### PAY ATTENTION HERE
  }
      # read all individual mediators
  l <- lapply(qhet_files, read_tsv)
  df_qhet <- l %>% purrr::reduce(bind_rows) 
    
  file_prefix <- paste0(mvmr_type, current_trait_category, "-to-BCAC_", bc_data, "_Qhet")
 
  write_tsv(df_qhet, paste0(results_path, current_trait_category, "/merged/merged_", file_prefix,"new.tsv"))
  print(paste0(" -> Merged and saved."))   
    
  
  }
  
  
}

```