R - Defining Your Own Color schemes for HeatMaps

This post is intended at those who are beginners at R, and is inspired by a small post in Martin's bioblog.

First, we plot a "correlation heatmap" using the same logic that Martin uses. In our example, let's use the Movies dataset that comes with ggplot2.

We take the 6 genre columns, and we can compute the correlation matrix for those 6 columns.

	library(ggplot2)
	library(reshape2)
	data(movies)
	movieGenres <- movies[c(18:23)] #subset to 6 genres
	cor(movieGenres) # 6x6 cor matrix
	#ggplot likes the data 'melted' one value per row
	m <-melt(cor(movieGenres))
	p <- ggplot(data=m, aes(x=Var1, y=Var2, fill=value)) + geom_tile()

view raw gistfile1.r hosted with ❤ by GitHub

Here's what the matrix looks like:

> cor(movieGenres) # 6x6 cor matrix
                  Action   Animation      Comedy        Drama
Action       1.000000000 -0.05443315 -0.08288728  0.007760094
Animation   -0.054433153  1.00000000  0.17967294 -0.179155441
Comedy      -0.082887284  0.17967294  1.00000000 -0.255784957
Drama        0.007760094 -0.17915544 -0.25578496  1.000000000
Documentary -0.069487718 -0.05204238 -0.14083580 -0.173443622
Romance     -0.023355368 -0.06637362  0.10986485  0.103545195
            Documentary     Romance
Action      -0.06948772 -0.02335537
Animation   -0.05204238 -0.06637362
Comedy      -0.14083580  0.10986485
Drama       -0.17344362  0.10354520
Documentary  1.00000000 -0.07157792
Romance     -0.07157792  1.00000000

When we plot with the default colors we get:

It is difficult to see the details in the tiles. Now, if you want to better control the colors, you can use the handy colorRampPalette() function and combine that with scale_fill_gradient2.
Let's say that we want "red" colors for negative correlations and "green" for positives.
(We can gray out the 1 along the diagonal.)

	#set up a coloring scheme using colorRampPalette
	red=rgb(1,0,0); green=rgb(0,1,0); blue=rgb(0,0,1); white=rgb(1,1,1)
	RtoWrange<-colorRampPalette(c(red, white ) )
	WtoGrange<-colorRampPalette(c(white, green) )
	p <- p + scale_fill_gradient2(low=RtoWrange(100), mid=WtoGrange(100), high="gray")

view raw gistfile1.r hosted with ❤ by GitHub

Doing this produces:

If there are values close to 1 or to -1, those will pop out visually. Values close to 0 are a lot more muted.

Hope that helps someone.

References: Using R: Correlation Heatmap with ggplot2

R - Simple Recursive XML Parsing

This is intended for those who are starting out in R and interested in parsing an XML document recursively. It uses DT Lang's XML package.

If you want to just read certain types of nodes, then XPATH is great. This document by DT Lang is perfect for that.

However, if you want to read the whole document, then you have to recursively visit every node. Here's the way I ended up doing it. The generic function visitNode could be useful if you are just starting out reading XML in R.

	library(XML)
	#Recursive Function to visit the XML tree (depth first)
	visitNode <- function(node) {
	if (is.null(node)) {
	#leaf node reached. Turn back
	return()
	}
	print(paste("Node: ", xmlName(node)))
	num.children = xmlSize(node)
	if(num.children == 0 ) {
	# Add your code to process the leaf node here
	print( paste(" ", xmlValue(node)))
	}
	#Go one level deeper
	for (i in 1 : num.children) {
	visitNode(node[[i]]) #the i-th child of node
	}
	}
	xmlfile <- "books.xml"
	#read the XML tree into memory
	xtree <- xmlInternalTreeParse(xmlfile)
	root <- xmlRoot(xtree)
	visitNode(root)

view raw visitNode.R hosted with ❤ by GitHub

The full code, along with a sample XML file to test it is here.

Simulating Allele Counts in a population using R

This post is inspired by the Week 7 lectures of the Coursera course "Introduction to Genetics and Evolution" (I highly recommend this course for anyone interested in genetics, BTW.) Professor Noor uses a Univ Washington software called AlleleA1 for trying out scenarios.

We can just as well use R to get an intuitive feel for how Alleles and Genotypes propagate or die out in populations.

Basic Scenario

There are N individuals in an isolated island. Say, we are interested in two specific Alleles (Big "A", or small "a"). This in turn means that they can have 3 types of genotypes: AA, Aa or aa. The individuals mate in pairs, and produce two offspring and die out. (Thus the total population remains the same generation after generation.)
The genotype of the offspring depends on those of the parents. A 'gamete' has only one parental allele, depending on what the parent's genotype was. AA type parent can only product gamete type A, aa parent can only produce gamete type a, but Aa can produce either type of gamete.

A Punnett square of parents gametes to offspring's genotypes. 

  | A  | a

 ---------------- 

A | AA | Aa 

a | Aa | aa 

With these simple rules, we can use R Simulation scripts to observe what happens to the Allele Frequencies over generations. (The goal here is to learn to use R for Monte Carlo simulations.)

Writing the R Script from scratch

 I toyed around with the idea of using character strings for the genotypes and the alleles. But then I realized that are only three types and I could just as easily represent them with the numbers 1, 2, 3 as a simple R vector.

	# Genotype AA = 1, Aa = 2 and aa = 3
	# Gamete A =1 and a = 0
	#Uniformly distribute AA, Aa and aa
	start.pop <- sample(1:3, kStartPop, replace=TRUE)

view raw gistfile1.r hosted with ❤ by GitHub

With that done, we can write very simple functions for the procreation process.

	#Given a parent individual, get one of their Alleles in the gamete
	getGamete <- function(indiv) {
	if (indiv == 1) return(1) #AA
	if (indiv == 3) return(0) #aa
	#if Parent is Aa, the gamete is one binomial trial with prob.big.A
	if (indiv == 2) return(rbinom(1, size=1, prob.big.A)) #Aa
	}
	#Two parental Gametes combine to form a zygote
	combineGametes <- function(mg,dg) {
	if ((mg == 1) && (dg == 1)) return(1) #AA
	else if ((mg == 0) && (dg == 0)) return(3) #aa
	else return(2) #Aa
	}
	#Given two individuals, get an offspring for next generation
	getOffspring <- function(mom, dad) {
	momGam <- getGamete(mom)
	dadGam <- getGamete(dad)
	return(combineGametes(momGam, dadGam))
	}

view raw gistfile1.r hosted with ❤ by GitHub

With these useful functions, we can take one generation and produce another, 2 offspring for each set of 2 parents.

	getNextGen <- function(x) {
	nextgen <- list()
	for(i in seq(1, kStartPop, by=2)) {
	firstborn <- getOffspring(x[i], x[i+1])
	secondborn <- getOffspring(x[i], x[i+1])
	nextgen[i] <- firstborn
	nextgen[i+1] <- secondborn
	}
	return(unlist(nextgen))
	}

view raw gistfile1.r hosted with ❤ by GitHub

Putting it all together to generate multiple trials:

	simulationOneTrial <- function(start.pop, knumGenerations, trial.index) {
	df.allele <- NULL
	df.gen <- NULL
	#Keep track of the individuals in each generation
	df.gen <- rbind(df.gen, start.pop)
	x <- start.pop
	for ( gen in 1:knumGenerations) {
	#print(unlist(nex))
	#print(paste("Gen:", gen))
	numA <- calcAllelesInGeneration(x)
	df.gen <- rbind(df.gen, x)
	df.allele <- rbind(df.allele, c(gen, trial.index, numA))
	nex <- getNextGen(x)
	x <- sample(nex) #shuffle the population for breeding
	}

view raw gistfile1.r hosted with ❤ by GitHub

We also need to compute the Allele counts for each generation, and for plotting I use ggplot.

	#Bookkeeping
	calcAllelesInGeneration <- function(x) {
	AA = sum(x==1)
	AB = sum(x==2)
	BB = sum(x==3)
	#print(unlist(x))
	#print(c(AA, AB, BB))
	num.A <- (2 * AA) + AB
	num.B <- (2 * BB) + AB
	return(c(num.A, num.B))
	}
	###plotting function
	plotAllelesWithTime <- function(df, num.trials) {
	colorRange<-colorRampPalette(c(rgb(0,0,1), rgb(1,0.7,0) ))
	p <- ggplot(df, aes(x= Generation, y= value, group=Trial, color=factor(Trial))) + geom_line()
	p <- p + scale_colour_manual(values = colorRange(num.trials),
	name="Trial")
	p <- p + labs(title = "Allele A Frequencies Across Generations (2011) ")
	p <- p + ylab("Number of \"A\" Allele in the Population")
	return(p)
	}

view raw gistfile1.r hosted with ❤ by GitHub

Using this simple Monte Carlo "toy" we can develop quite a bit of intuition.



For small starting populations, either the big A or the small a allele takes over the entire population fairly quickly. Given large enough number of generations, invariably one of the alleles gets wiped out.














As one example, we can see that even a small increase in the probability of Allele A to be 0.53 (up from 0.5) makes it take over quite dramatically.

Conversely, setting it to any value under 0.5 means that the Big A allele gets wiped out of the entire population.

The entire R script can be found here. You can download the code and try playing with various starting scenarios, changing the starting population counts, generations and probabilities.


References:

1. Coursera.org (Introduction to Genetics and Evolution by Md. Noor, Week 7 lectures) 
2. AlleleA1 software at Univ Washington