tests/testthat/test_2_fitOneMarker.R

context("tests of function fitOneMarker")

# setwd("tests/testthat") # for testing interactively

outdir <- "../../../fitPoly_testoutput"
#output directory is sister directory of package directory

cleanup <- function(outdir, create) {
  while (names(dev.cur()) != "null device") suppressWarnings( dev.off() )
  if (dir.exists(outdir)) unlink(outdir, recursive=TRUE) #remove outdir completely
  if (file.exists(outdir)) file.remove(outdir)
  if (file.exists(outdir)) return(FALSE)
  if (create) {
    dir.create(outdir)
    if (!dir.exists(outdir)) return(FALSE)
  }
  TRUE
}

if (!cleanup(outdir, create=TRUE))
  stop("Start of test2_fitOneMarker: cannot cleanup or create outdir")

data(fitPoly_data)

test_that("4x - multiple populations example Gerrit Gort", {
  #without population structure, ploidy=4:
  loadvars <- load("../testdata/SNP5815.RData")
  expect_equal(loadvars, "SNP5815")
  MarkerName <- rep("SNP5815", nrow(SNP5815))
  SampleName <- character(nrow(SNP5815))
  for (pp in unique(SNP5815$pop)) {
    n <- sum(SNP5815$pop == pp)
    SampleName[SNP5815$pop==pp] <- paste(pp, 1:n, sep="_")
  }
  dat <- data.frame(MarkerName=MarkerName, SampleName=SampleName,
                    ratio=SNP5815$ratio)
  fp <- fitOneMarker(ploidy=4, marker="SNP5815", data=dat)
  expect_equal(nrow(fp$allmodeldata), 24)
  failedmodels <- fp$allmodeldata$m[is.na(fp$allmodeldata$LL)]
  expect_equal(failedmodels, c(5, 13, 16, 43, 45, 48))
  expect_equal(fp$modeldata$LL, 2203.933, tolerance=5)
  #with all plots:
  suppressWarnings(unlink(paste(outdir, "*", sep="/")))
  expect_equal(length(dir(outdir)), 0)
  fp <- fitOneMarker(ploidy=4, marker=1, data=dat,
                     plot="all", plot.dir=outdir)
  expect_equal(length(dir(outdir)), 4)
  #with population structure, ploidy=4:
  pop.parents <- data.frame( #not pre-sorted
    popID = c("P2", "X", "F1", "P1"),
    Par1 = c("", "", "P1", ""),
    Par2 = c("", "", "P2", ""),
    stringsAsFactors = TRUE #must be converted to character in fitOneMarker
  )
  population <- data.frame(
    SampleName = dat$SampleName,
    Population = SNP5815[, 2]
  )
  fp <- fitOneMarker(ploidy=4, marker="SNP5815", data=dat,
                     population=population, pop.parents=pop.parents)
  expect_equal(nrow(fp$allmodeldata), 16)
  failedmodels <- fp$allmodeldata$m[is.na(fp$allmodeldata$LL)]
  expect_equal(failedmodels, c(5, 13, 16))
  expect_equal(fp$modeldata$LL, 2274.178, tolerance=5)
  #with population structure, ploidy=4, and all plots:
  unlink(paste(outdir, "*", sep="/"))
  expect_equal(length(dir(outdir)), 0)
  fp <- fitOneMarker(ploidy=4, marker=1, data=dat,
                     population=population, pop.parents=pop.parents,
                     plot="all", plot.dir=outdir)
  expect_equal(length(dir(outdir)), 5)
  expect_equal(fp$scores$geno[substring(fp$scores$SampleName, 1, 1) == "P"],
               c(1, 1, 3, 3)) #correct parental dosages
})

test_that("6x - multiple populations", {
  #loadvars <- load("../testdata/6x.RData")
  #expect_equal(loadvars, c("dat6x", "pop6x", "pop.par6x",
  #                         "parPriors6x", "sampPriors6x", "startmeans6x"))

  #with population structure, ploidy=6:
  unlink(paste(outdir, "*", sep="/"))
  expect_equal(length(dir(outdir)), 0)
  fp <- fitOneMarker(ploidy=6, marker=1, data=fitPoly_data$ploidy6$dat6x,
                     population=fitPoly_data$ploidy6$pop6x,
                     pop.parents=fitPoly_data$ploidy6$pop.par6x,
                     plot="all", plot.dir=outdir)
  expect_equal(length(dir(outdir)), 13)
  expect_equal(nrow(fp$allmodeldata), 24)
  failedmodels <- fp$allmodeldata$m[is.na(fp$allmodeldata$LL)]
  expect_equal(failedmodels, c(13, 14, 16, 41, 43, 45, 46, 47, 48))
  expect_equal(fp$modeldata$LL, 626.3242, tolerance=5)

  #with population structure and startmeans, ploidy=6:
  unlink(paste(outdir, "*", sep="/"))
  expect_equal(length(dir(outdir)), 0)
  fp <- fitOneMarker(ploidy=6, marker=1, data=fitPoly_data$ploidy6$dat6x,
                     population=fitPoly_data$ploidy6$pop6x,
                     pop.parents=fitPoly_data$ploidy6$pop.par6x,
                     startmeans=fitPoly_data$ploidy6$startmeans6x,
                     plot="fitted", plot.dir=outdir)
})

test_that("3x", {
  #loadvars <- load("../testdata/3x.RData")
  #expect_equal(loadvars, c("dat3x", "pop3x", "pop.par3x"))

  #without population structure, ploidy=3:
  fp <- fitOneMarker(ploidy=3, marker="mrk039", data=fitPoly_data$ploidy3$dat3x)
  expect_equal(nrow(fp$allmodeldata), 16)
  expect_equal(length(fp$allmodeldata), 53)
  failedmodels <- fp$allmodeldata$m[is.na(fp$allmodeldata$LL)]
  expect_equal(failedmodels, c(7, 12, 14, 15, 16))
  expect_equal(fp$modeldata$LL, 132.6224, tolerance=5)

  #with population structure, ploidy=3, no F1's, plot fitted:
  #specify as matrix/integer population nr instead of data.frame/ population names:
  nf <- length(dir(outdir))
  pop.parents <- matrix(c(NA,NA,NA,NA), ncol=2, byrow=TRUE) #2 panels
  rownames(pop.parents) <- c("panel1", "panel2")
  pop <- data.frame(
    SampleName = fitPoly_data$ploidy3$dat3x$SampleName[fitPoly_data$ploidy3$dat3x$MarkerName == "mrk039"],
    Population = c(rep(1, 90), rep(2, 62)))
  fp <- fitOneMarker(ploidy=3, marker="mrk039", data=fitPoly_data$ploidy3$dat3x,
                     population=pop, pop.parents=pop.parents,
                     plot="fitted", plot.dir=outdir)
  expect_equal(length(dir(outdir)) - nf, 1)
  expect_equal(nrow(fp$allmodeldata), 16)
  failedmodels <- fp$allmodeldata$m[is.na(fp$allmodeldata$LL)]
  expect_equal(failedmodels, c(7, 14, 15, 16))
  expect_equal(fp$modeldata$LL, 132.624, tolerance=5)
  expect_equal(fp$modeldata$mu0, 0.401, tolerance=1e-3)
  expect_equal(fp$modeldata$mu1, 0.531, tolerance=1e-3)

  #with population structure including F1,
  #not allowed for odd ploidy:
  pop.parents <- matrix(c(2,3,NA,NA,NA,NA), ncol=2, byrow=TRUE)
  rownames(pop.parents) <- c("F1", "P1", "P2")
  pop <- data.frame(
    SampleName = fitPoly_data$ploidy3$dat3x$SampleName[fitPoly_data$ploidy3$dat3x$MarkerName == "mrk039"],
    Population = c(2, 3, 3, rep(1, 149)))
  expect_error(
    fp <- fitOneMarker(ploidy=3, marker="mrk040", data=fitPoly_data$ploidy3$dat3x,
                       population=pop, pop.parents=pop.parents,
                       plot="fitted", plot.dir=outdir),
    "invalid population structure: F1 populations not allowed with odd ploidy"
  )
})

test_that("2x", {
  #loadvars <- load("../testdata/2x.RData")
  #expect_equal(loadvars, "dat2x")

  #without population structure, ploidy=2:
  fp <- fitOneMarker(ploidy=2, marker=1, data=fitPoly_data$ploidy2$dat2x)
  expect_equal(nrow(fp$allmodeldata), 16)
  expect_equal(length(fp$allmodeldata), 45)
  failedmodels <- fp$allmodeldata$m[is.na(fp$allmodeldata$LL)]
  expect_equal(failedmodels, c(3, 4, 7, 8, 11, 12, 15, 16))
  #            all models with nonlinear mu's fail
  expect_equal(fp$modeldata$LL, 213.3432, tolerance=5)
})

test_that("4x - with parental priors", {
  #loadvars <- load("../testdata/4x.RData")
  #expect_equal(loadvars, c("dat4x", "pop4x", "pop.par4x",
  #                         "parPriors4x", "sampPriors4x"))

  #without population structure, ploidy=4:
  unlink(paste(outdir, "*", sep="/"))
  expect_equal(length(dir(outdir)), 0)
  fp <- fitOneMarker(ploidy=4, marker=2, data=fitPoly_data$ploidy4$dat4x,
                     plot="fitted", plot.dir=outdir)
  expect_equal(nrow(fp$allmodeldata), 16)
  failedmodels <- fp$allmodeldata$m[is.na(fp$allmodeldata$LL)]
  expect_equal(failedmodels, c(10, 15, 16))
  expect_equal(fp$modeldata$LL, 467.1877, tolerance=5)

  #with population structure, no parental priors:
  fp <- fitOneMarker(ploidy=4, marker=2, data=fitPoly_data$ploidy4$dat4x,
                     population=fitPoly_data$ploidy4$pop4x,
                     pop.parents=fitPoly_data$ploidy4$pop.par4x,
                     plot="fitted", plot.dir=outdir)
  expect_equal(nrow(fp$allmodeldata), 16)
  failedmodels <- fp$allmodeldata$m[is.na(fp$allmodeldata$LL)]
  expect_equal(failedmodels, c(15, 16))
  expect_equal(fp$modeldata$LL, 590.1281, tolerance=5)

  #with population structure and parental and sample priors:
  unlink(paste(outdir, "*", sep="/"))
  expect_equal(length(dir(outdir)), 0)
  fp <- fitOneMarker(ploidy=4, marker=2, data=fitPoly_data$ploidy4$dat4x,
                     population=fitPoly_data$ploidy4$pop4x,
                     pop.parents=fitPoly_data$ploidy4$pop.par4x,
                     parentalPriors=fitPoly_data$ploidy4$parPriors4x,
                     samplePriors=fitPoly_data$ploidy4$sampPriors4x,
                     plot="all", plot.dir=outdir)
  expect_equal(length(dir(outdir)), 13)
  expect_equal(nrow(fp$allmodeldata), 24)
  failedmodels <- fp$allmodeldata$m[is.na(fp$allmodeldata$LL)]
  expect_equal(failedmodels, c(15, 16))
  expect_equal(fp$modeldata$LL, 590.1281, tolerance=5)
})

if (!cleanup(outdir, create=FALSE))
  stop("End of test2_fitOneMarker: cannot cleanup or remove outdir")

Try the fitPoly package in your browser

Any scripts or data that you put into this service are public.

fitPoly documentation built on May 2, 2019, 6:13 a.m.