inst/doc/categoryCompare_vignette.R

In categoryCompare: Meta-analysis of high-throughput experiments using feature annotations

## ----loadLibs, message=FALSE--------------------------------------------------
library("affy")
library("hgu95av2.db")
library("genefilter")
library("estrogen")
library("limma")

## ----loadMeta-----------------------------------------------------------------
datadir <- system.file("extdata", package = "estrogen")
pd <- read.AnnotatedDataFrame(file.path(datadir,"estrogen.txt"), 
  	header = TRUE, sep = "", row.names = 1)
pData(pd)

## ----loadAffy-----------------------------------------------------------------
currDir <- getwd()
setwd(datadir)
a <- ReadAffy(filenames=rownames(pData(pd)), phenoData = pd, verbose = TRUE)
setwd(currDir)

## ----normalizeAffy, message=FALSE---------------------------------------------
eData <- rma(a)

## ----edata10------------------------------------------------------------------
e10 <- eData[, eData$time.h == 10]
e10 <- nsFilter(e10, remove.dupEntrez=TRUE, var.filter=FALSE, 
        feature.exclude="^AFFX")$eset

e10$estrogen <- factor(e10$estrogen)
d10 <- model.matrix(~0 + e10$estrogen)
colnames(d10) <- unique(e10$estrogen)
fit10 <- lmFit(e10, d10)
c10 <- makeContrasts(present - absent, levels=d10)
fit10_2 <- contrasts.fit(fit10, c10)
eB10 <- eBayes(fit10_2)
table10 <- topTable(eB10, number=nrow(e10), p.value=1, adjust.method="BH")
table10$Entrez <- unlist(mget(rownames(table10), hgu95av2ENTREZID, ifnotfound=NA))

## ----edata48------------------------------------------------------------------
e48 <- eData[, eData$time.h == 48]
e48 <- nsFilter(e48, remove.dupEntrez=TRUE, var.filter=FALSE, 
        feature.exclude="^AFFX" )$eset

e48$estrogen <- factor(e48$estrogen)
d48 <- model.matrix(~0 + e48$estrogen)
colnames(d48) <- unique(e48$estrogen)
fit48 <- lmFit(e48, d48)
c48 <- makeContrasts(present - absent, levels=d48)
fit48_2 <- contrasts.fit(fit48, c48)
eB48 <- eBayes(fit48_2)
table48 <- topTable(eB48, number=nrow(e48), p.value=1, adjust.method="BH")
table48$Entrez <- unlist(mget(rownames(table48), hgu95av2ENTREZID, ifnotfound=NA))

## ----gUniverse----------------------------------------------------------------
gUniverse <- unique(union(table10$Entrez, table48$Entrez))

## ----createGeneList, message=FALSE--------------------------------------------
library("categoryCompare")
library("GO.db")
library("KEGG.db")

g10 <- unique(table10$Entrez[table10$adj.P.Val < 0.05])
g48 <- unique(table48$Entrez[table48$adj.P.Val < 0.05])

## ----viewLists----------------------------------------------------------------
list10 <- list(genes=g10, universe=gUniverse, annotation='org.Hs.eg.db')
list48 <- list(genes=g48, universe=gUniverse, annotation='org.Hs.eg.db')

geneLists <- list(T10=list10, T48=list48)
geneLists <- new("ccGeneList", geneLists, ccType=c("BP","KEGG"))
fdr(geneLists) <- 0 # this speeds up the calculations for demonstration
geneLists

## ----runEnrichment------------------------------------------------------------
enrichLists <- ccEnrich(geneLists)
enrichLists

## ----modPVal------------------------------------------------------------------
pvalueCutoff(enrichLists$BP) <- 0.001
enrichLists

## ----ccOptions----------------------------------------------------------------
ccOpts <- new("ccOptions", listNames=names(geneLists), outType='none')
ccOpts

## ----ccResults----------------------------------------------------------------
ccResults <- ccCompare(enrichLists,ccOpts)
ccResults

## ----setupRunningCytoscape, echo=FALSE, error=FALSE, results='asis'-----------
push_graph <- try(
  {
    # borrowed from createNetworkFromDataFrames
    nodes <- data.frame(id=c("node 0","node 1","node 2","node 3"),
           stringsAsFactors=FALSE)
    edges <- data.frame(source=c("node 0","node 0","node 0","node 2"),
           target=c("node 1","node 2","node 3","node 3"),
           stringsAsFactors=FALSE)

    RCy3::createNetworkFromDataFrames(nodes,edges)
  }
)

if (inherits(push_graph, "try-error")) {
  runCy <- FALSE
  print("No connection to Cytoscape available, subsequent visualizations were not run")
} else {
  runCy <- TRUE
}

## ----cwBP, eval=runCy---------------------------------------------------------
#  ccBP <- breakEdges(ccResults$BP, 0.2)
#  
#  cw.BP <- ccOutCyt(ccBP, ccOpts)

## ----breakHighBP, eval=runCy--------------------------------------------------
#  breakEdges(cw.BP,0.4)
#  breakEdges(cw.BP,0.6)
#  breakEdges(cw.BP,0.8)

## ----deleteBP, include=FALSE, eval=runCy--------------------------------------
#  RCy3::deleteNetwork(cw.BP)

## ----GOHierarchical-----------------------------------------------------------
graphType(enrichLists$BP) <- "hierarchical"
ccResultsBPHier <- ccCompare(enrichLists$BP, ccOpts)
ccResultsBPHier

## ----hier2Cytoscape, eval=runCy-----------------------------------------------
#  cw.BP2 <- ccOutCyt(ccResultsBPHier, ccOpts, "BPHier")

## ----deleteHier, include=FALSE, eval=runCy------------------------------------
#  RCy3::deleteNetwork(cw.BP2)

## ----cytoscapeKEGG, eval=runCy------------------------------------------------
#  cw.KEGG <- ccOutCyt(ccResults$KEGG,ccOpts)

## ----deleteKEGG, include=FALSE, eval=runCy------------------------------------
#  RCy3::deleteNetwork(cw.KEGG)

## ----sessionInfo--------------------------------------------------------------
sessionInfo()