Simple Limma Function

runlimma.R

library(limma)
library(edgeR)

# function to return the toptable from limma
# indata: raw count matrix for all kras comparisons
# subsetcols: column indices for columns that you want to subset
# firstannot: the name of the first annotation
# secondannot: the name of the second annotation
get_toptable <- function(indata, subsetcols, firstannot, secondannot){
  subset_gfp_wt <- indata[,subsetcols]
  design <- model.matrix(~factor(c(rep(firstannot,9),rep(secondannot,9))))
  subset_gfp_wt <- subset_gfp_wt[rowSums(subset_gfp_wt) != 0,]
  dge <- DGEList(counts=subset_gfp_wt)
  dge <- calcNormFactors(dge)
  logCPM <- cpm(dge, log=TRUE, prior.count=3)
  fit <- lmFit(logCPM, design)
  fit <- eBayes(fit, trend=TRUE)
  return(topTable(fit, coef=ncol(design), number = nrow(subset_gfp_wt)))
}