Entering edit mode
McGee, Monnie
▴
300
@mcgee-monnie-1108
Last seen 10.5 years ago
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) !allis.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