# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #01 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - library(aroma.affymetrix) cdf <- AffymetrixCdfFile$fromChipType("Mapping50K_Hind240") print(cdf) getPathname(cdf) getPath(cdf) getFilename(cdf) getFileSize(cdf) # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #02 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - library(aroma.affymetrix) cdf <- AffymetrixCdfFile$fromChipType("Mapping50K_Hind240") nbrOfUnits(cdf) nbrOfCells(cdf) nbrOfRows(cdf) nbrOfColumns(cdf) AffymetrixCdfFile unitNames <- getUnitNames(cdf) unitNames <- getUnitNames(cdf, units=101:120) unit <- indexOf(cdf, "SNP_A-1683825") data <- readUnits(cdf, units=121) names(data) names(data[[1]]$groups) idxs <- getCellIndices(cdf, units=121) # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #03 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", chipType="Mapping50K_Hind240") print(cs) getName(cs) getPath(cs) getNames(cs) getFullNames(cs) getFileSize(cs) cdf <- getCdf(cs) AffymetrixCelSet # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #04 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - as.list(cs) cf <- getFile(cs, 1) print(cf) # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #05 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - library(aroma.affymetrix) cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", chipType="Mapping50K_Hind240") plotDensity(cs, stratifyBy="pm") # ~2s/array # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #06 [Requires: EBImage] # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - library(aroma.affymetrix) cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", chipType="Mapping50K_Hind240") cf <- getFile(cs, 1) print(cf) plotImage(cf) # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #07 [Requires: EBImage] # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - library(aroma.affymetrix) cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", chipType="Mapping50K_Hind240") ae <- ArrayExplorer(cs) setColorMaps(ae, "sqrt,yellow") process(ae, verbose=-3) display(ae) # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #08 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - library(aroma.affymetrix) cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", chipType="Mapping50K_Hind240") qn <- QuantileNormalization(cs) print(qn) csQN <- process(qn, verbose=TRUE) plotDensity(csQN, stratifyBy="pm") # ~2s/array gc() # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #09 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - quit("no") # Restart R library(aroma.affymetrix) cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", chipType="Mapping50K_Hind240") qn <- QuantileNormalization(cs) csQN <- process(qn, verbose=TRUE) print(csQN) # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #10 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - library(aroma.affymetrix) cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", chipType="Mapping50K_Hind240") qn <- QuantileNormalization(cs) csQN <- process(qn, verbose=TRUE) plm <- RmaCnPlm(csQN, combineAlleles=TRUE, mergeStrands=TRUE) ces <- getChipEffectSet(plm) theta <- readUnits(ces, units=10) print(theta) fit(plm, units=10) theta <- readUnits(ces, units=10) print(theta) # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #11 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - fit(plm, units=1:1000, verbose=TRUE) fit(plm, units=1:1000, verbose=TRUE) # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #12 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - library(aroma.affymetrix) cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", chipType="Mapping50K_Hind240") qn <- QuantileNormalization(cs) csQN <- process(qn, verbose=TRUE) plm <- RmaCnPlm(csQN, combineAlleles=TRUE, mergeStrands=TRUE) fit(plm, verbose=TRUE) # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #13 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - cdf <- AffymetrixCdfFile$fromChipType("Mapping50K_Hind240") # dChip genome information files gi <- getGenomeInformation(cdf) print(gi) units <- getUnitsOnChromosome(gi, 16) x <- getPositions(gi, units) plot(density(x/1e6, from=0, adjust=0.2)) # dChip SNP information files si <- getSnpInformation(cdf) print(si) # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #14 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - library(aroma.affymetrix) cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", chipType="Mapping50K_Hind240") qn <- QuantileNormalization(cs) csQN <- process(qn, verbose=TRUE) ces <- getChipEffectSet(plm) theta <- extractMatrix(ces) thetaAvg <- rowMedians(theta) M <- log2(theta/thetaAvg) cdf <- getCdf(ces) si <- getSnpInformation(cdf) fl <- getFragmentLengths(si) layout(matrix(1:6, ncol=3, byrow=TRUE)) for (kk in seq(ces)) { smoothScatter(fl, M[,kk], ylim=c(-1.5,1.5), main=kk) lines(smooth.spline(fl, M[,kk]), col="red", lwd=3) } # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #15 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - fln <- FragmentLengthNormalization(ces) cesN <- process(fln, verbose=TRUE) theta <- extractMatrix(cesN) thetaAvg <- rowMedians(theta) M <- log2(theta/thetaAvg) layout(matrix(1:6, ncol=3, byrow=TRUE)) for (kk in seq(ces)) { smoothScatter(fl, M[,kk], ylim=c(-1.5,1.5), main=kk) lines(smooth.spline(fl, M[,kk]), col="orange", lwd=3) } # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #16 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - glad <- GladModel(cesN) fit(glad, arrays=1, chromosomes=19, verbose=TRUE) plot(glad, arrays=1, chromosomes=19) plot(glad, arrays=1, chromosomes=18) # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #17 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - glad <- GladModel(cesN) ce <- ChromosomeExplorer(glad) print(ce) process(ce, arrays=1, chromosomes=1:4, verbose=TRUE) display(ce) # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #18 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - library(aroma.affymetrix) cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", chipType="Mapping50K_Hind240") qn <- QuantileNormalization(cs) csQN <- process(qn, verbose=TRUE) plm <- RmaCnPlm(csQN, combineAlleles=TRUE, mergeStrands=TRUE) fit(plm, verbose=TRUE) ces <- getChipEffectSet(plm) fln <- FragmentLengthNormalization(ces) cesN <- process(fln, verbose=TRUE) glad <- GladModel(cesN) ce <- ChromosomeExplorer(glad) process(ce, verbose=-3) # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #19 # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - library(aroma.affymetrix) chipTypes <- c("Mapping50K_Hind240", "Mapping50K_Xba240") cesList <- list() for (chipType in chipTypes) { cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", chipType=chipType) print(cs) qn <- QuantileNormalization(cs) csQN <- process(qn, verbose=TRUE) plm <- RmaCnPlm(csQN, combineAlleles=TRUE, mergeStrands=TRUE) fit(plm, verbose=TRUE) ces <- getChipEffectSet(plm) fln <- FragmentLengthNormalization(ces) cesN <- process(fln, verbose=TRUE) cesList[[chipType]] <- cesN } glad <- GladModel(cesList) ce <- ChromosomeExplorer(glad) process(ce, verbose=-3) # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # Code chunk #20 [Requires: EBImage] # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - library(aroma.affymetrix) cs <- AffymetrixCelSet$byName("HapMap100K-bioc2007", chipType="Mapping50K_Hind240") qn <- QuantileNormalization(cs) csQN <- process(qn, verbose=TRUE) plm <- RmaPlm(csQN) fit(plm, verbose=TRUE) rs <- calculateResidualSet(plm, verbose=TRUE) # ~5 sec/array rf <- getFile(rs, 1) plotImage(rf) ae <- ArrayExplorer(rs) setColorMaps(ae, "log2,log2abs,rainbow") process(ae, verbose=-3)