tests/testthat/test-merlin.R

In pedprobr: Probability Computations on Pedigrees

if(!checkMerlin())
  skip("Merlin not installed")

lm = function(..., verbose = FALSE)
  likelihoodMerlin(..., verbose = verbose)

expect_signif = function(a, b)
  expect_equal(signif(a, 3), signif(b, 3))

test_that("MERLIN gives correct likelihood in selfing ped linked markers", {
  x = setMarkers(selfingPed(1), alleleMatrix = data.frame(m1 = c("1/2", "1/1"), m2 = c("1/2", "1/1")))
  p = q = 0.5
  r = 0.5
  expect_signif(lm(x), .5*p^2*q^2*(r^2 + (1-r)^2))

  # Complete linkage (requires that the markers are one the same chrom)
  chrom(x) = 1
  r = 0
  expect_signif(lm(x, rho = 0), .5*p^2*q^2*(r^2 + (1-r)^2))

  # Complete  linkage (requires that the markers are one the same chrom)
  r = 0.1
  expect_signif(lm(x, rho = 0.1), .5*p^2*q^2*(r^2 + (1-r)^2))
})


test_that("MERLIN agrees with likelihood2() with two linked markers", {
  x = linearPed(2)
  m = marker(x, geno = c("1/1", NA, "1/2", NA, "1/1"))
  x = setMarkers(x, list(m, m))

  rho = 0.1
  expect_signif(likelihood2(x, 1, 2, rho = rho),
               lm(x, 1:2, rho = rho, verbose = F))

  # With marker names and pre-set chroms
  y = x
  name(y) = c("m1", "m2")
  chrom(y) = 2
  expect_signif(likelihood2(y, 1, 2, rho = rho),
                lm(y, 1:2, rho = rho))

  # With pre-set map
  map = data.frame(CHROM = 1, MARKER = name(y), CM = 1:2)
  expect_signif(likelihood2(y, 1, 2, rho = haldane(cM = 1)),
                lm(y, 1:2, linkageMap = map))

})

test_that("MINX agrees with likelihood2() with two linked markers", {
  x = linearPed(2, sex = 2:1)
  m = marker(x, geno = c("1/1", NA, NA, "1/2", "1/1"), chrom = "X")
  x = setMarkers(x, list(m, m))

  rho = 0.5
  expect_signif(likelihood2(x, 1, 2, rho = rho),
               lm(x, 1:2, rho = rho))

  # With marker names and pre-set chroms
  y = x
  name(y) = c("m1", "m2")
  expect_signif(likelihood2(y, 1, 2, rho = rho),
                lm(y, 1:2, rho = rho))

  # With pre-set map
  map = data.frame(CHROM = "X", MARKER = name(y), CM = 1:2)
  expect_signif(likelihood2(y, 1, 2, rho = haldane(cM = 1)),
                lm(y, 1:2, linkageMap = map))
})

test_that("MINX gives correct likelihood with two markers", {
  m = data.frame(m1 = c("1/1", "1/2", "1/1"), m2 = c("1/1", "1/2", "1/1"), row.names = c(1,4,5))
  x = setMarkers(linearPed(2, sex=2:1), alleleMatrix = m)
  chrom(x) = "X"
  # plot(x, 1:2)
  expect_signif(prod(likelihood(x)), lm(x))
})