tests/testthat/test-in_silico_population.R

In oliviaAB/sismonr: Simulation of in Silico Multi-Omic Networks

context("Creating in silico populations")
library(sismonr)
library(XRJulia)

test_that("creation of variants works",{
  check_julia()
  indargs = insilicoindividualargs()
  genes = createGenes(insilicosystemargs(G = 5, PC.p = 0))
  variants = createVariants(genes, indargs)

  expect_equal(length(variants), nrow(genes))
  expect_equal(ncol(variants[[1]]), indargs$ngenevariants)
  expect_true(all(variants[[1]][c("qtlTLrate", "qtlPDrate", "qtlTLregbind", "qtlPDregrate", "qtlPTMregrate"),] == 0))
})

test_that("creation of in silico individuals works", {
  check_julia()
  mysystem = createInSilicoSystem(G = 3, ploidy = 4)
  indargs = insilicoindividualargs(ngenevariants = 2)
  genvariants = createVariants(mysystem$genes, indargs)
  genvariants.freq = list(`1` = c(0, 1),
                          `2` = rep(0.5, 2),
                          `3` = rep(0.5, 2))
  myind = createIndividual(mysystem, genvariants, genvariants.freq, indargs)
  qtlvalues = sapply(rownames(genvariants[[1]]), function(x){myind$QTLeffects$GCN1[[x]][1]})

  expect_equal(dim(myind$haplotype), c(3, 4))
  expect_false(any(myind$haplotype[1, ] == 1))
  expect_equal(qtlvalues, genvariants[[1]][, 2])
  expect_is(myind, "insilicoindividual")
})

test_that("creation of in silico individuals works - Initial abundance", {
  check_julia()
  mysystem = createInSilicoSystem(G = 3, ploidy = 2)
  indargs = insilicoindividualargs(ngenevariants = 1)
  genvariants = createVariants(mysystem$genes, indargs)
  genvariants.freq = list(`1` = c(1),
                          `2` = rep(1),
                          `3` = rep(1))
  myind1 = createIndividual(mysystem, genvariants, genvariants.freq, indargs, initialNoise = F)
  myind2 = createIndividual(mysystem, genvariants, genvariants.freq, indargs, initialNoise = F)

  myInitVar = list("R" = c(1, 1, 0), "P" = c(1, 2, 1))
  myind3 = createIndividual(mysystem, genvariants, genvariants.freq, indargs, initialNoise = F, InitVar = myInitVar)

  myind4 = createIndividual(mysystem, genvariants, genvariants.freq, indargs, initialNoise = T)

  expect_equal(myind1$InitAbundance$GCN1$R, myind1$InitAbundance$GCN2$R)
  expect_equal(myind1$InitAbundance$GCN1$P, myind1$InitAbundance$GCN2$P)

  expect_equal(myind1$InitAbundance$GCN1$R, myind2$InitAbundance$GCN1$R)
  expect_equal(myind1$InitAbundance$GCN1$P, myind2$InitAbundance$GCN1$P)

  expect_equal(myind3$InitAbundance$GCN1$R[3], 0)
  expect_equal(myind3$InitAbundance$GCN2$R[3], 0)

  expect_equal(abs(myind3$InitAbundance$GCN1$P[2] - 2 * myind1$InitAbundance$GCN1$P[2]) <= 1, TRUE)
  expect_equal(abs(myind3$InitAbundance$GCN2$P[2] - 2 * myind1$InitAbundance$GCN2$P[2]) <= 1, TRUE)

  ## expect_false(all(myind4$InitAbundance$GCN1$R == myind4$InitAbundance$GCN2$R)) removed just in case
  ##expect_false(all(myind4$InitAbundance$GCN1$P == myind4$InitAbundance$GCN2$P))

  ##expect_false(all(myind4$InitAbundance$GCN1$R == myind1$InitAbundance$GCN1$R))
  ##expect_false(all(myind4$InitAbundance$GCN1$P == myind1$InitAbundance$GCN1$P))
})

test_that("creation of in silico population works", {
  check_julia()
  mysystem = createInSilicoSystem(G = 5, ploidy = 4)
  mypop = createInSilicoPopulation(nInd = 3, mysystem, ngenevariants = 2)

  expect_equal(length(mypop$individualsList), 3)
  expect_equal(length(mypop$GenesVariants), 5)
  expect_equal(ncol(mypop$GenesVariants[[1]]), 2)
  expect_equal(dim(mypop$individualsList$Ind1$haplotype), c(5, 4))
})

#removeJuliaEvaluator(getJuliaEvaluator())