merlin: Pedigree likelihood computed by MERLIN

Description Usage Arguments Details Value Author(s) References Examples

View source: R/merlin.R

Description

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.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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()

Arguments

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

Details

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.

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 MERLIN is installed and available on the system path, and FALSE otherwise.

Author(s)

Magnus Dehli Vigeland

References

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

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
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)))
}

pedprobr documentation built on Nov. 13, 2020, 5:06 p.m.