```{r setup, include=FALSE}
knitr::opts_chunk$set(
  results = 'asis',
  echo = FALSE,
  warning = FALSE,
  message = FALSE,
  fig.align = 'center'
)

#Load libraries
library(tidyverse)
library(kableExtra)
library(lubridate)
library(cowplot)
#test
```


# Introduction

  I have been an RPO at CSIRO for the last 10 years. I was involved in several projects with the objective of achieving higher yields in wheat through transferring desirable phenotypic or physiological traits. My current project is examining the impact of transferring chromosomes from rye into wheat to determine which characteristics can contribute to increased yield. My daily routine usually involves nurturing my plants in the glasshouse, or monitoring how they survive in the paddock. Hail and covid have therefore been destructive to project goals in 2020!
  
  My background is 'classical' plant breeding and I have worked as a breeder on a number of crops in the last 30 years ranging from wheat to cotton. This has been at government organisations and an international breeding company. The cotton company bred coventional and transgenic varieties. The San Remo pasta you eat is likely one of my durum varieties or their derivatives.   
  
  Before I attended Data School FOCUS I was an R raw recruit. My data collection revolved around manhandling Excel. 
  
# My Project

My project was to examine if incorporating alternative novel plant height reduction genes into commercial varieties of wheat improved yield? Varieties had their present dwarfing gene, either Rht 1 or Rht 2, replaced with either Rht 3_c26, Rht 18 or Rht 13. Furthermore, lines with an additional genetic marker for the gene, Rht 8, were developed. A range of traits which contribute to grain yield were measured in the field. I proposed to determine if there was a yield response to the height gene/s incorporated? From the traits measured it was anticipated to identify those which contributed to yield. The following table demonstrates a sample of the data collected.     

```{r read data}

beauties2017 <- read_csv('data/beauties2017.csv') #loading your data file
# head(beauties2017)

# determine mean of reps

beauties_summary <- beauties2017 %>% 
  select(entry, 10:19) %>% 
  group_by(entry) %>% 
  summarise_all(mean, na.rm = TRUE) # remove NA to estimate mean of reps 

beauties_raw <- beauties2017 %>% 
  select(-c(10:19), -REP, -row, -range) %>% 
  left_join(beauties_summary, by="entry") %>% 
  distinct()

```



```{r mytable, out.width='100%', echo = F}
knitr::kable(head(beauties_summary, n = 5), format = "html", caption = "Example of data measured for each trait") %>% 
  kable_styling("striped")

```


## Preliminary results

The overall objective of the project was to see if different dwarfing genes (QTLs) had an impact on grain yield, and if so, what combination of genes was most desirable? Yield is the driving force for which varieties are adopted by farmers. 

```{r group gene combos for main varieties and traits}

#gather and spread
#tidy is a rectangular data - 2 dimensions i.e. rows and columns
#if, say, traits names are column names - wide format data - as you add
#new variable, the data becomes wider
#long format is in its simplicity a pair-wise data
#for your data you will have 72x9 number of rows
#spread makes data wider
#gather makes data long

#functions - anything that follows () is a function ##generalisation - CAREFUL

beauties_raw <- beauties_raw %>% 
  select(genotype, rhtmajor, rhtminor, 9:18) %>% #select variable columns with data
  group_by(genotype, rhtmajor, rhtminor) %>% # group rht gene combo
  summarise_all(mean) # calculate mean for each variety's gene combo
  
beauties_raw <- beauties_raw %>% 
  mutate(group = paste(rhtmajor, rhtminor, sep = "_")) #groups gene combo label

```


```{r yield by QTL}

yield <- beauties_raw %>% 
  filter(!genotype %in% c("C6", "Sceptre", "Shield")) %>% # drop these var
  ggplot() +
  geom_boxplot(aes(x = reorder(group, desc(yieldkgm2)), y = yieldkgm2))+ # reorder in descending yield
  #facet_wrap(~genotype, scales = "free_x", nrow = 1)+
   theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 8)) +
   labs(title="Plot yield response to height QTL combination", x = "QTL combination", y = "Grain Yield (kg)")

yield


```

```{r biomass by QTL not shewn}

biomass <- beauties_raw %>% 
  filter(!genotype %in% c("C6", "Sceptre", "Shield")) %>% # drop these var
  ggplot() +
  geom_boxplot(aes(x = reorder(group, desc(bm)), y = bm))+
  #facet_wrap(~genotype, scales = "free_x", nrow = 1)+
   theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 8)) +
   labs(title="Biomass Z89 response to height QTL combination", x = "QTL combination", y = "Biomass (g)")

# biomass


```



As plant development stage, particularly time of heading or flowering, influences yield I wanted a visual to see if different QTL combinations had an effect on growth rate for  varieties.

```{r plot, fig.height= 4, fig.width=8 }
# fig.cap="Heading time"

headplot <- beauties_raw %>% 
  filter(!genotype %in% c("C6", "Sceptre", "Shield")) %>% # drop these var
  pivot_longer(cols = c(Z20170920, Z20170926), names_to = "Z_date", values_to = "Z_score") %>% # make long 
  mutate(Z_date=str_sub(string = Z_date, 2, -1)) %>% # -1 is last character
  mutate(Z_date=ymd(Z_date)) %>% # make it a real date
  ggplot(aes(x = group, y = Z_score)) +
  geom_point(aes(colour= as.factor(Z_date)))+
  facet_wrap(~genotype, scales = "free_x", nrow = 1)+
   theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 8)) +
   labs(title="Heading score on two dates for varieties and their QTLs", x = "QTL combination", y = "Zadoks Score", colour = "Date") +
  theme(legend.position = "bottom")
# +
#   theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 5),
#         plot.title = (element_text(size=10, hjust = 0.5)),
#         panel.grid.minor.x  = element_blank(),
#         axis.text.y = element_text(size = 5),
#          axis.title = element_text(size = 8))

  

headplot

```



```{r yield for genotype by QTL}

yield_gen <- beauties_raw %>% 
  filter(!genotype %in% c("C6", "Sceptre", "Shield")) %>% # drop these var
  ggplot() +
  geom_point(aes(x = group, y = yieldkgm2))+
  facet_wrap(~genotype, scales = "free_x", nrow = 1)+
   theme_bw() +
  theme(axis.text.x = element_blank()) +
   labs(title="Plot yield", x="", y = "Grain Yield")

# yield_gen


```



```{r HI by QTL}

hi <- beauties_raw %>% 
  filter(!genotype %in% c("C6", "Sceptre", "Shield")) %>% # drop these var
  ggplot(aes(x = group, y = HI)) +
  geom_point()+
  facet_wrap(~genotype, scales = "free_x", nrow = 1)+
   theme_bw() +
  theme(axis.text.x = element_blank()) +
   labs(title="Harvest Index", x = "", y = "HI")

# hi


```


```{r}
tiller <- beauties_raw %>% 
  filter(!genotype %in% c("C6", "Sceptre", "Shield")) %>% # drop these var
  ggplot(aes(x = group, y = tillerperm2)) +
  geom_point()+
  facet_wrap(~genotype, scales = "free_x", nrow = 1)+
   theme_bw() +
  theme(axis.text.x = element_blank()) +
   labs(title="Tillering", x = "", y = "Tillers m2")

# tiller
```

```{r}

earl <- beauties_raw %>% 
  filter(!genotype %in% c("C6", "Sceptre", "Shield")) %>% # drop these var
  ggplot(aes(x = group, y = earlgthcm)) +
  geom_point()+
  facet_wrap(~genotype, scales = "free_x", nrow = 1)+
   theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 8)) +
   labs(title="Ear size", x = "QTL combination", y = "Ear length cm")

# earl
```

```{r}

rnn <- beauties_raw %>% 
  filter(!genotype %in% c("C6", "Sceptre", "Shield")) %>% # drop these var
  ggplot(aes(x = group, y = RNN)) +
  geom_point()+
  facet_wrap(~genotype, scales = "free_x", nrow = 1)+
   theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 8)) +
   labs(title="Rachis number", x = "QTL combination", y = "RNN")

# rnn


```

```{r plot8}

klwt <- beauties_raw %>% 
  filter(!genotype %in% c("C6", "Sceptre", "Shield")) %>% # drop these var
  ggplot(aes(x = group, y = kwt)) +
  geom_point()+
  facet_wrap(~genotype, scales = "free_x", nrow = 1)+
   theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 8)) +
   labs(title="Grain Weight", x = "QTL combination", y = "200 kwt")

# klwt


```

In assessing yield for the gene combinations in different varieties, it is of interest to examine the yield components to determine which, if any, have a major influence on yield, or whether it may be a cumulative response. If traits can be used to predict yield, they can be used as a selection tool in early generations when yield plots are not viable.


```{r to do cowplot}

plot_grid(yield_gen, hi, tiller, earl, rnn, klwt, rel_heights = c(.9,1))
          

```



```{r regression ear length and tillers}

corr_ear_till <- ggplot(beauties2017, aes(x= tillerperm2, y = earlgthcm)) +
  geom_point(aes(colour = as.factor(genotype), size=kwt)) +
  geom_smooth(method = "lm")+ # linear model, grey line is sd
  theme_classic() +
  labs(title="Ear length as response to tillers per m2", x= "Tillers per m2", y = "Ear length cm", colour = "Variety") + # add these in
   theme(legend.text = element_text(size=5), # reduce size of text in legend box
        legend.title = element_text(size = 7), # reduce size of legend title
        legend.key.size = unit(0.2,"cm"), # reduce size of legend box
        legend.box = "horizontal", # align boxes horizontally
        legend.position = "bottom",
       # legend.position = c(.65,0.8), # grid position within graph
        legend.direction = "horizontal") # key in horizontal direction

# corr_ear_till

```


```{r reg seed wt and yield by genotype}

reg_sdwt_yld <- ggplot(beauties2017, aes(x= kwt, y =  yieldkgm2)) +
  geom_point(aes(colour = as.factor(genotype), size = RNN)) +
  geom_smooth(method = "lm")+ # linear model, grey line is sd
  theme_classic() +
  labs(title="Kernel wt and yield", x = "Kernel weight g", y = "Yield (kgm2)", colour = "Variety") +
   theme(legend.text = element_text(size=5),
        legend.title = element_text(size = 7),
        legend.key.size = unit(0.2,"cm"),
        legend.box = "horizontal", # align boxes horizontally
        legend.position = "bottom",
        legend.direction = "horizontal") # key in horizontal direction



# reg_sdwt_yld

```

```{r reg biomass and yield by genotype}

reg_bm_yld <- ggplot(beauties_raw, aes(x= bm, y =  yieldkgm2)) +
  geom_point(aes(colour = as.factor(group), size = htcm)) +
  geom_smooth(method = "lm")+ # linear model, grey line is sd
  theme_classic() +
  labs(title="Biomass and yield", x = "Biomass at Z89 (g)", y = "Yield (kgm2)", colour = "QTL") +
  theme(legend.text = element_text(size=5),
        legend.title = element_text(size = 7),
        legend.key.size = unit(0.2,"cm"),
        legend.box = "horizontal", # align boxes horizontally
        legend.position = "bottom",
        legend.direction = "horizontal") # key in horizontal direction


# reg_bm_yld

```



```{r reg seed wt and yield by QTL}

corr_sdwt_yld <- ggplot(beauties_raw, aes(x= kwt, y =  yieldkgm2)) +
  geom_point(aes(colour = as.factor(group), size = RNN)) +
  geom_smooth(method = "lm")+ # linear model, grey line is sd
  theme_classic() +
  labs(title="Kernel wt and yield", x = "Kernel weight g", y = "Yield (kgm2)", colour = "QTL") +
   theme(legend.text = element_text(size=5),
        legend.title = element_text(size = 7),
        legend.key.size = unit(0.2,"cm"),
        legend.box = "horizontal", # align boxes horizontally
        legend.position = "bottom",
        legend.direction = "horizontal") # key in horizontal direction


# corr_sdwt_yld

```

Furthermore, where there are relationships between a trait and yield or between traits outliers may identify potential sources of genetic gain. Outliers can become subject to physiological investigations to identify what is contributing to their response.


```{r cowplot of regressions, out.width='100%', fig.align='center', fig.height= 4, fig.width=8}

plot_grid(reg_sdwt_yld, corr_sdwt_yld)

```

```{r, out.width='100%', fig.align='center', fig.height= 4, fig.width=8}

plot_grid(corr_ear_till, reg_bm_yld)

```



# My Digital Toolbox

![](resources/img/tidyverse.png){width=100px} ![](resources/img/lubridate.jpg){width=100px}
 ![](resources/img/Rmarkdown.png){width=100px}  ![](resources/img/ggplot.png){width=100px}

## Favourite tool

![](resources/img/cowplot.png){width=100px}

The cowplot Longhorn icon is iconic in agriculture, and the tool allowed me to incorporate a large number of graphs relating to different traits into a confined area. This meant I did not have to omit the graphs where I had taken the results into consideration.

# My time went ...

My initial time was spent trying to write code in the right structure using functions. But then I found I spent a larger proportion of my time playing with the aesthetics of plots.

# Next steps

I would add the data taken in subsequent years to identify interactions, and plot a PCA analysis.

# My Data School Experience

The program has been a great opportunity to expand my networking throughout the organisation. 

The experience has given me a foundation to coding and has encouraged me to continue using into the future. It has given me a second language for which I can converse at work.

It has also encouraged me to produce tidy data!


