###################################################################### # R commands Genomic Data Science # Lab 15:GSEA University of South Carolina ###################################################################### # This file contains the R commands for the lab. # # Lines beginning with the symbol '#' are comments in R. All other # lines contain code. # # In R for Windows, you may wish to open this file from the menu bar # (File:Display file); you can then copy commands into the command # window. (Use the mouse to highlight one or more lines; then # right-click and select "Paste to console".) ###################################################################### ################################################### ### chunk: load0 ################################################### ################################################### ### chunk: ALLdata ################################################### library("limma") library("annaffy") library("GOstats") ################# library(ALL) library("hgu95av2.db") 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: nsfilter ################################################### group<-as.numeric(ALLfilt_af4bcr$mol.biol=="BCR/ABL") design = cbind(mean = 1, diff =group) ###!!!!!!!!!! ################################################### ### chunk: limma ################################################### fit = lmFit(exprs(ALLfilt_af4bcr), design) fit = eBayes(fit) t2<-topTable(fit, coef=2, sort.by="p", p.value=0.05, number=nrow(exprs(ALLfilt_af4bcr))) #### Annotation ALLset1Syms = unlist(mget(featureNames(ALLfilt_af4bcr), env = hgu95av2SYMBOL)) ## Alternative Approach # ALLset1Syms<-as.character(hgu95av2SYMBOL[featureNames(ALLfilt_af4bcr)]) ##### fit$gene<-ALLset1Syms t3<-topTable(fit, coef=2, sort.by="p", p.value=0.05, genelist=fit$gene, number=nrow(exprs(ALLfilt_af4bcr))) ################################################### ### chunk: Get annotation ################################################### library("hgu95av2.db") #library("hgu133plus2.db") ALLset1Syms = unlist(mget(featureNames(ALLfilt_af4bcr), env = hgu95av2SYMBOL)) ################################################### ### chunk: ALLsub ################################################### selid<-rownames(t2) ## probe_id of DE genes sel<-match(selid, featureNames(ALLfilt_af4bcr)) ALLsub = ALLfilt_af4bcr[sel,] ################################################### ### chunk: heatmap ################################################### library(gplots) mycolor<-colorpanel(n=100, low="black", mid="white", high="darkblue") color.map <- function(mol.biol) { if (mol.biol=="NEG") { ans<-"blue" } else if (mol.biol=="ALL1/AF4") { ans="#008000" } else { ans="#7879FF" } return(ans) } patientcolors <- unlist(lapply(ALLsub$mol.bio, color.map)) mat<-exprs(ALLsub) myheat1<-heatmap.2(mat,cexRow=0.02,col=mycolor, trace="none", ColSideColors=patientcolors) #pdf("Heatmap1.pdf") myheat2<-heatmap.2(mat,cexRow=0.02, col=heat.colors(100), trace="none", ColSideColors=patientcolors) #dev.off() ################################################### ### chunk: EGdup1 ################################################### apropos("hgu95av2") EG = as.character(hgu95av2ENTREZID[featureNames(ALL)]) EGsub = as.character(hgu95av2ENTREZID[featureNames(ALLsub)]) ################################################### ### chunk: tabletable ################################################### tab<-table(EG) table(table(EG)) ################################################### ### chunk: figprofile ################################################### syms = as.character(hgu95av2SYMBOL[featureNames(ALLsub)]) #syms = unlist(mget(featureNames(ALLsub), env= hgu95av2SYMBOL)) 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: html simple ################################################### library("annotate") probeids<-featureNames(ALLsub) geneList<-list(probeids, EGsub) repository<-list("affy", "en") otherNames<-list(syms) head<-c("Probe ID", "EntrezID", "Symbol") filename<-"out1.html" htmlpage(genelist=geneList, filename=filename, title="ALL Interesting Genes", othernames=otherNames, table.head=head, repository=repository) ################################################### ### 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", widget = F) localURL = file.path("file:/", getwd(), "ALLsub.html") browseURL(localURL) ################################################### ### 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.db", ontology="BP", pvalueCutoff=0.001, conditional=FALSE, testDirection="over") hgOver = hyperGTest(params) ################################################### ### chunk: figmfhyper ################################################### hist(pvalues(hgOver), breaks=50, col="mistyrose") ################################################### ### chunk: sigCats ################################################### summ = GOstats::summary(hgOver, p=0.001) head(summ) df = summary(hgOver) names(df) df = summary(hgOver, pvalue=0.05, categorySize=350) nrow(df) ################################################### ### chunk:hgOver html ################################################### htmlReport(hgOver, file="ALL_hgo.html") #################################################### # Conditional HyperGeometric Test #################################################### hgCutoff = 0.001 params = new("GOHyperGParams", geneIds=EGsub, universeGeneIds=entrezUniverse, annotation="hgu95av2.db", ontology="BP", pvalueCutoff=hgCutoff, conditional=FALSE, testDirection="over") paramsCond = params conditional(paramsCond) = TRUE ################################################### ### chunk: condhgRun ################################################### hgCond = hyperGTest(paramsCond) ################################################### ### chunk: condhgResult ################################################### summary(hgCond) ################################################### ### chunk: KEGGHyperGparam ################################################### library("GSEABase") library("Category") library("org.Hs.eg.db") frame = toTable(org.Hs.egPATH) keggframeData = data.frame(frame$path_id, frame$gene_id) head(keggframeData) keggFrame = KEGGFrame(keggframeData, organism = "Homo sapiens") gsc <- GeneSetCollection(keggFrame, setType = KEGGCollection()) kparams <- GSEAKEGGHyperGParams(name = "My Custom GSEA based annot Params", geneSetCollection = gsc, geneIds = EGsub, universeGeneIds = entrezUniverse, pvalueCutoff = 0.1, testDirection = "over") kOver <- hyperGTest(kparams) summary(kOver) ################################################### ### chunk: ans1 ################################################### library("ALL") data("ALL") ################################################### ### chunk: bcrabl ################################################### bcell = grep("^B", as.character(ALL$BT)) moltyp = which(as.character(ALL$mol.biol) %in% c("NEG", "BCR/ABL")) ALL_bcrneg = ALL[, intersect(bcell, moltyp)] ALL_bcrneg$mol.biol = factor(ALL_bcrneg$mol.biol) ################################################### ### chunk: nsFilter ################################################### library("genefilter") ALLfilt_bcrneg = nsFilter(ALL_bcrneg, var.cutoff=0.5)$eset ################################################### ### chunk: ans1 ################################################### table(ALLfilt_bcrneg$mol.biol) ################################################### ### chunk: Enrichment Analysis ################################################### library("GSEABase") gsc = GeneSetCollection(ALLfilt_bcrneg, setType=KEGGCollection()) Am = incidence(gsc) dim(Am) Am[1:5, 1:5] ################################################### ### chunk: EsetfromKEGG ################################################### nsF = ALLfilt_bcrneg[colnames(Am),] ################################################### ### chunk: compttests ################################################### rtt = rowttests(nsF, "mol.biol") rttStat = rtt$statistic ################################################### ### chunk: reducetoInt ################################################### selectedRows = (rowSums(Am)>10) Am2 = Am[selectedRows, ] ################################################### ### chunk: pctests ################################################### tA = as.vector(Am2 %*% rttStat) tAadj = tA/sqrt(rowSums(Am2)) names(tA) = names(tAadj) = rownames(Am2) ################################################### ### chunk: ttperm ################################################### library(Category) set.seed(123) NPERM = 1000 pvals = gseattperm(nsF, nsF$mol.biol, Am2, NPERM) pvalCut = 0.025 lowC = names(which(pvals[, 1]<=pvalCut)) highC = names(which(pvals[, 2]<=pvalCut)) lowC highC xx<-as.list(KEGGPATHID2NAME) wlowC<-match(lowC, names(xx)) lowCName<-xx[wlowC] lowCName whighC<-match(highC, names(xx)) highCName<-xx[whighC] highCName ############################################################################################## ##### More detail about parsing out overlap genesets please check ###### Bioconductor Case Studies: Chapter 13 ################################################### ### chunk: Overlap genes of two terms ################################################### Ams = Am2[union(lowC, highC),] Amx = Ams %*% t(Ams) minS = outer(diag(Amx), diag(Amx), pmin) overlapIndex = Amx/minS ################################################### ### chunk: overlapSetsPlot ################################################### library("lattice") myFun = function(x) { dd.row = as.dendrogram(hclust(dist(x))) row.ord = order.dendrogram(dd.row) dd.col = as.dendrogram(hclust(dist(t(x)))) col.ord = order.dendrogram(dd.col) colnames(x) = sapply(getPathNames(colnames(x)), "[", 1L) levelplot(x[row.ord,col.ord], scales = list(x = list(rot = 90), cex=0.5), xlab="", ylab="", main="overlapIndex", col.regions=colorRampPalette(c("white", "darkblue"))(256)) } pdf("Overlap.pdf") print(myFun(overlapIndex)) dev.off() ################################################### ### chunk: lmfits ################################################### P04512 = Ams["04512",] P04510 = Ams["04510",] lm1 = lm(rttStat ~ P04512) summary(lm1)$coefficients lm2 = lm(rttStat ~ P04510) summary(lm2)$coefficients lm3 = lm(rttStat ~ P04510+P04512) summary(lm3)$coefficients ################################################### ### chunk: stopifnot ################################################### stopifnot( coefficients(summary(lm1))["P04512", "Pr(>|t|)"]<0.05, coefficients(summary(lm2))["P04510", "Pr(>|t|)"]<0.05, coefficients(summary(lm3))["P04512", "Pr(>|t|)"]>0.05, coefficients(summary(lm3))["P04510", "Pr(>|t|)"]<0.05) ################################################### ### chunk: threegroups ################################################### P04512.Only = ifelse(P04512 != 0 & P04510 == 0, 1, 0) P04510.Only = ifelse(P04512 == 0 & P04510 != 0, 1, 0) Both = ifelse(P04512 != 0 & P04510 != 0, 1, 0) lm4 = lm(rttStat ~ P04510.Only + P04512.Only + Both)