05-data-visualization.Rmd

---
layout: topic
title: Data Visualization
author: Data Carpentry contributors
minutes: 30
---

```{r, echo=FALSE}
metadata <- read.csv('./data/Ecoli_metadata.csv')
```

> ## Learning Objectives {.objectives}
>
> * Basic plots 
> * Advanced plots (introducing `ggplot`)
> * Writing images (and other things) to file

# Basic plots in R

The mathematician Richard Hamming once said, "The purpose of computing is insight, not numbers", and the best way to develop insight is often to visualize data. Visualization deserves an entire lecture (or course) of its own, but we can explore a few features of R's plotting packages.

When we are working with large sets of numbers it can be useful to display that information graphically. R has a number of built-in tools for basic graph types such as hisotgrams, scatter plots, bar charts, boxplots and much [more](http://www.statmethods.net/graphs/). We'll test a few of these out here on the `genome_size` vector from our metadata.


```{r simplestats}
genome_size <- metadata$genome_size

```

## Scatterplot
Let's start with a **scatterplot**. A scatter plot provides a graphical view of the relationship between two sets of numbers. We don't have a variable in our metadata that is a continous variable, so there is nothing to plot it against but we can plot the values against their index values just to demonstrate the function.

```{r scatter-plot1, fig.align='center'}
plot(genome_size)
```

Each point represents a clone and the value on the x-axis is the clone index in the file, where the values on the y-axis correspond to the genome size for the clone. For any plot you can customize many features of your graphs (fonts, colors, axes, titles) through [graphic options](http://www.statmethods.net/advgraphs/parameters.html)
For example, we can change the shape of the data point using `pch`.

```{r, fig.align='center'}
plot(genome_size, pch=8)
```

We can add a title to the plot by assigning a string to `main`:

```{r, fig.align='center'}
plot(genome_size, pch=8, main="Scatter plot of genome sizes")
```

## Histogram
Another way to visualize the distribution of genome sizes is to use a histogram, we can do this buy using the `hist` function:

```{r, fig.align='center'}
hist(genome_size)
```

##Boxplot

Using additional information from our metadata, we can use plots to compare values between the different citrate mutant status using a **boxplot**. A boxplot provides a graphical view of the median, quartiles, maximum, and minimum of a data set. 

```{r boxplot, fig.align='center'}
# Boxplot
boxplot(genome_size ~ cit, metadata)
```

Similar to the scatterplots above, we can pass in arguments to add in extras like plot title, axis labels and colors.

```{r, fig.align='center'}
boxplot(genome_size ~ cit, metadata,  col=c("pink","purple", "darkgrey"),
        main="Average expression differences between celltypes", ylab="Expression")
```


# Advanced figures (`ggplot2`)

More recently, R users have moved away from base graphic options and towards a plotting package called [`ggplot2`](http://docs.ggplot2.org/) that adds a lot of functionality to the basic plots seen above. The syntax takes some getting used to but it's extremely powerful and flexible. We can start by re-creating some of the above plots but using ggplot functions to get a feel for the syntax.

`ggplot` is best used on data in the `data.frame` form, so we will will work with `metadata` for the following figures. Let's start by loading the `ggplot2` library.

```{r}
library(ggplot2)
```

The `ggplot()` function is used to initialize the basic graph structure, then we add to it. The basic idea is that you specify different parts of the plot, and add them together using the `+` operator.

We will start with a blank plot and will find that you will get an error, because you need to add layers.

```{r, eval=FALSE}
ggplot(metadata) # note the error 
```

Geometric objects are the actual marks we put on a plot. Examples include:

* points (`geom_point`, for scatter plots, dot plots, etc)
* lines (`geom_line`, for time series, trend lines, etc)
* boxplot (`geom_boxplot`, for, well, boxplots!)

A plot **must have at least one geom**; there is no upper limit. You can add a geom to a plot using the + operator

```{r, eval=FALSE}
ggplot(metadata) +
  geom_point() # note what happens here
```

Each type of geom usually has a **required set of aesthetics** to be set, and usually accepts only a subset of all aesthetics --refer to the geom help pages to see what mappings each geom accepts. Aesthetic mappings are set with the aes() function. Examples include:

* position (i.e., on the x and y axes)
* color ("outside" color)
* fill ("inside" color) shape (of points)
* linetype
* size

To start, we will add position for the x- and y-axis since `geom_point` requires mappings for x and y, all others are optional.

```{r, fig.align='center'}
ggplot(metadata) +
  geom_point(aes(x = sample, y= genome_size))
```

The labels on the x-axis are quite hard to read. To do this we need to add an additional theme layer. The ggplot2 `theme` system handles non-data plot elements such as:

* Axis labels
* Plot background
* Facet label backround
* Legend appearance

There are built-in themes we can use, or we can adjust specific elements. For our figure we will change the x-axis labels to be plotted on a 45 degree angle with a small horizontal shift to avoid overlap. We will also add some additional aesthetics by mapping them to other variables in our dataframe. _For example, the color of the points will reflect the number of generations and the shape will reflect citrate mutant status._ The size of the points can be adjusted within the `geom_point` but does not need to be included in `aes()` since the value is not mapping to a variable.

```{r, fig.align='center'}
ggplot(metadata) +
  geom_point(aes(x = sample, y= genome_size, color = generation, shape = cit), size = rel(3.0)) +
  theme(axis.text.x = element_text(angle=45, hjust=1))
```

When learning ggplot, it is often useful to break things up into several steps and view the plot after each step.  For example we could create a plot without a custom theme, view it and then add the customizations to our incomplete plot.

```{r, fig.align='center',eval=FALSE}
# Create the plot object without actually displaying anything
p <- ggplot(metadata) +
  geom_point(aes(x = sample, y= genome_size, color = generation, shape = cit), size = rel(3.0)) 

p # Display the plot
```

```{r, fig.align='center',eval=FALSE}
p <- p + theme(axis.text.x = element_text(angle=45, hjust=1))
p
```


## Histogram

To plot a histogram we require another geometric object `geom_bar`, which requires a statistical transformation. Some plot types (such as scatterplots) do not require transformations, each point is plotted at x and y coordinates equal to the original value. Other plots, such as boxplots, histograms, prediction lines etc. need to be transformed, and usually has a default statistic that can be changed via the `stat_bin` argument. 

```{r, eval=FALSE}
ggplot(metadata) +
  geom_bar(aes(x = genome_size))
  
```

Try plotting with the default value and compare it to the plot using the binwidth values. How do they differ?

```{r, fig.align='center'}
ggplot(metadata) +
  geom_bar(aes(x = genome_size), stat = "bin", binwidth=0.05)
  
```

## Boxplot

Now that we have all the required information on let's try plotting a boxplot similar to what we had done using the base plot functions at the start of this lesson. We can add some additional layers to include a plot title and change the axis labels. Explore the code below and all the different layers that we have added to understand what each layer contributes to the final graphic.

```{r, fig.align='center'}
ggplot(metadata) +
  geom_boxplot(aes(x = cit, y = genome_size, fill = cit)) +
  ggtitle('Boxplot of genome size by citrate mutant type') +
  xlab('citrate mutant') +
  ylab('genome size') +
  theme(panel.grid.major = element_line(size = .5, color = "grey"),
          axis.text.x = element_text(angle=45, hjust=1),
          axis.title = element_text(size = rel(1.5)),
          axis.text = element_text(size = rel(1.25)))

```

## Raw data on the boxplot

If we want to see the raw data points on top of our boxes we need to add another layer.  Since we want points we use `geom_point()` but because the `x` axis is categorical we are likely to have points plotted over the top of one another.  To avoid this we set the position to `jitter` which adds a small random value to the `x` position.

```{r, fig.align='center'}
ggplot(metadata) +
  geom_boxplot(aes(x = cit, y = genome_size, fill = cit)) +
  ggtitle('Boxplot of genome size by citrate mutant type') +
  xlab('citrate mutant') +
  ylab('genome size') +
  theme(panel.grid.major = element_line(size = .5, color = "grey"),
          axis.text.x = element_text(angle=45, hjust=1),
          axis.title = element_text(size = rel(1.5)),
          axis.text = element_text(size = rel(1.25))) +
  geom_point(aes(x = cit, y= genome_size ),position="jitter")
```

Note that the `x` and `y` aesthetics for geom_boxplot and geom_point are the same in this plot.  We could have saved a bit of typing by adding these aesthetics to the base ggplot layer. 

```{r, fig.align='center',eval=FALSE}
ggplot(metadata,aes(x = cit, y = genome_size)) +
  geom_boxplot(aes(fill = cit)) +
  ggtitle('Boxplot of genome size by citrate mutant type') +
  xlab('citrate mutant') +
  ylab('genome size') +
  theme(panel.grid.major = element_line(size = .5, color = "grey"),
          axis.text.x = element_text(angle=45, hjust=1),
          axis.title = element_text(size = rel(1.5)),
          axis.text = element_text(size = rel(1.25))) +
  geom_point(position="jitter")
```



### Challenge 1

Try adding `+ facet_wrap(~clade)` to the previous plot.  It should produce a multi-facetted plot with data for each clade in a separate panel. Now use what you learned in the lesson on `dplyr` to pre-filter the data so that the `NA` and `unknown` panels don't appear in this plot.

### Challenge 2

Create a scatterplot with generation on the x axis, genome_size on the y axis and points coloured by `cit`. 



## Writing figures to file

There are two ways in which figures and plots can be output to a file (rather than simply displaying on screen). The first (and easiest) is to export directly from the RStudio 'Plots' panel, by clicking on `Export` when the image is plotted. This will give you the option of `png` or `pdf` and selecting the directory to which you wish to save it to. The second option is to use R functions in the console, allowing you the flexibility to specify parameters to dictate the size and resolution of the output image. Some of the more popular formats include `pdf()`, `png`.

Initialize a plot that will be written directly to a file using `pdf`, `png` etc. Within the function you will need to specify a name for your image, and the with and height (optional). Then create a plot using the usual functions in R. Finally, close the file using the `dev.off()` function. There are also `bmp`, `tiff`, and `jpeg` functions, though the jpeg function has proven less stable than the others.


```{r, eval=FALSE}
pdf("figure/boxplot.pdf")

ggplot(example_data) +
  geom_boxplot(aes(x = cit, y =....) +
  ggtitle(...) +
  xlab(...) +
  ylab(...) +
  theme(panel.grid.major = element_line(...),
          axis.text.x = element_text(...),
          axis.title = element_text(...),
          axis.text = element_text(...)

dev.off()
```

## Tidyr (Optional)

A great package that works very well together with `dplyr` and `ggplot` is `tidyr`.  Using these packages could fill up a course in itself.  If you're curious this is a quick illustration of how they all work together.

For this example you will need another dataset which you can download  [here](./data/squid.rda)

The data have been packaged into a special R data format.  Instead of reading this data in with the `read.csv` command we load it using `load` like this;

```{r}
load("data/squid.rda")
```

When you run this command it will automatically create an object in your environment called `squid`.  Take a look at the data using head

```{r}
head(squid)
```

There are several identifying columns `prot_id`, `protein_group`.  There is a column with protein expression measurements (from one particular tissue, Slime) and there are 5 columns with RNA expression measurements for different tissues `Arms`, `Brain`, `Slime`, `VM` and `DM`. (VM and DM are Ventral and Dorsal Mantle). This data needs tidying because we have got multiple observations per row.  Specifically we are interested in collapsing all of the RNA expression measurements into a single column with another column to code for the tissue they came from.  This kind of manipulation is exactly what the `tidyr` package is for.  

Load the `tidyr` package

```{r}
library(tidyr)
library(dplyr)
```

It also turns out that there is another issue with the data.  Important information is encoded in the rownames but we would prefer to have that information in its own column.  We deal with that using the `add_rownames()` function from `dplyr`. 

```{r}
squid %>% add_rownames("gene_id")
```

Now we gather the RNA expression information into a pair of columns, and we also use piping to combined the `tidyr` operation with the `add_rownames()` operation from `dplyr`.

```{r}
tidy_squid <- squid %>% add_rownames("gene_id") %>% 
  gather(Tissue , TPM, -prot_id, -protein_group,-iBAQ,-gene_id)
```

### Challenge 3

Use the `tidy_squid` data to create a scatterplot of `iBAQ` protein expression versus `RNA` expression. Colour the points to distinguish different tissue types.



### Challenge 4

This time use the `tidy_squid` data to create  a basic heatmap with tissues as columns, `gene_id` as rows and with heatmap squares coloured according to `TPM` value. You will need to use `geom_tile()` for this.



### Solutions

Challenge 1

<!--
```{r, fig.align='center',echo=FALSE}
library(dplyr)
library(tidyr)
```
-->

```{r, fig.align='center'}

metadata_c <- metadata %>% filter(!is.na(clade)) %>% filter(clade!='unknown')

ggplot(metadata_c) +
  geom_boxplot(aes(x = cit, y = genome_size, fill = cit)) +
  ggtitle('Boxplot of genome size by citrate mutant type') +
  xlab('citrate mutant') +
  ylab('genome size') +
  theme(panel.grid.major = element_line(size = .5, color = "grey"),
        axis.text.x = element_text(angle=45, hjust=1),
        axis.title = element_text(size = rel(1.5)),
        axis.text = element_text(size = rel(1.25))) + facet_wrap(~clade)
```

Challenge 2

```{r, fig.align='center'}
ggplot(data=metadata,aes(x=generation,y=genome_size)) + geom_point(aes(color=cit))
```

Challenge 3

```{r, fig.align='center'}
ggplot(tidy_squid,aes(x=iBAQ,y=TPM)) + geom_point(aes(color=Tissue))
```

Challenge 4

```{r, fig.align='center'}
ggplot(tidy_squid,aes(x=Tissue,y=gene_id)) + geom_tile(aes(fill=TPM))
```



Resources:
---------
We have only scratched the surface here. To learn more, see the [ggplot reference site](http://docs.ggplot2.org/), and Winston Chang's excellent [Cookbook for R](http://wiki.stdout.org/rcookbook/Graphs/) site. Though slightly out of date, [ggplot2: Elegant Graphics for Data Anaysis](http://www.amazon.com/ggplot2-Elegant-Graphics-Data-Analysis/dp/0387981403) is still the definative book on this subject. Much of the material here was adpapted from [Introduction to R graphics with ggplot2 Tutorial at IQSS](http://tutorials.iq.harvard.edu/R/Rgraphics/Rgraphics.html).