inst/testScripts/archive/system/chipTypes/Mapping250K_Nsp,Sty/31.ASCN,paired.R

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

library("aroma.affymetrix")
log <- verbose <- Arguments$getVerbose(-8, timestamp=TRUE)

dataSetName <- "Affymetrix_2006-TumorNormal"
tags <- "ACC,-XY,BPN,-XY,RMA,FLN,-XY"
chipType <- "Mapping250K_Nsp"

pairs <- matrix(c(
  "CRL-2325D", "CRL-2324D",
  "CRL-5957D", "CRL-5868D",
  "CCL-256.1D", "CCL-256D",
  "CRL-2319D", "CRL-2320D",
  "CRL-2362D", "CRL-2321D",
  "CRL-2337D", "CRL-2336D",
  "CRL-2339D", "CRL-2338D",
  "CRL-2341D", "CRL-2340D",
  "CRL-2346D", "CRL-2314D"
), ncol=2, byrow=TRUE)
colnames(pairs) <- c("normal", "tumor")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Setting up data set
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if (!exists("cesN")) {
  cdf <- AffymetrixCdfFile$byName(chipType)
  gi <- getGenomeInformation(cdf)

  cdfM <- getMonocellCdf(cdf)
  cesN <- SnpChipEffectSet$byName(dataSetName, tags=tags, cdf=cdfM)
  print(cesN)

  # Reorder arrays according to 'pairs' matrix
  cesN <- cesN[indexOf(cesN, pairs)]
  names <- getNames(cesN)
  dim(names) <- dim(pairs)
  print(names)


  Blim <- c(0,1)
  Clim <- c(0,4)
  Clab <- expression(C[T]/C[N])
  BTlab <- expression(beta[T])
  BNlab <- expression(beta[N])
  Flab <- expression(beta[T]-beta[N])
  Flim <- c(-1/2,1/2)
}

# Identify SNPs
chromosome <- 9
units <- getUnitsOnChromosome(gi, chromosome)
pos <- getPositions(gi, units=units)/1e6
o <- order(pos)
units <- units[o]
pos <- pos[o]
str(units)

isSNP <- (getUnitTypes(cdf, units=units) == 2)
units <- units[isSNP]
pos <- pos[isSNP]

data <- extractTotalAndFreqB(cesN, units=units)
dim(data) <- c(dim(data)[1:2], dim(names))
dimnames(data) <- list(NULL, c("total", "freqB"), NULL, colnames(pairs))
C <- data[,"total",,"tumor"] / data[,"total",,"normal"]
B <- data[,"freqB",,]
F <- B[,,"tumor"] - B[,,"normal"]

F2 <- 2*(1/2-abs(B[,,"normal"]-1/2))
dAA <- sqrt((B[,,"normal"]-0)^2)
dBB <- sqrt((B[,,"normal"]-1)^2)
#F[dAA < 0.2 | dBB < 0.2] <- NA


# Plot (x,C) and (x,F) along genome
devSet(2)
subplots(3, ncol=1)
par(mar=c(2,4,1,1)+0.1)
par(ask=TRUE)
for (kk in seq_len(ncol(C))) {
  plot(pos, C[,kk], ylim=Clim, ylab=Clab)
  abline(h=1, lty=3, col="red")
  stext(side=3, pos=0, paste(pairs[kk,], collapse="/"))
  stext(side=3, pos=1, sprintf("Chr %d", chromosome))
  plot(pos, B[,kk,"tumor"], ylim=Blim, ylab=BTlab)
#  smoothScatter(pos, B[,kk,"tumor"], bandwidth=c(1,0.0001), ylim=Blim, ylab=BTlab)
  plot(pos, F[,kk], cex=F2, ylim=Flim, ylab=Flab)
#  smoothScatter(pos, F[,kk], bandwidth=c(1,0.0001), ylim=Flim, ylab=Flab)
  abline(h=0, lty=3, col="red")
}


devSet(3)
BNlab <- expression(beta[N])
BTlab <- expression(beta[T])
unitNames <- getUnitNames(cdf, units=units)
col <- seq_len(ncol(C))
par(ask=TRUE)
for (kk in seq_len(nrow(C))) {
  plot(NA, xlim=Blim, ylim=Blim, xlab=BNlab, ylab=BTlab)
  stext(side=3, pos=0, unitNames[kk])
  points(B[kk,,], col=col, pch=19)
}