# Little plants in Sweden and Italy

_Arabidopsis_, also known as "mouse-ear cress" is the mouse of the plant world -- it is the model plant.  It is a "cousin" of mustanrd, and small annual, with a quick growing season.  In this notebook we will analyze data from an experiment that was designed to investigate the genetics of adaptation in _Arabidopsis_.

A set of recombinant inbred (RI) lines was derived from a cross of a strain native to Sweden and a strain native to Italy.  Italy and Sweden represent the northern and southern ends of the European range of _Arabidopsis_; these strains may have genetic elements that confer adaptation to their home range; they may also have genetic elements that help them adapt everywhere.

About 400 RI lines were grown in Sweden and Italy for three consecutive years 2009, 2010, and 2011), and a measure of fitness (number of fruits per plant) was recorded.  Genomewide genetic markers were measured.  In this analysis we will only look at the gene [FLC](https://en.wikipedia.org/wiki/Flowering_Locus_C) which is very well-studied and would make a nice candidate gene for conferring genetic adaptation to the environment.

## Question

Does fitness measured in previous years, in the other site, and FLC help predict fitness in Italy in 2011?  We will investigate this question using linear regression.

## A data analysis notebook

This notebook is a demonstration of what a data analysis report might look like.  It is also a live notebook where you can try out R by modifying what has been already done by someone else.  You can also tweak the analysis by trying out different things compared to the original author.

The code is in R; the commentary is in [Markdown](https://daringfireball.net/projects/markdown/), which is a lightweight text markup format.

## Data source

These data are from Ågren et. al. (2013). _PNAS_ 110:21077-21082.  More details are [here](https://bitbucket.org/linen/smalldata/src/47d88b298b2b1cb18d900e13aea8de8ec14de311/arabidopsis). The data is a spreadsheet CSV form, so we use the `read.csv` function.  The CSV file is a plain text file that can be read by virtually any program (Excel, SAS, R, Python, etc), is portable across operating systems, fairly space efficient, and therefore highly recommended.  Whenever you have some data, you should also include a plain text README file that describes the data source, and contents.  For elaborate datasets, you may need more than one such description file.

## Read data

We read data from an online source (Bitbucket) by specifying the URL.  Instead of a URL, you could have a filename in your own computer.

```{r}
# read data from bitbucket site in CSV (comma separated value) format
# the result is the object called "agren" which tabular data
agren <- read.csv("https://bitbucket.org/linen/smalldata/raw/47d88b298b2b1cb18d900e13aea8de8ec14de311/arabidopsis/agren2013.csv")
```

We now have the `agren` object which has the spreadsheet data.  It is called a `data.frame`.

```{r}
# the tabular data has a class called "data frame" used for tabular data
class(agren)
head(agren)
```

## Data checks

### How many columns and rows?

Let us get the names of the columns using the command `names`.

```{r}
# get the names of the columns
names(agren)
agren <- agren[,-7]
```

The first six columns have the fitness (average number of fruit per plant) for each RI line.  The abbreviations indicate the year and site: `it` for Italy and `sw` for Sweden.  The `id` column is for RI line ID, and `flc` is the genotype of the FLC gene mentioned earlier.

Let us now see how many rows of data we have using `nrow`.

```{r}
nrow(agren)
```

Each command manipulates objects, and sometimes creates new ones.  Let us see what we have in our workspace using `ls`.

```{r}
# list what we have in our workspace
ls()
```

As expected, we only have one object, the data we read in earlier.  Let us dig deeper.

### Looking at the first few rows

We can look at the first few rows using `head`.  You can customize how many rows to see and how many decimal points to display.

```{r}
# get the first few rows of the data
head(agren)
```

```{r}
# round off to two decimals and get the first 10 rows
round(head(agren,n=10),1)
```

To look at just one column, we have to refer to it by the `$` sign.

```{r}
round(head(agren$it09),2)
```

We can also equivalently use the column number as follows.

```{r}
round(head(agren[,1]),2)
```

## Visualizations

### Histogram

Now let's visualize the data.  We will first make a histogram of fitness in Italy for 2009.  We have to tweak the figure size using the `options` command.  The default size is bigger.

```{r}
# set plot size
options(repr.plot.width=4, repr.plot.height=4)
# make histogram using command "hist"
# set title using argument "main", x axis label using "xlab"
# and number of bins using "nbins"
hist(agren$it09,main="Histogram of fitness Italy 2009",xlab="Number of fruits",nclass=20)
```

### Pairwise scatterplots

Now let us plot every fitness measure in each site and year against every other site and year, to see how they are related.

```{r}
# make plot bigger
options(repr.plot.width=6, repr.plot.height=6)
# plot all data columns against each other
# pch is for plotting character
pairs(agren[,1:6],pch=".")
```

Notice the flaring in the scatterplots.  This means that there is more variation in fitness when fitness is high.  Often this can be remedied by a logarithmic transform.  We will use base 2, but any other base can be used.

```{r}
# make the same scatterplot, but taking logarithm (base 2)
pairs(log2(agren[,1:6]),pch=".")
```

The scatterplot clouds look more oval (or circle) shaped.  The log transform can also be justified on intuitive grounds since addition on the log scale is multiplication on the original scale.  So a value 1 log2 units higher is 2 times the original value.

For ease of use, we will make a new data frame with the log2 transformed data.  We will store it in a new object called `agren2`.

```{r}
# make new data frame with log2-transformed fitness data, and FLC
agren3 <- cbind(log2(agren[,1:6]),flc=agren$flc)
# check names of data frame
names(agren3)
# see what's in workspace now
ls()
```

## Linear regression

Now we will fit a linear regression model to predict fitness in 2011 using data from all other sites and years.  This can be done using the `lm` function.  The result is stored in the object called `out1`.  To see the results, we use the `summary` function.

```{r}
# fit linear model using data in data frame called "agren3"
# the response variable is it11 and all the rest are predictors
out1 <- lm(it11~it09+it10+sw11+sw09+sw10+flc,data=agren3)
# print summary of linear regression
summary(out1)
```

Most of the output is self-explanatory if you are familiar with linear regression.  Essentially, about 50% of the variation in fitness is explained by fitness in other years and site.  Fitness in Italy in the other two years seems to be most important.


## Conclusion

About 50% of the variation in fitness in Italy in 2011 is explained by fitness in Sweden in 2011 and fitness in both sites in prior years. Fitness in Italy in prior years is most predictive.  When fitness in Italy in prior years is considered, fitness in Sweden and FLC genotype do not add to predictive value.
