# Clustering genes in a time course experiment

I have seen a few queries on clustering in time series and specifically on clustering, but I don't think they answer my question.

Background: I want to cluster genes in a time course experiment in yeast. There are four time points say: t1 t2 t3 and t4 and total number of genes G. I have the data in form a matrix M in which the columns represent the treatments (or time points) t1 t2 t3 and t4 and the rows represent the genes. Therefore, M is a Gx4 matrix.

Problem: I want to cluster the genes which behave the same across all time points t1 t2 t3 and t4 as well as within a particular time point ti , where i is in {1, 2, 3, 4} (In case we cannot do both the clusterings together, the clustering within a time point is more important than clustering across time points). In addition to this, I also want to draw a heatmap.

My Solution: I use the R code below to obtain a heatmap as well as the clusters using hclust function in R (performs hierarchical clustering with euclidean distance)

    row.scaled.expr <- (expr.diff - rowMeans(expr.diff)) / rowSds(expr.diff)

breaks.expr <- c(quantile(row.scaled.expr[row.scaled.expr < 0],
seq(0,1,length=10)[-9]), 0,
quantile(row.scaled.expr[row.scaled.expr > 0],
seq(0,1,length=10))[-1] )

blue.red.expr <- maPalette(low = "blue", high = "red", mid = "white",
k=length(breaks.expr) - 1)

pdf("images/clust.pdf",
height=30,width=20,pointsize=20)
ht1 <- heatmap.2(row.scaled.expr, col = blue.red.expr, Colv = FALSE, key = FALSE,
dendrogram = "row", scale = "none", trace = "none",
cex=1.5, cexRow=1, cexCol=2,
density.info = "none", breaks = breaks.expr,
labCol = colnames(row.scaled.expr),
labRow="",
lmat=rbind( c(0, 3), c(2,1), c(0,4) ), lhei=c(0.25, 4, 0.25 ),
main=expression("Heat Map"),
ylab="Genes in the Microarray",
xlab="Treatments"
)
dev.off()


I recently discovered hopach package in Bioconductor which can be used to estimate the number of clusters. Previously, I was randomly assigning the number of bins for the heatmap and cutting the tree at an appropriate height to get a pre-specified number of clusters.

Possible Problems in my solution:

1. I may be not clustering the genes within a particular treatment and clustering genes only across treatments or vice versa.
2. There may be better ways of obtaining a heatmap for the pattern I want to see (similar genes within a treatment and across treatments).
3. There may be better visualization methods which I am not aware of.

Note:

1. csgillespie (moderator) has a more general document on his website in which he discusses all the aspects of time course analysis (including heatmaps and clustering). I would appreciate if you can point me to an articles which describe heatmaps and clustering in detail.

2. I have tried the pvclust package, but it complains that M is singular and then it crashes.

It seems you just want to make a fair standard analysis, so I am not a best person to answer your question; yet I would suggest you to dive deeper into Bioconductor; it has a lot of useful stuff, nevertheless finding what you want is painful. For instance Mfuzz package looks promising.

In complement to @mbq's response (Mfuzz looks fine), I'll just put some references (PDFs) about clustering of time-course gene expression data:

1. Futschik, ME and Charlisle, B (2005). Noise robust clustering of gene expression time-course data. Journal of Bioinformatics and Computational Biology, 3(4), 965-988.
2. Luan, Y and Li, H (2003). Clustering of time-course gene expression data using a mixed-effects model with B-splines. Bioinformatics, 19(4), 474-482.
3. Tai YC and Speed, TP (2006). A multivariate empirical Bayes statistic for replicated microarray time course data. The Annals of Statistics, 34, 2387–2412.
4. Schliep, A, Steinhoff, C, and Schönhuth, A (2004). Robust inference of groups in gene expression time-courses using mixtures of HMMs. Bioinformatics, 20(1), i283-i228.
5. Costa, IG, de Carvalho, F, and de Souto, MCP (2004). Comparative analysis of clustering methods for gene expression time course data. Genetics and Molecular Biology, 27(4), 623-631.
6. Inoue, LYT, Neira, M, Nelson, C, Gleave, M, and Etzioni, R (2006). Cluster-based network model for time-course gene expression data. Biostatistics, 8(3), 507-525.
7. Phang, TL, Neville, MC, Rudolph, M, and Hunter, L (2003). Trajectory Clustering: A Non-Parametric Method for Grouping Gene Expression Time Courses with Applications to Mammary Development. Pacific Symposium on Biocomputing, 8, 351-362.

Did you try the timecourse package (as suggested by @csgillespie in his handout)?

• The timecourse package isn't really for determining clusters, rather it's for calculating which genes are differentially expressed. – csgillespie Oct 5 '10 at 9:41
• @csgillespie (+1) Thanks. I thought it might be used to isolate genes with varying temporal profiles across biological conditions, or as a first step before using a clustering procedure (in fact, I was thinking of kml but I'm not really an expert in that domain). – chl Oct 5 '10 at 9:53
• You are correct in that you would tend use to isolate interesting genes before any clustering - basically thin down your list of genes. I suppose it does perform clustering of a sort, i.e.differentially expressed vs non-differentially expressed. – csgillespie Oct 5 '10 at 10:09

Just to add to the other answers (which look like they should solve your problem), did you try using standard clustering algorithms for you data when constructing your dendrogram? For example,

heatmap.2(dataset, <standard args>,
hclustfun = function(c){hclust(c, method= 'average')}
)


Instead of using the average distance for clustering, you can also use "ward", "single", "median", ... See ?hclust for a full list.

To extract clusters, use the hclust command directly and then use the cutree command. For example,

hc = hclust(dataset)
cutree(hc)


More details can be found at my webpage.

• . Yes, but one of the difficulties I had was in extracting the clusters from the heatmap.2 object. Is there an easy way of extracting the clusters? I am aware of the cutree command which can be used to extract clusters from the heatmap.2 object. – suncoolsu Oct 4 '10 at 17:33
• @suncoolsu: I've updated my answer. Does that help? – csgillespie Oct 4 '10 at 21:13