################################################### ### chunk number 1: set width ################################################### options(width=60) options(continue=" ") X11.options(type="xlib") ################################################### ### chunk number 2: load packages ################################################### library(GEOquery) library(eisa) library(affy) library(gcrma) library(GO.db) library(KEGG.db) library(xtable) ################################################### ### chunk number 3: the id of the data set ################################################### GEO <- "GSE9954" ################################################### ### chunk number 4: download files if needed ################################################### if (!file.exists(GEO)) { rawfiles <- getGEOSuppFiles(GEO) tarfile <- grep("\\.tar$", rownames(rawfiles), value=TRUE) untar(tarfile, exdir=GEO) } ################################################### ### chunk number 5: read in the CEL files ################################################### ABfile <- paste(GEO, sep="/", "AB.Rdata") if (!file.exists(ABfile)) { celfiles <- list.files(GEO, pattern="\\.CEL.gz", full.names=TRUE) AB <- ReadAffy(filenames=celfiles) save(AB, file=ABfile) } ################################################### ### chunk number 6: PA calls ################################################### PAfile <- paste(GEO, sep="/", "PA.Rdata") if (!file.exists(PAfile)) { load(ABfile) PA <- mas5calls(AB) save(PA, file=PAfile) } else { load(PAfile) } ################################################### ### chunk number 7: normalize ################################################### EXPfile <- paste(GEO, sep="/", "EXP.Rdata") if (!file.exists(EXPfile)) { tissueEXP <- gcrma(AB, fast=FALSE) save(tissueEXP, file=EXPfile) } else { load(EXPfile) } ################################################### ### chunk number 8: filter based on the PA calls ################################################### keep <- rowSums(exprs(PA) == "P") >= 3 EXP <- tissueEXP[keep,] ################################################### ### chunk number 9: load the required annotation package ################################################### annotation(EXP) annpackage <- paste(annotation(EXP), sep=".", "db") library(annpackage, character.only=TRUE) ################################################### ### chunk number 10: add sample annotation ################################################### if (is.null(EXP$tissue)) { annot <- getGEOfile(GEO, amount="quick") annotSample <- function(x) { x <- sub(".CEL.gz", "", x, fixed=TRUE) f <- getGEOfile(x, amount="quick") l <- readLines(f) s <- grep("Sample_source_name", l, value=TRUE)[1] strsplit(s, " = ")[[1]][2] } tissues <- sapply(sampleNames(EXP), annotSample) tissueEXP$tissue <- sub("Mouse ", "", tissues, fixed=TRUE) EXP$tissue <- tissueEXP$tissue save(tissueEXP, file=EXPfile) } ################################################### ### chunk number 11: ISA ################################################### NEXP <- ISANormalize(EXP) modulesfile <- paste(GEO, sep="/", "modules.Rdata") if (!file.exists(modulesfile)) { set.seed(123) modules <- ISA(EXP, flist=NA, thr.gene=3, thr.cond=c(1,2), no.seeds=1000) save(modules, file=modulesfile) } else { load(modulesfile) } ################################################### ### chunk number 12: check modules ################################################### length(modules) getNoFeatures(modules) getNoSamples(modules) ################################################### ### chunk number 13: pick a module ################################################### ## We choose a module that has opposite sample scores and ## at least two tissues. mymod <- which(sapply(getSampleScores(modules), function(x) { length(unique(sign(x)))==2 }) & sapply(getSamples(modules), function(x) { length(unique(EXP$tissue[x]))>=2 }) ) if (length(mymod)==0) mymod <- length(modules) else mymod <- mymod[1] ################################################### ### chunk number 14: cond plot parameters ################################################### sep <- c(which(!duplicated(NEXP$tissue))[-1]-1, ncol(NEXP)) names(sep) <- NEXP$tissue[!duplicated(NEXP$tissue)] colbar <- rainbow(length(sep)) col <- colbar[as.factor(NEXP$tissue)] ################################################### ### chunk number 15: condplot eval=FALSE ################################################### ## condPlot(modules, mymod, NEXP, sep=sep, col=col) ################################################### ### chunk number 16: condplot2 ################################################### condPlot(modules, mymod, NEXP, sep=sep, col=col) ################################################### ### chunk number 17: function to decide tissue selectivity ################################################### tspec <- function(mod) { t <- EXP$tissue[getSamples(modules, mod)[[1]]] s <- getSampleScores(modules, mod)[[1]] t[ s<0 ] <- paste(t[s<0], sep=" ", "(d)") unique(t) } spec <- sapply(seq_len(length(modules)), tspec) specc <- sapply(spec, paste, collapse=", ") specc ################################################### ### chunk number 18: condplot3 eval=FALSE ################################################### ## to.plot <- which(sapply(spec, length)==1)[1:4] ## layout(rbind(1:2,3:4)) ## for (p in to.plot) { ## condPlot(modules, p, NEXP, sep=sep, col=col) ## } ################################################### ### chunk number 19: condplot4 ################################################### to.plot <- which(sapply(spec, length)==1)[1:4] layout(rbind(1:2,3:4)) for (p in to.plot) { condPlot(modules, p, NEXP, sep=sep, col=col) } ################################################### ### chunk number 20: enrichment calculation ################################################### enrichFile <- paste(GEO, sep="/", "enrichment.Rdata") if (!file.exists(enrichFile)) { GO <- ISAGO(modules) KEGG <- ISAKEGG(modules) save(GO, KEGG, file=enrichFile) } else { load(enrichFile) } GO KEGG ################################################### ### chunk number 21: function to query GO terms ################################################### goterm <- function(x) { gt <- sapply(mget(na.omit(x), GOTERM), Term) res <- character(length(x)) res[!is.na(x)] <- gt res[is.na(x)] <- NA res } keggpath <- function(x) { kp <- unlist(mget(na.omit(x), KEGGPATHID2NAME)) res <- character(length(x)) res[!is.na(x)] <- kp res[is.na(x)] <- NA res } ################################################### ### chunk number 22: first function ################################################### first <- function(x) if (length(x) >= 1) x[[1]] else NA ################################################### ### chunk number 23: create GO ################################################### topgo <- sapply(seq_len(length(modules)), function(x) { p <- c(BP=summary(GO$BP, p=2)[[x]][1,]$Pvalue, CC=summary(GO$CC, p=2)[[x]][1,]$Pvalue, MF=summary(GO$MF, p=2)[[x]][1,]$Pvalue) mc <- names(which.min(p)) paste(goterm(sigCategories(GO[[mc]], p=2)[[x]][1]), sep=", p=", format(min(p), digits=3)) }) ################################################### ### chunk number 24: create KEGG ################################################### keggterms <- keggpath(sapply(sigCategories(KEGG, p=2), first)) keggp <- sapply(summary(KEGG, p=2), function(x) x$Pvalue[1]) ################################################### ### chunk number 25: put together table ################################################### eTable <- data.frame(Tissue=specc, GO=topgo, KEGG=paste(keggterms, sep=", p=", format(keggp, digits=3))) ################################################### ### chunk number 26: print the table eval=FALSE ################################################### ## caption <- paste("The most significantly enriched GO categories", ## "and KEGG pathways for the (first 20)", ## "transcription modules.", ## "Oppositely regulated tissues are marged with `d'.") ## if (nrow(eTable)>20) { eTable2 <- eTable[1:20,] } ## print(xtable(eTable2, caption=caption, ## align=c("r@{\\hspace{1em}}", ## "p{4.5cm}@{\\hspace{1em}}", ## "p{7.5cm}@{\\hspace{1em}}", "p{7.5cm}")), ## floating.environment="sidewaystable", ## tabular.environment="tabular") ################################################### ### chunk number 27: really print it ################################################### caption <- paste("The most significantly enriched GO categories", "and KEGG pathways for the (first 20)", "transcription modules.", "Oppositely regulated tissues are marged with `d'.") if (nrow(eTable)>20) { eTable2 <- eTable[1:20,] } print(xtable(eTable2, caption=caption, align=c("r@{\\hspace{1em}}", "p{4.5cm}@{\\hspace{1em}}", "p{7.5cm}@{\\hspace{1em}}", "p{7.5cm}")), floating.environment="sidewaystable", tabular.environment="tabular") ################################################### ### chunk number 28: sessioninfo ################################################### toLatex(sessionInfo())