1A. Creating a Project

1A. We will begin by getting organized. This entails creating a new RStudio Project, creating a new R script, and downloading the input data we will use For this hands-on session.

• Begin by opening RStudio and creating a new Project File > New Project > New Directory > New Project make sure you are working in the directory (a.k.a. folder!) where you want to keep all your work for this class organized. For example, for me this is a directory on my Desktop with the class name (see animated figure below). We will create our project as a sub-directory called class05 in this location.

The key step here is to name your project after this class session (i.e. “class05”) and make sure it is a sub-directory of where ever you are organizing all your work for this class (see image above for an example).

1B. Getting data to plot

1B. Your next task is to download the ZIP file containing all the input data for this lab to your computer and move it to your RStudio project.

• Download the ZIP file
• Note that on most computers this will typically unzip itself after downloading to create a new directory/folder named bimm143_05_rstats. Check your Downloads directory/folder. However, some computers may require you to double click on the ZIP file to begin the unzipping process.
• Once unzipped move the resulting bimm143_05_rstats folder into your R project directory. You can use your Finder window (on Mac) or File Explorer (on Windows) to do this.

2A. Line plot

2A. The file weight_chart.txt from the example data you downloaded above contains data for a growth chart for a typical baby over the first 9 months of its life.

Your first task is to get this data into R:

• Find the file weight_chart.txt in your RStudio Files panel. This should typically be on the right lower side of your RStudio window.
• Click on the weight_chart.txt file to get a preview of its contents in RStudio.
• How are the records in the file separated? Is there a comma, space, tab or other character between the data entries? Does the file have a “Header” line that contains the names of the variables (this is often the first line of a file)?

Now we know something about the file format we can decide on the R function to use for data reading into R.

• Look at the help page for the read.table() function and decide on the function and arguments to use for importing the weight_chart.txt file.
• What effect does setting and changing the argument header=FALSE to header=TRUE in your chosen read function?
• Make sure to examine the object you create and store
• Add your working code for reading the file to your R script and make sure to save your file.

weight <- read.table("bimm143_05_rstats/weight_chart.txt", header=TRUE)

Use the plot() function to draw this as a point and line graph with the following changes:

• Change the point character to be a filled square (pch=15)
• Change the plot point size to be 1.5x normal size (cex=1.5)
• Change the line thickness to be twice the default size (lwd=2)
• Change the y-axis to scale between 2 and 10kg (ylim=c(2,10))
• Change the x-axis title to be Age (months) (xlab="Age (months)")
• Change the y-axis title to be Weight (kg) (ylab="Weight (kg)")
• Add a suitable title to the top of the plot (main="Some title")

plot(weight$Age, weight$Weight, typ="o", 
     pch=15, cex=1.5, lwd=2, ylim=c(2,10), 
     xlab="Age (months)", ylab="Weight (kg)", 
     main="Baby weight with age") 

2B. Barplot

2B. The file feature_counts.txt contains a summary of the number of features of different types in the mouse GRCm38 genome.

• Read this data into R using the same overall procedure you used in section A above (i.e. check the format of the file to decide on the read function and arguments to use).
• What is the field separator and header format of this file?
• Can you plot this as a barplot. Note that the data you need to plot is contained within the Count column of the data frame, so you pass only that column as the data.

mouse <- read.table("bimm143_05_rstats/feature_counts.txt", sep="\t", header=TRUE)

barplot(mouse$Count)

Once you have the basic plot above make the following changes:

• The bars should be horizontal rather than vertical (horiz=TRUE).
• The count axis should be labelled (ylab="A title")
• The feature names should be added to the y axis. (set names.arg to the Feature column of the data frame)
• The plot should be given a suitable title (main="Some title")
• The text labels should all be horizontal (las=1) Note that you can pass this parameter either via par, or as an additional option to barplot.
• The margins should be adjusted to accommodate the labels (par mar parameter). You need to supply a 4 element vector for the bottom, left, top and right margin values.

Look at the value of par()$mar to see what the default values are so you know where to start. Note that you will have to redraw the barplot after making the changes to par.

par(mar=c(3.1, 11.1, 4.1, 2))
barplot(mouse$Count, names.arg=mouse$Feature, 
        horiz=TRUE, ylab="", 
        main="Number of features in the mouse GRCm38 genome", 
        las=1, xlim=c(0,80000))

2C. Histograms

[Extension if you have time] Use this hist() function to plot out the distribution of 10000 points sampled from a standard normal distribution (use the rnorm() function) along with another 10000 points sampled from the same distribution but with an offset of 4.

• Example: c(rnorm(10000),rnorm(10000)+4)
• Find a suitable number of breaks to make the plot look nicer (breaks=10 for example)

x <- c(rnorm(10000),rnorm(10000)+4)
hist(x, breaks=80)

3A. Providing color vectors

3A. The file male_female_counts.txt contains a time series split into male and female count values.

• Plot this as a barplot
  
• Make all bars different colors using the rainbow() function.
• Note that the rainbow() function takes a single argument, which is the number of colors to generate, eg rainbow(10). Try making the vector of colors separately before passing it as the col argument to barplot (col=rainbow(10)).
• Rather than hard coding the number of colors, think how you could use the nrow() function to automatically generate the correct number of colors for the size of dataset.
• Re-plot, and make the bars for the males a different color to those for the females. In this case the male and female samples alternate so you can just pass a 2 color vector to the col parameter to achieve this effect. (col=c("blue2","red2")) for example.

mf <- read.delim("bimm143_05_rstats/male_female_counts.txt")

barplot(mf$Count, names.arg=mf$Sample, col=rainbow(nrow(mf)), 
        las=2, ylab="Counts")

barplot(mf$Count, names.arg=mf$Sample, col=c("blue2","red2"), 
        las=2, ylab="Counts")

3B. Coloring by value

3B. The file up_down_expression.txt contains an expression comparison dataset, but has an extra column which classifies the rows into one of 3 groups (up, down or unchanging).

Plot this as a scatterplot (plot) with the up being red, the down being blue and the unchanging being grey.

• Read in the file to an object named genes

genes <- read.delim("bimm143_05_rstats/up_down_expression.txt")

• How many genes are detailed in this file (i.e. how many rows are there)?
• To determine how many genes are up, down and unchanging in their expression values between the two conditions we can inspect the genes$State column. A useful function for this is the table() function, e.g. try table(genes$State).

table(genes$State)

## 
##       down unchanging         up 
##         72       4997        127

• For graphing start by just plotting the Condition1 column against the Condition2 column in a plot
• Pass the State column as the col parameter (col=genes$State for example). This will set the color according to the state of each point, but the colors will be set automatically from the output of palette.

plot(genes$Condition1, genes$Condition2, col=genes$State, 
     xlab="Expression condition 1", ylab="Expression condition 2")

• Run palette() to see what colors are there initially and check that you can see how these relate to the colors you get in your plot.
• Run levels() on the State column and match this with what you saw in palette() to see how each color was selected. Work out what colors you would need to put into palette to get the color selection you actually want.
• Use the palette() function to set the corresponding colors you want to use (eg palette(c("red","green","blue")) – but using the correct colors in the correct order.
• Redraw the plot and check that the colors are now what you wanted.

palette(c("blue","gray","red"))
plot(genes$Condition1, genes$Condition2, col=genes$State, xlab="Expression condition 1", ylab="Expression condition 2")

3C. Dynamic use of color

The uses of color described above would be what you would use to do categorical coloring, but often we want to use color for a more quantitative purpose. In some graph types you can specify a color function and it will dynamically generate an appropriate color for each data point, but it can be useful to do this more manually in other plot types.

# Lets plot expresion vs gene regulation
meth <- read.delim("bimm143_05_rstats/expression_methylation.txt")

plot(meth$gene.meth, meth$expression)

dcols <- densCols(meth$gene.meth, meth$expression)

# Plot changing the plot character ('pch') to a solid circle
plot(meth$gene.meth, meth$expression, col = dcols, pch = 20)

# Find the indices of genes with above 0 expresion
inds <- meth$expression > 0

# Plot just these genes
plot(meth$gene.meth[inds], meth$expression[inds])

## Make a desnisty color vector for these genes and plot
dcols <- densCols(meth$gene.meth[inds], meth$expression[inds])

plot(meth$gene.meth[inds], meth$expression[inds], col = dcols, pch = 20)

dcols.custom <- densCols(meth$gene.meth[inds], meth$expression[inds],
                         colramp = colorRampPalette(c("blue2",
                                                      "green2",
                                                      "red2",
                                                      "yellow")) )

plot(meth$gene.meth[inds], meth$expression[inds], 
     col = dcols.custom, pch = 20)

4A. Mapping colors

Whilst color in general is not a great way to represent quantitative data it can be sometimes useful to use the color of points in a plot as an extra quantitative variable. This will then require a way to map colors from a palette to values in a quantitative set. There is no built in base function within R to do this so we have borrowed one from a fictional lab mate - see the file color_to_value_map.r that contains a function map.colors() that may help us do this.

#meth <- read.delim("bimm143_05_rstats/expression_methylation.txt")

plot(meth$promoter.meth, meth$gene.meth, ylab="Gene Methylation", xlab="Promoter Methylation")

# source the provided function so we can use it
source("bimm143_05_rstats/color_to_value_map.r")

mycols=map.colors(meth$expression, 
                  c(max(meth$expression), min(meth$expression)), 
                  colorRampPalette(c("blue","red"))(100))

plot(meth$promoter.meth, meth$gene.meth, 
     ylab="Gene Methylation", 
     xlab="Promoter Methylation", 
     col=mycols)