[BioC] Mean Plots and Heat Maps for Chromosome Bands

[McGee, Monnie]

Dear BioC Users,

I am following the code for Chapter 13 of "Bioconductor Case Studies" (dealing with GSEA.  The example of GSEA analysis in the text shows the use of KEGG pathways to determine gene sets.  Then the functions KEGGmnplot and KEGGheatmap are used to obtain a mean plot and a heat map for one of the gene sets (corresponding to one pathway) that has a very large (negative) t-statistic.

Section 13.2.4 suggests using gene sets defined by chromosome bands to do the same analysis.  I have found an "interesting" chromosome band (cytogenetic) location, and I would like to get a heatmap of all the genes from that location, as well as a mean plot.  However, there is no nice function "CHRLOCmnplot"  or "CHRLOCheatmap".  What is the most elegant way to get mean plots and heat maps of genes from a gene set associated with a particular cytogenetic location?   

The code that I used and the session information follow.

Thank you,
Monnie

> library(GSEABase)
> library(ALL)
> data(ALL)
> 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)
> library(genefilter)
> ALLfilt = nsFilter(ALL_bcrneg,var.cutoff=0.5)$eset
> fnames = featureNames(ALLfilt)
> if(is(hgu95av2MAP,"environment")){
+ chrLocs = mget(fnames,hgu95av2MAP)
+ mapping = names(chrLocs[sapply(chrLocs,function(x) !all(is.na(x)))])
+ }else{
+ mapping = toTable(hgu95av2MAP[fnames])$probe_id
+ }
> psWithMAP = unique(mapping)
> nsF2 = ALLfilt[psWithMAP,]
> chrMat1 = chrMat[rowSums(chrMat) >= 5,]
> rtt = rowttests(nsF2,"mol.biol")
> rttStat = rtt$statistic
> tA = as.vector(chrMat1 %*% rttStat)
> tAadj = tA/sqrt(rowSums(chrMat1))
> names(tA) = names(tAadj) = rownames(chrMat1)
> smBand = tAadj[tAadj < (-5)]
> smBand
  19q13.3 
-5.741334 
> EGtable = toTable(hgu95av2ENTREZID[featureNames(nsF2)])
> entrezUniv = unique(EGtable$gene_id)
> chrMat = MAPAmat("hgu95av2", univ = entrezUniv)
> EGlist = mget(featureNames(nsF2),hgu95av2ENTREZID)
> EGIDs = sapply(EGlist , "[",1)
> idx = match(EGIDs,colnames(chrMat1))
> chrInc = chrMat1[, idx]
## Now - can I get a heatmap and mean plot of the genes in the gene set "smBand" in an elegant way with the information I have?

> sessionInfo()
R version 2.7.2 (2008-08-25) 
i386-apple-darwin8.11.1 

locale:
C

attached base packages:
[1] splines   tools     stats     graphics  grDevices utils     datasets  methods  
[9] base     

other attached packages:
 [1] humanCHRLOC_2.1.4   Category_2.6.0      graph_1.18.1        KEGG.db_2.2.0      
 [5] hgu95av2.db_2.2.0   genefilter_1.20.0   survival_2.34-1     ALL_1.4.4          
 [9] GSEABase_1.2.2      annotate_1.18.0     xtable_1.5-3        AnnotationDbi_1.2.2
[13] RSQLite_0.7-0       DBI_0.2-4           Biobase_2.0.1      

loaded via a namespace (and not attached):
[1] RBGL_1.16.0     XML_1.96-0      cluster_1.11.11