R/heatmapsM.R

In maigesPack: Functions to handle cDNA microarray data, including several methods of data analysis

## Function to plot heatmaps separating groups generated by SOM or k-means
##
## Parameters: data      -> numeric matrix to be ploted
##             groups    -> numeric or char vector (with length equal
##                          to columns or rows of data) with identification of
##                          the groups to be separated
##             sampleT   -> list with 2 vectors. First specifying the
##                          samples (first character of each sample) and second
##                          specifying the colors for then
##             doHier    -> plot hierarquical branch at other dimension?
##             distfun   -> function for distance calculation
##             hclustfun -> function for hierarquical cluster calculation
##             ...       -> additional parameters for image
##
## Gustavo H. Esteves
## (adapted from Elier B. Cristo)
## Data: 15/05/07
##
##


heatmapsM <- function (data, groups, sampleT=NULL, doHier=FALSE, distfun=dist,
hclustfun=hclust, ...) {
    
    
    ## Doing a basic test
    if (length(di <- dim(data)) != 2 || !is.numeric(data))
        stop("data must be a numeric matrix")
    if (di[1] <= 1 || di[2] <= 1)
        stop("data must have at least 2 rows and 2 columns")
    
    
    ## Defining the range of values
    tmp <- range(data)
    zLim <- c(-max(abs(tmp)), max(abs(tmp)))
    
    
    ## Defining the length of the figures
    if(doHier) {
        tam1 <- 0.5
        tmpMat <- matrix(c(0:(2*length(unique(groups))+1), 0, 0), 2,
        length(unique(groups))+2)
        
    }
    else {
        tam1 <- 0.1
        tmpMat <- matrix(c(0,0:(2*length(unique(groups))), 0, 0), 2,
        length(unique(groups))+2)
        
    }
    for(i in unique(groups))
        tam1 <- c(tam1, sum(groups == i)*0.06)
    
    tam1 <- c(tam1, 0.3)
    tam2 <- nrow(data)*0.07
    
    
    ## Defining the layout of the figure
    if(length(groups) == dim(data)[2])
        layout(tmpMat, widths=tam1, heights=c(0.15, 0.85))
    else if(length(groups) == dim(data)[1])
        layout(t(tmpMat), widths=c(0.15, 0.85), heights=tam1)
    else
        stop("'groups' must be a vector of length nrow(data) or ncol(data)!")
    
    
    ## ordering genes (or samples) by hierarquical cluster
    if(doHier){
        if(length(groups) == dim(data)[2]) {
            tmp <- rowMeans(data)
            hcr <- hclustfun(distfun(data))
            ddr <- as.dendrogram(hcr)
            ddr <- reorder(ddr, tmp)
            rowInd <- order.dendrogram(ddr)
        }
        else {
            tmp <- colMeans(data)
            hcr <- hclustfun(distfun(t(data)))
            ddr <- as.dendrogram(hcr)
            ddr <- reorder(ddr, tmp)
            colInd <- order.dendrogram(ddr)
        }
        
        if(length(groups) == dim(data)[2]) {
            par(mar=c(8, 1, 0, 0))
            plot(ddr, horiz=TRUE, axes=FALSE, yaxs="i", leaflab="none")
        }
        else {
            par(mar=c(0, 0, 1, 8))
            plot(ddr, axes=FALSE, xaxs="i", leaflab="none")
        }
    }
    else {
        if(length(groups) == dim(data)[2])
            rowInd <- 1:dim(data)[1]
        else
            colInd <- 1:dim(data)[2]
    }
    
    
    ## Doing the heatmaps for each group
    for(i in unique(groups)) {
        
        if(length(groups) == dim(data)[2])
            x <- data[, groups == i]
        else
            x <- data[groups == i, ]
        
        if(length(groups) == dim(data)[2]) {
            tmp <- colMeans(x)
            hcc <- hclustfun(distfun(t(x)))
            ddc <- as.dendrogram(hcc)
            ddc <- reorder(ddc, tmp)
            colInd <- order.dendrogram(ddc)
        }
        else {
            tmp <- rowMeans(x)
            hcc <- hclustfun(distfun(x))
            ddc <- as.dendrogram(hcc)
            ddc <- reorder(ddc, tmp)
            rowInd <- order.dendrogram(ddc)
        }
        
        x <- x[rowInd, colInd]
        
        if(length(groups) == dim(data)[2]) {
            par(mar=c(0, 0, 1, 1))
            plot(ddc, axes=FALSE, xaxs="i", leaflab="none")
        }
        else {
            par(mar=c(1, 1, 0, 0))
            plot(ddc, horiz=TRUE, axes=FALSE, yaxs="i", leaflab="none")
        }
        
        if(length(groups) == dim(data)[2])
            par(mar=c(8, 0, 0, 1))
        else
            par(mar=c(1, 0, 0, 8))

        image(1:ncol(x), 1:nrow(x), t(x), axes=FALSE, xlim=c(0.5,ncol(x)+0.5),
        ylim=c(0.5,nrow(x)+0.5), xlab="", ylab="", zlim=zLim, ...)
        
        
        if(length(groups) == dim(data)[2]) {
            ## putting color at sample descritions
            nlab <- ncol(x)
            if(!is.null(colnames(x)) & !is.null(sampleT)) {
                for (ilab in 1:nlab) {
                    samp <- (colnames(x))[ilab]
                    sampTmp <- substr(samp, 1, 1)
                    sampC <- sampleT[[2]][sampleT[[1]] == sampTmp]

                    axis(1, ilab, las=2, line=-0.5, tick=0, labels=samp,
                    col.axis=sampC)
                    
                }
            }
            else if(is.null(sampleT) & !is.null(colnames(x)))
                axis(1, 1:ncol(x), las=2, line=-0.5, tick=0, labels=colnames(x))
            else
                axis(1, 1:ncol(x), las=2, line=-0.5, tick=0,
                labels=(1:ncol(x))[colInd])
            
        }
        else { ## Putting gene descriptions
            if(is.null(rownames(x)))
                tmp <- (1:nrow(x))[rowInd]
            else
                tmp <- rownames(x)
            
            axis(4, 1:nrow(x), las=2, line=-0.5, tick=0, labels=tmp)
        }
        
    }
    
    
    if(length(groups) == dim(data)[2]) {
        ## Puting gene descriptions
        if(is.null(rownames(x)))
            tmp <- (1:nrow(x))[rowInd]
        else
            tmp <- rownames(x)
        axis(4, 1:nrow(x), las=2, line=-0.5, tick=0, labels=tmp)
    }
    else {
        ## putting color at sample descritions
        nlab <- ncol(x)
        if(!is.null(colnames(x)) & !is.null(sampleT)) {
            for (ilab in 1:nlab) {
                samp <- (colnames(x))[ilab]
                sampTmp <- substr(samp, 1, 1)
                sampC <- sampleT[[2]][sampleT[[1]] == sampTmp]

                axis(1, ilab, las=2, line=-0.5, tick=0, labels=samp,
                col.axis=sampC)
                
            }
        }
        else if(is.null(sampleT) & !is.null(colnames(x)))
            axis(1, 1:ncol(x), las=2, line=-0.5, tick=0, labels=colnames(x))
        
        else
            axis(1, 1:ncol(x), las=2, line=-0.5, tick=0,
            labels=(1:ncol(x))[colInd])
        
    }
    
}