---
title: "Tidy Data for Data Management Modules Wheat CAP 2018"
author: 
  name: "Jean-Luc Jannink"
  affiliation: "USDA-ARS"
date: "May 31, 2018"
output:
  html_document:
    fig_height: 6
    fig_width: 12
    highlight: haddock
    theme: cosmo
    toc: yes
    toc_float:
      collapsed: no
      toc_depth: 2
editor_options: 
  chunk_output_type: console
---

<style>
pre {
  overflow-x: auto;
}
pre code {
  word-wrap: normal;
  white-space: pre;
}
</style>

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Pre-preparation

If you rated your familiarity with R as "low" please go to [Data Carpentry: R for Genomics](http://www.datacarpentry.org/R-genomics/) *in advance* of the first class and do the following: 

* Work through the material between *at minimum* chapters 1 ("Before We Start") and 4 ("Data frames"). If you have time and can go through the whole thing, that would be excellent. It should not take too long (one hour or so).
* This website includes instructions for installing R and Rstudio and *basic* instruction on how to interact with R. 

**Please install R and Rstudio in advance**

For quick reference about a number of the functions discussed here, go to the [R Cheat Sheet page](https://www.rstudio.com/resources/cheatsheets/)

## Learning objectives
1. Start from the **raw** data
  - Data straight out of T3 or whatever database.  No manipulation of the data in software that does not leave a record of the manipulations.
2. Write a **script**
  - No matter what problem you discover downstream, efficiently correct it by correcting the script.
3. All manipulations available in XL are scriptable.  This module will teach some necessary commands.
4. Tidy data here are 2D matrices.
  - Distinguish experimental variables having to do with the design and that can explain variation (i.e., by used to model variation), and measured variables that are the result of the experiment.
  - There should be just one column for each measured variable.
  - Other columns should contain levels of factors that enable you to model measured values.
  
```{r html output}
# The script code is in gray boxes.  One line of script code per line of html.  Slide to the right if it's  too long.
print("The output of the code is in white boxes and each line starts with `##`")
```
  
## Set up needed packages
The main package we will be using here is called "tidyverse"
```{r packages}
ip <- installed.packages()
tidyverse_installed <- "tidyverse" %in% rownames(ip)
if (!tidyverse_installed) stop("ERROR: you need to install the tidyverse package")
if (tidyverse_installed) library(tidyverse)
readxl_installed <- "readxl" %in% rownames(ip)
if (!readxl_installed) stop("ERROR: you need to install the readxl package")
if (readxl_installed) library(readxl)
```
## Show traits.txt
This is what the text file that we are working with looks like
```{r show traits, echo=FALSE}
trait_file <- read_lines("data/traits.txt"); trait_file[7] <- "";
for (i in 1:7)
  cat(paste(substr(trait_file[i], 1, 75), "..."), "\n")
old_readr.num_columns <- options(readr.num_columns=0) # Prevent a long message from printing
```

## Read traits.txt
The file provided by T3 has **TWO** header rows, one row showing the trait and one row showing the trial where it was measured. R doesn't really know what to do with that.  So we read the trait names in first, then the trial names, then we concatenate (paste) the trait and trial names together, then we read in all the data and apply single column names to the data.
```{r read traits}
trait_names <- scan("data/traits.txt", what=character(), nlines=1, sep="\t", quiet=T)
trial_names <- scan("data/traits.txt", what=character(), nlines=1, sep="\t", quiet=T, skip=1)
pasted_names <- paste(trait_names, trial_names, sep="@")
pasted_names[1] <- "line_name"
phenoByTrial <- read_tsv("data/traits.txt", skip=2, col_names=FALSE)
colnames(phenoByTrial) <- pasted_names
print(phenoByTrial)
```
```{r restore default, echo=FALSE}
options(readr.num_columns=old_readr.num_columns[[1]])
```

## Explanatory vs. measured variables
The problem with `phenoByTrial` is that there is information about the values embedded in the column **names**. That is, the column names say what the trait is AND what trial it was measured in.  For R analyses to act on that information, it needs to be coded into its own column. Then, each row will have a single value in one column and the levels of all the independent variables needed to model that value in other columns.  The `gather` function of tidyr can do that.
```{r gather columns}
pheno_one_col <- gather(phenoByTrial, key="trait_trial", value, -line_name)
print(pheno_one_col)
```

## Parse trait and trial names
The result isn't very pretty because trait and trial names are glued together.  That is, the independent variable levels we need are pasted together in the `trait_trial` variable and are not accessible to a modeling function.  So we want to parse them.  Fortunately, we put in the unique `@` character for parsing.
```{r parse trait trial}
pheno_one_col <- separate(pheno_one_col, col=trait_trial, into=c("trait", "trial"), sep="@")
print(pheno_one_col)
```
### Exercise Parse Trial Names
Can you parse the trial name to extract both the year and the location of the trial?

## Simple lm model
Now this data is ready to analyze by a simple linear model using `lm`
```{r simple lm}
mod_simple <- lm(value ~ (line_name + trial + trait)^2, data=pheno_one_col)
anova(mod_simple)
```
The interpretation of the output here is not entirely simple. Focusing on the line_name:trait interaction, it has the smallest F value of any effect.  This small value says that there is not much "genotype by trait" interaction, meaning that flag leaf senescence and heading date change in parallel across genotypes.  We can verify that by estimating their genetic correlation. Before that, let's look at these traits side by side.  The function `spread` from tidyr does that.
Script note: `spread` orders the rows by lines rather than by trials. We can rearrange by trials using the `arrange` function of dyplr.  To feed in to that `arrange` function, we use the "pipe" `%>%`.  Pipes are very handy.  The result from the function on the left-hand side of the pipe is given as the first parameter of the function on the right-hand side of the pipe.
```{r two column matrix}
pheno_two_col <- spread(pheno_one_col, key=trait, value=value) %>%
  arrange(trial)
```
Without the pipe, we would have had to write:
```{r no pipe}
pheno_two_col <- spread(pheno_one_col, key=trait, value=value)
pheno_two_col <- arrange(pheno_two_col, trial)
```
Which is not so bad, but the first version is (imho) clearer because it has less text, making it easier to understand when looking over the code.
```{r stealth print, echo=FALSE}
print(pheno_two_col)
```

## Compare trials
```{r plot across trials}
ggplot(pheno_two_col) +
  geom_point(mapping=aes(x=days_to_flag_leaf_senescence, y=heading_date, color=trial))
```
The interpretation here is pretty clear: the two trials from Huntley did not allow good evaluation of maturity traits compared to the two trials from Bozeman.  So in any downstream analysis we should remove those trials.  Note: the `|` character is the logical "or".  The statement says to keep 2012 or 2013 Bozeman trials.

## Remove data from Huntley
```{r remove Huntley}
pheno_Bozeman <- filter(pheno_two_col, trial == "SW-AMPanel_2012_Bozeman" | trial == "SW-AMPanel_2013_Bozeman")
```

### Exercise Remove Huntley data
If you were successful in the "Exercise Parse Trial Names", you could do this more easily.  How?

## Manipulate variables
```{r manipulate variables}
pheno_new_var <- mutate(pheno_Bozeman, stdHD=scale(heading_date))
pheno_replace_var <- mutate_if(pheno_Bozeman, is.numeric, scale)
print(pheno_replace_var)
```

### Exercise Square root transform
Write a line of script to square root transform numeric variables in `pheno_Bozeman`

## Read genetic characters
These genetic characters of lines are accessed from T3 by showing line information and then they come in an xls table.  So we have to read .xls into R.
```{r read genetic characters}
gen_char <- read_xls("data/Line_Details.xls") %>%
  rename(line_name=Name)
print(gen_char)
```

### Exercise Filter lines
Write a statement to remove lines for which awn information is missing.  Hint: `is.na(Awned)` returns `TRUE` when awn information is missing.

## Add line info to the phenotypes
```{r line info to pheno}
pheno_Bozeman <- select(gen_char, line_name, Program, Hardness, Height, Color, Awned) %>%
  right_join(pheno_Bozeman)
print(pheno_Bozeman)
```

## Group by breeding programs
```{r group by program}
prog_means <- group_by(pheno_Bozeman, Program) %>%
  summarize(mean_leaf=mean(days_to_flag_leaf_senescence), mean_head=mean(heading_date))
print(prog_means)
ggplot(prog_means, aes(mean_head, mean_leaf)) +
  geom_point() +
  geom_text(mapping=aes(label=Program), nudge_x=0.14) +
  geom_smooth(method="lm", formula = y ~ x)
```

### Exercise Group by height
And calculate the mean and standard deviations of heading date in those groups
What are the group counts?