################################################### ### chunk number 1: ################################################### library(aCGH) datadir <- system.file(package = "aCGH") datadir <- paste(datadir, "/examples", sep="") clones.info <- read.table(file = file.path(datadir, "clones.info.ex.txt"), header = T, sep = "\t", quote="", comment.char="") log2.ratios <- read.table(file = file.path(datadir, "log2.ratios.ex.txt"), header = T, sep = "\t", quote="", comment.char="") pheno.type <- read.table(file = file.path(datadir, "pheno.type.ex.txt"), header = T, sep = "\t", quote="", comment.char="") ex.acgh <- create.aCGH(log2.ratios, clones.info, pheno.type) ################################################### ### chunk number 2: ################################################### ex.acgh <- aCGH.process(ex.acgh, chrom.remove.threshold = 23, prop.missing = .25, sample.quality.threshold = .4, unmapScreen=TRUE, dupRemove = FALSE) ################################################### ### chunk number 3: ################################################### log2.ratios.imputed(ex.acgh) <- impute.lowess(ex.acgh, maxChrom=24) ################################################### ### chunk number 4: ################################################### data(colorectal) colorectal summary(colorectal) ################################################### ### chunk number 5: ################################################### plot(colorectal) ################################################### ### chunk number 6: ################################################### sample.names(colorectal) phenotype(colorectal)[1:4,] ################################################### ### chunk number 7: ################################################### datadir <- system.file("examples", package = "aCGH") latest.mapping.file <- file.path(datadir, "human.clones.info.Jul03.txt") ex.acgh <- aCGH.read.Sprocs(dir(path = datadir,pattern = "sproc", full.names = TRUE), latest.mapping.file, chrom.remove.threshold = 23) ex.acgh ################################################### ### chunk number 8: ################################################### plot(ex.acgh) ################################################### ### chunk number 9: ################################################### cr <- colorectal[ ,1:3] ################################################### ### chunk number 10: ################################################### plotGenome(ex.acgh, samples=2, Y = FALSE) ################################################### ### chunk number 11: ################################################### ## Determining hmm states of the clones. In the interest of time, ##we have commented this step out and used pre-computed results. ##hmm(ex.acgh) <- find.hmm.states(ex.acgh) hmm(ex.acgh) <- ex.acgh.hmm ## Merging hmm states hmm.merged(ex.acgh) <- mergeHmmStates(ex.acgh, model.use = 1, minDiff = .25) ## Calculating the standard deviations for each array. Standard error is ##calculated for each region and then averaged across regions. The final ##SDs for each samples are contained in sd.samples(exa.acgh)$madGenome. sd.samples(ex.acgh) <- computeSD.Samples(ex.acgh) ## Finding the genomic events associated with each sample using ##results of the partitioning into the states. genomic.events(ex.acgh) <- find.genomic.events(ex.acgh) ## Plotting and printing the hmm states either to the screen or into the ##postscript file. Each chromosome for each sample is plotted on a separate ##page ##postscript("hmm.states.temp.ps");plotHmmStates(ex.acgh, sample.ind=1);dev.off() ################################################### ### chunk number 12: ################################################### plotHmmStates(colorectal, sample.ind = 1, chr = 1) ################################################### ### chunk number 13: ################################################### plotFreqStat(colorectal, all = T) ################################################### ### chunk number 14: ################################################### summarize.clones(colorectal)[1:10 ,] ################################################### ### chunk number 15: ################################################### factor <- 3 tbl <- threshold.func(log2.ratios(colorectal), posThres=factor*(sd.samples(colorectal)$madGenome)) rownames(tbl) <- clone.names(colorectal) colnames(tbl) <- sample.names(colorectal) tbl[1:5,1:5] ################################################### ### chunk number 16: ################################################### col.fga <- fga.func(colorectal, factor=3,chrominfo=human.chrom.info.Jul03) cbind(gainP=col.fga$gainP,lossP=col.fga$lossP)[1:5,] ################################################### ### chunk number 17: ################################################### colnames(phenotype(colorectal)) sex <- phenotype(colorectal)$sex sex.na <- !is.na(sex) index.clones.use <- which(clones.info(colorectal)$Chrom < 23) colorectal.na <- colorectal[ index.clones.use,sex.na , keep=TRUE] dat <- log2.ratios.imputed(colorectal.na) resT.sex <- mt.maxT(dat, sex[sex.na], test = "t.equalvar", B = 1000) ################################################### ### chunk number 18: ################################################### plotFreqStat(colorectal.na, resT.sex, sex[sex.na], factor=3, titles = c("Female", "Male"), X = FALSE, Y = FALSE) ################################################### ### chunk number 19: ################################################### plotSummaryProfile(colorectal, response = sex, titles = c("Female", "Male"), X = FALSE, Y = FALSE, maxChrom = 22) ################################################### ### chunk number 20: ################################################### factor <- 3 minChanged <- 0.1 gainloss <- gainLoss(log2.ratios(colorectal)[,sex.na], cols=1:length(which(sex.na)), thres=(factor*(sd.samples(colorectal)$madGenome))[sex.na]) ind.clones.use <- which(gainloss$gainP >= minChanged | gainloss$lossP>= minChanged & clones.info(colorectal)$Chrom < 23) colorectal.na <- colorectal[ind.clones.use,sex.na, keep=TRUE] dat <- log2.ratios.imputed(colorectal.na) resT.sex <- mt.maxT(dat, sex[sex.na],test = "t.equalvar", B = 1000) ################################################### ### chunk number 21: ################################################### plotFreqStat(colorectal.na, resT.sex, sex[sex.na],factor=factor,titles = c("Male", "Female")) ################################################### ### chunk number 22: ################################################### time <- rexp(ncol(colorectal), rate = 1 / 12) events <- rbinom(ncol(colorectal), size = 1, prob = .5) surv.obj <- Surv(time, events) surv.obj stat.coxph <- aCGH.test(colorectal, surv.obj, test = "coxph", p.adjust.method = "fdr") stat.coxph[1:10 ,] ################################################### ### chunk number 23: ################################################### plotFreqStat(colorectal, stat.coxph, events, titles = c("Survived/Censored", "Dead"), X = FALSE, Y = FALSE) ################################################### ### chunk number 24: ################################################### age <- phenotype(colorectal)$age age.na <- which(!is.na(age)) age <- age[age.na] colorectal.na <- colorectal[, age.na] stat.age <- aCGH.test(colorectal.na, age, test = "linear.regression", p.adjust.method = "fdr") stat.age[1:10 ,] ################################################### ### chunk number 25: ################################################### plotFreqStat(colorectal.na, stat.age, ifelse(age < 70, 0, 1), titles = c("Young", "Old"), X = FALSE, Y = FALSE) ################################################### ### chunk number 26: ################################################### sex <- phenotype(colorectal)$sex sex.na <- !is.na(sex) index.clones.use <- which(clones.info(colorectal.na)$Chrom < 23) colorectal.na <- colorectal[ index.clones.use,sex.na , keep=TRUE] dat <- log2.ratios.imputed(colorectal.na) resT.sex <- mt.maxT(dat, sex[sex.na], test = "t.equalvar", B = 1000) sex.tbl <- summarize.clones(colorectal.na, resT.sex, sex[sex.na], titles = c("Male", "Female")) sex.tbl[1:5,] ################################################### ### chunk number 27: ################################################### par(mfrow=c(2,1)) clusterGenome(colorectal.na, response = sex[sex.na], titles = c("Female", "Male"), byclass = FALSE, showaber = TRUE, vecchrom = c(4,8,9), dendPlot = FALSE, imp = FALSE) clusterGenome(colorectal.na, response = sex[sex.na], titles = c("Female", "Male"), byclass = TRUE, showaber = TRUE, vecchrom = c(4,8,9), dendPlot = FALSE, imp = FALSE)