pedprobr (version 0.9.5)

merlin: Pedigree likelihoods computed by MERLIN

Description

These functions enable users to call MERLIN (Abecasis et al., 2002) from within R.

Usage

merlin(
  x,
  options,
  markers = NULL,
  linkageMap = NULL,
  verbose = TRUE,
  generateFiles = TRUE,
  cleanup = TRUE,
  dir = tempdir(),
  logfile = NULL,
  merlinpath = NULL,
  checkpath = TRUE
)
likelihoodMerlin(
  x,
  markers = NULL,
  linkageMap = NULL,
  rho = NULL,
  logbase = NULL,
  perChrom = FALSE,
  options = "--likelihood --bits:100 --megabytes:4000 --quiet",
  ...
)
checkMerlin(program = NULL, version = TRUE, error = FALSE)

Value

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.

Arguments

x

A pedtools::ped() object.

options

A single string containing all arguments to merlin except for the input file indications.

markers

A vector of names or indices of markers attached to x. (Default: all markers).

linkageMap

A data frame with three columns (chromosome; marker name; centiMorgan position) to be used as the marker map by MERLIN.

verbose

A logical.

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.

cleanup

A logical. If TRUE (default), the MERLIN input files are deleted after the call to MERLIN.

dir

The name of the directory where input files should be written.

logfile

A character. If this is given, the MERLIN screen output will be dumped to a file with this name.

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).

checkpath

A logical indicating whether to check that the merlin executable is found.

rho

A vector of length one less than the number of markers, specifying the recombination rate between each consecutive pair.

logbase

Either NULL (default) or a positive number indicating the basis for logarithmic output. Typical values are exp(1) and 10.

perChrom

A logical; if TRUE, likelihoods are reported per chromosome.

...

Further arguments passed on to merlin().

program

A character containing "merlin", "minx" or both (default), optionally including full paths.

version

A logical. If TRUE (default), it is checked that running program produces a printout starting with "MERLIN 1.1.2".

error

A logical, indicating if an error should be raised if program is not found. Default: FALSE.

Author

Magnus Dehli Vigeland

Details

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():

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.

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)))
}