Run GAGE analysis against wikipathways (modify as neccessary)

gage_wikipath.R

library(gage)
library(org.Hs.eg.db)
library(pathview)
library(rWikiPathways)
library(DESeq2)

## Get annotation from csv file (obtained via biomaRt)
annot <- read.table("../annotation/annotation.csv", sep = ",", header = T,
                    quote = "\"", fill = T, stringsAsFactors = F)

## Get wikipathways
hs.pathways <- listPathways('Homo sapiens')
hs.pathways[[1]]
ids <- sapply(hs.pathways, `[[`, 1) # get pathway ids
names <- sapply(hs.pathways, `[[`, 3) # get pathway names
path_names <- paste(ids, names)
gene_lists <- lapply(ids, function(x) getXrefList(x,'L'))
names(gene_lists) <- path_names
gene_lists

# Function to run uni-directional GAGE analysis
runGAGE <- function(exp, ct_cols, tr_cols, annotation, prefix, gene.gs) 
{
  # Change gene symbols to entrez ids
  symb <- annot[,c(2,3)]
  symb <- symb[!is.na(symb$entrezgene_id),]
  exp_eg <- merge(exp, symb, by.x=0, by.y="hgnc_symbol")
  exp_eg <- exp_eg[!duplicated(exp_eg$entrezgene_id),]
  rownames(exp_eg) <- exp_eg$entrezgene_id
  exp_eg <- exp_eg[,-c(1, length(names(exp_eg)))]
  
  # Perform gene set analysis against kegg pathways and go categories
  wikipath.p <- gage(exp_eg, gsets = gene.gs, ref = ct_cols, samp = tr_cols, 
                 set.size=c(10, 500), same.dir=TRUE, compare="unpaired", 
                 use.fold=T, rank.test=T)
  
  # Output data
  write.table(rbind(wikipath.p$greater, wikipath.p$less), 
              file = paste(prefix, "_wikipath.p.csv", sep = ""), sep=",")
  
  # Plot heatmaps
  sigGeneSet(wikipath.p, cutoff=0.05, qpval="q.val", heatmap=TRUE, pdf.size=c(14, 14), 
             cexRow=0.6, cexCol = 0.6, 
             out=paste(prefix, "_wikipath.p.RData", sep=""))
  
}


# Run GAGE analysis against wikipathways for each of the experiments
# load pre-saved DESeq object 'dds'
dds
colData(dds)
vsd <- varianceStabilizingTransformation(dds)
runGAGE(exp=assay(vsd), ct_cols=c(1:3), tr_cols=c(4:6), 
                    annotation=annot, prefix="Series1_NHBE", gene.gs=gene_lists)