```{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(cowplot)
library(ggrepel)

```


# Introduction
I have been working as a molecular biologist for over 25 years on both animal and plant systems. The last 10 years I have worked at CSIRO and have been involved with several projects investigating different aspects of wheat physiology, growth and development. This has included carrying out laboratory, glasshouse and field experiments and previously data generated has been entered and analysed using Excel. Prior to Data School I did not know how to code or use R. Data School FOCUS has opened my eyes to a whole new world of data analysis and I hope to use this as a starting point to learn and apply new and exciting methodologies. 

# My Project
I work on a project within the GRDC's National Phenology Initiative which aims to predict the flowering time of wheat cultivars in many different growing regions within Australia.

Wheat cultivars have an optimal flowering window. Crops that flower too early can have lower yield due to insufficient biomass accummulation or exposure to cold or frost events. Conversely, crops that flower too late risk being exposed to water stress or heat events which can negatively impact yield. Growers require accurate information to select the correct cultivar and sowing date for their farming conditions. For new cultivars it takes several years conducting field trials to accumulate data to predict flowering time. With new cultivars being released it is essential that growers have accurate information by the time of release.

The major environmental factors influencing flowering time are thermal time, photoperiod and vernalisation. APSIM models have been developed that use parameters based on these factors to model cultivar flowering times across the many cereal growing regions within Australia.  This project aims to improve and modify the existing APSIM models of wheat through derivation of parameters in controlled conditions and association with genomic data to more accurately predict flowering time in the field. This study involved four controlled temperature experiments being carried out on 54 Australian wheat cultivars and 15 Wheat NILs (Near Isogenic Lines). 



## Preliminary results

The controlled experiments were carried out under four environmental conditions, SN (short days, no vernalisation), LN (long days, no vernalisation), SV (short days plus vernalisation) and LV (long days plus vernalisation) to determine the influence of photoperiod and vernalisation on flowering times. Vernalisation was carried out by imbiding seeds at 4<sup>o</sup>C for 8 weeks prior to planting. The temperature for all experiments was set at 22<sup>o</sup>C and measured every 30 minutes using a TinyTag data logger.  Traits measured include emergence date, flowering date, heading date, final leaf number and spikelet number. In addition Haun stage, a measure of developmental growth stage based on leaf emergence, was recorded every third day from emergence to flowering. The data was initially entered into several sheets within one Excel file. I have used the tidyverse package to clean and arrange the data. I have separated the data into two files, Haun_temp, containing the haun stage scores and a second file all_data_wide, containing all other traits.  I used the temperature data to convert dates to accumulated thermal time (degree-days). Tables 1 and 2 show a representation of the two tidied data files.


**Tables**

``````{r mytable, out.width='80%', echo = FALSE}

#read files

all_data_wide <- read_csv("data/analysed_data/all_data_wide.csv") # all data wide

cult_desc <- read_csv("data/rawdata/cultivar_descriptions.csv") # cultivar description

haun_data <- read_csv("data/analysed_data/Haun_temp.csv") # haun data

flw_hd_data <- read_csv("data/analysed_data/fwr_hd_tt.csv") # long form of flowering  and heading data

#join files

haun_data2 <- cult_desc %>% 
  right_join(haun_data, by=c("genotype")) %>% 
rename(Genotype="genotype", Environment="environment", 
         Rep="reps", 
         'degree-days'="tt_haun",
         'Haun stage'="haun",
       Type= "type", 
       Maturing="speed") %>%  
  filter(Genotype=="EMU_ROCK",Rep== 1) %>% 
  mutate_if(is.numeric, round, digits=1)

knitr::kable(head(haun_data2, n = 5), format = "html", caption = " Controlled Environment - Haun Stage Data") %>% 
  kable_styling("striped", full_width = F)

```

```{r secondtable, out.width='80%', echo = FALSE}

  all_data_wide2 <- cult_desc %>% 
  right_join(all_data_wide, by=c("genotype")) %>% 
  rename(Genotype="genotype", Environment="environment",
         Rep="reps", Type= "type", Maturing="speed") %>% 
    mutate_if(is.numeric, round, digits=0)

knitr::kable(head(all_data_wide2, n = 5), format = "html", caption = " Controlled Environment - Flowering Time Data") %>% 
  kable_styling("striped")%>%
  kableExtra::scroll_box(width = "100%")
```




**Plots**

**Figure 1 : Flowering time for plants grown in four different environments**

```{r standard-plot, echo=FALSE, fig.align='center', fig.height=4, fig.width=6, out.width='80%'}


cult_desc %>% 
  right_join(flw_hd_data, by=c("genotype")) %>% 
  rename(Genotype="genotype", Environment="environment",
         Rep="reps", Type= "type", Maturing="speed")  %>% 
  unite(mat_type, c(Maturing, Type), sep = " ", remove = FALSE) %>% 
  mutate( (mat_type = as.factor(mat_type))) %>% 
  mutate(mat_type=fct_relevel(mat_type, c("Slow Spring", "Mid Spring", "Fast Spring",
                                          "Slow Winter", "Mid Winter", 
                                          "Fast Winter"))) %>% 
  filter(obs_type=="half_fwr") %>% 
  group_by(Genotype, Environment, mat_type) %>% 
  summarise(mean_fwr = mean(thermaltime)) %>% 
  filter(mean_fwr != "NA") %>% 
ggplot(mapping=aes(x=Genotype, y=mean_fwr, colour=Environment))+
  geom_point(size=1)+ 
  facet_wrap(~fct_relevel(mat_type), scales="free_x", nrow=1)+
  scale_colour_manual(values=c("SV" = "#e41a1c", "SN" = "#984ea3", "LN" = "#4daf4a", "LV" = "#377eb8"), aesthetics = c("colour", "fill"))+
  labs(title="", x = "Genotypes", y = "Mean Flowering Time (degree-days)", 
       colour="Environment", size = 4)+
  theme_bw()+
  theme(axis.text.x = element_blank(),
        plot.title = element_text(size=10),
        legend.position="bottom",
        panel.grid.minor.x  = element_blank(),
        panel.grid.major.x  = element_blank())+
  ylim(0,5000)

```
Figure 1. Flowering time results for all four environments.

Plants usually flower quickest in the LV environment, where all vernalisation requirements are met and the day length is long. Shortening the day length or removing vernalisation can delay flowering time, but there is a genotypic effect. The flowering times of spring lines are least affected by vernalisation, whilst slow winter lines are most affected. The flowering times of fast lines are least affected by photoperiod, whilst slow lines are most affected. 
    
    
    
**Figure 2 : Vernalisation effect on flowering time**

```{r second_plot, echo=FALSE, fig.align='center', fig.height=4, fig.width=6, out.width='80%'}



#vern effect

flw_hd_data_vern <-  flw_hd_data %>% 
  filter(obs_type=="half_fwr") %>% 
  filter(environment=="LV"|environment=="LN" ) %>% 
  group_by(genotype, environment) %>%
  summarise(mean_tt=mean(thermaltime)) %>% 
  filter(mean_tt!="NA") %>% 
  spread(key=environment, value=mean_tt)
  
cult_sel <- flw_hd_data_vern %>% 
  filter(genotype=="AXE" |genotype=="FORREST"|genotype=="LONGSWORD"|genotype=="KITTYHAWK")

ggplot()+
  geom_point(data=flw_hd_data_vern,mapping=aes(x=reorder(genotype, LN), y=LN,colour="#4daf4a"),size=1)+
  geom_point(data=flw_hd_data_vern,mapping=aes(x=reorder(genotype, LN),y=LV, colour="#377eb8"),size=1)+
  geom_label_repel(data=cult_sel,mapping=aes(x=reorder(genotype, LN),y=LV,  label=genotype), colour="black", angle=90, vjust=3,
             size=1.5)+
    labs(y = "Mean Flowering Time, degree-days", colour= "Environment", size = 4)+
  scale_color_identity(name = "Environment",
                       breaks = c("#4daf4a", "#377eb8"),
                       labels = c("LN", "LV"),
                       guide = "legend")+
  theme_bw()+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 5),
          plot.title = element_blank(),
          axis.title.x= element_blank(),
          panel.grid.minor.x  = element_blank(),
          panel.grid.major.x  = element_blank(),
          legend.position = "bottom",
          axis.text.y = element_text(size = 5))+
  ylim(0,4000)



   
```
Figure 2. The effect of vernalisation on flowering time of cultivars.

Spring wheats, Axe and Forrest, have no vernalisation requirement and therefore flower at approximately the same time in both environments. Winter wheats, Longsword and Kittyhawk, require vernalisation, and without vernalisation flowering is delayed.



**Figure 3 : Photoperiod effect on flowering time of cultivars**

```{r third_plot, echo=FALSE, fig.align='center', fig.height=4, fig.width=6, out.width='80%'}

#Pp effect

flw_hd_data_pp <-  flw_hd_data %>% 
  filter(obs_type=="half_fwr") %>% 
  filter(environment=="SV"|environment=="LV" ) %>% 
  group_by(genotype, environment) %>%
  summarise(mean_tt=mean(thermaltime)) %>% 
  filter(mean_tt!="NA") %>% 
  spread(key=environment,value=mean_tt) 
  
ggplot()+
    geom_point(data=flw_hd_data_pp,mapping=aes(x=reorder(genotype, SV), y=LV,colour="#377eb8"),size=1)+
    geom_point(data=flw_hd_data_pp,mapping=aes(x=reorder(genotype, SV),y=SV, colour="#e41a1c"),size=1)+
    geom_label_repel(data=cult_sel,mapping=aes(x=reorder(genotype, LN),y=LV,  label=genotype), colour="black",
             size=1.5, angle=90, vjust=3)+
    labs(title="Photoperiod effect on Flowering time", 
         y = "Mean Flowering Time, degree-days", colour= "Environment", size = 6)+
    scale_color_identity(name = "Environment",
                         breaks = c("#377eb8", "#e41a1c"),
                         labels = c("LV", "SV"),
                         guide = "legend")+
    theme_bw()+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 5),
          plot.title = (element_text(size=10, hjust = 0.5)),
          axis.title.x= element_blank(),
          panel.grid.minor.x  = element_blank(),
          panel.grid.major.x  = element_blank(),
          legend.position = "bottom",
          axis.text.y = element_text(size = 5))+
  ylim(0,4000)


```

Figure 3. The effect of photoperiod on cultivars. 

There is a genotyptic response in daylength. Fast maturing cultivars, Axe and Longsword, are photoperiod insensitive, the effect of changing day length on flowering time is small. Conversely, slow cultivars, Forrest and Kittyhawk, are photoperiod sensitive and shorter day length delays flowering. Photoperiod sensitivity and vernalisation requirements have separate genetic controls. The Australian cultivars used in this experiment display a range of  vernalisation and photoperiod sensitivities. 



**Figure 4 : Vernalisation requirement and photoperiod sensitivity of wheat cultivars**

```{r r fourth_plot, echo=FALSE, fig.align='center', fig.height=4, fig.width=6, out.width='80%'}


  
     flw_data2 <- flw_hd_data %>% 
  filter(obs_type=="half_fwr") %>% 
  spread(key=environment,value=thermaltime) %>% 
    group_by(genotype) %>% 
    summarise(vern_effect=(mean(LN-LV)), pp_effect=(mean(SV-LV))) %>%  
left_join(cult_desc, by=c("genotype"))  %>% 
  unite(mat_type, c(speed, type), sep = " ", remove = FALSE) 
  
cult_sel2 <- flw_data2 %>% 
  filter(genotype=="AXE" |genotype=="FORREST"|genotype=="LONGSWORD"|genotype=="KITTYHAWK") 

  ggplot(data=flw_data2, mapping=aes(x=vern_effect, y=pp_effect, shape=type, colour=speed))+
  geom_point(size=1.5)+
  geom_label_repel(data=cult_sel2, mapping=aes(x=vern_effect, y=pp_effect, label=genotype), colour="black",
                   size=1.5, angle=90, vjust=-2, fill="white")+
  labs(title="Vernalisation requirement and photoperiod sensitivity of Wheat cultivars", x="Vernalisation effect", y="Photoperiod effect", colour= "Flowering rate", shape= "Cultivar Type", size = 4)+
  scale_shape_manual(values=c(16,17)) +
  theme_bw()+
  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))

```

Figure 4. The vernalisation requirement and photoperiod sensitivity of the cultivars studied. 

Vernalisation effect is calculated as the difference between mean flowering time in LV and LN. Photoperiod sensitivity is calculated as the difference between mean flowering time in LV and SV.
 
 
 
 **Figure 5: Correlation of flowering time with the traits measured**
 
```{r standard-plot4, echo=FALSE, fig.align='center', fig.height=4, fig.width=6, out.width='80%'}

#Correlations with flowering

plot1 <- ggplot(all_data_wide, mapping=aes(x=half_fwr, y=tt_final_leaf))+
  geom_point(alpha=0.4, colour="grey49")+
  geom_smooth(method=lm, colour="blue")+
  labs(title="Relationship of flowering time and final leaf emergence", x = "Flowering (degree-days)", y = "Final leaf emergence (degree-days)", size = 2)+
  theme_bw()+
  theme(plot.title = element_text(size=6, hjust = 0.5),
        panel.grid.minor.x  = element_blank(),
        axis.title = element_text(size = 6),
        axis.text = element_text(size = 6))

plot2 <- ggplot(all_data_wide, mapping=aes(x=half_fwr, y=half_hd))+
  geom_point(alpha=0.4, colour="grey49")+
  geom_smooth(method=lm, colour="blue")+
  labs(title="Relationship of flowering time and heading time", x = "Flowering (degree-days)", y = "Heading (degree-days)", size = 2)+
  theme_bw()+
  theme(plot.title = element_text(size=6, hjust = 0.5),
        panel.grid.minor.x  = element_blank(),
        axis.title = element_text(size = 6),
        axis.text = element_text(size = 6))

plot3 <- all_data_wide %>% 
  ggplot( mapping=aes(x=half_fwr, y=spikelet_no))+
  geom_point(alpha=0.4, colour="grey49")+
  geom_smooth(method=lm, colour="blue")+
  labs(title="Relationship of flowering time versus spikelet number", x = "Flowering (degree-days)", y = "Spikelet number", size = 2)+
  theme_bw()+
  theme(plot.title = element_text(size=6, hjust = 0.5),
        panel.grid.minor.x  = element_blank(),
        panel.grid.major.x  = element_blank(),
        axis.title = element_text(size = 6),
        axis.text = element_text(size = 6))

plot4 <- all_data_wide  %>%
  ggplot( mapping=aes(x=half_fwr, y=final_leaf))+
  geom_point(alpha=0.4, colour="grey49")+
  geom_smooth(method=lm, colour="blue")+
  labs(title="Relationship of flowering time versus final leaf number", x = "Flowering (degree-days)", y = "Final leaf number", size = 2)+
  theme_bw()+
  theme(plot.title = element_text(size=6, hjust = 0.5),
        panel.grid.minor.x  = element_blank(),
        axis.title = element_text(size = 6),
        axis.text = element_text(size = 6))


plot_grid(plot1, plot4, plot3, plot2)

```

Figure 5. Correlation of flowering time with other traits measured


**Calculating APSIM Parameters **

Seven parameters were calculated from the original data. Descriptions of the seven parameters are shown in Table 3. The parameters were calculated using R scripts and a representation of the output is in Table 4.

```{r thirdtable, out.width='100%', echo = FALSE}

#Parameters
parameter_list <- data.frame("Parameter"=c("Parameter_1","Parameter_2","Parameter_3","Parameter_4","Parameter_5","Parameter_6","Parameter_7"),
                             "Name"=c("Minimal_Leaf_Number","Pp_Sensitivity","Vrn_Sensitivity","Base_Phyllochron",
                                      "Phyllochron_Photoperiod_effect","Early_Reproductive_Long_Day_Base", "Early_Reproductive_Pp_Sensitivity"),
                             "Calculation"=c("LV mean final leaf number","SV mean final leaf number minus LV mean final leaf number",
                                             "LN mean final leaf number minus LV mean final leaf number",
                                             "Slope of accummulated thermal time against mean leaf number between 3 and 7 in LV",
                                             "Ratio of Phyllochron SV to Phyllochron LV","Accumulated thermal time at flowering divided by BasePhyllochron",
                                             "Difference of accumulated time between LV and SV treatment divided by BasePhyllochron"))

knitr::kable(head(parameter_list, n = 7), format = "html", caption = " APSIM Parameters description") %>% 
  kable_styling("striped")

```


  
```{r parameter_calculations, out.width='100%', echo = FALSE }


#Upload files 
final_leaf <- read_csv("data/analysed_data/final_leaf.csv")
haun_tt <- read_csv("data/analysed_data/Haun_temp.csv")
spikelet_data <- read_csv("data/analysed_data/spikelet_data.csv")
flw_hd_data <- read_csv("data/analysed_data/fwr_hd_tt.csv")
all_data <- read_csv("data/analysed_data/all_data_wide.csv")
mean_final_leaf <- read_csv("data/analysed_data/mean_final_leaf.csv")
final_leaf_tt <- read_csv("data/analysed_data/final_leaf_tt.csv")





#####################Parameter 1 #######################

# MinimumLeafNumber defined as final leaf number LV Environment

Parameter_1 <- filter(final_leaf, environment=="LV") %>% 
  group_by(genotype) %>% 
  summarise(Minimum_leaf_number=mean(final_leaf))


###############Parameter 2#################################  

#PPsensitivity defined as difference in leaf number between LV and SV

Parameter_2 <- filter(final_leaf,environment=="LV"  |environment=="SV")%>% 
  spread(environment,value=final_leaf) %>% 
  rename("LV_final_leaf"="LV", "SV_final_leaf"="SV") %>% 
  filter(LV_final_leaf!="NA") %>% 
  filter(SV_final_leaf!="NA") %>% 
  mutate(Ppsensitivity=SV_final_leaf - LV_final_leaf)%>%
  group_by(genotype) %>%
  summarise(Ppsensitivity=mean(Ppsensitivity)) %>% 
  mutate(Ppsensitivity=ifelse(Ppsensitivity>0,Ppsensitivity,0))


######################Parameter 3############################

#VrnSensitivity defined as differnece between final leaf number LV and LN


Parameter_3 <- filter(final_leaf,environment=="LV"  |environment=="LN") %>% 
  spread(environment,value=final_leaf) %>% 
  rename("LV_final_leaf"="LV", "LN_final_leaf"="LN") %>% 
  filter(LN_final_leaf!="NA") %>%  
  mutate(vrnsensitivity= LN_final_leaf - LV_final_leaf)%>%
  group_by(genotype) %>% 
  summarise(Vrn_sensitivity=mean(vrnsensitivity)) 



##########################Parameter_4##############################################



Phyllochron <-  haun_tt %>% filter(haun>3, haun<7) %>% 
  group_by(genotype, environment) %>% 
  do((lm(tt_haun ~ haun, data= .)) %>% 
       coef %>% 
       as_data_frame()) %>% 
  mutate(sel=1:n()) %>%
  filter(sel==2) %>%
  rename(phyllochron=value) %>% 
  select(1:3)



Parameter_4 <- Phyllochron %>% 
  filter(environment=="LV") %>% 
  select(1,3) %>% 
  rename(Base_Phyllochron="phyllochron")
  


  ############Parameter 5###################
  
  Parameter_5 <- Phyllochron %>% 
  filter(environment=="LV" | environment=="SV") %>% 
  spread(environment, phyllochron) %>% 
  summarise(Phyllochron_PhotoPeriod_Effect=SV/LV)


##########Parameter 6######################


  flw_half <- flw_hd_data %>% 
        filter(obs_type=="half_fwr")
  
tt_flag_flower <- flw_half %>% 
  full_join(final_leaf_tt,by=c("genotype", "environment", "reps")) %>% 
    filter(environment=="LV") %>% 
  mutate(tt_flag_flower=thermaltime -tt_final_leaf) %>% 
  group_by(genotype) %>% 
  summarise(mean_tt_ff=mean(tt_flag_flower))
 
  
  Parameter_6 <- tt_flag_flower%>% 
    full_join(Parameter_4, by=c("genotype")) %>% 
   mutate(Early_Reproductive_Long_Day_Base=mean_tt_ff / Base_Phyllochron) %>% 
    select(1,5)
  


###Parameter7#############

tt_flag_flower_pp <- flw_half %>% 
  full_join(final_leaf_tt,by=c("genotype", "environment", "reps")) %>% 
  filter(environment=="LV"| environment=="SV") %>% 
  mutate(tt_flag_flower=thermaltime -tt_final_leaf) %>%
  select(1,2,3,8) %>% 
  spread(environment,tt_flag_flower) %>% 
  mutate(tt_flag_flower_pp=SV -LV)%>% 
    filter(tt_flag_flower_pp!="NA") %>% 
  group_by(genotype) %>% 
  summarise(mean_tt_ff_pp=mean(tt_flag_flower_pp))
  
Parameter_7 <- tt_flag_flower_pp%>% 
  full_join(Parameter_4, by=c("genotype")) %>% 
  mutate(EarlyReproductivePpSensitivity=mean_tt_ff_pp / Base_Phyllochron) %>% 
  select(1,5)

```

```{r fourthtable, out.width='100%', echo = FALSE}
#One parameter file

Parameters <- Parameter_1 %>% 
  full_join(Parameter_2) %>% 
  full_join(Parameter_3) %>% 
  full_join(Parameter_4) %>% 
  full_join(Parameter_5) %>% 
  full_join(Parameter_6) %>% 
  full_join(Parameter_7) %>% 
  rename('Parameter 1'=Minimum_leaf_number,
         'Parmeter 2'=Ppsensitivity,
         'Parameter 3'= Vrn_sensitivity,
         'Parameter 4'= Base_Phyllochron,
         'Parameter 5'= Phyllochron_PhotoPeriod_Effect,
         'Parameter 6' = Early_Reproductive_Long_Day_Base,
         'Parameter 7' = EarlyReproductivePpSensitivity) %>% 
  select(1:4, 6:9) %>% 
  mutate_if(is.numeric, round, digits=1)

knitr::kable(head(Parameters, n = 5), format = "html", caption = " APSIM Parameter results") %>% 
  kable_styling("striped", full_width = F)
```


# My Digital Toolbox


![](resources/img/tidyverse.png){width=100px}
![](https://raw.githubusercontent.com/tidyverse/ggplot2/master/man/figures/logo.png){width=90px}
![](resources/img/cowplot.png){width=90px}
![](resources/img/lubridate.jpg){width=110px}
![](resources/img/rmarkdown.png){width=100px}

## Favourite tool

Tidyverse is my favorite tool as it is so incredibly useful! I also like ggplot but find it frustrating because there are too many options.


# My time went ...

My time went on understanding and tidying the data so I could attempt to work out what it meant. I also spent a large amount of time working out scripts to calculate the APSIM parameters, especially Parameter 4, BasePhyllochron, which required use of a lm function (Thanks Aswin). I also spent quite a lot of time trying to make the ggplots and improve their appearance.


# Next steps
It would be interesting to examine the data further. One possibility is to see how the data looks when considering photo-thermal time, which takes into account daylength in addition to temperature.
I am planning to learn how to use Tassel to carry out GWAS (Genome-wide Association Studies) to identify the molecular markers linked to the parameters I have calculated.  Recently a version of Tassel based on R, aptly named RTassel, has been released and it would be cool to explore the possibilities of using this for GWAS analysis. 
I would love to learn more and find more uses for R. Hopefully this is only the start of my data journey.

# My Data School Experience

I have learned so much from Data School Focus. I think one of the most important things I have learnt is to keep rawdata and manipulated data separate and to use R and Git to track any changes I make. The other thing I have learnt is to be consistent as to how data is entered, using consistent cultivar names, only having one type of data in each column and having notes and comments in separate columns or files. I have become a convert to using R to manage my data.
