-
Notifications
You must be signed in to change notification settings - Fork 0
/
GetGainEst.Rmd
78 lines (70 loc) · 2.54 KB
/
GetGainEst.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
---
title: "Estimate Genetic Gain"
author: "wolfemd"
date: "2020-2-13"
output: workflowr::wflow_html
editor_options:
chunk_output_type: inline
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = F, tidy = T)
```
# Objective
Given a selection index and the GEBV previously obtained [previously](IITA_StageII_Predict_C4.html), estimate genetic gain.
```{r}
library(tidyverse); library(magrittr)
gebvs<-read.csv("output/GEBV_IITA_OutliersRemovedTRUE_73119.csv", stringsAsFactors = F) %>%
select(-contains("_"))
unique(gebvs$GeneticGroup)
gebvs$GeneticGroup <-factor(gebvs$GeneticGroup,levels=c("GGetc","TMS13F","TMS14F","TMS15F","TMS18F"))
traits<-colnames(gebvs) %>% grep("GID|GeneticGroup",.,value = T, invert = T)
```
## Boxplot of GEBVs
```{r}
boxplotGenGain<-function(gebvs,traits){
# Input
# traits: vector of columns in df containing gebvs
# gebvs: an input dataframe with columns containing gebvs
## gebvs$GeneticGroup: one column should be named GeneticGroup...
## and be a factor with levels ordered
## in sequence by breeding cycle
gebvs_long<-gebvs %>%
tidyr::pivot_longer(cols=traits,names_to = "Trait",values_to = "GEBV")
gebvs_long %>%
ggplot2::ggplot(.,aes(x=GeneticGroup,y=GEBV, fill=GeneticGroup)) +
geom_boxplot() +
facet_wrap(~Trait,scales='free_y') +
theme(axis.text.x = element_text(face = 'bold',angle = 90),
legend.position = 'none') }
```
```{r}
boxplotGenGain(gebvs = gebvs,traits = traits)
```
## Barplot (Mean + SE) GEBVs
```{r}
barplotGenGain<-function(gebvs,traits){
# Input
# traits: vector of columns in df containing gebvs
# gebvs: an input dataframe with columns containing gebvs
## gebvs$GeneticGroup: one column should be named GeneticGroup...
## and be a factor with levels ordered
## in sequence by breeding cycle
gebvs_long<-gebvs %>%
tidyr::pivot_longer(cols=traits,names_to = "Trait",values_to = "GEBV")
gebvs_long %>%
group_by(Trait,GeneticGroup) %>%
summarize(meanGEBV=mean(GEBV),
stdErr=sd(GEBV)/sqrt(n()),
upperSE=meanGEBV+stdErr,
lowerSE=meanGEBV-stdErr) %>%
ggplot(.,aes(x=GeneticGroup,y=meanGEBV,fill=Trait)) +
geom_bar(stat = 'identity') +
geom_linerange(aes(ymax=upperSE,
ymin=lowerSE)) +
facet_wrap(~Trait,scales='free_y') +
theme(axis.text.x = element_text(face = 'bold',angle = 90),
legend.position = 'none') }
```
```{r}
barplotGenGain(gebvs = gebvs,traits = traits)
```