inst/testScripts/replication/chipTypes/GenomeWideSNP_6/11.justSNPRMA_vs_oligo.R

In HenrikBengtsson/aroma.affymetrix: Analysis of Large Affymetrix Microarray Data Sets

###########################################################################
# Replication test
#
# Description:
# This test verifies that aroma.affymetrix can reproduce the SNPRMA
# chip-effect estimates as estimated by oligo.
# It verifies that they give the same results whether or not one
# is normalizing towards the HapMap reference (as defined by oligo).
#
# Author: Henrik Bengtsson
# Created: 2008-12-04
# Last modified: 2010-09-01
###########################################################################
library("aroma.affymetrix")
library("oligo")
verbose <- Arguments$getVerbose(-8, timestamp=TRUE)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Local functions
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
compareESets <- function(eSet1, eSet2, tolerance=1e-4) {
  # Default
  res <- all.equal(eSet1, eSet2, tolerance=tolerance)
  if (!isTRUE(res)) return(res)

  # Compare assayData()
  data1 <- assayData(eSet1)
  data2 <- assayData(eSet2)
  data1 <- as.list(data1)
  data2 <- as.list(data2)
  res <- all.equal(data1, data2, tolerance=tolerance)
  if (!isTRUE(res)) return(res)
  TRUE
} # compareESets()


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Setup
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
dataSet <- "GSE13372,testset"
chipType <- "GenomeWideSNP_6"

cdf <- AffymetrixCdfFile$byChipType(chipType)
csR <- AffymetrixCelSet$byName(dataSet, cdf=cdf)

# Process only a subset of the arrays.  Since this data set
# contains many replicates (cf. GEO), they need to be for
# different samples.
sampleNamesMap <- c(
  GSM337641="HCC1143_GLEYS_A02",
#  GSM337646="HCC1143_TRIBE_H11",
  GSM337662="HCC1143BL_GLEYS_A01",
#  GSM337666="HCC1143BL_TRIBE_D02",
#  GSM337668="HCC1143BL_GHATS_H04",
#  GSM337674="HCC1143BL_TRIGS_G07",
  GSM337683="HCC1954_GLEYS_B02",
#  GSM337688="HCC1954_TRIBE_G12",
  GSM337696="HCC1954BL_GLEYS_B01",
#  GSM337700="HCC1954BL_TRIBE_B01",
#  GSM337702="HCC1954BL_GHATS_G10",
#  GSM337703="HCC1954BL_TRIGS_G11",
  GSM337707="NCI-H2347",
  GSM337708="NCI-H2347BL"
)
sampleNames <- names(sampleNamesMap)

sampleNames <- rev(sampleNames)
csR <- csR[sampleNames]
setFullName(csR, sprintf("%s,crlmmSubset", dataSet))
print(csR)

# Assert that the file header of the first CEL file in
# truly a GenomeWideSNP_6 (and not GenomeWideEx_6 because
# oligo::justSNPRMA() cannot handle chip type aliases)
stopifnot(getHeader(csR[[1]])$chiptype == chipType)

# Assert that the correct Platform Design package is installed
pdPkgName <- cleanPlatformName(chipType)
library(pdPkgName, character.only=TRUE)

for (normalizeToHapmap in c(FALSE, TRUE)) {
  verbose && enter(verbose, sprintf("justSNPRMA(..., normalizeToHapMap=%s)", normalizeToHapmap))

  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  # SNPRMA according to aroma.affymetrix
  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  eSet <- justSNPRMA(csR, normalizeToHapmap=normalizeToHapmap, verbose=verbose)
  print(eSet)


  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  # SNPRMA according to oligo
  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  eSet0 <- justSNPRMA(getPathnames(csR), normalizeToHapmap=normalizeToHapmap, verbose=as.logical(verbose))
  print(eSet0)

  if (!normalizeToHapmap) {
    # CLEAN UP: justSNPRMA() stores a target distribution file
    # in the working directory that we don't need
    filename <- sprintf("%s.quantileReference.rda", pdPkgName)
    if (isFile(filename)) {
      file.remove(filename)
    }
  }


  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  # Compare
  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  res <- compareESets(eSet, eSet0, tolerance=0.02)
  verbose && print(verbose, res)

  # Sanity check
  stopifnot(res)

  verbose && exit(verbose)
} # for (normalizeToHapmap ...)