R/DEm-7-wrapper.R

In th789/ffl: Identifies Feed Forward Gene Networks

#file: DEm-7-wrapper-fn-v2-percent-denom-upordown-genes-not-all-DE-genes
#removed: code for MM subgroups
#removed: setwd("/Users/than/Dropbox/research/miRNA/ffl/DE-method")
#removed: loading files

# wrapper function --------------------------------------------------------

# wrapper function --------------------------------------------------------

#ffl_analysis: wraps all the ffl analysis steps into one function
#' @title ffl_analysis
#' @description wraps FFL analysis steps into one function
#'
#' @param mirna_expr miRNA expression matrix (miRNAs x samples):
#' matrix where rows are miRNAs and columns are samples
#' (row names are miRNA names in hsa format and column names are sampleIDs)
#' @param mrna_expr mRNA expression matrix (genes x samples):
#' matrix where rows are genes and columns are samples
#' (row names are gene names and column names are sampleIDs)
#' @param pheno phenotype matrix (samples x 1):
#' matrix where rows are samples and there is column ("group") that indicates group membership of samples
#' (row names are sampleIDs, column name is "group)
#' @param ffl_fdr_cutoff FDR cutoff for FFL gene enrichment analysis
#'
#' @return list with all analysis results (length 11: data, arm1, arm2, arm3a, arm3b,
#' ffl1, ffl2, ge1, ge2, plot_ffl1, plot_ffl2)
#' @export

ffl_analysis <- function(mirna_expr, mrna_expr, pheno, ffl_fdr_cutoff){
  #####arm1
  print("step 1/7 -- arm1: find miRNA-gene pairs")
  arm1 <- DEm_arm1_miRNA_gene(mirna_expr = mirna_expr, mrna_expr = mrna_expr, pheno = pheno)

  #####arm2
  print("step 2/7 -- arm2: find TF-gene pairs")
  arm2 <- DEm_arm2_tf_gene(mrna_DE = arm1$mrna_DE)

  #####arm3a
  print("step 3/7 -- arm3a: find miRNA-TF pairs")
  arm3a <- DEm_arm3a_miRNA_tf(miRNA_gene_pairs = arm1$db_hits)

  #####arm3b
  print("step 4/7 -- arm3b: find TF-miRNA pairs")
  arm3b <- DEm_arm3b_tf_miRNA(mrna_DE = arm1$mrna_DE, mirna_DE = arm1$mirna_DE)

  #####assemble ffl
  print("step 5/7 -- assemble FFLs (FFL1 & FFL2)")
  #ffl1
  ffl1 <- DEm_assemble_ffl(miRNA_gene_pairs = arm1$db_hits, tf_gene_pairs = arm2$db_hits, miRNA_tf_pairs = arm3a, ffl_type = "ffl1")
  #ffl2
  ffl2 <- DEm_assemble_ffl(miRNA_gene_pairs = arm1$db_hits, tf_gene_pairs = arm2$db_hits, miRNA_tf_pairs = arm3b$db_hits, ffl_type = "ffl2")

  #####gene enrichment analysis & plotting
  print("step 6/7 -- gene enrichment analysis (FFL1 & FFL2)")
  print("step 7/7 -- plot results (FFL1 & FFL2)")

  #ffl1 -- case 1: no ffl1
  if(dim(ffl1)[1] == 0){
    ge1 <- "no ffl1"
    plot_ffl1 <- "no ffl1"
    print("no ffl1")}
  #ffl1 -- case 2: there is at least one ffl1
  if(dim(ffl1)[1] >= 1){
    ge1 <- DEm_ge_analysis(arm1 = arm1, arm2 = arm2, ffls = ffl1)
    plot_ffl1 <- DEm_summary_plot(arm1 = arm1, arm2 = arm2, ge = ge1, ffl_fdr_cutoff = ffl_fdr_cutoff)
  }
  #ffl2 -- case 1: no ffl2
  if(dim(ffl2)[1] == 0){
    ge2 <- "no ffl2"
    plot_ffl2 <- "no ffl2"
    print("no ffl2")}
  #ffl2 -- case 2: there is at least one ffl2
  if(dim(ffl2)[1] >= 1){
    ge2 <- DEm_ge_analysis(arm1 = arm1, arm2 = arm2, ffls = ffl2)
    plot_ffl2 <- DEm_summary_plot(arm1 = arm1, arm2 = arm2, ge = ge2, ffl_fdr_cutoff = ffl_fdr_cutoff)
  }

  print("***completed all analyses***")

  return(list(arm1 = arm1, arm2 = arm2, arm3a = arm3a, arm3b = arm3b,
              ffl1 = ffl1, ffl2 = ffl2, ge1 = ge1, ge2 = ge2,
              plot_ffl1 = plot_ffl1, plot_ffl2 = plot_ffl2))
}