R/fimo.R

In snystrom/memes: motif matching, comparison, and de novo discovery using the MEME Suite

#' Find instances of motifs using FIMO
#'
#' FIMO scans input sequences to identify the positions of matches to each input
#' motif. FIMO has no sequence length or motif number restrictions.
#'
#' @param sequences path to fasta file, or stringset input.
#' @param motifs path to .meme format file, or universalmotif/universalmotif list input.
#' @param bfile path to background file, or special values: "motif" to use
#'   0-order frequencies contained in the motif, or "uniform" to use a uniform
#'   letter distribution. (default: "motif")
#' @param outdir output directory location. Only used if text = FALSE. Default:
#'   "auto" to autogenerate directory name. Note: if not using a fasta path as
#'   input, this will be a temporary location unless explicity set.
#' @param parse_genomic_coord `logical(1)` whether to parse genomic position
#'   from fasta headers. Fasta headers must be UCSC format positions (ie
#'   "chr:start-end"), but base 1 indexed (GRanges format). If names of fasta
#'   entries are genomic coordinates and parse_genomic_coord == TRUE, results
#'   will contain genomic coordinates of motif matches, otherwise FIMO will return
#'   relative coordinates (i.e. positions from 1 to length of the fasta entry).
#' @param skip_matched_sequence `logical(1)` whether or not to include the DNA
#'   sequence of the match. Default: `FALSE`. Note: jobs will complete faster if
#'   set to `TRUE`. `add_sequence()` can be used to lookup the sequence after data import if
#'   `parse_genomic_coord` is `TRUE`, so setting this flag is not strictly needed.
#' @param max_strand if match is found on both strands, only report strand with
#'   best match (default: TRUE).
#' @param text `logical(1)` (default: `TRUE`). No output files will be created
#'   on the filesystem. The results are unsorted and no q-values are computed.
#'   This setting allows fast searches on very large inputs. When set to `FALSE`
#'   FIMO will discard 50% of the lower significance matches if >100,000 matches are
#'   detected. `text = FALSE` will also incur a performance penalty because it
#'   must first read a file to disk, then read it into memory. For these reasons,
#'   I suggest keeping `text = TRUE`.
#' @param meme_path path to `meme/bin/` (optional). Defaut: `NULL`, searches
#'   "MEME_PATH" environment variable or "meme_path" option for path to "meme/bin/".
#' @param silent `logical(1)` whether to suppress stdout/stderr printing to
#'   console (default: TRUE). If the command is failing or giving unexpected
#'   output, setting `silent = FALSE` can aid troubleshooting.
#' @param ... additional commandline arguments to fimo. See the FIMO Flag table below.
#'
#' @return GRanges object containing positions of each match. Note: if
#'   `parse_genomic_coord = FALSE`, each `seqnames` entry will be the full fasta
#'   header, and start/end will be the relative position within that sequence of the
#'   match. The GRanges object has the following additional `mcols`:
#'     * motif_id = primary name of matched motif
#'     * motif_alt_id = alternate name of matched motif
#'     * score = score of match (higher score is a better match)
#'     * pvalue = pvalue of the match
#'     * qvalue = qvalue of the match
#'     * matched_sequence = sequence that matches the motif
#' @export
#'
#' @details Additional arguments passed to `...`. See: [Fimo web manual](http://meme-suite.org/doc/fimo.html?man_type=web)
#'   for a complete description of FIMO flags.
#'
#'
#' | FIMO Flag         | allowed values | default | description                |
#' |:-----------------:|:--------------:|:-------:|:---------------------------|
#' | alpha             | `numeric(1)`   | 1       | alpha for calculating position-specific priors. Represents fraction of sites that are binding sites of TF of interest. Used in conjunction with `psp` |
#' | bfile             | "motif", "motif-file", "uniform", file path, | "motif" | If "motif" or "motif-file", use 0-order letter frequencies from motif. "uniform" sets uniform letter frequencies. |
#' | max_stored_scores | `integer(1)`   | NULL    | maximum number of scores to be stored for computing q-values. used when `text = FALSE`, see FIMO webpage for details |
#' | motif_pseudo      | `numeric(1)`   | 0.1     | pseudocount added to motif matrix |
#' | no_qvalue         | `logical(1)`   | FALSE   | only needed when `text = FALSE`, do not compute q-value for each p-value |
#' | norc              | `logical(1)`   | FALSE   | Do not score reverse complement strand |
#' | prior_dist        | file path      | NULL    | file containing binned distribution of priors |
#' | psp               | file path      | NULL    | file containing position specific priors. Requires `prior_dist` |
#' | qv_thresh         | `logical(1)`   | FALSE   | use q-values for the output threshold |
#' | thresh            | `numeric(1)`   | `1e-4`  | output threshold for returning a match, only matches with values less than `thresh` are returned. |
#'
#'
#' @details # Citation
#' If you use `runFimo()` in your analysis, please cite:
#'
#' Charles E. Grant, Timothy L. Bailey, and William Stafford Noble, "FIMO:
#' Scanning for occurrences of a given motif", Bioinformatics, 27(7):1017-1018,
#' 2011. [full text](http://bioinformatics.oxfordjournals.org/content/early/2011/02/16/bioinformatics.btr064.full)
#'
#' @details ## Licensing
#' The MEME Suite is free for non-profit use, but for-profit users should purchase a
#' license. See the [MEME Suite Copyright Page](http://meme-suite.org/doc/copyright.html) for details.
#'
#' @md
#'
#' @examples
#' if (meme_is_installed()){
#' # Generate some example input sequences
#' seq <- universalmotif::create_sequences()
#' # sequences must have values in their fasta headers
#' names(seq) <- seq_along(seq)
#' # Create random example motif to search for
#' motif <- universalmotif::create_motif()
#'
#' # Search for motif in sequences
#' # parse_genomic_coord set to FALSE since fasta headers aren't in "chr:start-end" format.
#' runFimo(seq, motif, parse_genomic_coord = FALSE)
#' }
runFimo <- function(
  sequences,
  motifs,
  bfile = "motif",
  outdir = "auto",
  parse_genomic_coord = TRUE,
  skip_matched_sequence = FALSE,
  max_strand = TRUE,
  text = TRUE,
  meme_path = NULL,
  silent = TRUE,
  ...
) {
  sequences <- sequence_input(sequences)
  motifs <- motif_input(motifs)

  if (outdir == "auto")
    outdir <- paste0(outdir_name(sequences, motifs$path), "_fimo")

  user_flags <- prepareFimoFlags(
    bfile = bfile,
    parse_genomic_coord = parse_genomic_coord,
    skip_matched_sequence = skip_matched_sequence,
    max_strand = max_strand,
    text = text,
    outdir = outdir,
    ...
  )

  flags <- c(user_flags, motifs$path, sequences)

  command <- search_meme_path(path = meme_path, util = "fimo")

  ps_out <- processx::run(command, flags, error_on_status = FALSE)

  ps_out %>%
    process_check_error(
      help_fun = ~ {
        fimoHelp(command)
      },
      user_flags = cmdfun::cmd_help_parse_flags(user_flags) %>%
        # filter out special inputs to bfile
        grep("--$", ., invert = TRUE, value = TRUE),
      flags_fun = ~ {
        gsub("-", "_", .x)
      },
      default_help_fun = TRUE
    )

  print_process_stdout(ps_out, silent = silent)
  print_process_stderr(ps_out, silent = silent)

  if (!is.na(text)) {
    if (text) {
      fimo_res <- ps_out$stdout %>%
        I() %>%
        parseFimo()
      return(fimo_res)
    }
  }

  fimo_out <- cmdfun::cmd_file_combn("fimo", "tsv", outdir = outdir)

  #print(fimo_out$tsv)
  fimo_out$tsv %>%
    parseFimo()
}

#' Parse flags for FIMO input
#'
#' @param bfile background file
#' @param parse_genomic_coord parse genomic coordinates from name
#' @param skip_matched_sequence skip matched
#' @param max_strand max strand
#' @param text text
#' @param outdir outdir
#' @param ... ...
#'
#' @return
#'
#' @noRd
prepareFimoFlags <- function(
  bfile,
  parse_genomic_coord,
  skip_matched_sequence,
  max_strand,
  text,
  outdir = outdir,
  ...
) {
  argsDict <- c("outdir" = "oc")

  flags <- cmdfun::cmd_args_all() %>%
    cmdfun::cmd_list_interp(argsDict) %>%
    purrr::set_names(
      ~ {
        gsub("_", "-", .x)
      }
    )

  if (!is.null(bfile)) {
    if (file.exists(bfile) & bfile %in% c("motif", "uniform")) {
      message(paste0(
        "Working directory contains a file named: ",
        bfile,
        " that will be used as the `bfile` argument.",
        "To force use of the FIMO keyword, pass argument as: `bfile = --",
        bfile,
        "--`"
      ))
    }
    if (!file.exists(bfile) & bfile %in% c("motif", "uniform")) {
      # if bfile isn't a path, user is probably inputting special keywords which
      # get wrapped in '--', but first drop all "-" in bfile input in case user
      # added '--' already.
      flags$bfile %<>% gsub("-", "", .) %>% gsub("(.+)", "--\\1--", .)
    }
  }

  flags %>%
    cmdfun::cmd_list_to_flags(prefix = "--")
}

#' Returns help string for FIMO
#'
#' @param command
#'
#' @return
#'
#' @noRd
fimoHelp <- function(command) {
  processx::run(command, error_on_status = FALSE)$stderr
}

#' Import fimo matches as GRanges object.
#'
#' @param fimo_tsv path to fimo.tsv output
#'
#' @return GenomicRanges object for each match position. Note: if
#'   parse_genomic_coord == FALSE, each `seqnames` entry will be the fasta
#'   header, and start/end will be the position within that sequence of the
#'   match.
#'
#' @noRd
parseFimo <- function(fimo_tsv) {
  fimo_lines <- tryCatch(
    readr::read_tsv(
      fimo_tsv,
      comment = "#",
      col_types = c(
        "motif_id" = "c",
        "motif_alt_id" = "c",
        "sequence_name" = "c",
        "start" = "i",
        "stop" = "i",
        "strand" = "c",
        "score" = "d",
        "p-value" = "d",
        "q-value" = "d",
        "matched_sequence" = "c"
      )
    ),
    error = function(e) {
      # move the empty string check down
      if (fimo_tsv == "") {
        return(NULL)
      }
      stop(e)
    },
    warning = function(w) {
      # If the above file import fails w/ warning
      # (usually because fimo.tsv is empty)
      # Double check that the file is actually empty
      # (no other lines except for comments & empty lines)
      # Then return NULL if that's the case.
      lines <- readr::read_lines(fimo_tsv) %>%
        grep("^#", ., invert = TRUE, value = TRUE)
      line_lengths <- vapply(lines, nchar, integer(1))

      if (any(line_lengths) > 0) {
        stop(paste("Error reading file:", fimo_tsv))
      }

      return(NULL)
    }
  )

  # NULL is only returned above if fimo_tsv is empty, therefore no matches
  if (is.null(fimo_lines)) {
    message("No matches were detected")
    return(NULL)
  }

  fimo_matches <- fimo_lines %>%
    dplyr::rename_all(
      ~ {
        gsub("-", "", .)
      }
    ) %>%
    dplyr::rename("seqnames" = "sequence_name") %>%
    # NOTE: FIMO uses 1-based coordinates, so no need to shift for GRanges conversion
    GenomicRanges::GRanges()
  # Compute q-value?
  #dplyr::group_by(motif_alt_id) %>%
  #dplyr::mutate(q.value = p.value/n())
}

#' Import FIMO results
#'
#' @param fimo_tsv path to fimo.tsv output file
#'
#' @return GenomicRanges object for each match position. Note unless coordinates
#'   are genomic positions, each `seqnames` entry will be the fasta header, and
#'   start/end will be the position within that sequence of the match.
#' @export
#'
#' @examples
#' fimo_tsv <- system.file("extdata", "fimo.tsv", package = "memes")
#' importFimo(fimo_tsv)
importFimo <- function(fimo_tsv) {
  parseFimo(fimo_tsv)
}