FA_drawing.Rmd

---
title: "GBLUP_h2"
author: "Ye Bi"
date: "`r Sys.Date()`"
output: html_document
---

```{r}
library(tidyverse)
library(ggplot2)
library(plyr)
library(ggpubr)
library(ggrepel)
library(tidyverse)
library(stargazer)
library(patchwork)
library(ggplot2)
```

## Drawing plots for h2, varE, varG
```{r}
source("../Functions/functions.R")
met_names = met_name_func(mode = 'name')
```

##Loading h2 from local files.
```{r}
load_h2 = function(prepro, gblup){
  load(paste0("../../temp/", prepro,"/", gblup, "/Control_",gblup,"_h2.rda" ))
  h2varcon <- as.data.frame(do.call(rbind, h2L))
  # h2varcon = h2varcon[-30,] #remove d3 from met. 
  load(paste0("../../temp/", prepro,"/", gblup, "/Stress_",gblup,"_h2.rda" ))
  h2vartrt <- as.data.frame(do.call(rbind, h2L))
  # h2vartrt = h2vartrt[-30,]
  h2 = c(h2varcon$h2, h2vartrt$h2)
  varE = c(h2varcon$varE, h2vartrt$varE)
  varG = c(h2varcon$varG, h2vartrt$varG)
  df = data.frame(Met=met_names, 
                  h2=h2, 
                  varE = varE, 
                  varG = varG, 
                  Treatment = rep(c("Control", "Stress"), each=66))
  
  df.long <- reshape2::melt(df, id.vars = c('Met', 'Treatment'))
  
  
  h2.df.long <- df.long %>% filter(variable == 'h2') %>% droplevels()
  colnames(h2.df.long) <- c("Met", "Treatment", "variable", "value")
  h2.df.long$Treatment <- as.factor(h2.df.long$Treatment)
  h2.df.long$Met <- factor(h2.df.long$Met, levels=unique(h2.df.long$Met))
  h2.df.long$Kernel <- paste0(prepro, "-G-", gblup)
  return(h2.df.long)
}

h22L = list()
m= 1
for(i in c("log_lmer")){
  for(j in c("sommer")){
        h22L[[m]] = load_h2(i, j)
        m = m+1
  }
        }
```



## h2 density plot --- FA
```{r}
####functions.R
FA_h2_den=FA_density_plot_generator(df=h22L[[1]],
                          my_colors = RColorBrewer::brewer.pal(11, 'Spectral')[c(11,10,9,8,7,5,4,3,2,1)],
                          limits_range = c(0,1),
                          breaks_range = seq(0,1,0.2),
                          x_name = "Heritability", 
                          ylim = NULL)
FA_h2_den
# dev.print(pdf, file=paste0("../temp/FA_h2_density.pdf"), height=4, width=8)
```