########################################################################################################
# The Broad Institute
# SOFTWARE COPYRIGHT NOTICE AGREEMENT
# This software and its documentation are copyright 2007 by the
# Broad Institute/Massachusetts Institute of Technology.
# All rights are reserved.
#
# This software is supplied without any warranty or guaranteed support
# whatsoever. Neither the Broad Institute nor MIT can be responsible for
# its use, misuse, or functionality.
########################################################################################################
#
# Revlimid Analysis -- Library of auxiliary functions

MSIG.Gct2Frame <- function(filename = "NULL") { 
#
# Reads a gene expression dataset in GCT format and converts it into an R data frame
#
# The Broad Institute
# SOFTWARE COPYRIGHT NOTICE AGREEMENT
# This software and its documentation are copyright 2003 by the
# Broad Institute/Massachusetts Institute of Technology.
# All rights are reserved.
#
# This software is supplied without any warranty or guaranteed support
# whatsoever. Neither the Broad Institute nor MIT can be responsible for
# its use, misuse, or functionality.

   ds <- read.delim(filename, header=T, sep="\t", skip=2, row.names=1, blank.lines.skip=T, comment.char="", as.is=T, na.strings = "")
   descs <- ds[,1]
   ds <- ds[-1]
   row.names <- row.names(ds)
   names <- names(ds)
   return(list(ds = ds, row.names = row.names, descs = descs, names = names))
}

MSIG.ReadClsFile <- function(file = "NULL") { 
#
# Reads a class vector CLS file and defines phenotype and class labels vectors (numeric and character) for the samples in a gene expression file (RES or GCT format)
#
# The Broad Institute
# SOFTWARE COPYRIGHT NOTICE AGREEMENT
# This software and its documentation are copyright 2003 by the
# Broad Institute/Massachusetts Institute of Technology.
# All rights are reserved.
#
# This software is supplied without any warranty or guaranteed support
# whatsoever. Neither the Broad Institute nor MIT can be responsible for
# its use, misuse, or functionality.

      cls.cont <- readLines(file)
      num.lines <- length(cls.cont)
      class.list <- unlist(strsplit(cls.cont[[3]], " "))
      s <- length(class.list)
      t <- table(class.list)
      l <- length(t)
      phen <- vector(length=l, mode="character")
      class.v <- vector(length=s, mode="numeric")
     
      current.label <- class.list[1]
      current.number <- 1
      class.v[1] <- current.number
      phen[1] <- current.label
      phen.count <- 1

      if (length(class.list) > 1) {
         for (i in 2:s) {
             if (class.list[i] == current.label) {
                  class.v[i] <- current.number
             } else {
                  phen.count <- phen.count + 1
                  current.number <- current.number + 1
                  current.label <- class.list[i]
                  phen[phen.count] <- current.label
                  class.v[i] <- current.number
             }
        }
       }
     return(list(phen = phen, class.v = class.v, class.list = class.list))
}
MSIG.HeatMapPlot.3 <- function(
V, 
row.names = "NA", 
col.labels = "NA", 
col.classes = "NA", 
phen.cmap = "NA", 
col.names = "NA", 
main = " ", 
sub = " ", 
xlab=" ", 
ylab=" ",
row.norm = TRUE,
char.rescale = 1.0,                               
cmap.type = 1,   # 1 = vintage pinkogram, 2 = scale of blues, 3 = high-resolution pinkogram for probabilities [0, 1], 4 = high-resolution pinkogram for general values, 5 = color map for normalized enrichment scores, 6 = GenePattern color map
max.v = "NA")
{
#
# Plots a heatmap "pinkogram" of a gene expression matrix including phenotype vector and gene, sample and phenotype labels
#
# The Broad Institute
# SOFTWARE COPYRIGHT NOTICE AGREEMENT
# This software and its documentation are copyright 2003 by the
# Broad Institute/Massachusetts Institute of Technology.
# All rights are reserved.
#
# This software is supplied without any warranty or guaranteed support
# whatsoever. Neither the Broad Institute nor MIT can be responsible for
# its use, misuse, or functionality.

       n.rows <- length(V[,1])
       n.cols <- length(V[1,])

       if ((cmap.type == 3) | (cmap.type == 5)) {
          row.norm <- F
       }

       if (row.norm == TRUE) {
          row.mean <- apply(V, MARGIN=1, FUN=mean)
          row.sd <- apply(V, MARGIN=1, FUN=sd)
          row.n <- length(V[,1])
          for (i in 1:n.rows) {
	     if (row.sd[i] == 0) {
    	         V[i,] <- 0
             } else {
	         V[i,] <- (V[i,] - row.mean[i])/(0.333 * row.sd[i])
             }
             V[i,] <- ifelse(V[i,] < -6, -6, V[i,])
             V[i,] <- ifelse(V[i,] > 6, 6, V[i,])
          }
        }

        if (cmap.type == 1) { 
             mycol <- c("#0000FF", "#4040FF", "#7070FF", "#8888FF", "#A9A9FF", "#D5D5FF", "#EEE5EE", "#FFAADA", "#FF9DB0", "#FF7080", 
                        "#FF5A5A", "#FF4040", "#FF0D1D") # blue-pinkogram colors. This is the 1998-vintage, pre-gene cluster, original pinkogram color map
        } else if (cmap.type == 2) {
             mycol <- c(
                        "#FCFBFD","#F4F2F8","#F8F7FB","#EFEDF5","#E1E1EF","#E8E7F2","#DADAEB","#C6C7E1","#D0D1E6","#BCBDDC","#A8A6CF",
                        "#B2B2D6","#9E9AC8","#8A87BF","#9491C4","#807DBA","#7260AB","#796FB3","#6A51A3","#5C3596","#63439D","#54278F","#460D83","#4D1A89","#3F007D")
        } else if ((cmap.type == 3) | (cmap.type == 4) | (cmap.type == 5)) {
            mycol <- vector(length=512, mode = "numeric")

            for (k in 1:256) {
               mycol[k] <- rgb(255, k - 1, k - 1, maxColorValue=255)
            }
            for (k in 257:512) {
               mycol[k] <- rgb(511 - (k - 1), 511 - (k - 1), 255, maxColorValue=255)
            }
            mycol <- rev(mycol)
        } else if (cmap.type == 6) {

            mycol <- c("#4500AD", "#2700D1", "#6B58EF", "#8888FF", "#C7C1FF", "#D5D5FF", "#FFC0E5", "#FF8989", "#FF7080", "#FF5A5A", "#EF4040", "#D60C00") # blue-pinkogram colors. This is the GenePattern ComparativeMarker Selection pinkogram color map
        }

        ncolors <- length(mycol)
        mycol <- c(mycol, phen.cmap[1:length(col.classes)])

       if (cmap.type == 5) {
           if (max.v == "NA") {
              max.v <- max(max(V), -min(V))
            }
           V <- ceiling(ncolors * (V - (- max.v))/(1.001*(max.v - (- max.v))))
       } else {
           V <- ceiling(ncolors * (V - min(V))/(1.001*(max(V) - min(V))))
       }

       heatm <- matrix(0, nrow = n.rows + 1, ncol = n.cols)
       heatm[1:n.rows,] <- V[seq(n.rows, 1, -1),]
       heatm[n.rows + 1,] <- ncolors + col.labels

       if (cmap.type == 2) {
           par(mar = c(3, 7, 3, 1))
       } else {
           par(mar = c(4, 15, 4, 1))
       }

        print(c("range=", range(V)))
        if (cmap.type == 5) {
           image(1:n.cols, 1:(n.rows + 1), t(heatm), zlim = c(0, ncolors + max(col.labels)), col=mycol, axes=FALSE, main=main, sub = sub, xlab= xlab, ylab=ylab)
         } else {
           image(1:n.cols, 1:(n.rows + 1), t(heatm), col=mycol, axes=FALSE, main=main, sub = sub, xlab= xlab, ylab=ylab)
         }

       if (row.names[1] != "NA") {
            numC <- nchar(row.names)
            size.row.char <- char.rescale*30/(n.rows + 15)
            size.col.char <- char.rescale*30/(n.cols + 15)
            for (i in 1:n.rows) {
               row.names[i] <- substr(row.names[i], 1, 30)
            }
            row.names <- c(row.names[seq(n.rows, 1, -1)], "Class")
            axis(2, at=1:(n.rows + 1), labels=row.names, adj= 0.5, tick=FALSE, las = 1, cex.axis=size.row.char, font.axis=2, line=-1)
        }

        if (col.names[1] != "NA") {
           axis(1, at=1:n.cols, labels=col.names, tick=FALSE, las = 3, cex.axis=size.col.char, font.axis=2, line=-1)
        }

	return()

     }

write.gct <- function (gct.data.frame, descs = "", filename) 
{
    f <- file(filename, "w")
    cat("#1.2", "\n", file = f, append = TRUE, sep = "")
    cat(dim(gct.data.frame)[1], "\t", dim(gct.data.frame)[2], "\n", file = f, append = TRUE, sep = "")
    cat("Name", "\t", file = f, append = TRUE, sep = "")
    cat("Description", file = f, append = TRUE, sep = "")

    names <- names(gct.data.frame)
    cat("\t", names[1], file = f, append = TRUE, sep = "")

    if (length(names) > 1) {
       for (j in 2:length(names)) {
           cat("\t", names[j], file = f, append = TRUE, sep = "")
       }
     }
    cat("\n", file = f, append = TRUE, sep = "\t")

    oldWarn <- options(warn = -1)
    m <- matrix(nrow = dim(gct.data.frame)[1], ncol = dim(gct.data.frame)[2] +  2)
    m[, 1] <- row.names(gct.data.frame)
    if (length(descs) > 1) {
        m[, 2] <- descs
    } else {
        m[, 2] <- row.names(gct.data.frame)
    }
    index <- 3
    for (i in 1:dim(gct.data.frame)[2]) {
        m[, index] <- gct.data.frame[, i]
        index <- index + 1
    }
    write.table(m, file = f, append = TRUE, quote = FALSE, sep = "\t", eol = "\n", col.names = FALSE, row.names = FALSE)
    close(f)
    options(warn = 0)

}
