19.9 Exercise 12: ggplot2

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

correction

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

19.9.1 Exercise 12a- Scatter plot

1- Load ggplot2 package

correction

library(ggplot2)

2- Download the data we will use for plotting:

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

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

About this file:

• It is comma separated (csv format).
• The first row is the header.
• Take the row names from the first column.

correction

project1 <- read.table("ex12_normalized_intensities.csv", 
  sep=",", 
  header=TRUE, 
  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.

correction

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”.

correction

# 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”

correction

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”.

correction

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:

correction

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

19.9.2 Exercise 12b- 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.

correction

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)

correction

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.

correction

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 !!

correction

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

5- Finally, add a title of your choice to the plot.

correction

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")

19.9.3 Exercise 12c- Bar plot

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

correction

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”

correction

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”.

correction

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 !**

correction

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

19.9.4 Exercise 12d- Histogram

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

correction

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.

correction

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.

correction

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”.

correction

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.

correction

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

7- “Zoom in” the plots: set the x-axis limits from from 6 to 13.
Add the xlim() layer.

correction

hist4 <- hist3 + xlim(6, 13)

8- Change the default theme to theme_minimal()

correction

hist5 <- hist4 + theme_minimal()

9- Save that last plot to a file (format of your choice) with ggsave()

correction

ggsave(filename="myhistograms.png", plot=hist5, device="png", width=20)