Scripts used for running SPIA on a data set.

entrez.r

#This script converts the Agilent prob ids to entrez ids

library(gdata)

library(RCurl)

getentrez <- function(hgnc)
{
  content = getURLContent( paste0( "http://skimlab.tgen.org:3000/convert?key=", hgnc ) ) #Query conversion server
  id = strsplit(content, "\t", fixed=TRUE)[[1]][1] #Get first id returned
  
  return(id) #Return that
}

getentrezlist <- function(list )
{
  returnval <- vector()
  
  len = length(list)
  
  for(i in 1:len )
  {
    returnval <- append( returnval, getentrez(list[i] ) )
  }
  
  return(returnval)
}

print("Converting Agilent Chip...")
all_genes <- as.vector( read.delim("agilent.txt", header=FALSE)[,1] ) #Read list of all genes on microarray
all_entrez <- getentrezlist( all_genes )

#Convert Gene Symbols to Entrez ID's
print("Converting Classical...")
classicalresults$ENTREZ <- getentrezlist(classicalresults$ID )

print("Converting Mesenchymal...")
mesenchymalresults$ENTREZ <- getentrezlist( mesenchymalresults$ID )

print("Converting Neural...")
neuralresults$ENTREZ <- getentrezlist( neuralresults$ID )

print("Converting Proneural...")
proneuralresults$ENTREZ <- getentrezlist( proneuralresults$ID )

#Write output to "out folder"
write.table(classicalresults, "out/classicalresults.tab", sep = "\t", row.names = FALSE)
write.table(mesenchymalresults, "out/mesenchymalresults.tab", sep = "\t", row.names = FALSE)
write.table(neuralresults, "out/neuralresults.tab", sep = "\t", row.names = FALSE)
write.table(proneuralresults, "out/proneuralresults.tab", sep = "\t", row.names = FALSE)

exec.r

#This script creates loads the data matrix and uses limma to find the list of differentially expressed genes for each of the subtypes

setwd("C:/Ivy/SPIA_core") #Set working directory
core <- read.delim("output.tab", row.names = 1) #Read in data matrix

library(limma)

classicalcolumn = c( rep(1, 38), rep(0, 56), rep(0, 26), rep(0, 53) )
mesenchymalcolumn = c( rep(0, 38), rep(1, 56), rep(0, 26), rep(0, 53) )
neuralcolumn = c( rep(0, 38), rep(0, 56), rep(1, 26), rep(0, 53) )
proneuralcolumn = c( rep(0, 38), rep(0, 56), rep(0, 26), rep(1, 53) )

design <- cbind( rep(1,173), rep(0,173) )
colnames(design) <- c("WT", "MU")
rownames(design) <- colnames(core)

classicaldesign = design
classicaldesign[, 2] <- classicalcolumn

mesenchymaldesign = design
mesenchymaldesign[, 2] <- mesenchymalcolumn

neuraldesign = design
neuraldesign[, 2] <- neuralcolumn

proneuraldesign = design
proneuraldesign[, 2] <- proneuralcolumn

#P-Value cutoff
p_value <- 0.05

classicalfit <- lmFit(core, classicaldesign)
classicalfit <- eBayes(classicalfit)
classicalresults = topTable(classicalfit, coef="MU", adjust="BH", number=20000, p.value=p_value, lfc=1)

mesenchymalfit <- lmFit(core, mesenchymaldesign)
mesenchymalfit <- eBayes(mesenchymalfit)
mesenchymalresults = topTable(mesenchymalfit, coef="MU", adjust="BH", number=20000, p.value=p_value, lfc=1)

neuralfit <- lmFit(core, neuraldesign)
neuralfit <- eBayes(neuralfit)
neuralresults = topTable(neuralfit, coef="MU", adjust="BH", number=20000, p.value=p_value, lfc=1)

proneuralfit <- lmFit(core, proneuraldesign)
proneuralfit <- eBayes(proneuralfit)
proneuralresults = topTable(proneuralfit, coef="MU", adjust="BH", number=20000, p.value=p_value, lfc=1)

#Write output to "out folder"
write.table(classicalresults, "out/classicalresults.tab", sep = "\t", row.names = FALSE)
write.table(mesenchymalresults, "out/mesenchymalresults.tab", sep = "\t", row.names = FALSE)
write.table(neuralresults, "out/neuralresults.tab", sep = "\t", row.names = FALSE)
write.table(proneuralresults, "out/proneuralresults.tab", sep = "\t", row.names = FALSE)

spia.r

# This script actually runs the DEG's through SPIA, and saves the result to a file

library(SPIA)

DE_Classical = classicalresults$logFC
DE_Mesenchymal = mesenchymalresults$logFC
DE_Neural = neuralresults$logFC
DE_Proneural = proneuralresults$logFC

names( DE_Classical ) <- as.vector(classicalresults$ENTREZ)
names( DE_Mesenchymal ) <- as.vector(mesenchymalresults$ENTREZ)
names( DE_Neural ) <- as.vector(neuralresults$ENTREZ)
names( DE_Proneural ) <- as.vector(proneuralresults$ENTREZ)

res_classical = spia(de=DE_Classical,all=all_entrez,organism="hsa",nB=2000,plots=FALSE,beta=NULL,combine="fisher",verbose=TRUE)
res_mesenchymal = spia(de=DE_Mesenchymal,all=all_entrez,organism="hsa",nB=2000,plots=FALSE,beta=NULL,combine="fisher",verbose=TRUE)
res_neural = spia(de=DE_Neural,all=all_entrez,organism="hsa",nB=2000,plots=FALSE,beta=NULL,combine="fisher",verbose=TRUE)
res_proneural = spia(de=DE_Proneural,all=all_entrez,organism="hsa",nB=2000,plots=FALSE,beta=NULL,combine="fisher",verbose=TRUE)

#Write out pathway rankings
write.table(res_classical, "out/pathways_classical.tab", sep = "\t", row.names = FALSE)
write.table(res_mesenchymal, "out/pathways_mesenchymal.tab", sep = "\t", row.names = FALSE)
write.table(res_neural, "out/pathways_neural.tab", sep = "\t", row.names = FALSE)
write.table(res_proneural, "out/pathways_proneural.tab", sep = "\t", row.names = FALSE)