## 16.8 Exercise 9: ggplot2

Create the script “exercise9.R” and save it to the “Rcourse/Module3” directory: you will save all the commands of exercise 9 in that script.
Remember you can comment the code using #.

``````getwd()
setwd("~/Rcourse/Module3")``````

### 16.8.1 Exercise 9a- Scatter plot

``library(ggplot2)``

``download.file("https://raw.githubusercontent.com/biocorecrg/CRG_RIntroduction/master/ex9_normalized_intensities.csv", "ex9_normalized_intensities.csv", method="curl")``

3- Read file into object “project1” (remember the input/output tutorial!)

• It is comma separated (csv format).
• The first row is the header.
• Take the row names from the first column.
``````project1 <- read.table("ex9_normalized_intensities.csv",
sep=",",
row.names = 1)``````

4- Using ggplot, create a simple scatter plot representing gene expression of “sampleB” on the x-axis and “sampleH” on the y-axis.

``````ggplot(data=project1, mapping=aes(x=sampleB, y=sampleH)) +
geom_point()``````

5- Add a column to the data frame “project1” (call this column “expr_limits”), that will be filled the following way:

• if the expression of a gene is > 13 in both sampleB and sampleH, set to the value in “expr_limits” to “high”
• if the expression of a gene is < 6 in both sampleB and sampleH, set it to “low”
• if different, set it to “normal.”
``````# Initialize all values to "normal"
project1\$expr_limits <- "normal"
# "high" if project1\$sampleB > 13 and project1\$sampleH > 13
project1\$expr_limits[project1\$sampleB > 13 & project1\$sampleH > 13] <- "high"
# "low" if project1\$sampleB < 6 and project1\$sampleH < 6
project1\$expr_limits[project1\$sampleB < 6 & project1\$sampleH < 6] <- "low"

## more complicated version, using a for loop and if statement
# initialize column "expr_limits" with "normal"
project1\$expr_limits <- "normal"
# loop around each row of "project1"
for(i in 1:nrow(project1)){
# create an object that contains only row "i" (the row will be different at every iteration)
rowi <- project1[i,]
# test values in rowi: assign expr_limits accordingly
if(rowi\$sampleB > 13 & rowi\$sampleH > 13){
project1\$expr_limits[i] <- "high"
}else if(rowi\$sampleB < 6 & rowi\$sampleH < 6){
project1\$expr_limits[i] <- "low"
}
}``````

6- Color the points of the scatter plot according to the newly created column “expr_limits.” Save that plot in the object “p”

``````p <- ggplot(data=project1, mapping=aes(x=sampleB, y=sampleH, color=expr_limits)) +
geom_point()``````

7- Add a layer to “p” in order to change the points colors to blue (for low), grey (for normal) and red (for high). Save this plot in the object “p2.”

``p2 <- p + scale_color_manual(values=c("red", "blue", "grey"))``

8- Save p2 in a jpeg file. a. Try with RStudio Plots window (Export)
b. Try in the console:

``````jpeg("myscatterggplot.jpg")
p2
dev.off()``````

### 16.8.2 Exercise 9b- Box plot

1- Convert “project1” from a wide format to a long format: save in the object “project_long” Note: remember melt function from reshape2 package.

``````library(reshape2)
project_long <- melt(data=project1)``````

2- Produce a boxplot of the expression of all samples (i.e. each sample is represented by a box)

``````ggplot(data=project_long, mapping=aes(x=variable, y=value)) +
geom_boxplot()``````

3- Modify the previous boxplot so as to obtain 3 “sub-boxplots” per sample, each representing the expression of either “low,” “normal” or “high” genes.

``````ggplot(data=project_long, mapping=aes(x=variable, y=value, color=expr_limits)) +
geom_boxplot()``````

4- Rotate the x-axis labels (90 degrees angle).
This is new ! Google it !!

``````ggplot(data=project_long, mapping=aes(x=variable, y=value, color=expr_limits)) +
geom_boxplot() +
theme(axis.text.x = element_text(angle = 90))``````

``````ggplot(data=project_long, mapping=aes(x=variable, y=value, color=expr_limits)) +
geom_boxplot() +
theme(axis.text.x = element_text(angle = 90)) +
ggtitle("My boxplots")``````

### 16.8.3 Exercise 9c- Bar plot

1- Produce a bar plot of how many low/normal/high genes are in the column “expr_limits” of “project1.”

``````ggplot(data=project1, mapping=aes(x=expr_limits)) +
geom_bar()``````

2- Add an horizontal line at counts 250 (y-axis). Save the plot in the object “bar”

``````bar <- ggplot(data=project1, mapping=aes(x=expr_limits)) +
geom_bar() +
geom_hline(yintercept=250)``````

3- Swap the x and y axis. Save in object “bar2.”

``bar2 <- bar + coord_flip()``

4- Save “bar” and “bar2” plots in a “png” file, using the png()** function, in the console: use grid.arrange (from the gridExtra package) to organize both plots in one page !**

``````png("mybarplots.png", width=1000)
grid.arrange(bar, bar2, nrow=1, ncol=2)
dev.off()``````

### 16.8.4 Exercise 9d- Histogram

1- Create a simple histogram using project_long (column “value”).

``````ggplot(data=project_long, mapping=aes(x=value)) +
geom_histogram()``````

2- Notice that you get the following warning message" stat_bin() using `bins = 30`. Pick better value with `binwidth`.
Set “bins”" parameter of geom_histogram() to 50.

``````ggplot(data=project_long, mapping=aes(x=value)) +
geom_histogram(bins=50)``````

3- The histogram plots the expression values for All samples.
Change the plot so as to obtain one histograms per sample.

``````ggplot(data=project_long, mapping=aes(x=value, fill=variable)) +
geom_histogram(bins=50)``````

4- By default, geom_histogram produces a stacked histogram.
Change argument “position” to “dodge.”

``````hist1 <- ggplot(data=project_long, mapping=aes(x=value, fill=variable)) +
geom_histogram(position="dodge")``````

5- A bit messy ?? Run the following:

``````hist2 <- ggplot(data=project_long, mapping=aes(x=value, fill=variable)) +
geom_histogram(bins=50) +
facet_grid(~variable)``````

facet_grid() is another easy way to split the views!

6- Change the default colors with scale_fill_manual().
You can try the rainbow() function for coloring.

``hist3 <- hist2 + scale_fill_manual(values=rainbow(8))``

7- “Zoom in” the plots: set the x-axis limits from from 6 to 13.
``hist4 <- hist3 + xlim(6, 13)``
``hist5 <- hist4 + theme_minimal()``
``ggsave(filename="myhistograms.png", plot=hist5, device="png", width=20)``