###---------------------------------------------------------------------
## processNimblegenData.R:  Example script for reading in and
## processing data from Nimblegen chips.
## This script used to process chips using the XVGA (768 x 1024)
## 1:2 feature density format, (probes at every 2nd location, approx
## 380000 data values per sample per chip).
##
## 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.
###---------------------------------------------------------------------


## Load required R/BioConductor libraries
library(snapCGH)
library(limma)
library(IDPmisc)


## Assuming that NimblegenUtilityFunctions.R is in working
## directory.  Add appropriate file path otherwise.
source("NimblegenUtilityFunctions.R")

## NimRootDir should point to the directory containing
## all the NimbleGen data material as delivered on CD
## or DVD.  This root directory will contain a sample key
## file, e.g. ORDER_ID3051_SampleKey.txt, and a subdirectory
## "Raw_data_files" containing the 'pair' data files with names
## such as "63392_532.pair".
## You may need to specify a different path here, if your
## NimbleGen data is not in the subdirectory "Raw_data_files"
## of the R working directory.
NimRootDir <- getwd()

## readSampleKey will by default look for a file whose name
## contains "SampleKey.txt", e.g. ORDER_ID3051_SampleKey.txt.
## If your NimbleGen order comes split over several CDs, you may have
## to manually create one SampleKey file containing all the Chip
## information for the chip data you want to process.
samplekey <- readSampleKey(path = NimRootDir,  row.names = NULL)

## NimblegenTargetFiles builds the names of the data files
## based on the information from the SampleKey file.
data.files <-
  NimblegenTargetFiles(samplekey,
                       path = paste(NimRootDir,
                         "/Raw_Data_Files",
                         sep=""))

## NimblegenTargets retrieves sample descriptive information, e.g.
## sample type and fluorescent labeling.
targets <- as.data.frame(NimblegenTargets(samplekey, path = NimRootDir))

## read.nimblegen() reads in all the 'pair' raw data intensity files.
RG <- read.nimblegen(files = data.files)

## Some of our order was customized, and the returned data files had
## extra quality control spot intensities in the raw data files.
## We want to examine only the spots that map to chromosomes.
## Data for these spots have PROBE_ID labels that start with "CHR",
## so we keep these and filter out the rest.
if ( length(chrIDxs <- grep("^CHR", as.character(RG$genes$PROBE_ID))) < nrow(RG) ) {
  RG <- RG[chrIDxs, ]
}
if ( exists("chrIDxs") ) rm(chrIDxs)

## Our data is from female samples, so we filter out Y chromosome
## data as well.
if( length(chrYIDxs <- grep("^CHRY", as.character(RG$genes$PROBE_ID))) ) {
  RG <- RG[-chrYIDxs, ]
}
if ( exists("chrYIDxs") ) rm(chrYIDxs)


## After filtering, we need to get rid of unnecessary factor classes
## for several of the variables.
RG$genes$SEQ_ID <- as.factor(as.character(RG$genes$SEQ_ID))
RG$genes$GENE_EXPR_OPTION <- as.factor(as.character(RG$genes$GENE_EXPR_OPTION))
RG$genes$PROBE_ID <- as.character(RG$genes$PROBE_ID)


## Produce heatmaps of the raw data.
NimblegenHeatmaps(RG)

## Produce spotmaps of the raw data to examine spot locations
## and chip features.
NimblegenHeatmaps(input = RG, filesuffix = "spotmap.pdf", ncolors = 1)

## Now we are ready to start processing the raw data.
## readPositionalInfo() sets up probe position information.
RG <- readPositionalInfo(RG, source = "nimblegen")

## Now we are ready to normalize.  NimbleGen uses qspline, and we confirmed
## that no background subtraction plus within-chip q-spline normalization
## works well for our situation.
## Set up a matrix to hold normalization data.
pmi.normmtx <- matrix(NA, nrow = nrow(RG$G), ncol = 2)

for (iptr in seq(ncol(RG$G))) {
  ## Perform q-spline normalization for each chip.
  pmi.normmtx[] <- normalize.qspline(as.matrix(cbind(RG$G[, iptr], RG$R[, iptr])))
  RG$G[, iptr] <- pmi.normmtx[, 1]
  RG$R[, iptr] <- pmi.normmtx[, 2]
}
rm(pmi.normmtx)

## Set up the 'design' component, as discussed in snapCGH
## documentation.  In our case, CY5 (Red) is the reference sample.
RG$design <- rep(-1, nrow(targets))

## Get MA form for data
MA <- MA.RG(RG)


## MA plots
## plotMA():  limma library standard MA plot routine.
## Our chips have about 380000 data values, so standard
## MA plots are mostly black.  Graphs and output files
## are large.  Adobe reader chokes on these standard plotMA plots.
pdf(file = "plotMA.pdf")
lapply(seq(ncol(MA)),
       function(x) { plotMA(MA, array = x)
                     mtext("MA after q-spline normalization") })
dev.off()

## I found MA 'heatmap' scatterplots to work well for this
## plot conumdrum.  I've set up alternative MA plot routines
## using ixyplot() from package IDPmisc.
## iplot() produces an ixyplot with heat legend.  I'd like
## to have vertical and horizontal lines in the plot, just haven't
## had the time to work them in.
pdf(file = "iplotMA.pdf")
lapply(seq(ncol(MA)), function(x) iplotMA(MA, array = x))
dev.off()

## To get e.g. horizontal line at 0, I use ImageplotMA for now.
pdf(file = "ImageplotMA.pdf")
lapply(seq(ncol(MA)), function(x) {ImageplotMA(MA, array = x); abline(h = 0)})
dev.off()


## Set up the data (e.g. put the clone data in order along the chromosomes).
MA.process <- processCGH(MA,
                         maxChromThreshold = 23, minChromThreshold = 1,
                         method.of.averaging = NULL, ID = "PROBE_ID")


## Run segmentation analysis.
SegInfo.HomHMM <- runHomHMM(MA.process, criteria = "AIC")

###---------------------------------------------------------------------
## End  processNimblegenData.R
###---------------------------------------------------------------------
