inst/doc/ChromBand.R

In Category: Category Analysis

### R code from vignette source 'ChromBand.Rnw'

###################################################
### code chunk number 1: setup1
###################################################
library("Category")
library("ALL")
library("hgu95av2.db")
library("annotate")
library("genefilter")
##library("SNPchip")
library("karyoploteR")
library("geneplotter")
library("limma")
library("lattice")
library("graph")


###################################################
### code chunk number 2: chr12ideogram
###################################################
build <- "hg19" ## or hg18
kp <- plotKaryotype(genome=build, chromosomes="chr12", main="Human chromosome 12")
## cyt2 <- getCytoband(build=build)
## cyt2$gieStain <- "foo"
## c12p12_idx <- intersect(grep("^q21", cyt2$name),
##                         which(cyt2$chrom == "12"))
## cyt2[c12p12_idx, "gieStain"] <- rep(c("gpos50", "gpos75"),
##                                     length=length(c12p12_idx))
##
##
## plotCytoband2(chromosome="12", build=build, cytoband=cyt2, outer=FALSE,  ## From SNPchip
##     cex.axis=0.6, main="Human chromosome 12")
## plotCytoband2(chromosome="12", build=build, cytoband=cyt2, outer=FALSE,
##               cex.axis=0.6,
##               main="Human chromosome 12")


###################################################
### code chunk number 3: bcrAblOrNegSubset
###################################################
data(ALL, package="ALL")

subsetType <- "BCR/ABL"
Bcell <- grep("^B", as.character(ALL$BT))
bcrAblOrNegIdx <- which(as.character(ALL$mol.biol) %in% c("NEG", subsetType))

bcrAblOrNeg <- ALL[, intersect(Bcell, bcrAblOrNegIdx)]
bcrAblOrNeg$mol.biol <- factor(bcrAblOrNeg$mol.biol)


###################################################
### code chunk number 4: annMaps
###################################################
annType <- c("db", "env")
entrezMap <- getAnnMap("ENTREZID", annotation(bcrAblOrNeg),
                       type=annType, load=TRUE)
symbolMap <- getAnnMap("SYMBOL", annotation(bcrAblOrNeg),
                       type=annType, load=TRUE)
bandMap <- getAnnMap("MAP", annotation(bcrAblOrNeg),
                     type=annType, load=TRUE)


###################################################
### code chunk number 5: nsFiltering
###################################################
filterAns <- nsFilter(bcrAblOrNeg,
                      require.entrez = TRUE, 
                      remove.dupEntrez = TRUE, 
                      var.func = IQR, var.cutoff = 0.5)
nsFiltered <- filterAns$eset


###################################################
### code chunk number 6: ChromBand.Rnw:208-213
###################################################
hasSYM <- sapply(mget(featureNames(nsFiltered), symbolMap, ifnotfound=NA),
                 function(x) length(x) > 0 && !is.na(x[1]))
hasMAP <- sapply(mget(featureNames(nsFiltered), bandMap, ifnotfound=NA),
                 function(x) length(x) > 0 && !is.na(x[1]))
nsFiltered <- nsFiltered[hasSYM & hasMAP, ]


###################################################
### code chunk number 7: defineGeneUniverse
###################################################
affyUniverse <- featureNames(nsFiltered)
entrezUniverse <- unlist(mget(affyUniverse, entrezMap))
names(affyUniverse) <- entrezUniverse
if (any(duplicated(entrezUniverse)))
    stop("error in gene universe: can't have duplicate Entrez Gene Ids")


###################################################
### code chunk number 8: parametric1
###################################################
design <- model.matrix(~ 0 + nsFiltered$mol.biol)
colnames(design) <- c("BCR/ABL", "NEG")
contr <- c(1, -1) ## NOTE: we thus have BCR/ABL w.r.t NEG
fm1 <- lmFit(nsFiltered, design)
fm2 <- contrasts.fit(fm1, contr)
fm3 <- eBayes(fm2)
ttestLimma <- topTable(fm3, number = nrow(fm3), adjust.method = "none")
ttestLimma <- ttestLimma[featureNames(nsFiltered), ]

tstats <- ttestLimma$t
names(tstats) <- entrezUniverse[rownames(ttestLimma)]
##


###################################################
### code chunk number 9: selectedSubset
###################################################
ttestCutoff <- 0.01
smPV  <- ttestLimma$P.Value < ttestCutoff
pvalFiltered <- nsFiltered[smPV, ]
selectedEntrezIds <- unlist(mget(featureNames(pvalFiltered), entrezMap))
##


###################################################
### code chunk number 10: ChromBand.Rnw:418-434
###################################################

chrSortOrder <- function(df) {
    chrs <- sub("([^pq]+).*$", "\\1", rownames(df))
    xyIdx <- chrs %in% c("X", "Y")
    xydf <- NULL
    if (any(xyIdx)) {
        chrs <- chrs[!xyIdx]
        xydf <- df[xyIdx, ]
        df <- df[!xyIdx, ]
    }
    ord <- order(as.integer(chrs), rownames(df))
    df <- df[ord, ]
    if (!is.null(xydf))
      df <- rbind(df, xydf)
    df
}


###################################################
### code chunk number 11: ChromBand.Rnw:439-454
###################################################

gseaTstatStripplot <- function(bands, g, ..., include.all = FALSE)
{
    chroms <- c(1:22, "X", "Y")
    chromArms <- c(paste(chroms, "p", sep=""), paste(chroms, "q", sep=""))
    egid <- lapply(nodeData(g, bands), "[[", "geneIds")
    if (include.all) {
        egid$All <- 
            unique(unlist(lapply(nodeData(g)[chromArms], "[[", "geneIds")))
    }
    tdf <- do.call(make.groups, lapply(egid, function(x) tstats[x]))
    stripplot(which ~ data, tdf, jitter = TRUE, ...)
}




###################################################
### code chunk number 12: ChromBand.Rnw:459-506
###################################################

esetBWPlot <- function(tmpSet, ..., layout=c(1, nrow(emat)))
{
    emat <- exprs(tmpSet)
    pd <- pData(tmpSet)
    probes <- rownames(emat)
    syms <- 
        sapply(mget(probes, hgu95av2SYMBOL, ifnotfound=NA),
               function(x) if (all(is.na(x))) "NA" else as.character(x)[1])
    selectedAffy <- 
        probes %in% affyUniverse[selectedEntrezIds]
    symsSelected <- syms[selectedAffy]
    symsWithStatus <- 
        paste(syms, 
              ifelse(selectedAffy, "*", ""), 
              sep = "")
    pdat <- 
        cbind(exprs=as.vector(emat),
              genes=factor(probes, levels = probes, labels = syms),
              pd[rep(seq_len(nrow(pd)), each=nrow(emat)), ])
    pdat <- transform(pdat, genes = reorder(genes, exprs))
    panels.to.shade <- levels(pdat$genes) %in% symsSelected
    bwplot(mol.biol ~ exprs | genes, data=pdat, 
           layout = layout,
           auto.key=TRUE,
           scales=list(x=list(log=2L)),
           xlab="Log2 Expression",
           panels.to.shade = panels.to.shade,
           panel = function(..., panels.to.shade) {
               if (panels.to.shade[packet.number()]) 
                   panel.fill(col = "lightgrey")
               panel.bwplot(...)
           },
           strip=FALSE,
           strip.left=TRUE, ...)
}

g1 <- makeChrBandGraph(annotation(nsFiltered), univ=entrezUniverse)
ct <- ChrBandTreeFromGraph(g1)

subsetByBand <- function(eset, ct, band) {
    egIDs <- unlist(nodeData(ct@toChildGraph, n=band, 
                             attr="geneIds"), use.names=FALSE)
    wantedProbes <- affyUniverse[as.character(egIDs)]
    eset[intersect(wantedProbes, featureNames(eset)), ]
}



###################################################
### code chunk number 13: basicParams
###################################################
params <- new("ChrMapHyperGParams",
              conditional=FALSE,
              testDirection="over",
              universeGeneIds=entrezUniverse,
              geneIds=selectedEntrezIds,
              annotation="hgu95av2",
              pvalueCutoff=0.05)

paramsCond <- params
paramsCond@conditional <- TRUE


###################################################
### code chunk number 14: basicTest
###################################################
hgans <- hyperGTest(params)
hgansCond <- hyperGTest(paramsCond)


###################################################
### code chunk number 15: result1
###################################################
sumUn <- summary(hgans, categorySize=1)
chrSortOrder(sumUn)

sumCond <- summary(hgansCond, categorySize=1)
chrSortOrder(sumCond)


###################################################
### code chunk number 16: ChromBand.Rnw:597-612
###################################################
gseaTstatStripplot(c("12q21.1", "12q21", "12q2", "12q"),
                   include.all = TRUE, 
                   g = g1,
                   xlab = "Per-gene t-statistics", 
                   panel = function(...) {
                       require(grid, quietly = TRUE)
                       grid.rect(y = unit(2, "native"), 
                                 height = unit(1, "native"),
                                 gp = 
                                 gpar(fill = "lightgrey", 
                                      col = "transparent"))
                       panel.grid(v = -1, h = 0)
                       panel.stripplot(...)
                       panel.average(..., fun = mean, lwd = 3)
                   })


###################################################
### code chunk number 17: gseaStripplot
###################################################
plot(trellis.last.object())


###################################################
### code chunk number 18: LMtestSetup
###################################################
params <- new("ChrMapLinearMParams",
              conditional = FALSE,
              testDirection = "up",
              universeGeneIds = entrezUniverse,
              geneStats = tstats,
              annotation = "hgu95av2",
              pvalueCutoff = 0.01, 
              minSize = 4L)
params@graph <- makeChrBandGraph(params@annotation, params@universeGeneIds)
params@gsc <- makeChrBandGSC(params@graph)
paramsCond <- params
paramsCond@conditional <- TRUE


###################################################
### code chunk number 19: LMtest
###################################################

lmans <- linearMTest(params)
lmansCond <- linearMTest(paramsCond)

chrSortOrder(summary(lmans))
chrSortOrder(summary(lmansCond))

##


###################################################
### code chunk number 20: ChromBand.Rnw:716-719
###################################################
tmpSet <- subsetByBand(nsFiltered, ct, "1p36.2")
esetBWPlot(tmpSet, ylab="1p36.2", layout = c(2, 8),
           par.strip.text = list(cex = 0.8))


###################################################
### code chunk number 21: dotplot1p362
###################################################
plot(trellis.last.object())