In rmflight/ccPaper: work for categoryCompare paper revisions

cssFile <- system.file("extdata", "style.css", package="ccPaper")
options(markdown.HTML.stylesheet = cssFile)

library(ccPaper)

Muscle Skin Data processing

How do we process the data for muscle and skin to something we can start running through categoryCompare?

options(stringsAsFactors=FALSE)
library(affy)

Skin

Get the raw data in (this should actually be on GEO, but doesn't seem to be there yet)

currLoc <- getwd()
dataLoc <- "/mlab/data/rmflight/Documents/projects/work/Petruska/drg_sc_skin_denervation/Skin/"
setwd(dataLoc)

attFile <- "Sample information file.txt"
dataAttributes <- read.table(attFile, sep="\t", header=T)
dataAttributes$timePoint <- 0
dataAttributes$timePoint[(dataAttributes$Time == "7 day")] <- 7
dataAttributes$timePoint[(dataAttributes$Time == "14 day")] <- 14

celData <- ReadAffy(phenoData=dataAttributes)
skinData <- rma(celData)

setwd(currLoc)
.sessionInfo <- sessionInfo()
.timeDate <- Sys.time()
save(skinData, .sessionInfo, .timeDate, file="inst/data/skin_rma_data.RData")

Muscle

library(GEOquery)
muscleData <- getGEO(GEO="GSE4411")[[1]]

controlSamples <- grep("Control", pData(muscleData)$title)
muscleData <- muscleData[,controlSamples]

tmpP <- pData(muscleData)
tmpP$innervation <- "innervated"
denLoc <- grep("Denervated", tmpP$title)
tmpP$innervation[denLoc] <- "denervated"
pData(muscleData) <- tmpP
.sessionInfo <- sessionInfo()
.timeDate <- Sys.time()
save(muscleData, .sessionInfo, .timeDate, file="inst/data/muscle_data.RData")

Gene Lists

For both the skin and muscle data sets, we will calculate the p-values using limma. We will probably attempt to do both hypergeometric enrichment and GSEA.

Do not forget that these are from different organisms. Skin: rat, muscle: mouse.

Skin Gene Lists and Enrichment

library(limma)
library(rat2302.db)
data(skin_rma_data)

skinComps <- c("T7 - T0", "T14 - T0")
geneID <- unlist(mget(featureNames(skinData), rat2302ENTREZID))

skinExpr <- exprs(skinData)
skinCollapse <- collapseProbes(skinExpr, geneID) # collapse to single genes using median of expression
skinCharacter <- pData(skinData)
skinCharacter$timePoint2 <- paste("T", skinCharacter$timePoint, sep="")

skinFC <- rankGenes(skinCollapse, skinCharacter$timePoint2, skinComps, doAggregation=FALSE, aggregateIndex=NA)
names(skinFC) <- c("T7", "T14")

skinDiff <- lapply(skinFC, getDiffGenes, id="id")

library(GO.db)
library(org.Rn.eg.db)
library(limma)
rnGO <- as.list(org.Rn.egGO2ALLEGS)
rnGO <- rnGO[(Ontology(names(rnGO)) == "BP")]
rnGO <- lapply(rnGO, unique)

skinSet <- symbols2indices(rnGO, rownames(skinCollapse))

skinStatus <- skinCharacter$timePoint2
f <- factor(skinStatus)
skinDesign <- model.matrix(~0 + f)
colnames(skinDesign) <- levels(f)

skinContrast <- makeContrasts(contrasts=skinComps, levels=skinDesign)

save(list=ls(), file="runSkinRomer.RData")

load("runSkinRomer.RData")
options(mc.cores=12) # at least if we are on hera
t1 <- Sys.time()
skinRomer <- multicontrastRomer(skinSet, skinCollapse, skinDesign, skinContrast, nrot=10000)
names(skinRomer) <- c("T7", "T14")
save(skinRomer, file="inst/data/skinRomer.RData")
t2 <- Sys.time()
difftime(t2, t1)

Run using geneSetTest

load("runSkinRomer.RData")
options(mc.cores=12) # at least if we are on hera
t1 <- Sys.time()
skinGeneSetTest <- lapply(skinFC, function(x){
  mmGeneSetTest(skinSet, x$t)
})

t2 <- Sys.time()
difftime(t2, t1)
save(skinGeneSetTest, t1, t2, file="inst/data/skinGeneSetTest.RData")

Muscle Gene Lists and Enrichment

library(limma)
library(mouse4302.db)
data(muscle_data)

muscleComps <- c("denervated - innervated")
geneID <- unlist(mget(featureNames(muscleData), mouse4302ENTREZID))

muscleExpr <- exprs(muscleData)
muscleCollapse <- collapseProbes(muscleExpr, geneID) # collapse to single genes using median of expression
muscleCharacter <- pData(muscleData)

muscleFC <- rankGenes(muscleCollapse, muscleCharacter$innervation, muscleComps, doAggregation=FALSE, aggregateIndex=NA)
names(muscleFC) <- "denervation"

muscleDiff <- lapply(muscleFC, getDiffGenes, id="id")

library(GO.db)
library(org.Mm.eg.db)
library(limma)
mmGO <- as.list(org.Mm.egGO2ALLEGS)
mmGO <- mmGO[(Ontology(names(mmGO)) == "BP")]
mmGO <- lapply(mmGO, unique)

muscleSet <- symbols2indices(mmGO, rownames(muscleCollapse))

muscleStatus <- muscleCharacter$innervation
f <- factor(muscleStatus)
muscleDesign <- model.matrix(~0 + f)
colnames(muscleDesign) <- levels(f)

muscleContrast <- makeContrasts(contrasts=muscleComps, levels=muscleDesign)

save(muscleSet, muscleCollapse, muscleDesign, muscleContrast, mmGO, muscleFC, file="runMuscleRomer.RData")

load("runMuscleRomer.RData")
options(mc.cores=12) # at least if we are on hera
t1 <- Sys.time()
muscleRomer <- multicontrastRomer(muscleSet, muscleCollapse, muscleDesign, muscleContrast, nrot=10000)
names(muscleRomer) <- "denervation"
save(muscleRomer, file="inst/data/muscleRomer.RData")
t2 <- Sys.time()
difftime(t2, t1)

Compare Skin and Muscle

data(muscleRomer)
data(skinRomer)

bothGO <- union(names(rnGO), names(mmGO))
combMapping <- lapply(bothGO, function(inName){
  isMM <- inName %in% names(mmGO)
  isRN <- inName %in% names(rnGO)

  tmpGene <- character(0)

  if (isMM){
    tmpGene <- c(tmpGene, mmGO[[inName]])
  }
  if (isRN){
    tmpGene <- c(tmpGene, rnGO[[inName]])
  }
  tmpGene <- unique(tmpGene)
  return(tmpGene)
})

names(combMapping) <- bothGO

geneAnnMapping <- new("namedList", .Data=combMapping, names=names(combMapping))


getSigID <- function(inRomer, pCut=0.05, whichCol=c("Up", "Down")){
  sigID <- lapply(whichCol, function(inCol){
    rownames(inRomer[(inRomer[, inCol] <= pCut),])
  })
  names(sigID) <- whichCol
  return(sigID)
}

muscleSigRomer <- lapply(muscleRomer, getSigID, pCut=0.01)
muscleSigRomer <- unlist(muscleSigRomer, recursive=FALSE)
names(muscleSigRomer) <- c("muscle.Up", "muscle.Down")

skinSigRomer <- lapply(skinRomer, getSigID, pCut = 0.01)
skinSigRomer <- unlist(skinSigRomer, recursive=FALSE)

allSigRomer <- c(skinSigRomer, muscleSigRomer)

genCCSigList <- function(inSig){
  tmp <- new("ccSigList", sigID=inSig)
}

allCCSig <- lapply(allSigRomer, genCCSigList)

allCCRomer <- new("GENccEnrichResult", allCCSig, categoryName="GSEAGO", geneAnnMapping=geneAnnMapping, overlapType="overlap", annDescription=Term(names(geneAnnMapping)))

allRomerOpts <- new("ccOptions", listNames=names(allCCRomer), colorType="pie")
compareSMRomer <- ccCompare(allCCRomer, allRomerOpts)
compareSMRomer

save(compareSMRomer, allRomerOpts, file="skinMuscleCompare.RData")

compareSMRomer <- breakEdges(compareSMRomer, 0.8)
cwSM <- ccOutCyt(compareSMRomer, allRomerOpts, postText="GSEAGO", rpcPort=9001)

.sessionInfo <- sessionInfo()
.timeDate <- Sys.time()
save(cwSM, compareSMRomer, allRomerOpts, .sessionInfo, .timeDate, file="inst/data/skinMuscleCompareOutput.RData")

Get out the groups of related GO terms:

smNodes <- cytOutNodes("typeII hypersensitivity - T7.Up", cwSM)
smNodes <- cytOutNodes("pH and lysosome regulation - Skin.T7.Down,Skin.T14.Up,Muscle.Up", cwSM, smNodes)
smNodes <- cytOutNodes("Neurotransmitter secretion and transport - Skin.T7.Down,Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("entry into host and movement - Skin.T7.Up", cwSM, smNodes)
smNodes <- cytOutNodes("steriod biosynthesis - Skin.T7.Down", cwSM, smNodes)
smNodes <- cytOutNodes("microtubule organization - Skin.T7.Up,Skin.T14.Down", cwSM, smNodes)
smNodes <- cytOutNodes("cellular respiration - Skin.T14.Down,Muscle.Down", cwSM, smNodes)
smNodes <- cytOutNodes("epithelial tube branching - Skin.T7.Down,Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("rRNA processing - Muscle.Up", cwSM, smNodes)
smNodes <- cytOutNodes("digestion - Skin.T14.Up", cwWM, smNodes)
smNodes <- cytOutNodes("digestion - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("response to virus - Skin.T7.Up", cwSM, smNodes)
smNodes <- cytOutNodes("purine metabolism - Muscle.Down", cwSM, smNodes)
smNodes <- cytOutNodes("protein localization in mitochondrion - Skin.T14.Down", cwSM, smNodes)
smNodes <- cytOutNodes("synapse assembly and nervous system development - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("folic acid compound metabolism - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("retinoic acid biosynthesis - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("mitotic spindle assembly - Skin.T7.Up", cwSM, smNodes)
smNodes <- cytOutNodes("regulation of phospholipase C - Muscle.Down", cwSM, smNodes)
smNodes <- cytOutNodes("acetyl CoA biosynthesis - Skin.T14.Down,Muscle.Down", cwSM, smNodes)
smNodes <- cytOutNodes("response to glucose - Skin.T7.Down, Muscle.Down", cwSM, smNodes)
smNodes <- cytOutNodes("response to VEGF - Muscle.Down", cwSM, smNodes)
smNodes <- cytOutNodes("mitotic cell cycle checkpoint - Skin.T7.Down,Skin.T14.Down,Muscle.Up", cwSM, smNodes)
smNodes <- cytOutNodes("urea cycle - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("response to prostaglandin - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("keratinocyte migration - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("atp biosynthesis - Skin.T7.Down,Skin.T14.Down,Muscle.Down", cwSM, smNodes)
smNodes <- cytOutNodes("membrane repolarization - Muscle.Down", cwSM, smNodes)
smNOdes <- cytOutNodes("muscle contraction and development - Skin.T7.Down,Skin.T14.Down,Muscle.Down", cwSM, smNodes)
smNodes <- cytOutNodes("muscle contraction and development - Skin.T7.Down,Skin.T14.Down,Muscle.Down", cwSM, smNodes)
smNodes <- cytOutNodes("muscle contraction and development - Skin.T7.Down,Skin.T14.Down,Muscle.Down", cwSM, smNodes)
smNodes <- cytOutNodes("blood vessel endothelial cell differentiation - Skin.T7.Up,Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("endocrine regulation of blood pressure - Skin.T7.Down,Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("tube development - Skin.T7.Down,Muscle.Up", cwSM, smNodes)
smNodes <- cytOutNodes("amino acid transport - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("cation channel activity - Skin.T14.Down", cwSM, smNodes)
smNodes <- cytOutNodes("interferon gamma response - Skin.T14.Down", cwSM, smNodes)
smNodes <- cytOutNodes("lung cell differentiation - Skin.T7.Down", cwSM, smNodes)
smNodes <- cytOutNodes("neg regulation of astrocyte diff. - Skin.T7.Down", cwSM, smNodes)
smNodes <- cytOutNodes("fear response - Skin.T7.Down", cwSM, smNodes)
smNodes <- cytOutNodes("N acetylglucosamine metabolism - Skin.T7.Up", cwSM, smNodes)
smNodes <- cytOutNodes("proteoglycan biosynthesis - Skin.T7.Up", cwSM, smNodes)
smNodes <- cytOutNodes("carbohydrate catabolism - Muscle.Down", cwSM, smNodes)
smNodes <- cytOutNodes("response to muscle activity - Skin.T14.Down,Muscle.Down", cwSM, smNodes)
smNodes <- cytOutNodes("response to indole 3 methanol - Muscle.Up", cwSM, smNodes)
smNodes <- cytOutNodes("negative regulation of response to granulocyte / myeloid cell diff. - Skin.T7.Up,Muscle.Up", cwSM, smNodes)
smNodes <- cytOutNodes("vitamin metabolism - Skin.T14.Up,Muscle.Up", cwSM, smNodes)
smNodes <- cytOutNodes("DNA break repair - Muscle.Up", cwSM, smNodes)
smNodes <- cytOutNodes("calcineurin NFAT signaling cascade - Skin.T14.Down", cwSM, smNodes)
smNodes <- cytOutNodes("histone demethylation - Skin.T14.Down", cwSM, smNodes)
smNodes <- cytOutNodes("negative regulation of Ras signal transduction - Skin.T14.Down", cwSM, smNodes)
smNodes <- cytOutNodes("positive regulation of metalloenzyme activity - Skin.T14.Down", cwSM, smNodes)
smNodes <- cytOutNodes("glucocorticoid receptor signaling pathway - Skin.T14.Down", cwSM, smNodes)
smNodes <- cytOutNodes("mesenchymal cell diff in kidney / renal system dev. - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("axon extension - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("lung lobe dev. and morphogenesis - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("embryonic digestive tract dev. and morphogenesis - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("negative regulation of protein transport - Skin.T7.Up,Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("phosphate ion transport - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("fluid transport - Skin.T7.Down,Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("fluid secretion - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("interleukin 13 production - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("glia guided migration - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("central nervous system maturation - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("bile acid transport - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("positive regulation of muscle cell apoptosis - Skin.T14.Up", cwSM, smNodes)
smNodes <- cytOutNodes("atrioventricular valve dev. and morphogenesis - Skin.T14.Up,Muscle.Down", cwSM, smNodes)
smNodes <- cytOutNodes("dopamine transport - Skin.T7.Down", cwSM, smNodes)
smNodes <- cytOutNodes("collateral sprouting - Skin.T7.Down", cwSM, smNodes)
smNodes <- cytOutNodes("response to osmotic stress - Skin.T7.Down", cwSM, smNodes)
smNodes <- cytOutNodes("positive regulation of epidermal growth factor signaling - Skin.T7.Down", cwSM, smNodes)
smNodes <- cytOutNodes("response to leptin - Skin.T7.Down", cwSM, smNodes)
smNodes <- cytOutNodes("catecholamin biosynthesis - Skin.T7.Down", cwSM, smNodes)
smNodes <- cytOutNodes("negative regulation of cell junction assembly - Skin.T7.Up", cwSM, smNodes)
smNodes <- cytOutNodes("endocardial cell differentiation - Skin.T7.Up", cwSM, smNodes)
smNodes <- cytOutNodes("sequestering of actin monomers - Skin.T7.Up", cwSM, smNodes)
smNodes <- cytOutNodes("beta amyloid formation - Skin.T7.Up", cwSM, smNodes)
smNodes <- cytOutNodes("response to platelet derived growth factor stimulus - Skin.T7.Up", cwSM, smNodes)
smNodes <- cytOutNodes("lipopolysaccharid biosynthesis - Skin.T7.Up", cwSM, smNodes)
smNodes <- cytOutNodes("aminoglycan metabolism - Skin.T7.Up", cwSM, smNodes)
smNodes <- cytOutNodes("epithelial to mesenchymal transition, endocardial cushion formation - Skin.T7.Up", cwSM, smNodes)
smNodes <- cytOutNodes("pulmonary valve dev. and morphogenesis - Skin.T7.Up", cwSM, smNodes)
.sessionInfo <- sessionInfo()
.timeDate <- Sys.time()
save(smNodes, cwSM, allCCRomer, allRomerOpts, compareSMRomer, file="inst/data/skinMuscleCCOutput.RData")

Lets generate a couple of different tables for output.

data(skinMuscleCCOutput)

allDescStrings <- sapply(smNodes, function(x){x$descStr})
string2List <- strsplit(allDescStrings, " - ", fixed=TRUE)
justDesc <- sapply(string2List, function(x){x[1]})
listMem <- sapply(string2List, function(x){x[2]})
listMem[1] <- "Skin.T7.Up" # fix an error in the initial setup
listMemSplit <- strsplit(listMem, ",", fixed=TRUE)

descMembershipTable <- matrix("", nrow=length(allDescStrings), ncol=7)
colnames(descMembershipTable) <- c("Description", "Skin.T7.Down", "Skin.T7.Up", "Skin.T14.Down", "Skin.T14.Up", "Muscle.Down", "Muscle.Up")

descMembershipTable[,"Description"] <- justDesc

for (inRow in seq(1, nrow(descMembershipTable))){
  useSplit <- listMemSplit[[inRow]]
  trimSplit <- gsub(" ", "", useSplit)
  descMembershipTable[inRow, trimSplit] <- "X"
}

orderBy <- c("Skin.T7.Down,Skin.T14.Down,Muscle.Down",
             "Skin.T7.Down,Skin.T14.Down,Muscle.Up",
             "Skin.T7.Down,Skin.T14.Down",
             "Skin.T7.Down,Skin.T14.Up,Muscle.Up",
             "Skin.T7.Down,Skin.T14.Up",
             "Skin.T7.Down,Muscle.Up",
             "Skin.T7.Down,Muscle.Down",
             "Skin.T7.Down",
             "Skin.T7.Up,Skin.T14.Down",
             "Skin.T7.Up,Skin.T14.Up",
             "Skin.T7.Up,Muscle.Up",
             "Skin.T7.Up",
             "Skin.T14.Down,Muscle.Down",
             "Skin.T14.Down",
             "Skin.T14.Up,Muscle.Down",
             "Skin.T14.Up,Muscle.Up",
             "Skin.T14.Up",
             "Muscle.Down",
             "Muscle.Up")

listMem <- gsub(" ", "", listMem)
newOrder <- unlist(lapply(orderBy, function(x){
  which(listMem %in% x)
}))

descMembershipTable <- descMembershipTable[newOrder,]
require(xtable)

# add an html link to each entry in the table
useLink <- paste('<a href="#loc', seq(1, nrow(descMembershipTable)), '">', descMembershipTable[,"Description"], '</a>', sep="")
descMembershipTable[, "Description"] <- useLink

Print out the table.

cat('<a name="tableLink"></a>\n') # a link back to the table
print(xtable(descMembershipTable), type = 'html', html.table.attributes = 'style="border-spacing:20px 5px;"', include.rownames=FALSE, sanitize.text.function=function(x){x})

We also want the set of GO terms, their descriptions, and list memberships so that the inquisitive mind can see how we generated the groups.

smNodes <- smNodes[newOrder]

for (iNode in seq(1:length(smNodes))){
  nodeSet <- smNodes[[iNode]]
  toLink <- paste('<a name="loc', iNode, '"></a>', sep="")
  cat(toLink, "\n")
  cat("####", nodeSet$descStr, sep=" ")
  cat("\n")

  nGO <- length(nodeSet$nodes)
  outMatrix <- matrix("", nrow=nGO, ncol=8)
  colnames(outMatrix) <- c("GOID", "Description", "T7.Down", "T7.Up", "T14.Down", "T14.Up", "muscle.Down", "muscle.Up")
  rownames(outMatrix) <- nodeSet$nodes
  outMatrix[,"GOID"] <- nodeSet$nodes

  for (useGO in nodeSet$nodes){
    outMatrix[useGO, "Description"] <- nodeSet$nodeData[[useGO]]$Desc
    listMems <- strsplit(nodeSet$nodeData[[useGO]]$listMembership, ",")[[1]]
    listMems <- listMems[(nchar(listMems) > 0)]
    outMatrix[useGO, listMems] <- "X"
  }
  allMems <- sapply(nodeSet$nodeData, function(x){x$listMembership})
  reOrder <- order(allMems, decreasing=TRUE)
  outMatrix <- outMatrix[reOrder,]
  print(xtable(outMatrix), type = 'html', html.table.attributes = 'style="border-spacing:20px 5px;"', include.rownames=FALSE)
  cat('\n<a href="#tableLink">back</a>\n')
}

Collateral Sprouting

One of our groups above seems to be related to collateral sprouting, which is the process the skin data was actually seeking to measure, but it is in an unexpected direction. We should look at the genes involved and see what is actually going on.

allDescStrings <- sapply(smNodes, function(x){x$descStr})
isCollateral <- grep("collateral sprouting", allDescStrings, ignore.case=TRUE)

collateralGO <- smNodes[[isCollateral]]$nodes

skinRomer$T7[collateralGO,]
skinRomer$T14[collateralGO,]

rnGenes <- unique(unlist(rnGO[collateralGO]))

rnFCGenes <- rownames(skinFC$T7)
rnGenes <- rnGenes[(rnGenes %in% rnFCGenes)]
rnExtra <- cbind(unlist(mget(rnGenes, org.Rn.egSYMBOL)), unlist(mget(rnGenes, org.Rn.egGENENAME)))

rnData <- lapply(skinFC, function(inFC){
  tmp <- cbind(inFC[rnGenes,], rnExtra)
  names(tmp) <- c(names(tmp)[1:7], "Symbol", "Name")
  return(tmp)
})
lapply(rnData, dim)

Lets actually generate tables for T7 and T14.

for (iTable in 1:length(rnData)){
  tmpData <- rnData[[iTable]]
  negEntries <- which(tmpData$logFC <= 0)
  htmlOut <- capture.output(print(xtable(tmpData), type = 'html', html.table.attributes = 'style="border-spacing:20px 5px;"', include.rownames=FALSE))
  negEntries <- negEntries + 4
  for (iNeg in negEntries){
    htmlOut[iNeg] <- sub("<TR>", '<TR style="color: red;">', htmlOut[iNeg])
  }

  cat('###', names(rnData)[iTable], sep=' ')
  cat('\n\n')
  cat(htmlOut, sep="\n")

  cat('\n')
}

Run using geneSetTest

load("runMuscleRomer.RData")
options(mc.cores=12) # at least if we are on hera
t1 <- Sys.time()
muscleGeneSetTest <- lapply(muscleFC, function(x){
  mmGeneSetTest(muscleSet, x$t)
})

t2 <- Sys.time()
difftime(t2, t1)
save(muscleGeneSetTest, t1, t2, file="inst/data/muscleGeneSetTest.RData")