tests/testthat/test_pac.R

In cgmisc-team/cgmisc: Various functions useful in computational genetics, e.g. in GWAS

library(cgmisc)
library(GenABEL)
#library(visualTest)
library(RSQLite)
library(DBI)

#setwd('tests/testthat/')
data.qc1 <- readRDS('data.qc1.rds')
pop <- readRDS('pop.rds')
mm <- readRDS('mm.rds')

test_that("pop allele counts work", {
  pac <- pop.allele.counts(data = data.qc1[ ,data.qc1@gtdata@chromosome == 2], pops = pop, progress=F)
  expect_equal_to_reference(pac, 'pac.rds')
  
  # Test that the image generated by plot.pac is as it should be
  #tmp <- tempfile(fileext = ".png")
  #png(filename = tmp)
  #plot.pac(data = data.qc1[ ,data.qc1@gtdata@chromosome == 2], allele.cnt = pac, plot.LD = T, legend.pos='topleft')
  #dev.off()
  #expect_true(isSimilar(file = tmp, fingerprint = "E9379273B462B1C1"))
})

test_that("F_ST computations work", {
  fst <- compute.fstats(data = data.qc1[ ,data.qc1@gtdata@chromosome == 2], pops = pop)
  expect_equal_to_reference(fst, 'fst.rds')
  pop2 <- pop
  pop2[2] <- 3
  expect_error(compute.fstats(data = data.qc1[ ,data.qc1@gtdata@chromosome == 2], pops = pop2))
})

test_that("get adjacent markers works", {
  adjacent <- get.adjacent.markers(data = data.qc1,
                                 marker = 'BICF2S2365880',
                                 size.bp = 1e4)
  expect_equal(dim(adjacent), c(205, 2))
  expect_equal(colnames(adjacent), c("BICF2P425207", "BICF2S2365880"))
  expect_equal(adjacent[1,1], 1)
  expect_equal(adjacent['dog234', 'BICF2S2365880'], 0)
  expect_equal(adjacent['dog235', 'BICF2P425207'], 2)
})

#test_that("check LD plot", {
#  tmp <- tempfile(fileext = ".png")
#  png(filename = tmp)
#  plot.manhattan.LD(data = data.qc1, gwas.result = mm, chr = 2, region = c(37256927, 39256927), index.snp = 'BICF2S2365880', legend.pos = 'topleft')
#  dev.off()
#  expect_true(isSimilar(file = tmp, "EB352A3471E1E2D2"))
#})

test_that('selection of top snps', {
  top_snps <- choose.top.snps(data = data.qc1, chr = 2, region = c(37256927, 39256927), index.snp = 'BICF2S2365880')
  expect_equal(dim(top_snps), c(26,2))
  #expect_equal(top_snps[1, 2], 'INDEX SNP')
  expect_equal(top_snps[1, 'coord'], 38256927)
})

test_that('clumping works', {
  clumps <- clump.markers(data = data, gwas.result = mm,
                          chr = 2)
  expect_equal_to_reference(clumps, 'clumps.rds')
})

test_that('get overlapping windows', {
  my.LW <- get.overlap.windows(data = data.qc1, chr = 2, size = 125e3, overlap = 2500)
  het.windows <- het.overlap.wind(data = data.qc1, LW = my.LW, progress = F)
  roh1 <- get.roh(data = data.qc1, chr = 2,
          LW = my.LW,
          hetero.zyg = het.windows, threshold = 0.30,
          strict = TRUE)
  roh2 <- get.roh(data = data.qc1, chr = 2,
                  LW = my.LW,
                  hetero.zyg = het.windows, threshold = 0.25,
                  strict = FALSE)
  expect_equal_to_reference(my.LW, 'LW.windows.rds')
  expect_equal(het.windows[c(349, 350), 2], c(0.3915125, NA), tolerance=1e-3)
  expect_equal(dim(het.windows), c(694, 2))
  expect_equal_to_reference(roh1, 'roh1.rds')
  expect_equal_to_reference(roh2, 'roh2.rds')
})

test_that('getting chromosome midpoints', {
  midpoints <- get.chr.midpoints(data.qc1)
  expect_equal(midpoints[c(1, 27)], c(3255, 100594))
  data.tmp <- data.qc1
  data.tmp@gtdata@chromosome <- as.factor(NULL)
  expect_error(get.chr.midpoints(data.tmp))
})

test_that('get.erv test', {
  ervs <- get.erv(chr = 'chr2', coords = c(10e6, 40e6))
  expect_equal(ervs[ervs$id == 2738, 'length'], 14763)
})

# test_that('plot.qq works', {
#   autosomal <- which(data.qc1@gtdata@chromosome != 39) # 39 is the canine X chromosome
#   data.qc1.gkin <- ibs(data.qc1, snpsubset = autosomal, weight = 'freq')
#   h2h <- polygenic(formula = ct ~ sex, kinship.matrix = data.qc1.gkin,
#                  data = data.qc1, llfun = 'polylik')
#   tmp <- tempfile(fileext = ".png")
#   png(filename = tmp)
#   plot.qq(data = data.qc1, obs = mm[,'P1df'],
#         emp = h2h,
#         N = 10,
#         plot.emp = T, step = 100)
#   dev.off()
#   qq_fp <- visualTest::getFingerprint(file = 'plot_qq.png')
#   testthat::expect_true(isSimilar(file = tmp, qq_fp))
# })