merlin {pedprobr} | R Documentation |
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)
Arguments
x |
A |
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 |
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 |
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 |
perChrom |
A logical; if TRUE, likelihoods are reported per chromosome. |
... |
Further arguments passed on to |
program |
A character containing "merlin", "minx" or both (default), optionally including full paths. |
version |
A logical. If TRUE (default), it is checked that running
|
error |
A logical, indicating if an error should be raised if |
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.
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.
Author(s)
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)))
}