merlin: Pedigree likelihood computed by MERLIN
In pedprobr: Probability Computations on Pedigrees

Description Usage Arguments Details Value Author(s) References Examples

View source: R/merlin.R

For these functions to work, the program MERLIN (see References below) must be installed 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.

merlin(
  x,
  options,
  markers = NULL,
  linkageMap = NULL,
  verbose = TRUE,
  generateFiles = TRUE,
  cleanup = TRUE,
  dir = tempdir(),
  logfile = NULL,
  merlinpath = NULL
)

likelihoodMerlin(
  x,
  markers = NULL,
  linkageMap = NULL,
  rho = NULL,
  logbase = NULL,
  options = "--likelihood --bits:100 --megabytes:4000 --quiet",
  ...
)

checkMerlin()

`x`	a `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).
`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.
`...`	Further arguments passed on to `merlin()`.

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

1 2	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.

merlin() returns the screen output of MERLIN invisibly.

likelihoodMerlin() returns a single number; the total likelihood using all indicated markers.

checkMerlin() returns TRUE if MERLIN is installed and available on the system path, and FALSE otherwise.

Magnus Dehli Vigeland

http://csg.sph.umich.edu/abecasis/Merlin/

if(checkMerlin()) {

### Trivial example for validation
x = nuclearPed(1)
m1 = marker(x, "1" = 1:2)           # likelihood = 1/2
m2 = marker(x, "1" = 1, "3" = 1:2)  # likelihood = 1/8
x = setMarkers(x, list(m1, m2))

# 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 = list(singleton(1), singleton(2))
y = setMarkers(y, locus = list(alleles = 1:2))
genotype(y[[1]], marker = 1, id = '1') = 1:2
genotype(y[[2]], marker = 1, id = '2') = 1
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)))
}