[This article was first published on Learning Omics » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

---

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

ls()
rm(list=ls())
ls()

getwd() 
setwd("c:/your_folder/project_1") 
getwd()

dir()

dat1 <- read.csv("OTU.csv", header=TRUE,row.names=1, sep=",") 
meta <- read.csv("metadata.csv", header=TRUE,row.names=1, sep=",")

taxa.names <- dat$taxonomy 
dat2 <- as.matrix(dat[,-dim(dat)[2]])

s_abundances<-apply(dat2,2,sum)

bads<-dat2[,s_abundances<1000]
goods<-dat2[,s_abundances>1000]

ncol(dat2)-ncol(goods)
ncol(bads)

dat2<-goods

write.table(badm, file="output/bads.csv", col.names=NA,row.names=TRUE, sep=",", quote=FALSE)
write.table(goodm, file="output/goods.csv", col.names=NA,row.names=TRUE, sep=",", quote=FALSE)

dat2 <- scale(dat2, center=F, scale=colSums(dat2))

To begin, you will need transpose the OTU table so that the OTU ID numbers become the top row rather than the first column.

dat2 <-t(dat2)

Next, take a closer look at the way taxonomical designations are written in your OTU table. Are there six or seven levels?  If it looks like

• “Root;Bacteria;Firmicutes;”"Clostridia”";Clostridiales;”"Lachnospiraceae”";”"Lachnospiraceae Incertae Sedis”
• “k__Bacteria; p__Firmicutes; c__Clostridia; o__Clostridiales; f__Lachnospiraceae; g__Coprococcus; s__Coprococcus eutactus”

there are seven levels designated by the Root/Kingdom, Phylum, Class, Order, Family, Genus and Specie. You will need to start numbering the levels at level 2 rather than 1 in the extraction functions. However, if the beginning of your taxonomical designations are phyla then start numbering at level 1; or if phyla are not listed until the third spot in the taxonomical designations, start numbering the levels at 3.

d.phylum = otu2taxonomy(dat4,level=2,taxa=taxa.names)
d.class = otu2taxonomy(dat4,level=3,taxa=taxa.names)
d.order = otu2taxonomy(dat4,level=4,taxa=taxa.names)
d.family = otu2taxonomy(dat4,level=5,taxa=taxa.names)
d.genus = otu2taxonomy(dat4,level=6,taxa=taxa.names)
d.species = otu2taxonomy(dat4,level=7,taxa=taxa.names)

To look at the output, transpose and export.

phylum2 <-t(d.phylum)
class2 <-t(d.class)
order2 <-t(d.order)
family2 <-t(d.family)
genus2 <-t(d.genus)
species2 <-t(d.species)
write.table(phylum2, file="output/phyla.csv", col.names=NA,row.names=TRUE, sep=",", quote=FALSE) 
write.table(class2, file="output/classes.csv", col.names=NA,row.names=TRUE, sep=",", quote=FALSE)
write.table(order2, file="output/orders.csv", col.names=NA,row.names=TRUE, sep=",", quote=FALSE)
write.table(family2, file="output/families.csv", col.names=NA,row.names=TRUE, sep=",", quote=FALSE)
write.table(genus2, file="output/genera.csv", col.names=NA,row.names=TRUE, sep=",", quote=FALSE)
write.table(species2, file="output/species.csv", col.names=NA,row.names=TRUE, sep=",", quote=FALSE)

These steps allow you to analyze your OTU data at various taxonomical levels within and outside of R.

Merging metadata and OTU tables

Next we will merge the metadata and OTU table together. Originally, I tried to figure out a way to avoid this step but R was not being so cooperative with reordering. We are merging the two tables so that when we rearrange the samples according to a metadata variable, the OTU table rearranges as well.

First choose which taxonomical level you want to look at. For simplicity’s sake, let’s look at phyla first. This first script merges the two tables by matching up the Sample ID/Names. The “all.x=FALSE” portion makes sure that metadata information is not included for samples that were dropped earlier from the data processing (if you skipped that step, this should not affect the output).

mergephylum<-merge(meta, d.phylum, by="row.names",all.x=FALSE)

There are multiple ways to double check the merge worked. The easiest is just by looking at the dimensions of each object involved.

dim(meta)
dim(d.phylum)
dim(mergephylum)ncol

The first number reflects the number of rows while the second refers to the amount of columns. The number of rows of the final merge should match the number of rows in the phyla file (or whichever level you are looking at).  The number of columns of the final merge should be the sum of the columns in d.phylum and meta +1. Again, here is the code for output.

write.table(mergephylum, file="output/mergephylum.csv", col.names=NA,row.names=TRUE, sep=",", quote=FALSE)

Reordering by metadata variable

If you would like your heatmap in order by sample or order does not matter, skip to the next step (making the actual matrix for the heatmap).

You can reorder your samples based on a variable in your metadata. This can be anything from diet type to sample site. The variable that will direct the way your table will be reordered is not limited to just being discrete (as in group A, group B, etc.); the variable can also be continuous (blood pressure, height, pH).

For this rearrangement code, you will need the exact name of the column of the variable. Since names might have changed since importation into R, it’s best to double check. This lists the names of all the columns of your metadata.

colnames(meta)

Now use the name of the variable in question and replace “VARIABLE” in the following code.

reordermerge <- mergephylum[order(mergephylum$VARIABLE, mergephylum$Row.names),]

This will cause the table to be rearranged by increasing value of the VARIABLE (or alphabetical order) and then by the Row.names or Sample ID/Names.

Making the actual matrix for the heatmap

This next set of  script splits OTU table and metadata again so that the heatmap will only display OTU data rather than metadata as well. It involves removing columns, renaming the rows and changing the format of the object to a matrix.

(OTUcol1<-ncol(meta)+2)
(OTUcol2<-ncol(reordermerge))
justOTU<-reordermerge[,OTUcol1:OTUcol2]
justOTU[1:5,1:5]
rownames(justOTU[1:10,])
rownames(justOTU)<-reordermerge$Row.names
rownames(justOTU[1:10,])
justOTU2<-as.matrix(t(justOTU))
justOTU2[1:5,1:5]

Heatmap!!!

Hurray! After all those steps we are finally to the steps of making the actual heatmap.

First install the package that contains the codes to make the heatmap. (Refresh your memory on packages here.)

install.packages("gplots")
library(gplots)

Now print out your heatmap!

heatmap.2(justOTU2,Rowv=TRUE, Colv=FALSE, scale="column", trace="none", col=redgreen, xlab="sample", ylab="phylum", margins=c(10,15))

If you are happy with this image, skip to “Export the Heatmap.“

This tutorial will be covering some of the arguments within this function for customization. In case you want to know about all the possibilities, read this documentation on heatmap.2.

Looking at some of the arguments already present:

• justOTU2-This is where you put the object name of the heatmap you would like to visualize.
• Colv and Rowv-Hold your horses and read about these in the next subsection.
• scale-For the scale argument you can input “none”, “column” or “row”. R normalizes the rows or columns based on the Z-score or standard score (not to be confused with Z-factor) which indicates how many standard deviations away the value is from the mean. In other words, the scale indicates which values are outliers and possibly significant.
• trace-You should probably set this to “none” because if you put “column” “row” or “both”, you get some funky lines separating the boxes in the heatmap. Unless your data set is pretty small, then you would benefit from adding the following code with a different color.

  • tracecol="cyan"

• col-This controls the color scheme. Read about changing the color scheme in the second following subsection.
• xlab and ylab-This is where you put labels for your X and Y axis.
• margins-Changing the numbers in this argument changes the amount of white space surrounding your heatmap. I highly recommend fiddling with these numbers to make sure the visualization space fits the heatmap of your data.

Optional: Dendrograms

A dendrogram is a tree placed on right and/or top sides of the heatmap. Instead of clustering phylogenetically similar samples or phyla, these trees cluster columns or rows that have similar values. This clustering causes the rows and columns to be reordered when you input “TRUE” for Rowv or Colv. If you input Colv=TRUE, then you will lose any customized ordering you did earlier when you reordered the metadata.

Optional: Change the color scale

Although many default heatmaps use the red-green color scheme, I would suggest choosing a different palette to A) better match the color scheme of your paper and more seriously, B) allow color blind people the chance to understand your graphic.

I will not go much into the specifics of choosing good color palettes, but the most important aspect of designing graphics is that people can actually read and understand them. The color scheme must assist in this goal. So yellow text on a white background is a huge no-no (regardless if the audience is color blind or not) and remember, there are people who can be blue-green color blind. In case you want to learn more about colors in R in general, I stumbled across this presentation that has excellent information but is displayed in the most horrendous manner graphically.

There are already a couple color schemes present in R. You can find them by entering

?cm.colors

The easiest way to change the color scheme to a customized scheme is by using the package “rcolorbrewer.”

install.packages("RColorBrewer")
library(RColorBrewer)

The information in your heatmap would probably benefit best from a monochromatic color scheme (white, light blue, dark blue, etc.) or a diverging color scheme (dark red, pink, white, light blue, dark blue, etc.) Creating both utilize a similar function from the RColorBrewer: colorRampPalette. (For more information about RColorBrewer, look at an explanation here and and another one  by Dr. Stewart MacArthur).

For a monochromatic scale, change the color names within the quotation markers:

scaleblue<-colorRampPalette(colors=c("white","blue"))(100)

For a diverging scale, input three colors.

scaleyellowblue<-colorRampPalette(colors=c("yellow","white","blue"))(100)

The number at the end corresponds to the number of colors generated in the scale. Thus if you inputted

scaleyellowblue<-colorRampPalette(colors=c("yellow","white","blue"))(3)

You would only get the 3 original colors in the scale.

To apply the color scale, input the name of your custom scale in the argument for col.

heatmap.2(justOTU2,Rowv=TRUE, Colv=FALSE, scale="column", trace="none", col=COLORSCALE, xlab="sample", ylab="phylum", margins=c(10,15))

For example, the diverging yellow-white-blue color scheme:

heatmap.2(justOTU2,Rowv=TRUE, Colv=FALSE, scale="column", trace="none", col=scaleyellowblue, xlab="sample", ylab="phylum", margins=c(10,15))

If you want to get fancy, try choosing colors from:

For example, the Science paper I mentioned earlier had a really cool protein heatmap with a customized color scheme. To replicate that scheme:

colors in R

• http://simplystatistics.org/2011/10/17/colors-in-r/
• http://www.compbiome.com/2010/12/r-using-rcolorbrewer-to-colour-your.html

Optional: Colored Side Bar

In the several papers, you can find some authors put a customized colored side bar along the top or left hand side of the heatmap to denote different groupings of the samples. You would most likely want to create a customized colored side bar based on the variable you chose to base your data reordering. The following function is based upon a code created by Dr. Peter Cock found here.

First, you will need to create a function that assigns a color to each category. This function will uses a whole bunch of “if” functions. Only the first category will be preceded by “if” and following categories will be “else if”. Add a line for each category. This example has 5 categories.

color.map<-function(diet.type){
 if(diet.type=="0") "red" 
 else if(diet.type=="1") "green" 
 else if(diet.type=="2") "yellow"
 else if(diet.type=="3") "blue"
 else if(diet.type=="4") "purple"
 }

The text in between quotation marks after “VARIABLE==” should be the exact text of the category designation within the column. So if your categories are “dry”, “muddy” and “sandy”, you should replace the numbers.

color.map<-function(CATEGORY){
 if(CATEGORY=="dry") "red" 
 else if(CATEGORY=="muddy") "green" 
 else if(CATEGORY=="sandy") "yellow"
 }

If the variable is continuous(such as pH), the color scheme would probably be better suited by a monochromatic scale (look at previous section for more information). Install rcolorbrewer and create a custom palette.You can change the colors and number of tones the function will produce.

install.packages("RColorBrewer")
library(RColorBrewer)
colorRampPalette(colors=c("yellow","red"))(5)

Use the hex codes outputted the last code in your color map function.

color.map<-function(CATEGORY){
 if(CATEGORY=="0") "#FFFF00" 
 else if(CATEGORY=="1") "#FFBF00" 
 else if(CATEGORY=="2") "#FF7F00"
 else if(CATEGORY=="3") "#FF3F00"
 else if(CATEGORY=="4") "#FF0000"
 }

Now that you have a function linking category and color, apply the color map to the whole column of data.

sidebarcolors <- unlist(lapply(reordermerge$VARIABLE, color.map))

Now add the argument ColSideColors=sidebarcolors to your heatmap code.

heatmap.2(justOTU2,Rowv=TRUE, Colv=FALSE, scale="column", trace="none", col=redgreen, xlab="sample", ylab="phylum", margins=c(10,15), ColSideColors=sidebarcolors)

There is actually a way to put in multiple side bars, but I have not looked much into the code. You can find an explanation here.

Here is an example I made with a dendrogram for rows, a customized color scale and customized sidebar colors. The color scale was created using a 9 tone palette created by rcolorbrewer.

Export the Heatmap

Once you have found the perfect way to graphically display your information as a heatmap, you will probably want to save the beauty. The export process of images is different from tables. It is similar to turning a camera on, placing an object in front of the camera, and then taking the picture. When you turn on your camera, you need to choose the size of the photo, name, and by extension, folder of the photo will be written to. The “placing the object” part is rewriting whatever final code you have come up with for your heatmap.

png(filename="output/phylum.png", width=1200, height=800)
heatmap.2(justOTU2,Rowv=TRUE, Colv=FALSE, scale="column", trace="none", 
 col=redgreen, xlab="diet type", ylab="phylum", 
 ColSideColors=sidebarcolors, margins=c(15,15)) 
dev.off()

This is it! In case your curiousity about heatmaps has not been filled, there are many other aspects of heatmaps that can still be twiddled. For example, this tutorial shows how to overlay numbers on the colored boxes. If you have any suggestions to make this tutorial better, please comment below!

---

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

---

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.