06_create_forest_plots.Rmd

---
title: "Generate figures"
author: "Marina Vabistsevits"
date: "`r Sys.Date()`"
output: html_document

---
```{r setup, include=FALSE}
library(cowplot)
library(data.table)
library(dplyr)
library(ggplot2)
library(readr)
library(tidyr)
library(TwoSampleMR)
library(wesanderson)
library(cowplot)
knitr::opts_chunk$set(echo = TRUE)
```

```{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")) 
bcac_lookup<-read_csv(paste0("metadata/data_lookup_BCAC.csv")) %>% filter(!grepl("BRCA1", trait))

source("functions.R")

# collecting data for Source Data file here:
source_data<- list()

```



# Figure S1

Collecting data and plotting 4 separate plots:

1. uniMR: BMIs -> MD (unadjusted)
2. uniMR: BMIs -> MD (adjusted)

3.  mvMR: BMIs -> MD (unadjusted)
4.  mvMR: BMIs -> MD (adjusted)


## prep: md_unadj

```{r}
current_trait_category <- "md_unadj"
pal<-c(unname(yarrr::piratepal("pony")))[c(8,2)]

# uniMR - generated in 02
dat<- read_tsv(paste0(results_path, current_trait_category, "/merged_mr_BMI-to-", current_trait_category,".tsv")) %>% 
      filter(method == "Inverse variance weighted") 

unadj_unimr<-
  ggplot(dat, aes(y=exposure, x=b, label=outcome, colour=exposure)) +
  geom_errorbarh(aes(xmin=lo_ci, xmax=up_ci), height=0.2) +
  geom_point(size=2)+
  xlim(-1.2,1.2)+
  scale_color_manual(values=pal)+
  scale_shape_manual(values = c(19)) +
  geom_vline(xintercept=0, linetype=3,) +
  theme_minimal_hgrid(9, rel_small = 1) +
  facet_wrap(~outcome, ncol=1,switch = 'x')+
  labs(color = "",y = "", x = "(total) effect size")+
  theme(legend.position = "none", plot.title.position  = "plot")

# mvmr
dat <- read_tsv(paste0(results_path, current_trait_category, 
                       "/mvmr/merged/merged_mvmr_BMIs-", current_trait_category,"-to-BCAC.tsv"))

unadj_mvmr<-
  ggplot(dat, aes(y=exposure, x=b, label=outcome, colour=exposure)) +
  geom_errorbarh(aes(xmin=lo_ci, xmax=up_ci), height=0.2) +
  geom_point(size=2)+
  xlim(-1.2,1.2)+
  scale_color_manual(values=pal)+
  scale_shape_manual(values = c(19)) +
  geom_vline(xintercept=0, linetype=3,) +
  theme_minimal_hgrid(9, rel_small = 1) +
  facet_wrap(~outcome, ncol=1, switch = "x")+
  labs(color = "",y = "", x = "direct (effect size)" )+
  theme(legend.position = "none", plot.title.position  = "plot")

```

## prep: md (adjusted)

```{r}
current_trait_category <- "md"
pal<-c(unname(yarrr::piratepal("pony")))[c(8,2)]

# uniMR - generated in 02
dat<- read_tsv(paste0(results_path, current_trait_category, "/merged_mr_BMI-to-", current_trait_category,".tsv")) %>% 
      filter(method == "Inverse variance weighted") 

adj_unimr<-
  ggplot(dat, aes(y=exposure, x=b, label=outcome, colour=exposure)) +
  geom_errorbarh(aes(xmin=lo_ci, xmax=up_ci), height=0.2) +
  geom_point(size=2)+
  xlim(-1.2,1.2)+
  scale_color_manual(values=pal)+
  scale_shape_manual(values = c(19)) +
  geom_vline(xintercept=0, linetype=3,) +
  theme_minimal_hgrid(9, rel_small = 1) +
  facet_wrap(~outcome, ncol=1,switch = 'x')+
  labs(color = "",y = "", x = "(total) effect size")+
  theme(legend.position = "none", plot.title.position  = "plot")

#mvmr
dat <- read_tsv(paste0(results_path, current_trait_category, 
                       "/mvmr/merged/merged_mvmr_BMIs-", current_trait_category,"-to-BCAC.tsv"))
adj_mvmr<-
  ggplot(dat, aes(y=exposure, x=b, label=outcome, colour=exposure)) +
  geom_errorbarh(aes(xmin=lo_ci, xmax=up_ci), height=0.2) +
  geom_point(size=2)+
  xlim(-1.2,1.2)+
  scale_color_manual(values=pal)+
  scale_shape_manual(values = c(19)) +
  geom_vline(xintercept=0, linetype=3,) +
  theme_minimal_hgrid(9, rel_small = 1) +
  facet_wrap(~outcome, ncol=1, switch = "x")+
  labs(color = "",y = "", x = "direct (effect size)" )+
  theme(legend.position = "none", plot.title.position  = "plot")
```

## complete figure

```{r}
# collect plots for supl figure here

#Figure S1
pl<-cowplot::plot_grid(
              unadj_unimr, NULL, unadj_mvmr,
              NULL,        NULL, NULL,
              adj_unimr,   NULL, adj_mvmr, 
            
              rel_widths = c(1, 0.1, 1,
                             1, 0.1, 1),
              rel_heights = c(1, 0.1, 1), 
               labels = c("a", "", "b",
                            "", "","",
                            "c", "", "d" ), label_size = 12,
                #labels= c("Unadjusted MD GWASs, MR", "Unadjusted MD GWASs, MVMR",
                #          "Adjusted for BMI MD GWASs, MR", "Adjusted for BMI MD GWASs, MVMR"),
                align='vh', vjust=1, nrow = 3, ncol=3)


```




# Figure 2

##  prep: total and direcrt effect of childhood body size and age at menarche

```{r}
current_trait_category <- "md_unadj"
# cBMI -> MD total effect
dat_total1<- read_tsv(paste0(results_path, current_trait_category, "/merged_mr_BMI-to-", current_trait_category,".tsv")) %>% 
      filter(method == "Inverse variance weighted") %>% 
      filter(exposure == "Childhood body size") %>% 
      mutate(effect = "Total effect")

# ageatM -> MD total effect
dat_total2<- read_tsv( paste0(results_path, current_trait_category,"/", current_trait_category,"_w_traits_biMR_results.tsv")) %>% 
      filter(method == "Inverse variance weighted") %>% 
      filter(dir == 1) %>% 
      filter(exposure == "Age when periods started (menarche)") %>% # ukb
      mutate(exposure = ifelse(exposure == "Age when periods started (menarche)", "Age at menarche", exposure)) %>%
      mutate(effect = "Total effect")


# mvmr results
dat_direct<- read_tsv(paste0(results_path, "md_w_traits/","mvmr_results_w_cBMI_",current_trait_category,".tsv")) %>% 
      mutate(exposure = ifelse(exposure == "age_at_menarche_ukb", "Age at menarche", exposure)) %>%
      filter(pair == "age_at_menarche_ukb") %>% 
      mutate(effect = "Direct effect") 


dat_merged <- bind_rows(dat_total1, dat_total2, dat_direct) %>%
              select(effect,  everything(), -method) %>% 
              mutate(effect=factor(effect, levels = c("Total effect", "Direct effect"))) %>% 
              mutate(exposure=factor(exposure, levels = c("Childhood body size",  "Age at menarche")),
                     exposure = factor(exposure, levels = rev(levels(exposure)))) %>% 
              arrange(exposure) %>% 
              mutate(beta_CI = paste0(round(b,2), " [",round(lo_ci,2) ,":",round(up_ci,2), "]")) %>% 
              mutate(effect_direction = ifelse(or_lci95 > 1 & or_uci95 >= 1, 'ok',
                                        ifelse(or_lci95 < 1 & or_uci95 <= 1, 'ok', 'overlaps null'))) 


pal<-c(unname(yarrr::piratepal("pony")))[c(1,2)]
p<-ggplot(dat_merged, aes(y=exposure, x=b, label=outcome, shape=effect_direction, colour=exposure)) +
  geom_errorbarh(aes(xmin=lo_ci, xmax=up_ci), height=0.2) +
  geom_point(size=2.5)+
  #xlim(-1.2,0.5)+
  scale_color_manual(values=pal)+
  scale_shape_manual(values=c(16,1))+
  geom_vline(xintercept=0, linetype=3) +
  geom_text(aes(label=beta_CI),hjust=0.5, vjust=-1, size =2.9, color = '#3B3B3B')+
  theme_minimal_hgrid(10, rel_small = 1) +
  panel_border() +
  facet_grid(rows=vars(outcome), cols=vars(effect)) +
  scale_y_discrete(position = "right")+
  labs(color = "Exposure",y = "Outcome", x = "effect size"
       #title= paste0("Univariate MR, 95% CI; direction 1: trait -> MD") )+
       #title= paste0("MVMR direct effect, accounted for Childhood body size, 95% CI"
                      )+
  theme(legend.position = "bottom", legend.justification = "center",
        plot.title.position  = "plot",
        strip.text.x = element_text(face = "bold"),
        strip.text.y = element_text(face = "bold"),
        axis.text.y=element_blank())+
  guides(shape="none", colour=guide_legend(reverse = T))


```



## prep:  total and direcrt effect of childhood body size and adult body size

```{r}
current_trait_category <- "md_unadj"
# cBMI or aBMI-> MD total effect
dat_total<- read_tsv(paste0(results_path, current_trait_category, "/merged_mr_BMI-to-", current_trait_category,".tsv")) %>% 
      filter(method == "Inverse variance weighted") %>% 
      #filter(exposure == "Childhood body size") %>% 
      mutate(effect = "Total effect")


# mvmr results
dat_direct<- read_tsv(paste0(results_path, current_trait_category, "/mvmr/merged/merged_mvmr_BMIs-", current_trait_category,"-to-BCAC.tsv")) %>% 
  mutate(effect = "Direct effect") 


dat_merged_BMIS <- bind_rows(dat_total, dat_direct) %>%
              select(effect,  everything(), -method) %>% 
              mutate(effect=factor(effect, levels = c("Total effect", "Direct effect"))) %>% 
              mutate(exposure=factor(exposure, levels = c("Childhood body size",  "Adult body size")),
                     exposure = factor(exposure, levels = rev(levels(exposure)))) %>% 
              arrange(exposure) %>% 
              mutate(beta_CI = paste0(round(b,2), " [",round(lo_ci,2) ,":",round(up_ci,2), "]")) %>% 
              mutate(effect_direction = ifelse(or_lci95 > 1 & or_uci95 >= 1, 'ok',
                                        ifelse(or_lci95 < 1 & or_uci95 <= 1, 'ok', 'overlaps null'))) 


pal<-c(unname(yarrr::piratepal("pony")))[c(1,2)]
p<-ggplot(dat_merged_BMIS, aes(y=exposure, x=b, label=outcome, shape=effect_direction, colour=exposure)) +
  geom_errorbarh(aes(xmin=lo_ci, xmax=up_ci), height=0.2) +
  geom_point(size=2.5)+
  #xlim(-1.2,0.5)+
  scale_color_manual(values=pal)+
  scale_shape_manual(values=c(16,1))+
  geom_vline(xintercept=0, linetype=3) +
  geom_text(aes(label=beta_CI),hjust=0.5, vjust=-1, size =2.9, color = '#3B3B3B')+
  theme_minimal_hgrid(10, rel_small = 1) +
  panel_border() +
  facet_grid(rows=vars(outcome), cols=vars(effect), switch = "y") +
  scale_y_discrete(position = "left")+
  labs(color = "Exposure",y = "Outcome", x = "effect size"
       #title= paste0("Univariate MR, 95% CI; direction 1: trait -> MD") )+
       #title= paste0("MVMR direct effect, accounted for Childhood body size, 95% CI"
                      )+
  theme(legend.position = "bottom", legend.justification = "center",
        plot.title.position  = "plot",
        strip.text.x = element_text(face = "bold"),
        strip.text.y = element_text(face = "bold"),
        axis.text.y=element_blank())+
  guides(shape="none", colour=guide_legend(reverse = T))


```


## complete figure - combined plot of cBMI, aBMI, aaM - total and direct effects

```{r}

dat_all <- bind_rows(dat_merged %>% mutate(cov = "Direct effect (MVMR) "),
                     dat_merged_BMIS %>% mutate(cov = "Direct effect (MVMR)")) %>% 
  select(-id.exposure, -id.outcome, -dir, -pair, -outcome.full ) %>% 
  mutate(cov = ifelse(effect == "Total effect", "Total effect (MR)", cov)) %>% 
    distinct() %>% 
  mutate(cov = factor(cov, levels = c("Total effect (MR)",
                                      "Direct effect (MVMR)", 
                                      "Direct effect (MVMR) "))) %>% 
  mutate(exposure = factor(exposure, levels = rev(c("Childhood body size", 
                                                "Adult body size",
                                                "Age at menarche"))))

pal<-c(unname(yarrr::piratepal("pony")))[c(1,8,2)]
## adhoc change outcome label
dat_all <- dat_all %>% ####### PICK ONE!!!
  # long names
  mutate(outcome_short = case_when(outcome == "Dense area" ~ "Dense \narea",
                                  outcome == "Non-dense area" ~ "Non-dense \narea",
                                  outcome == "Percent density" ~ "Percent \ndensity")) #%>% 
  # ALT short names - like for poster!
  #mutate(outcome_short = case_when(outcome == "Dense area" ~ "DA",
  #                                  outcome == "Non-dense area" ~ "NDA",
  #                                  outcome == "Percent density" ~ "PD"))

pX<-ggplot(dat_all, aes(y=exposure, x=b, label=outcome_short, shape=effect_direction, colour=exposure)) +
  geom_errorbarh(aes(xmin=lo_ci, xmax=up_ci), height=0.2) +
  geom_point(size=2)+
  xlim(-1.2,1.2)+
  scale_color_manual(values=pal)+
  scale_shape_manual(values=c(16,1))+
  geom_vline(xintercept=0, linetype=3, color = '#3B3B3B') +
  #geom_text(aes(label=beta_CI),hjust=0.28, vjust=-1, size =2.9, color = '#636363')+
  theme_minimal_hgrid(8, rel_small = 1) +
  #panel_border() +
  facet_grid(rows=vars(outcome_short), cols=vars(cov))+#, switch = "y") +
  scale_y_discrete(position = "right")+
  labs(color = "Exposure",y = "Outcome",
       x = "\nbeta coefficient (1 SD change in MD phenotype per body size category or age SD change)\n")+
  theme(legend.position = "bottom", legend.justification = "left", legend.title = element_text(face = "bold"),
        plot.title.position  = "plot",
        strip.text.x = element_text(face = "bold"),
        #strip.text.y = element_text(face = "bold"),
        axis.title.y = element_text(face="bold"),
        axis.text.y=element_blank(),
        axis.text.x =element_text(size = 8),
        panel.spacing = unit(0.5, "lines"))+
  guides(shape="none", colour=guide_legend(reverse = T, nrow = 1,byrow=TRUE))

source_data[["fig2"]]<-dat_all %>% select(-starts_with("or"))

ggsave(paste0("figures_manuscript/figure2_poster.png"), # uses ABBS
       plot=pX, scale=1.2, 
       width=15, height=5,
       units=c("cm"), dpi=300, limitsize=F)

ggsave(paste0("figures_manuscript/figure2_paper.svg"),
       plot=pX, #scale=1, 
       width=180, height=79,
       units=c("mm"), dpi=300, limitsize=F)
```




## extra: combined plot of total and direcrt effect of childhood body size and IGF

```{r}
current_trait_category <- "md_unadj"
# cBMI -> MD total effect
dat_total1<- read_tsv(paste0(results_path, current_trait_category, "/merged_mr_BMI-to-", current_trait_category,".tsv")) %>% 
      filter(method == "Inverse variance weighted") %>% 
      filter(exposure == "Childhood body size") %>% 
      mutate(effect = "Total effect")

# IGF -> MD total effect
dat_total2<- read_tsv( paste0(results_path, current_trait_category,"/", current_trait_category,"_w_traits_biMR_results.tsv")) %>% 
      filter(method == "Inverse variance weighted") %>% 
      filter(dir == 1) %>% 
      filter(exposure == "IGF") %>% # ukb
      mutate(effect = "Total effect")


# mvmr results
dat_direct<- read_tsv(paste0(results_path, "md_w_traits/","mvmr_results_w_cBMI_",current_trait_category,".tsv")) %>% 
      filter(pair == "igf") %>% 
      mutate(effect = "Direct effect") 


dat_merged <- bind_rows(dat_total1, dat_total2, dat_direct) %>%
              select(effect,  everything(), -method) %>% 
              mutate(effect=factor(effect, levels = c("Total effect", "Direct effect"))) %>% 
              mutate(exposure=factor(exposure, levels = c("Childhood body size",  "IGF")),
                     exposure = factor(exposure, levels = rev(levels(exposure)))) %>% 
              arrange(exposure) %>% 
              mutate(beta_CI = paste0(round(b,2), " [",round(lo_ci,2) ,":",round(up_ci,2), "]")) %>% 
              mutate(effect_direction = ifelse(or_lci95 > 1 & or_uci95 >= 1, 'ok',
                                        ifelse(or_lci95 < 1 & or_uci95 <= 1, 'ok', 'overlaps null'))) 


pal<-c(unname(yarrr::piratepal("pony")))[c(9,2)]
p<-ggplot(dat_merged, aes(y=exposure, x=b, label=outcome, shape=effect_direction, colour=exposure)) +
  geom_errorbarh(aes(xmin=lo_ci, xmax=up_ci), height=0.2) +
  geom_point(size=2.5)+
  #xlim(-1.2,0.5)+
  scale_color_manual(values=pal)+
  scale_shape_manual(values=c(16,1))+
  geom_vline(xintercept=0, linetype=3) +
  geom_text(aes(label=beta_CI),hjust=0.5, vjust=-1, size =2.9, color = '#3B3B3B')+
  theme_minimal_hgrid(10, rel_small = 1) +
  panel_border() +
  facet_grid(rows=vars(outcome), cols=vars(effect)) +
  scale_y_discrete(position = "right")+
  labs(color = "Exposure",y = "Outcome", x = "effect size"
       #title= paste0("Univariate MR, 95% CI; direction 1: trait -> MD") )+
       #title= paste0("MVMR direct effect, accounted for Childhood body size, 95% CI"
                      )+
  theme(legend.position = "bottom", legend.justification = "center",
        plot.title.position  = "plot",
        strip.text.x = element_text(face = "bold"),
        strip.text.y = element_text(face = "bold"),
        axis.text.y=element_blank())+
  guides(shape="none", colour=guide_legend(reverse = T))


```

# Figure S2 

## complete figure - plot for aaM to MD, comparing md and md_unadj
```{r}
current_trait_category <- "md"
adj_md<- read_tsv(paste0(results_path, current_trait_category,"/", current_trait_category,"_w_traits_biMR_results.tsv")) %>% 
  filter(exposure == "Age when periods started (menarche)") %>% 
  filter(method == "Inverse variance weighted") %>% 
  mutate(analysis = "Adjusted for BMI MD GWAS")

current_trait_category <- "md_unadj"
unadj_md<- read_tsv(paste0(results_path, current_trait_category,"/", current_trait_category,"_w_traits_biMR_results.tsv"))%>% 
  filter(exposure == "Age when periods started (menarche)") %>% 
  filter(method == "Inverse variance weighted") %>% 
  mutate(analysis = "Unadjusted MD GWAS")

merged <-bind_rows(adj_md,unadj_md ) %>%
        mutate(exposure = "Age at menarche") %>% 
        mutate(analysis = factor(analysis, levels = c("Unadjusted MD GWAS", "Adjusted for BMI MD GWAS"))) %>% 
        mutate(effect_direction = ifelse(or_lci95 > 1 & or_uci95 >= 1, 'ok',
                                        ifelse(or_lci95 < 1 & or_uci95 <= 1, 'ok', 'overlaps null'))) 


pal<-c(unname(yarrr::piratepal("pony")))[c(1)]
pX<-ggplot(merged, aes(y=exposure, x=b, label=outcome, shape=effect_direction, colour=exposure)) +
  geom_errorbarh(aes(xmin=lo_ci, xmax=up_ci), height=0.2) +
  geom_point(size=2.5)+
  xlim(-0.5,0.5)+
  scale_color_manual(values=pal)+
  scale_shape_manual(values=c(16,1))+
  geom_vline(xintercept=0, linetype=3, color = '#3B3B3B') +
  #geom_text(aes(label=beta_CI),hjust=0.28, vjust=-1, size =2.9, color = '#636363')+
  theme_minimal_hgrid(10, rel_small = 1) +
  #panel_border() +
  facet_grid(rows=vars(outcome), cols=vars(analysis), switch = "y") +
  scale_y_discrete(position = "right")+
  labs(color = "Exposure",y = "Outcome",
       x = "\nbeta coefficient (1 SD change in MD phenotype per age SD change)\n")+
  theme(legend.position = "bottom", legend.justification = "center", ,
        plot.title.position  = "plot",
        strip.text.x = element_text(face = "bold"),
        strip.text.y = element_text(face = "bold"),
        axis.text.y=element_blank(),
        axis.text.x =element_text(size = 8),
        strip.text.y.left = element_text(angle = 0),
        panel.spacing = unit(1, "lines"))+
  guides(shape="none", colour=guide_legend(reverse = T))

```



# Figures 3 and S3

## total effects only - for building S3 and loading data for 2
```{r message=F}
current_trait_category <- "md_unadj"

dat_all<-tibble()
for (bc_data in bcac_lookup$trait_file_name[-10]){
  dat <- read_tsv(paste0(results_path, current_trait_category,
           "/","merged_mr_mediators_", current_trait_category,"-to-BCAC_", bc_data,".tsv")) 
  dat_all<-bind_rows(dat_all,dat)
}
dat_all<- dat_all %>%
  mutate(OR_CI =  paste0(round(or,2), " [",round(or_lci95,2) ,":",round(or_uci95,2), "]")) %>%  
  filter(method %in% c("Inverse variance weighted", "Wald ratio")) %>% 
  mutate(outcome= case_when(outcome == 'Breast cancer' & id.outcome == 'ieu-a-1126' ~  'Breast cancer 2017',
                            outcome == 'Breast cancer' & id.outcome != 'ieu-a-1126' ~  'Breast cancer 2020',
                            TRUE ~ outcome)) %>% 
  mutate(exposure = case_when(grepl("non-dense", exposure, ignore.case = T) ~ 'Non-dense area',
                             grepl("Percent", exposure) ~ 'Percent density', 
                             TRUE ~ 'Dense area')) 

# rearrange to be in the required order
dat_all<-arrange(dat_all, exposure) %>% mutate(exposure=as.factor(exposure))
dat_all$exposure <- forcats::fct_inorder(dat_all$exposure) %>% forcats::fct_rev() 

dat_all2 <- dat_all %>% select(exposure, outcome, or, or_lci95, or_uci95, OR_CI) %>% 
  mutate(outcome = case_when(outcome =="Breast cancer 2020" ~ "BCAC 2020",
                             outcome =="Breast cancer 2017" ~ "BCAC 2017",
                             outcome =="ER+ Breast cancer" ~ "ER+",
                             outcome =="ER- Breast cancer" ~ "ER-",
                                           outcome =="Luminal A"    ~ "Luminal A",   
                                           outcome =="Luminal B"   ~ "Luminal B1",  
                                           outcome =="Luminal B-HER2 negative" ~ "Luminal B2" ,
                                           outcome =="HER2 enriched"  ~ "HER2-enriched" ,
                                           outcome =="TNBC Basal-like"  ~ "TNBC" )) %>% 
  filter(outcome != "BCAC 2020") %>% 
   mutate(outcome = case_when(outcome == "BCAC 2017" ~ "Overall breast cancer",
                              outcome == "ER+" ~ "ER+ breast cancer",
                              outcome == "ER-" ~ "ER- breast cancer",
                              outcome == "TNBC" ~ "Triple-negative subtype",
                              outcome == "Luminal A" ~ "Luminal A subtype",
                              outcome == "Luminal B1" ~ "Luminal B1 subtype",
                              outcome == "Luminal B2" ~ "Luminal B2 subtype",
                              outcome == "HER2-enriched" ~ "HER2-enriched subtype" )) %>% 
   mutate(outcome = factor(outcome, levels = rev(c("Overall breast cancer",  "ER+ breast cancer","ER- breast cancer"  ,
                                                "Luminal A subtype", "Luminal B1 subtype" , "Luminal B2 subtype",
                                               "HER2-enriched subtype" , "Triple-negative subtype"))))%>% 
   mutate(effect_direction = ifelse(or_lci95 > 1 & or_uci95 >= 1, 'ok',
                              ifelse(or_lci95 < 1 & or_uci95 <= 1, 'ok', 'overlaps null')))


pal<- c("#EB5291FF", "#FBBB68FF" ,"#F5BACFFF", "#9DDAF5FF", "#6351A0FF" ,"#FFEA5E", "#1794CEFF", "#972C8DFF") # yarrr::piratepal("pony")
p<-ggplot(dat_all2,
           aes(y=exposure, x=or,  colour=outcome, shape = effect_direction)) +
  geom_point(size=3, 
             position = ggstance::position_dodgev(height = 0.7))+
  geom_errorbarh(aes(xmin=or_lci95, xmax=or_uci95), height=.1, 
            position = ggstance::position_dodgev(height = 0.7)) +
  scale_x_log10(breaks = c(0,1,2,3))+

  theme_minimal_vgrid(11)+
  geom_vline(xintercept=1, linetype='longdash') +
  scale_color_manual(values=pal)+
  scale_shape_manual(values=c(16,1))+
  theme(strip.text = element_text(face = 'bold'),
        strip.text.y = element_text(face = 'bold'),
        panel.grid.major.y = element_blank(),
        panel.grid.major.x = element_line(),
        )+
  panel_border()+
  labs(color = "",y = "", x = "Odds ratio", shape="")+
  theme(legend.position = "none", legend.direction = "vertical", 
        plot.title.position = "plot")+
  guides(colour = guide_legend(reverse = TRUE), shape="none")


```


## Figure S3
```{r}
# figure S3 - produce  two objects with the right data
library(cowplot)
# save "md" plot
p_adj <- p
# save "md_unadj"
p_unadj <- p

x<- plot_grid( p_unadj, p_adj, labels = c("Unadjusted", "Adjusted for BMI"), hjust = -0.1,
               ncol = 2, label_size = 12, rel_widths  = c(0.45, 0.55))

```


## complete figure
```{r}
### same as main plot but also with  direct estimates: 
dat_direct <- read_tsv(paste0(results_path, current_trait_category, "/mvmr/merged/merged_mvmr_earlyBMI-", current_trait_category,"-to-BCAC.tsv")) %>%  # made in 04
              filter(exposure != "Childhood body size") %>% 
              filter(outcome!= "BRCA1_TN") %>% 
              select(exposure, outcome, or, or_lci95, or_uci95) %>% 
              mutate(effect_direction = ifelse(or_lci95 > 1 & or_uci95 >= 1, 'ok',
                              ifelse(or_lci95 < 1 & or_uci95 <= 1, 'ok', 'overlaps null'))) %>% 
    mutate(outcome = case_when(outcome =="Breast cancer (new)" ~ "BCAC 2020",
                                         outcome =="Breast cancer" ~ "BCAC 2017",
                                         outcome =="ER+ breast cancer" ~ "ER+",
                                         outcome =="ER- breast cancer" ~ "ER-",
                                                       outcome =="Luminal A"    ~ "Luminal A",   
                                                       outcome =="Luminal B, HER2+"   ~ "Luminal B1",  
                                                       outcome =="Luminal B, HER2-" ~ "Luminal B2" ,
                                                       outcome =="HER2-enriched"  ~ "HER2-enriched" ,
                                                       outcome =="TN Basal"  ~ "TNBC" ))  %>% 
  filter(outcome != "BCAC 2020") %>% 
   mutate(outcome = case_when(outcome == "BCAC 2017" ~ "Overall breast cancer",
                              outcome == "ER+" ~ "ER+ breast cancer",
                              outcome == "ER-" ~ "ER- breast cancer",
                              outcome == "TNBC" ~ "Triple-negative subtype",
                              outcome == "Luminal A" ~ "Luminal A subtype",
                              outcome == "Luminal B1" ~ "Luminal B1 subtype",
                              outcome == "Luminal B2" ~ "Luminal B2 subtype",
                              outcome == "HER2-enriched" ~ "HER2-enriched subtype" )) %>% 
   mutate(outcome = factor(outcome, levels = rev(c("Overall breast cancer",  "ER+ breast cancer","ER- breast cancer",
                                                "Luminal A subtype", "Luminal B1 subtype" , "Luminal B2 subtype",
                                               "HER2-enriched subtype" , "Triple-negative subtype")))) %>% 
  mutate(effect = "Direct effect \n (accounted for\n childhood body size)")
  
dat_total <- dat_all2 %>% 
              mutate(effect = "\nTotal effect") %>% distinct()

dat_merged <- bind_rows(dat_total, dat_direct) %>%   
  mutate(exposure = factor(exposure, levels= rev(c("Dense area", "Non-dense area", "Percent density")))) %>% 
  mutate(effect = factor(effect, levels= c(  "\nTotal effect", "Direct effect \n (accounted for\n childhood body size)")))


# having to manually add labels as this is not a gris, but a facet plot
ab_labels <- tibble(
   label = c("a", "b"),
  x = c(.2, .1),
  y = c(.8, .4),
  #color = c("black",  "black" ),
  size = c(4,6))
  

pal<- c("#EB5291FF", "#FBBB68FF" ,"#F5BACFFF", "#9DDAF5FF", "#6351A0FF" ,"#FFEA5E", "#1794CEFF", "#972C8DFF") # yarrr::piratepal("pony")
p<-ggplot(dat_merged,
           aes(y=exposure, x=or,  colour=outcome, shape = effect_direction)) +
  geom_point(aes(size = effect), 
             position = ggstance::position_dodgev(height = 0.7))+
  geom_errorbarh(aes(xmin=or_lci95, xmax=or_uci95), height=0.1, 
            position = ggstance::position_dodgev(height = 0.7)) +
  theme_minimal_vgrid(8)+
  panel_border()+
  scale_color_manual(values=pal)+
  scale_shape_manual(values=c(16,1))+
  scale_size_manual(values=c(2.7, 2))+
  scale_x_log10(breaks = c(0,1,2,3))+
  theme(#strip.text.y = element_text(face = 'bold'),
        strip.text.x = element_text(face = "bold", size=8),
        axis.text.y = element_text(angle = 90, vjust=1, hjust = 0.5),
        panel.grid.major.y = element_blank(),
        panel.grid.major.x = element_line(),
        axis.title.y = element_text(face="bold"),
         panel.spacing = unit(1, "lines"),
        )+
  geom_vline(xintercept=1, linetype=3) +
  facet_grid(cols=vars(effect))+
  labs(color = "Outcome",y = "Exposure\n", x = "Odds ratio", shape="")+
  theme(legend.position = "bottom", legend.direction = "vertical", 
        plot.title.position = "plot",
        legend.title = element_text(face = "bold"),
        strip.text = element_text(size = 7))+
  guides(colour = guide_legend(reverse = TRUE, ncol=2), shape="none",  size="none")


source_data[["fig3"]] <- dat_merged


ggsave(paste0("figures_manuscript/figure3_poster.png"),
       plot=p, scale=1, 
       width=14, height=15,
       units=c("cm"), dpi=300, limitsize=F)

ggsave(paste0("figures_manuscript/figure3_paper.svg"),
       plot=p, #scale=1.4, 
       width=88, height=140,
       units=c("mm"), dpi=300, limitsize=F)


```




# Figure 4

This figure is in 07_instruments_heterogeneity_review script

# Figure 5

This figure is in 07_instruments_heterogeneity_review script



# Figure 6

```{r}
current_trait_category = "md_unadj"
# read cBMI direct effct after accounting for all MDs
direct_data <- 
  read_tsv( paste0(results_path, current_trait_category, "/mvmr/merged/merged_mvmr_earlyBMI-md_unadj-to-BCAC.tsv")) %>% 
  mutate(mediator = case_when(mediator=="dense_area_unadj" ~ "Dense area",
                              mediator=="nondense_area_unadj" ~ "Non-dense area",
                              mediator=="percent_density_unadj" ~ "Percent density")) %>% 
  filter(exposure == "Childhood body size") %>% 
  filter(outcome != "BRCA1_TN") %>% 
  mutate(exposure = paste0(exposure, "\n accounted for\n ", mediator)) %>% 
  select(exposure, outcome,  b, se, pval, starts_with("or", ignore.case = F)) %>% 
  mutate(outcome = case_when(outcome =="Breast cancer" ~ "BCAC 2017",
                             outcome =="Breast cancer (new)" ~ "BCAC 2020",
                             outcome =="ER+ breast cancer" ~ "ER+",
                             outcome =="ER- breast cancer" ~ "ER-",
                                           outcome =="Luminal A"    ~ "Luminal A",   
                                           outcome =="Luminal B, HER2+"   ~ "Luminal B1",  
                                           outcome =="Luminal B, HER2-" ~ "Luminal B2" ,
                                           outcome =="HER2-enriched"  ~ "HER2-enriched" ,
                                           outcome =="TN Basal"  ~ "TNBC" ))


# read total effects of  cBMI on all BC outcomes
total_data<- 
  read_csv(paste0(results_path, "univariableMR_BMI_BCACall.tsv")) %>% 
  filter(exposure == 'Childhood BMI') %>% 
  select(exposure, outcome, pval, b, se, starts_with("or", ignore.case = F)) %>% 
  mutate(exposure = "Childhood body size") %>% 
  filter(outcome != "TNBC_CIMBA_BRCA1") %>% 
  mutate(outcome = case_when(outcome =="Breast cancer" ~ "BCAC 2017",
                             outcome =="Breast cancer (overall BCAC 2020)" ~ "BCAC 2020",
                             outcome =="ER+ Breast cancer" ~ "ER+",
                             outcome =="ER- Breast cancer" ~ "ER-",
                                           outcome =="Luminal A"    ~ "Luminal A",   
                                           outcome =="Luminal B"   ~ "Luminal B1",  
                                           outcome =="Luminal B-HER2 negative" ~ "Luminal B2" ,
                                           outcome =="HER2 enriched"  ~ "HER2-enriched" ,
                                           outcome =="TNBC Basal-like"  ~ "TNBC" ))
  
              
all_data <- bind_rows(total_data, direct_data) %>% 
  filter(outcome != "BCAC 2020") %>% 
   mutate(outcome = case_when(outcome == "BCAC 2017" ~ "Overall breast cancer",
                              outcome == "ER+" ~ "ER+ breast cancer",
                              outcome == "ER-" ~ "ER- breast cancer",
                              outcome == "TNBC" ~ "Triple-negative subtype",
                              outcome == "Luminal A" ~ "Luminal A subtype",
                              outcome == "Luminal B1" ~ "Luminal B1 subtype",
                              outcome == "Luminal B2" ~ "Luminal B2 subtype",
                              outcome == "HER2-enriched" ~ "HER2-enriched subtype" )) %>% 
   mutate(outcome = factor(outcome, levels = rev(c("Overall breast cancer",  "ER+ breast cancer","ER- breast cancer"  ,
                                                "Luminal A subtype", "Luminal B1 subtype" , "Luminal B2 subtype",
                                               "HER2-enriched subtype" , "Triple-negative subtype")))) %>% 
  mutate(effect = ifelse(exposure == "Childhood body size", "total", "direct")) %>% 
  mutate(exposure = ifelse(exposure == "Childhood body size", "Childhood body size \n(total effect)", exposure)) %>% 
  mutate(exposure = factor(exposure, levels= rev(c(
    "Childhood body size \n(total effect)",   "Childhood body size\n accounted for\n Dense area",   
    "Childhood body size\n accounted for\n Non-dense area",  "Childhood body size\n accounted for\n Percent density")))) %>% 
  mutate(effect_direction = ifelse(or_lci95 > 1 & or_uci95 >= 1, 'ok',
                              ifelse(or_lci95 < 1 & or_uci95 <= 1, 'ok', 'overlaps null'))) 



pal<- c("#EB5291FF", "#FBBB68FF" ,"#F5BACFFF", "#9DDAF5FF", "#6351A0FF" ,"#FFEA5E", "#1794CEFF", "#972C8DFF") # yarrr::piratepal("pony")
p<-ggplot(all_data,
           aes(y=exposure, x=or,  colour=outcome, shape = effect_direction)) +
  geom_point(size=2.2,
             position = ggstance::position_dodgev(height = 0.7))+
  geom_errorbarh(aes(xmin=or_lci95, xmax=or_uci95), height=.1, 
            position = ggstance::position_dodgev(height = 0.7)) +
  scale_x_log10()+

  theme_minimal_vgrid(8)+
  geom_vline(xintercept=1, linetype='longdash') +
  scale_color_manual(values=pal)+
  scale_shape_manual(values=c(16,1))+
  theme(strip.text = element_text(face = 'bold'),
        strip.text.y = element_text(face = 'bold'),
        panel.grid.major.y = element_blank(),
        panel.grid.major.x = element_line(),
        axis.title.y = element_text(face="bold")
        )+
  panel_border()+
  labs(color = "Outcome",y = "Exposure\n", x = "Odds ratio", shape="")+
  theme(legend.position = "bottom", legend.direction = "vertical", 
        plot.title.position = "plot",   legend.title = element_text(face = "bold"))+
  guides(colour = guide_legend(reverse = TRUE, ncol=2), shape="none")

source_data[["fig6"]] <- all_data

ggsave(paste0("figures_manuscript/figure6_poster.png"),
       plot=p, scale=1, 
       width=10, height=12,
       units=c("cm"), dpi=300, limitsize=F)

ggsave(paste0("figures_manuscript/figure6_paper.svg"),
       plot=p, #scale=1.3, 
       width=88, height=140,
       units=c("mm"), dpi=300, limitsize=F)

```

```{r}
# save source data
names(source_data)
openxlsx::write.xlsx(source_data, file = 'figures_manuscript/source_data1.xlsx')

```


# extra

Comparing total and direct of cBMI + aaM


```{r}
library(rio)
paper <- import_list("../../Mini-project2/05_Writing/manuscript/final_version/Supplementary_Data_Vabistsevits2022.xlsx")

aaM_total<- paper$`Suppl Data 2` %>%
            filter(exposure == "Age at menarche (UKB)") %>% 
            filter(method == "Inverse variance weighted") %>% 
            mutate(effect = "Univariable MR")

cBMI_total <- paper$`Suppl Data 6` %>%
            filter(method == "Inverse variance weighted") %>% 
            filter(exposure == "Childhood body size") %>% 
            mutate(effect = "Univariable MR")

cBMI_direct <- paper$`Suppl Data 5` %>%
            filter(trait_category == "reproductive_traits") %>% 
            filter(accounted_for =="Age at menarche (UKB)") %>% 
            mutate(effect = "Multivariable MR")

aaM_direct <- paper$`Suppl Data 5` %>%
            filter(trait_category == "reproductive_traits") %>% 
            filter(exposure =="Age at menarche (UKB)") %>% 
            mutate(effect = "Multivariable MR")


res <- bind_rows(aaM_total, cBMI_total,cBMI_direct, aaM_direct) %>% 
  select(effect, exposure, outcome, or, or_lci95, or_uci95) %>% 
  mutate(effect = factor(effect, levels= c("Univariable MR", "Multivariable MR"))) %>% 
  mutate(OR_CI = paste0(round(or,2), " [", round(or_lci95,2), ":", round(or_uci95,2), "]")) 





 p<-  ggplot(res, aes(y=exposure, x=or, label=outcome, colour=effect)) +
  geom_errorbarh(aes(xmin=or_lci95, xmax=or_uci95), height=0.15) +
  geom_point(size=2)+
  #xlim(-1.2,1.2)+
  scale_color_manual(values=c("steelblue", "tomato"))+
  #scale_shape_manual(values = c(19)) +
  geom_text(aes(label=OR_CI),hjust=0.3, vjust=-1.2, size =2.5, color = '#333232')+
  geom_vline(xintercept=1, linetype=3,) +
  theme_minimal_hgrid(10, rel_small = 1) +
  facet_wrap(~effect, ncol=1)+
  labs(color = "",y = "", x = "")+
  theme(legend.position = "none", plot.title.position  = "plot")



```