###---------------------------------------------------------------------
## NimblegenUtilityFunctions.R
## Utility functions for reading and processing NimbleGen data
## from NimbleGen pair files.
##
## Steven McKinney
## Molecular Oncology and Breast Cancer Program
## British Columbia Cancer Research Centre
## Oct 19, 2006
## smckinney@bccrc.ca
##
## Most credit goes to
## Natalie Thorne (snapCGH developer, currently in Dept of Oncology,
##                 University of Cambridge))
## Todd Richmond (Manager of Research Informatics, NimbleGen)
## Natalie Thorne and Todd Richmond provided most of the functions
## to allow us to get up and running in R for NimbleGen data.
## Many thanks to them.
##
## Note:  NimbleGen's business model in part aims to stay flexible, to
## allow customization of chip features when necessary.  NimbleGen
## therefore does not publish data or software specifications, and
## reserves the right to change file and data formats.  It is thus
## difficult to produce code similar to that available for the
## affy chips, as NimbleGen specifications can change.
## You should thus carefully check that your data formats
## work with this code, and modify these scripts and functions
## accordingly.
## These scripts and functions come with no guarantee whatsoever.
## However, help is available - various Bioconductor and NimbleGen
## personnel certainly bent over backwards to help us out, for
## which we are most appreciative.
###---------------------------------------------------------------------


"readSampleKey" <-
  function(file = "SampleKey.txt",
           path = NULL,
           sep = "\t",
           row.names = "SAMPLE_LABEL",
           quote = "\"",
           ...)
{

  ## Find the sample key file.
  ## Nimblegen file naming convention is of the form
  ## ORDER_ID3051_SampleKey.txt
  if (is.null(path)) path <- getwd()
  ## Get a list of files in the <path> directory.
  filenames.vec <- list.files(path)
  ## Use <file> as a mask
  SampleKeyPosn <- grep(file, filenames.vec)
  if (length(SampleKeyPosn) == 0) {
    msg <- paste("\nNo sample key file found containing string '",
                 file,
                 "'",
                 "\nin directory\n",
                 path,
                 "\n\n",
                 sep = "")
    stop(msg)
  }
  if (length(SampleKeyPosn) > 1) {
    msg <- paste("\nMultiple sample key files found containing string '",
                 file,
                 "'",
                 "\nin directory\n",
                 path,
                 "\n\n",
                 sep = "")
    stop(msg)
  }
  file <- filenames.vec[SampleKeyPosn]
  
  if (!is.null(path)) file <- file.path(path, file)
  SampleKeyHeader <- readSampleKeyHeader(file)

  if(is.null(row.names)) {
    ## Need to handle different experiment scenario where
    ## same reference sample is used in copy number comparisons.
    tab <- read.table(file, header = TRUE, as.is = TRUE, sep = sep, 
                      quote = quote, fill = TRUE, ...)

  } else {
    ## Natalie Thorne's original read table for her style data.
    row.num <- which(row.names == SampleKeyHeader)
    
    tab <- read.table(file, header = TRUE, as.is = TRUE, sep = sep, 
                      quote = quote, fill = TRUE, row.names = row.num,...)
    ##row.names(tab) <- removeExt(tab[,row.num])
  }
  tab
}


"readSampleKeyHeader" <- function(file) 
{
    con <- file(file, "r")
    on.exit(close(con))
    out <- list()
    i <- 0
    txt <- readLines(con, n = 1)
    if (txt == "") 
      stop("Not SampleKey data file: input stopped at blank line")
    as.vector(unlist(strsplit(txt,split="\t")))        
}


"NimblegenTargetFiles" <- function(samplekey, path = NULL) {

  chipid <- unique(samplekey[,which("CHIP_ID" == colnames(samplekey))])
  cy3 <- paste(chipid,"532",sep="_")
  cy3 <- paste(cy3,"pair",sep=".")
  cy5 <- paste(chipid,"635",sep="_")
  cy5 <- paste(cy5,"pair",sep=".")
  if(!is.null(path)){
    cy3 <- paste(path,cy3,sep="/")
    cy5 <- paste(path,cy5,sep="/")
  }
  cbind(cy3,cy5)
}

"NimblegenTargets" <-
  function(samplekey,
           path = NULL)
{
  ## Get sample information from sample key file.
  
  ## CHIP_ID: All should occur in pairs.
  unique.chipids <- unique(samplekey[, "CHIP_ID"])
  if ( !all(table(samplekey[, "CHIP_ID"]) == 2) ) {
    msg <- paste("The following chip IDs do not occur in pairs:\n",
                 labels(which(table(samplekey[, "CHIP_ID"]) != 2)),
                 "\n",
                 sep = "")
    stop(msg)
  }
  
  n.chips <- length(unique.chipids)
  out <- matrix("", ncol = 5, nrow = n.chips)
  colnames(out) <- c(
                     "Array",
                     "OID",
                     "Tissue",
                     "Cy3",
                     "Cy5"
                     )
  
  for ( i in seq(along = unique.chipids)) {
    chipid.i <- unique.chipids[i]
    chipid.i.rows <- which(samplekey[, "CHIP_ID"] == chipid.i)
    out[i, "Array"] <- chipid.i
    out[i, "OID"] <- samplekey[chipid.i.rows[1], "ORD_ID"]
    out[i, "Tissue"] <-
      unlist(strsplit(samplekey[chipid.i.rows[1], "SAMPLE_DESCRIPTION"], "-"))[1]
    out[i, "Cy3"] <-
      samplekey[which(samplekey[, "CHIP_ID"] == chipid.i &
                      tolower(samplekey[, "DYE"]) == "cy3"),
                "SAMPLE_DESCRIPTION"]
    out[i, "Cy5"] <-
      samplekey[which(samplekey[, "CHIP_ID"] == chipid.i &
                      tolower(samplekey[, "DYE"]) == "cy5"),
                "SAMPLE_DESCRIPTION"]
  }
  return(out)
}


"read.nimblegen" <-
  function (files,
            path = NULL,
            ext = NULL,
            names = NULL,
            columns = NULL, 
            wt.fun = NULL,
            verbose = TRUE,
            sep = "\t",
            quote = "\"", 
            ...) 
{
  if (is.null(dim(files))) {
    if (length(files)%%2 == 0)
      files <- matrix(files, ncol = 2, byrow = TRUE)
    else stop("Odd number of files: should be two data files for each array")
  }
  else {
    files <- as.matrix(files)
    if (ncol(files) != 2) 
      stop("Need a two column matrix of file names")
  }
  if (!is.null(ext)) 
    files <- array(paste(files, ext, sep = "."), dim(files))
  narrays <- nrow(files)
  if (is.null(columns)) 
    columns <- list(f = "PM")
  if (is.null(names)) 
    names <- removeExt(basename(files[, 1]))
  fullname <- files[1, 1]
  if (!is.null(path)) 
    fullname <- file.path(path, fullname)
  skip <- 1
  obj <- read.table(fullname, skip = skip, header = TRUE, 
                    sep = sep, quote = quote, check.names = FALSE, comment.char = "", 
                    fill = TRUE,...)
  nspots <- dim(obj)[1]
  Y <- matrix(0, nspots, narrays)
  colnames(Y) <- names
  RG <- list(R = Y, G = Y, Rb = Y, Gb = Y)
  if (!is.null(wt.fun)) 
    RG$weights <- Y
  for (i in 1:narrays) {
    fullname <- files[i, 1]
    if (!is.null(path)) 
      fullname <- file.path(path, fullname)
    obj <- read.table(fullname, skip = skip, header = TRUE, 
                      sep = sep, quote = quote, check.names = FALSE, comment.char = "", 
                      fill = TRUE,...)
    if (verbose) 
      cat(paste("Read", fullname, "\n"))
    if (i == 1) 
      RG$genes <- obj[, c("IMAGE_ID", "GENE_EXPR_OPTION", "SEQ_ID", 
                          "PROBE_ID", "POSITION", "X","Y","MATCH_INDEX","SEQ_URL")]
    ## Ensure that MATCH_INDEX alignment occurs for all files.
    if(!all.equal(RG$genes$MATCH_INDEX, obj[, "MATCH_INDEX"]))
      stop(paste("MATCH_INDEX out of sequence for file\n", fullname))
    RG$G[, i] <- obj[, columns$f]
    fullname <- files[i, 2]
    if (!is.null(path)) 
      fullname <- file.path(path, fullname)
    obj <- read.table(fullname, skip = skip, header = TRUE, 
                      sep = sep, quote = quote, check.names = FALSE, comment.char = "", 
                      fill = TRUE, ...)
    if (verbose) 
      cat(paste("Read", fullname, "\n"))
    ## Ensure that MATCH_INDEX alignment occurs for all files.
    if(!all.equal(RG$genes$MATCH_INDEX, obj[, "MATCH_INDEX"]))
      stop(paste("MATCH_INDEX out of sequence for file\n", fullname))
    RG$R[, i] <- obj[, columns$f]
    if (!is.null(wt.fun)) 
      RG$weights[, i] <- wt.fun(obj)
  }
  new("RGList", RG)
}

"NimblegenHeatmaps" <-
  function(input,
           filesuffix = "heatmap.pdf",
           ngAspectRatio = 1.024,
           ncolors = 123
           )
{
  ## NimblegenHeatmaps:  Produce heat map for each R and G component
  ##                     of input.
  ## Usage:  NimblegenHeatmaps(input = RG, filesuffix = "heatmap.pdf", ngAspectRatio = 1.024)
  ##         NimblegenHeatmaps(RG)
  ## Args:
  ##  input:  RGList data object
  ##  filesuffix:  file name suffix for output pdf files.
  ##  ngAspectRatio:  aspect ratio to produce plots with individual spots
  ##                  appearing as squares.  

  ## Verify input:
  if ( class(input) != "RGList") stop("input must be an RGList object.")
  
  ## Re-orient the data to match NimbleGen orientation
  idxmat <- cbind(input$genes$X, max(input$genes$Y) + 1 - input$genes$Y)
  
  ## Set up zmat with a border of 2 empty rows/cols
  ## so heatmap does not print on border box.
  zmat <- matrix(as.numeric(NA),
                 nrow = max(input$genes$X) + 4,
                 ncol = max(input$genes$Y) + 4)
  ## Raw data heat maps
  ## in NimbleGen orientation

  on.exit(cat("\n"))
  cat("\nProcessing heat maps...\n")
  for (iptr in seq(ncol(input$G))) {
    chipID <- dimnames(input$G)[[2]][iptr]
    filename <- paste("ChipID_", chipID, "_G_", filesuffix, sep = "")
    pdf(file = filename)
    cat("G:ChipNo = ", iptr, "(ChipID = ", chipID, ")  ")
    zmat[idxmat + 2] <-  input$G[, iptr]
    ## colors
    if(ncolors == 1) {
      low <- col2rgb("green")/255
    } else {
      low <- col2rgb("white")/255
    }
    high <- col2rgb("green")/255
    col <- rgb(seq(low[1], high[1], len = ncolors),
               seq(low[2], high[2], len = ncolors),
               seq(low[3], high[3], len = ncolors))
    image(x = (-1:770), y = (-1:1026), z = zmat,
          col = col, xlab = "X", ylab = "Y", asp = ngAspectRatio)
    title(main = paste("G Sample = ", chipID))
    dev.off()
  }
  cat("\n")
  for (iptr in seq(ncol(input$R))) {
    chipID <- dimnames(input$R)[[2]][iptr]
    filename <- paste("ChipID_", chipID, "_R_", filesuffix, sep = "")
    pdf(file = filename)
    cat("R:ChipNo = ", iptr, "(ChipID = ", chipID, ")  ")
    zmat[idxmat + 2] <-  input$R[, iptr]
    ## colors
    if(ncolors == 1) {
      low <- col2rgb("red")/255
    } else {
      low <- col2rgb("white")/255
    }
    high <- col2rgb("red")/255
    col <- rgb(seq(low[1], high[1], len = ncolors),
               seq(low[2], high[2], len = ncolors),
               seq(low[3], high[3], len = ncolors))
    image(x = (-1:770), y = (-1:1026), z = zmat,
          col = col, xlab = "X", ylab = "Y", asp = ngAspectRatio)
    title(main = paste("R Sample = ", chipID))
    dev.off()
  }
  
  
}
    

###----------------------------------------------


iplotMA <-
  function (MA,
            array = 1,
            xlab = "A",
            ylab = "M",
            main = colnames(MA)[array], 
            xlim = NULL,
            ylim = NULL,
            status,
            values,
            pch,
            col,
            cex, 
            legend = TRUE,
            zero.weights = FALSE,
            ...) 
{
  ## Plot function for MA plots with many rows.
  ## limma library plotMA function yields largely black
  ## graphs that are difficult to render in pdf readers,
  ## printers and documents when hundreds of thousands
  ## of data values are plotted.
  
  ## Need ixyplot() from package IDPmisc.
  require(IDPmisc)
  switch(class(MA),
         RGList = {
           MA <- MA.RG(MA[, array])
           array <- 1
           x <- MA$A
           y <- MA$M
           w <- MA$w
         },
         MAList = {
           x <- as.matrix(MA$A)[, array]
           y <- as.matrix(MA$M)[, array]
           if (is.null(MA$weights)) 
             w <- NULL
           else w <- as.matrix(MA$weights)[, array]
         },
         list = {
           if (is.null(MA$A) || is.null(MA$M)) 
             stop("No data to plot")
           x <- as.matrix(MA$A)[, array]
           y <- as.matrix(MA$M)[, array]
           if (is.null(MA$weights)) 
             w <- NULL
           else w <- as.matrix(MA$weights)[, array]
         },
         MArrayLM = {
           x <- MA$Amean
           y <- as.matrix(MA$coefficients)[, array]
           if (is.null(MA$weights)) 
             w <- NULL
           else w <- as.matrix(MA$weights)[, array]
         },
         matrix = {
           narrays <- ncol(MA)
           if (narrays < 2) 
             stop("Need at least two arrays")
           if (narrays > 5) 
             x <- apply(MA, 1, median, na.rm = TRUE)
           else x <- rowMeans(MA, na.rm = TRUE)
           y <- MA[, array] - x
           w <- NULL
         },
         exprSet = {
           narrays <- ncol(MA@exprs)
           if (narrays < 2) 
             stop("Need at least two arrays")
           if (narrays > 5) 
             x <- apply(MA@exprs, 1, median, na.rm = TRUE)
           else x <- rowMeans(MA@exprs, na.rm = TRUE)
           y <- MA@exprs[, array] - x
           w <- NULL
           if (missing(main)) 
             main <- colnames(MA@exprs)[array]
         },
         stop("MA is invalid object")
         )
  if (!is.null(w) && !zero.weights) {
    i <- is.na(w) | (w <= 0)
    y[i] <- NA
  }
  if (is.null(xlim)) 
    xlim <- range(x, na.rm = TRUE)
  if (is.null(ylim)) 
    ylim <- range(y, na.rm = TRUE)
  if (missing(status)) 
    status <- MA$genes$Status
  ixyplot(x, y, xlab = xlab, ylab = ylab, main = main)
  ## abline(h = 0) ## Will have to put this in ixyplot before key legend plot

  invisible()
}



ImageplotMA <- 
function (MA, array = 1, xlab = "A", ylab = "M", main = colnames(MA)[array], 
    xlim = NULL, ylim = NULL, status, values, pch, col, cex, 
    legend = TRUE, zero.weights = FALSE, ...) 
{
  ## Plot function for MA plots with many rows.
  ## limma library plotMA function yields largely black
  ## graphs that are difficult to render in pdf readers,
  ## printers and documents when hundreds of thousands
  ## of data values are plotted.
 
  ## Need ixyplot() from package IDPmisc.
  require(IDPmisc)

    switch(class(MA), RGList = {
        MA <- MA.RG(MA[, array])
        array <- 1
        x <- MA$A
        y <- MA$M
        w <- MA$w
    }, MAList = {
        x <- as.matrix(MA$A)[, array]
        y <- as.matrix(MA$M)[, array]
        if (is.null(MA$weights)) 
            w <- NULL
        else w <- as.matrix(MA$weights)[, array]
    }, list = {
        if (is.null(MA$A) || is.null(MA$M)) 
            stop("No data to plot")
        x <- as.matrix(MA$A)[, array]
        y <- as.matrix(MA$M)[, array]
        if (is.null(MA$weights)) 
            w <- NULL
        else w <- as.matrix(MA$weights)[, array]
    }, MArrayLM = {
        x <- MA$Amean
        y <- as.matrix(MA$coefficients)[, array]
        if (is.null(MA$weights)) 
            w <- NULL
        else w <- as.matrix(MA$weights)[, array]
    }, matrix = {
        narrays <- ncol(MA)
        if (narrays < 2) 
            stop("Need at least two arrays")
        if (narrays > 5) 
            x <- apply(MA, 1, median, na.rm = TRUE)
        else x <- rowMeans(MA, na.rm = TRUE)
        y <- MA[, array] - x
        w <- NULL
    }, exprSet = {
        narrays <- ncol(MA@exprs)
        if (narrays < 2) 
            stop("Need at least two arrays")
        if (narrays > 5) 
            x <- apply(MA@exprs, 1, median, na.rm = TRUE)
        else x <- rowMeans(MA@exprs, na.rm = TRUE)
        y <- MA@exprs[, array] - x
        w <- NULL
        if (missing(main)) 
            main <- colnames(MA@exprs)[array]
    }, stop("MA is invalid object"))
    if (!is.null(w) && !zero.weights) {
        i <- is.na(w) | (w <= 0)
        y[i] <- NA
    }
    if (is.null(xlim)) 
        xlim <- range(x, na.rm = TRUE)
    if (is.null(ylim)) 
        ylim <- range(y, na.rm = TRUE)
    if (missing(status)) 
        status <- MA$genes$Status
    plot(x, y, xlab = xlab, ylab = ylab, main = main, xlim = xlim, 
        ylim = ylim, type = "n", ...)
    if (is.null(status) || all(is.na(status))) {
        if (missing(pch)) 
            pch = 16
        if (missing(cex)) 
            cex = 0.3
        Image(x, y)
    }
    else {
        if (missing(values)) {
            if (is.null(attr(status, "values"))) 
                values <- names(sort(table(status), decreasing = TRUE))
            else values <- attr(status, "values")
        }
        sel <- !(status %in% values)
        nonhi <- any(sel)
        if (nonhi) 
            Image(x[sel], y[sel])
        nvalues <- length(values)
        if (missing(pch)) {
            if (is.null(attr(status, "pch"))) 
                pch <- rep(16, nvalues)
            else pch <- attr(status, "pch")
        }
        if (missing(cex)) {
            if (is.null(attr(status, "cex"))) {
                cex <- rep(1, nvalues)
                if (!nonhi) 
                  cex[1] <- 0.3
            }
            else cex <- attr(status, "cex")
        }
        if (missing(col)) {
            if (is.null(attr(status, "col"))) {
                col <- nonhi + 1:nvalues
            }
            else col <- attr(status, "col")
        }
        pch <- rep(pch, length = nvalues)
        col <- rep(col, length = nvalues)
        cex <- rep(cex, length = nvalues)
        for (i in 1:nvalues) {
            sel <- status == values[i]
            Image(x[sel], y[sel])
        }
        if (legend) {
            if (is.list(pch)) 
                legend(x = xlim[1], y = ylim[2], legend = values, 
                  fill = col, col = col, cex = 0.9)
            else legend(x = xlim[1], y = ylim[2], legend = values, 
                pch = pch, , col = col, cex = 0.9)
        }
    }
    invisible()
}
######################################################

###---------------------------------------------------------------------
## End  NimblegenUtilityFunctions.R
###---------------------------------------------------------------------
