################################################### ### chunk: load0 ################################################### library("BiocCaseStudies") ################################################### ### chunk: ALLdata ################################################### library("ALL") data("ALL") types = c("ALL1/AF4", "BCR/ABL") bcell = grep("^B", as.character(ALL$BT)) ALL_af4bcr = ALL[, intersect(bcell, which(ALL$mol.biol %in% types))] ALL_af4bcr$mol.biol = factor(ALL_af4bcr$mol.biol) ################################################### ### chunk: groupsize ################################################### table(ALL_af4bcr$mol.biol) ################################################### ### chunk: nsfilter ################################################### qrange <- function(x) diff(quantile(x, c(0.1, 0.9))) library("genefilter") filt_af4bcr = nsFilter(ALL_af4bcr, require.entrez=TRUE, require.GOBP=TRUE, var.func=qrange, var.cutoff=0.5) ALLfilt_af4bcr = filt_af4bcr$eset ################################################### ### chunk: loadlibs ################################################### library("Biobase") library("annotate") library("hgu95av2.db") ################################################### ### chunk: rt ################################################### rt = rowttests(ALLfilt_af4bcr, "mol.biol") names(rt) ################################################### ### chunk: figrtsol1 ################################################### hist(rt$statistic, breaks=100, col="skyblue", main="", xlab="t-statistic") ################################################### ### chunk: figrtsol2 ################################################### hist(rt$p.value, breaks=100, col="mistyrose", main="", xlab="p-value") ################################################### ### chunk: ALLsub ################################################### sel = order(rt$p.value)[1:400] ALLsub = ALLfilt_af4bcr[sel,] ################################################### ### chunk: EGdup1 ################################################### EG = as.character(hgu95av2ENTREZID[featureNames(ALL)]) EGsub = as.character(hgu95av2ENTREZID[featureNames(ALLsub)]) ################################################### ### chunk: tabletable ################################################### table(table(EG)) table(table(EGsub)) ################################################### ### chunk: figprofile ################################################### syms = as.character(hgu95av2SYMBOL[featureNames(ALLsub)]) whFeat = names(which(syms =="CD44")) ordSamp = order(ALLsub$mol.biol) CD44 = ALLsub[whFeat, ordSamp] plot(as.vector(exprs(CD44)), main=whFeat, col=c("sienna", "tomato")[CD44$mol.biol], pch=c(15, 16)[CD44$mol.biol], ylab="expression") ################################################### ### chunk: chrtab ################################################### z = toTable(hgu95av2CHR[featureNames(ALLsub)]) chrtab = table(z$chromosome) chrtab ################################################### ### chunk: figchrtab ################################################### chridx = sub("X", "23", names(chrtab)) chridx = sub("Y", "24", chridx) barplot(chrtab[order(as.integer(chridx))]) ################################################### ### chunk: annaffy1 ################################################### library("annaffy") anncols = aaf.handler(chip="hgu95av2.db")[c(1:3, 8:9, 11:13)] anntable = aafTableAnn(featureNames(ALLsub), "hgu95av2.db", anncols) saveHTML(anntable, "ALLsub.html", title="The Features in ALLsub") localURL = file.path("file:/", getwd(), "ALLsub.html") ################################################### ### chunk: annaffy3 ################################################### browseURL(localURL) ################################################### ### chunk: figduplo0 ################################################### probeSetsPerGene = split(names(EG), EG) j = probeSetsPerGene$"7013" j ################################################### ### chunk: figduplo1 ################################################### plot(t(exprs(ALL_af4bcr)[j[c(1,7)], ]), asp=1, pch=16, col=ifelse(ALL_af4bcr$mol.biol=="ALL1/AF4", "black", "grey")) ################################################### ### chunk: figduplo2 ################################################### library("lattice") mat = exprs(ALL_af4bcr)[j,] mat = mat - rowMedians(mat) ro = order.dendrogram(as.dendrogram(hclust(dist(mat)))) co = order.dendrogram(as.dendrogram(hclust(dist(t(mat))))) at = seq(-1, 1, length=21) * max(abs(mat)) lp = levelplot(t(mat[ro, co]), aspect = "fill", at = at, scales = list(x = list(rot = 90)), colorkey = list(space = "left")) print(lp) ################################################### ### chunk: chr1 ################################################### ps_chr = toTable(hgu95av2CHR) ps_eg = toTable(hgu95av2ENTREZID) chr = merge(ps_chr, ps_eg) chr = unique(chr[, colnames(chr)!="probe_id"]) head(chr) ################################################### ### chunk: chr2 ################################################### table(table(chr$gene_id)) ################################################### ### chunk: chr4 ################################################### chr = chr[!duplicated(chr$gene_id), ] ################################################### ### chunk: EGsub ################################################### isdiff = chr$gene_id %in% EGsub tab = table(isdiff, chr$chromosome) tab fisher.test(tab, simulate.p.value=TRUE) chisq.test(tab) ################################################### ### chunk: chrloc1 ################################################### chrloc = toTable(hgu95av2CHRLOC[featureNames(ALLsub)]) head(chrloc) ################################################### ### chunk: chrloc2 ################################################### table(table(chrloc$probe_id)) ################################################### ### chunk: chrloc3 ################################################### strds = with(chrloc, unique(cbind(probe_id, sign(start_location)))) table(strds[,2]) ################################################### ### chunk: getMFchildren ################################################### library("GO.db") as.list(GOMFCHILDREN["GO:0008094"]) ################################################### ### chunk: getMFoffspring ################################################### as.list(GOMFOFFSPRING["GO:0008094"]) ################################################### ### chunk: loadGOstats ################################################### library("GOstats") ################################################### ### chunk: GOHyperGparam ################################################### affyUniverse = featureNames(ALLfilt_af4bcr) uniId = hgu95av2ENTREZID[affyUniverse] entrezUniverse = unique(as.character(uniId)) params = new("GOHyperGParams", geneIds=EGsub, universeGeneIds=entrezUniverse, annotation="hgu95av2", ontology="BP", pvalueCutoff=0.001, conditional=FALSE, testDirection="over") ################################################### ### chunk: mfhyper ################################################### mfhyper = hyperGTest(params) ################################################### ### chunk: figmfhyper ################################################### hist(pvalues(mfhyper), breaks=50, col="mistyrose") ################################################### ### chunk: sigCats ################################################### sum = summary(mfhyper, p=0.001) head(sum) ################################################### ### chunk: GOTERM ################################################### GOTERM[["GO:0032945"]] ################################################### ### chunk: bmCon eval=FALSE ################################################### ## library("biomaRt") ## head(listMarts()) ################################################### ### chunk: choseMart eval=FALSE ################################################### ## mart = useMart("ensembl") ################################################### ### chunk: bmDataset eval=FALSE ################################################### ## head(listDatasets(mart)) ################################################### ### chunk: bmChoseset eval=FALSE ################################################### ## ensembl = useDataset("hsapiens_gene_ensembl", ## mart=mart) ################################################### ### chunk: bmutrDo eval=FALSE ################################################### ## utr = getSequence(id=EGsub, seqType="3utr", ## mart=ensembl, type="entrezgene") ## utr[1,] ################################################### ### chunk: bmfilt eval=FALSE ################################################### ## head(listFilters(ensembl)) ################################################### ### chunk: bmatt eval=FALSE ################################################### ## head(listAttributes(ensembl)) ################################################### ### chunk: bmDom eval=FALSE ################################################### ## domains = getBM(attributes=c("entrezgene", "pfam", ## "prosite", "interpro"), filters="entrezgene", ## value=EGsub, mart=ensembl) ## interpro = split(domains$interpro, domains$entrezgene) ## interpro[1] ################################################### ### chunk: loadDBanno ################################################### library("hgu133a.db") dbc = hgu133a_dbconn() ################################################### ### chunk: accessEnvs ################################################### get("201473_at", hgu133aSYMBOL) mget(c("201473_at","201476_s_at"), hgu133aSYMBOL) hgu133aSYMBOL$"201473_at" hgu133aSYMBOL[["201473_at"]] ################################################### ### chunk: GOcharacteristics1 ################################################### goCats = unlist(eapply(GOTERM, Ontology)) gCnums = table(goCats)[c("BP","CC", "MF")] ################################################### ### chunk: xtable eval=FALSE ################################################### ## library("xtable") ## xtable(as.matrix(gCnums), display=c("d", "d"), ## caption="Number of GO terms per ontology.", ## label="ta:GOprops") ################################################### ### chunk: xtabledo ################################################### library("xtable") cat(print(xtable(as.matrix(gCnums), display=c("d", "d"), caption="Number of \\indexTerm[Gene Ontology (GO)]{GO} terms per ontology.", label="ta:GOprops"), file="table1.tex")) ################################################### ### chunk: GOcharacteristics2 ################################################### query = "select ontology from go_term" goCats = dbGetQuery(GO_dbconn(), query) gCnums2 = table(goCats)[c("BP","CC", "MF")] identical(gCnums, gCnums2) ################################################### ### chunk: GOterms2 ################################################### query = paste("select term from go_term where term", "like '%chromosome%'") chrTerms = dbGetQuery(GO_dbconn(), query) nrow(chrTerms) head(chrTerms) ################################################### ### chunk: tfb ################################################### query = paste("select go_id from go_term where", "term = 'transcription factor binding'") tfb = dbGetQuery(GO_dbconn(), query) tfbps = hgu133aGO2ALLPROBES[[tfb$go_id]] table(names(tfbps)) ################################################### ### chunk: tfsol ################################################### query = paste("select term from go_term where term", "like '%transcription factor%'") tf = dbGetQuery(GO_dbconn(), query) nrow(tf) head(tf) ################################################### ### chunk: SYM2EG ################################################### queryAlias = function(x) { it = paste("('", paste(x, collapse="', '"), "'", sep="") paste("select _id, alias_symbol from alias", "where alias_symbol in", it, ");") } queryGeneinfo = function(x) { it = paste("('", paste(x, collapse="', '"), "'", sep="") paste("select _id, symbol from gene_info where", "symbol in", it, ");") } queryGenes = function(x) { it = paste("('", paste(x, collapse="', '"), "'", sep="") paste("select * from genes where _id in", it, ");") } findEGs = function(dbcon, symbols) { rs = dbSendQuery(dbcon, queryGeneinfo(symbols)) a1 = fetch(rs, n=-1) stillLeft = setdiff(symbols, a1[,2]) if( length(stillLeft)>0 ) { rs = dbSendQuery(dbcon, queryAlias(stillLeft)) a2 = fetch(rs, n=-1) names(a2) = names(a1) a1 = rbind(a1, a2) } rs = dbSendQuery(dbcon, queryGenes(a1[,1])) merge(a1, fetch(rs, n=-1)) } ################################################### ### chunk: findEGs ################################################### findEGs(dbc, c("ALL1", "AF4", "BCR", "ABL")) ################################################### ### chunk: revMSym ################################################### s1 = revmap(hgu133aSYMBOL) s1$BCR ################################################### ### chunk: toTable ################################################### toTable(hgu133aGO["201473_at"])