Function to calculate SPIA scores for DESeq2 results table

SPIAcalc.R

resFile <- "DESeq2_result_name" # has to have entrez ids as annotations
compName <- "comparison" # Name of the comparison

calcSPIA <- function(resFile, compName) 
{
  res <- read.table(resFile, header = T, sep = "\t", fill=T, quote = "\"", stringsAsFactors = F)
  
  # Select differentially expressed genes (adjusted p-value < 0.05)
  # This is a named vector, where names are entrez ids, and values are log2 fold changes
  ids <- res[res$padj < 0.05,]$entrezgene_id
  logFC <- res[res$padj < 0.05,]$log2FoldChange
  DE_genes <- logFC
  names(DE_genes) <- ids
  DE_genes <- DE_genes[!is.na(names(DE_genes))]
  
  # Get all expressed genes entrez ids
  ALL_genes <- res$entrezgene_id
  ALL_genes <- ALL_genes[!is.na(ALL_genes)]
  res <- spia(de=DE_genes, all=ALL_genes, organism="hsa", nB=2000, plots=FALSE, beta=NULL, combine="fisher", verbose=FALSE)
  write.table(res, file=paste(compName, "_SPIA_result.txt", sep = ""), sep = "\t", col.names = T, row.names = F, quote = F)
  
  tiff(file=paste(compName, "_SPIA_result.tiff", sep = ""), width=500, height=500)
  plotP(res, threshold = 0.05)
  dev.off()
}