R/merlin.R

In pedprobr: Probability Computations on Pedigrees

#' Pedigree likelihoods computed by MERLIN
#'
#' These functions enable users to call MERLIN (Abecasis et al., 2002) from
#' within R.
#'
#' For these functions to work, the program MERLIN must be installed (see link
#' in the Reference section below) and correctly pointed to in the PATH
#' variable. The `merlin()` function is a general wrapper which runs MERLIN with
#' the indicated options, after creating the appropriate input files. For
#' convenience, MERLIN's "--likelihood" functionality is wrapped in a separate
#' function.
#'
#' The `merlin()` function creates input files "_merlin.ped", "_merlin.dat",
#' "_merlin.map" and "_merlin.freq" in the `dir` directory, and then runs the
#' following command through a call to [system()]:
#'
#' \preformatted{merlin -p _merlin.ped -d _merlin.dat -m _merlin.map -f
#' _merlin.freq  <options> }
#'
#' `likelihoodMerlin()` first runs `merlin()` with `options = "--likelihood
#' --bits:100 --megabytes:4000 --quiet"`, and then extracts the likelihood
#' values from the MERLIN output. Note that the output is the *total* likelihood
#' including all markers.
#'
#' For likelihood computations with linked markers, the argument `rho` should
#' indicate the recombination fractions between each consecutive pair of markers
#' (i.e., `rho[i]` is the recombination rate between markers `i-1` and `i`).
#' These will be converted to centiMorgan distances using Haldane's map
#' function, and used to create genetic marker map in a MERLIN-friendly format.
#'
#' @param x A [pedtools::ped()] object.
#' @param options A single string containing all arguments to merlin except for
#'   the input file indications.
#' @param markers A vector of names or indices of markers attached to `x`.
#'   (Default: all markers).
#' @param linkageMap A data frame with three columns (chromosome; marker name;
#'   centiMorgan position) to be used as the marker map by MERLIN.
#' @param verbose A logical.
#' @param generateFiles A logical. If TRUE (default), input files to MERLIN
#'   named '_merlin.ped', '_merlin.dat', '_merlin.map', and '_merlin.freq' are
#'   created in the directory indicated by `dir`. If FALSE, no files are
#'   created.
#' @param cleanup A logical. If TRUE (default), the MERLIN input files are
#'   deleted after the call to MERLIN.
#' @param dir The name of the directory where input files should be written.
#' @param logfile A character. If this is given, the MERLIN screen output will
#'   be dumped to a file with this name.
#' @param merlinpath The path to the folder containing the merlin executables.
#'   If the executables are on the system's search path, this can be left as
#'   NULL (default).
#' @param checkpath A logical indicating whether to check that the merlin
#'   executable is found.
#' @param rho A vector of length one less than the number of markers, specifying
#'   the recombination rate between each consecutive pair.
#' @param logbase Either NULL (default) or a positive number indicating the
#'   basis for logarithmic output. Typical values are `exp(1)` and 10.
#' @param program A character containing "merlin", "minx" or both (default),
#'   optionally including full paths.
#' @param version A logical. If TRUE (default), it is checked that running
#'   `program` produces a printout starting with "MERLIN 1.1.2".
#' @param error A logical, indicating if an error should be raised if `program`
#'   is not found. Default: FALSE.
#'
#' @return `merlin()` returns the screen output of MERLIN invisibly.
#'
#'   `likelihoodMerlin()` returns a single number; the total likelihood using
#'   all indicated markers.
#'
#'   `checkMerlin()` returns TRUE if the MERLIN executable indicated by
#'   `program` is found on the system. Otherwise FALSE, or (if `error = TRUE`)
#'   an error is raised.
#'
#' @author Magnus Dehli Vigeland
#' @references Abecasis et al. (2002) Nat Gen 30:97-101.
#'   <https://csg.sph.umich.edu/abecasis/merlin/>.
#'
#' @examples
#'
#' if(checkMerlin()) {
#'
#' ### Trivial example for validation
#' x = nuclearPed(1) |>
#'   addMarker("1" = "1/2") |>            # likelihood = 1/2
#'   addMarker("1" = "1/1", "3" = "1/2")  # likelihood = 1/8
#'
#' # MERLIN likelihoods
#' lik1 = likelihoodMerlin(x, markers = 1, verbose = FALSE)
#' lik2 = likelihoodMerlin(x, markers = 2, verbose = FALSE)
#' likTot = likelihoodMerlin(x, verbose = FALSE)
#' stopifnot(all.equal(
#'   round(c(lik1, lik2, likTot), c(3,3,4)), c(1/2, 1/8, 1/16)))
#'
#' # Example with ped lists
#' y = singletons(1:2) |>
#'   addMarker(`1` = "1/2", `2` = "1/1", alleles = 1:2)
#' lik = likelihoodMerlin(y, verbose = FALSE)
#' stopifnot(all.equal(round(lik, 3), 1/8))
#'
#' ### Linked markers
#' z = nuclearPed(2)
#' m = marker(z, geno = c("1/1", "1/2", "1/2", "1/2"))
#' z = setMarkers(z, list(m, m))
#'
#' # By MERLIN...
#' L1 = likelihoodMerlin(z, markers = 1:2, rho = 0.25, verbose = FALSE)
#'
#' # ...and by pedprobr
#' L2 = likelihood2(z, marker1 = 1, marker2 = 2, rho = 0.25)
#'
#' # stopifnot(all.equal(signif(L1, 3), signif(L2, 3)))
#' }
#'
#' @export
merlin = function(x, options, markers = NULL, linkageMap = NULL, verbose = TRUE,
                  generateFiles = TRUE, cleanup = TRUE, dir = tempdir(),
                  logfile = NULL, merlinpath = NULL, checkpath = TRUE) {

  # Select markers
  if (!hasMarkers(x))
    stop2("Pedigree has no attached markers")

  # Set linkage map if provided
  if(!is.null(linkageMap))
    x = setMap(x, map = linkageMap, matchNames = NA)

  if(is.null(markers))
    markers = seq_len(nMarkers(x))
  x = selectMarkers(x, markers)

  # MERLIN or MINX?
  xchrom = isXmarker(x)
  if(all(xchrom)) {
    if(verbose) cat("All markers are X-linked; calling MINX\n")
    program = "minx"
  }
  else if(all(!xchrom)) {
    program = "merlin"
  }
  else
    stop2("Both autosomal and X-linked markers are selected\n",
          "Please use the `markers` argument to run these in separate calls")

  if(!is.null(merlinpath))
    program = file.path(merlinpath, program)

  if(checkpath)
    checkMerlin(program, version = FALSE, error = TRUE)

  prefix = file.path(dir, "_merlin")

  # Generate input files to MERLIN/MINX
  if (generateFiles) {
    files = writePed(x, prefix = prefix, merlin = TRUE, verbose = verbose)

    if(cleanup)
      on.exit({unlink(files); if (verbose) cat("MERLIN input files removed\n")})
  }

  commandArgs = c(program,
                  sprintf("-p %s.ped", prefix),
                  sprintf("-d %s.dat", prefix),
                  sprintf("-m %s.map", prefix),
                  sprintf("-f %s.freq", prefix),
                  options)
  command = paste(commandArgs, collapse = " ")

  if (verbose)
    cat("\nExecuting the following command:\n", paste0(commandArgs, collapse = "\n "), "\n", sep = "")

  # Run MERLIN and store output
  mout = suppressWarnings(system(command, intern = TRUE))

  # Write logfile if indicated
  if (!is.null(logfile))
    write(mout, logfile)

  err = NULL
  if (any(fatal <- substr(mout, 1, 11) == "FATAL ERROR"))
    err = mout[which(fatal)[1]:length(mout)]
  else if (any(warn <- substr(mout, 2, 8) == "WARNING"))
    err = mout[which(warn)[1] + 0:5]

  if(!is.null(err))
    warn2(err)
  else if(verbose)
    cat("\nMERLIN run completed\n")

  invisible(mout)
}

#' @param perChrom A logical; if TRUE, likelihoods are reported per chromosome.
#' @param ... Further arguments passed on to `merlin()`.
#'
#' @rdname merlin
#' @export
likelihoodMerlin = function(x, markers = NULL, linkageMap = NULL, rho = NULL,
                            logbase = NULL, perChrom = FALSE,
                            options = "--likelihood --bits:100 --megabytes:4000 --quiet",
                            ...) {

  # Select markers
  x = selectMarkers(x, markers %||% seq_len(nMarkers(x)))

  # If rho given, replace map
  if(!is.null(rho)) {

    if(!is.null(linkageMap))
      stop2("At least one of `rho` and `linkageMap` must be NULL")

    if(length(rho) != nMarkers(x) - 1)
      stop2("Argument `rho` must have length one less than the number of markers")

    # Avoid infinities
    rho[rho == 0.5] = haldane(cM = 500)

    # If no chromosome info given, place all markers on chrom 1
    chr = chrom(x)
    if(all(is.na(chr)))
      chr = rep_len(1, nMarkers(x))

    # Convert to centiMorgan positions
    cm = c(0, haldane(rho = rho))

    linkageMap = data.frame(CHROM = chr, MARKER = name(x), CM = cm)
  }

  # Run MERLIN
  mout = merlin(x, linkageMap = linkageMap, options = options, ...)

  if(any(skipped <- substr(mout, 3, 9) == "SKIPPED"))
    stop2(paste(mout[sort(c(which(skipped) - 1, which(skipped)))], collapse = "\n"))

  nFam = if(is.pedList(x)) length(x) else 1
  likPat = sprintf("^lnLikelihood for %d families", nFam)

  parse_lnlik = function(s) {
    m = regmatches(s, regexec("=\\s*([-+]?\\d*\\.?\\d+(?:[eE][-+]?\\d+)?)", s))[[1]]
    if (length(m) < 2) NA_real_ else as.numeric(m[2])
  }

  chromLines = which(grepl("^Analysing Chromosome", mout))

  # Single chromosome (no chrom info in file)
  if (length(chromLines) == 0) {
    likLine = which(grepl(likPat, mout))
    res0 = if (length(likLine) == 0) -Inf else parse_lnlik(mout[likLine[1]])
    res = fixMerlinLog(res0, logbase = logbase)

    badLines = mout[grepl("^\\s*Skipping Marker\\b.*\\[BAD INHERITANCE\\]", mout)]

    if(length(badLines)) {
      marker = sub("^\\s*Skipping Marker\\s+(\\S+).*", "\\1", badLines)
      if(isTRUE(ncol(linkageMap) > 2)) # safeguard
        chrom = linkageMap[,1][match(marker, linkageMap[,2])]
      else
        chrom = rep(NA_character_, length(marker))
      diagn = data.frame(chrom = chrom, marker = marker, status = "bad_inheritance")
      attr(res, "bad") = diagn
      warn2("BAD INHERITANCE: ", diagn$marker)
    }
    return(res)
  }

  chromIds = trimws(sub("^Analysing Chromosome\\s*", "", mout[chromLines]))
  endLines = c(chromLines[-1] - 1, length(mout))

  lnliks = rep(NA_real_, length(chromLines))
  badMarkers = vector("list", length(chromLines))

  for (i in seq_along(chromLines)) {
    block = mout[chromLines[i]:endLines[i]]

    badLines = block[grepl("^\\s*Skipping Marker\\b.*\\[BAD INHERITANCE\\]", block)]
    if (length(badLines))
      badMarkers[[i]] = sub("^\\s*Skipping Marker\\s+(\\S+).*", "\\1", badLines)

    likIdx = which(grepl(likPat, block))
    if (length(likIdx))
      lnliks[i] = parse_lnlik(block[likIdx[1]])
  }

  names(lnliks) = chromIds

  names(badMarkers) = chromIds
  badMarkers = badMarkers[lengths(badMarkers) > 0]

  diagn = NULL
  if(length(badMarkers))
    diagn = data.frame(chrom = rep(names(badMarkers), lengths(badMarkers)),
                       marker = unlist(badMarkers, use.names = FALSE),
                       status = "bad_inheritance")

  skippedChrom = chromIds[is.na(lnliks)]
  if(length(skippedChrom))
    warn2("MERLIN produced no likelihood for chromosome(s):", skippedChrom)
  if(!is.null(diagn))
    warn2("BAD INHERITANCE: ", diagn$marker)

  # Returns
  if(perChrom)
    res = fixMerlinLog(lnliks, logbase = logbase)
  else {
    res0 = if(anyNA(lnliks)) -Inf else sum(lnliks)
    res = fixMerlinLog(res0, logbase = logbase)
  }

  if(!is.null(diagn)) attr(res, "bad") = diagn
  res
}


#' @rdname merlin
#' @export
checkMerlin = function(program = NULL, version = TRUE, error = FALSE) {

  # Default NULL for back compatibility
  if(is.null(program))
    return(checkMerlin("merlin", version = version, error = error) &&
             checkMerlin("minx", version = version, error = error))

  if(!is.character(program) || length(program) != 1)
    stop2("`program` must be a character of length 1")

  if(!nzchar(pth <- Sys.which(program))) {
    if(error)
      stop2("Executable not found. Use `merlinpath` to supply the path to the MERLIN folder")
    return(FALSE)
  }

  if(version) {
    out = tryCatch(
      suppressWarnings(system(program, intern = TRUE)),
      error = function(e) if(error) stop2(e) else return(FALSE))

    if(!isTRUE(startsWith(out[1], "MERLIN 1.1.2"))) {
      if(error)
        stop2(pth, "\nExpected printout to start with 'MERLIN 1.1.2', but received:\n", paste(out, collapse = "\n"))
      return(FALSE)
    }
  }

  TRUE
}


fixMerlinLog = function(a, logbase = NULL) {
  if(!is.null(logbase)) {
    if(length(logbase) != 1 || !is.numeric(logbase) || logbase <= 0)
      stop2("`logbase` must be a positive number: ", logbase)

    if(logbase == exp(1))
      res = a
    else
      res = round(a/log(logbase),3)
  }
  else {
    res = signif(exp(a), digits = 3)
    uflow = res == 0 & a > -Inf
    if(any(uflow))
      warning("Underflow!\nSome lnLikelihoods reported by MERLIN are too small to be exp'ed.",
              immediate. = TRUE, call. = FALSE)
    oflow = res == Inf & a < Inf
    if(any(oflow))
      warning("Overflow!\nSome lnLikelihoods reported by MERLIN are too large to be exp'ed.",
              immediate. = TRUE, call. = FALSE)
  }
  res
}