tests/testthat/test_snpClust.R

In adjclust: Adjacency-Constrained Clustering of a Block-Diagonal Similarity Matrix

context("Consistency of the results of 'snpClust' across various input formats")

# check function kept for possible future usage
check_snpStat_data <- function() {
  sf <- system.file("data/ld.example.RData", package="snpStats")
  expected <- "497fcd532b5c2bcb082a0dad7ca0d44d"
  if (!(tools::md5sum(sf) == expected)) {
      skip("Different version of data('ld.example', package = 'snpStats')")
  }
}

test_that("'snpClust' gives identical results regardless of data input format", {
  skip_if_not_installed("snpStats")
  check_snpStat_data()
  
  data("ld.example", package = "snpStats")
  h <- 100
  ld.ceph <- snpStats::ld(ceph.1mb, depth = h, stats = "R.squared")
  p <- ncol(ceph.1mb)
  nSamples <- nrow(ceph.1mb)
  h <- 100
  ceph.1mb[4,286]@.Data[1,1] <- as.raw(3) ## to avoid NaNs
    
  # case0: Input belongs to class Matrix::dgCMatrix generated by snpStats::ld function
  # should throw error because input is not symmetric
  ld.ceph <- snpStats::ld(ceph.1mb, depth = h, stats = "R.squared")
  expect_false(isSymmetric(ld.ceph))
  expect_warning(expect_error(snpClust(ld.ceph, h = 100)), 
                 "Forcing the LD similarity to be smaller than or equal to 1")
  
  # case1: Input belongs to class Matrix::dsCMatrix generated by snpStats::ld function
  # with 'symmetric=TRUE'
  ## diagonal elements are 0
  ld.ceph <- snpStats::ld(ceph.1mb, depth = h, stats = "R.squared", symmetric = TRUE)
  # ld.ceph <- round(ld.ceph, digits = 10)
  
  expect_identical(unname(diag(ld.ceph)), rep(0, p))
  ld.ceph[ld.ceph > 1] <- 1
  expect_message(snpClust(ld.ceph, h = 100), 
                "Note: forcing the diagonal of the LD similarity matrix to be 1",
                all = FALSE)
  fit1 <- snpClust(ld.ceph, h = 100)

  # LD values less than 0 or larger than 1
  ld1 <- ld.ceph
  ld1[1,2] <- 1.1
  expect_warning(snpClust(ld1, h = 100))
  ld1[1,2] <- -0.1
  ld1[2,1] <- -0.1
  expect_warning(snpClust(ld1, h = 100))
  rm(ld1)
  
  #case2: Input belongs to class snpStats::SnpMatrix
  expect_warning(fit2 <- snpClust(ceph.1mb, h = 100, stats = "R.squared"), 
                 "Forcing the LD similarity to be smaller than or equal to 1")
  expect_equal(fit2$merge, fit1$merge)
  expect_equal(fit2$height, fit1$height)  

  expect_error(snpClust(ceph.1mb, h = ncol(ceph.1mb), stats = "R.squared"), 
                 "h should be strictly less than p")
      
  #case3: Input belongs class base::matrix
  ceph.1mb <- as.matrix(ceph.1mb)
  fit3 <- expect_warning(snpClust(ceph.1mb, h = 100, stats = "R.squared"), 
                         "Forcing the LD similarity to be smaller than or equal to 1")
  expect_equal(fit3$merge, fit1$merge)
  expect_equal(fit3$height, fit1$height)
  
  # increase test coverage
  ceph.1mb_nonames <- as.matrix(ceph.1mb)
  colnames(ceph.1mb_nonames) <- NULL
  rownames(ceph.1mb_nonames) <- NULL
  expect_warning(snpClust(ceph.1mb_nonames, h = 100, stats = "R.squared"), 
                         "Forcing the LD similarity to be smaller than or equal to 1")
  
  #case4: default h
  ld.ceph.2 <- snpStats::ld(ceph.1mb, depth = ncol(ceph.1mb) - 1, stats = "R.squared", symmetric = TRUE)
  fit4 <- suppressWarnings({ snpClust(ld.ceph.2, ncol(ceph.1mb) - 1) })
  fit5 <- suppressWarnings({ snpClust(ld.ceph.2) })
  fit6 <- expect_warning(snpClust(ceph.1mb, stats = "R.squared"), 
                         "Forcing the LD similarity to be smaller than or equal to 1")
  expect_equal(fit4$merge, fit5$merge)
  expect_equal(fit4$height, fit5$height)
  expect_equal(fit4$merge, fit6$merge) ## identical heights but different merges
  expect_equal(fit4$height, fit6$height)
  
  # test that hicClust methods returns expected 'calls'
  expect_identical(as.list(fit1$call)[[1]], as.symbol("snpClust"))
  expect_identical(as.list(fit2$call)[[1]], as.symbol("snpClust"))
  expect_identical(as.list(fit3$call)[[1]], as.symbol("snpClust"))
  expect_identical(as.list(fit4$call)[[1]], as.symbol("snpClust"))
  
})