R/limma_functions.R

#' create an ExpressionSet for a series matrix file
#' @param gse a valid GSE value, e.g.: "GSE8921"
#' @param path2geo_data full or relative where GEO data for this GSE is stored
#' or will be stored. The default value is "GEOtemp"
#' @export
get_eset <- function(gse, path2geo_data){
  eset <- getGEO(gse, destdir = path2geo_data)[[1]]
  exprs <- exprs(eset)
  return(exprs)
}


#' cretes an ExpressionSet from a customized series matrix file
#' @param path2table relative or full path to the table
#' @param delim delimiter of the table file. The default is "," for csv files.
#' This should be changed to "\t" for tsv files or to whatever delimiter used
#' in the target file.
get_eset_from_table <- function(path2table, delim = ","){
  exprs <- as.matrix(read.table(path2table,
                                header = TRUE,
                                sep = delim,
                                row.names = 1,
                                as.is = TRUE))
  exprs <- ExpressionSet(assayData = exprs)
  return(exprs)
}


#' create limma design object based on groups in a table called samples
#' @param groups input file based on sample_review.txt file. The default value
#' is samples$Group
#' @export
get_design <- function(groups = samples$Group){
  groups <- factor(na.omit(groups))
  design <- model.matrix(~0 +groups)
  colnames(design) <- levels(groups)
  return(design)
}

#' create a contrast.matrix to be used with limma's fit function
#' @param path2contrast_file the default is contrast.txt which is a table
#' that contains the contrasts. It should have a column named Formula in which
#' the contrast based on groups are given and anoter column named Contrast
#' which gives a title for each contrast.
#' @param design this is a design object created with teh get_design() function
#' its default value is design
#' @param lmFit_out output of lmFit function
#' @export
get_contrast.matrix <- function(path2contrast_file = "contrasts.txt", my_design = design, lmFit_out = fit){
  contrasts_table <- fread(path2contrast_file)
  contrast.matrix <- limma::makeContrasts(contrasts = contrasts_table$Formula,
                                   levels=lmFit_out$design)
  colnames(contrast.matrix) <- contrasts_table$Contrast
  rownames(contrast.matrix) <- colnames(my_design)
  return(contrast.matrix)
}

#' runs limma's lmFit, contrast.fit, and eBayes in succession
#' @param exprs a ExpressionSet default value is exprs
#' @param design this is a design object created with teh get_design() function
#' its default value is design
#' @param path2contrast_file the default is contrast.txt which is a table
#' that contains the contrasts. It should have a column named Formula in which
#' the contrast based on groups are given and anoter column named Contrast
#' which gives a title for each contrast.
fitting <- function(my_eset = exprs, design = my_design, path2contrast_file = "contrast.txt"){
  print("running lmFit")
  fit <- lmFit(my_eset, design)
  print("getting contrasts")
  #contrasts <- get_contrast.matrix()
  contrasts_table <- fread(path2contrast_file)
  contrast.matrix <- limma::makeContrasts(contrasts = contrasts_table$Formula,
                                          levels = fit$design)
  colnames(contrast.matrix) <- contrasts_table$Contrast
  rownames(contrast.matrix) <- colnames(my_design)

  print("running contrast.fit")
  fit2 <- contrasts.fit(fit, contrasts)
  print("running eBayes")
  fit3 <- eBayes(fit2)
  print("all done")
  return(fit3)
}

#' writes the result of the fitting() function to a text file.
#' @param fit a fit object created by the fitting() function. The default value
#' is fit
#' @param out_path full or relativ path to the desired location for the output
#' file. The default value is "Results.fit.txt
write_fit_file <- function(fit = fit, out_path = "Results/fit.txt"){
  write.fit(fit, results=NULL, file=out_path, digits=10, adjust="fdr", method="separate",F.adjust="none", sep="\t")
}
axelmuller/geometric2 documentation built on May 25, 2019, 6:26 p.m.