tests/testthat/test_ihh2ihs.R

In rehh: Searching for Footprints of Selection using 'Extended Haplotype Homozygosity' Based Tests

#compares ihh2ihs (standardization) with values obtained by this function in release 2.0.4
context("ihh2ihs")

test_that("checked_ihh2ihs", {
  library(rehh.data)
  data(wgscan.cgu)
  ihs.cgu <- ihh2ihs(wgscan.cgu, verbose = FALSE)
  expectedFreqClasses <- data.frame(
    c(
      1090,
      907,
      787,
      837,
      824,
      886,
      802,
      690,
      843,
      824,
      991,
      823,
      732,
      874,
      917,
      1038,
      966,
      816,
      954,
      1005,
      1160,
      970,
      1036,
      1049,
      897,
      1346,
      1129,
      1177,
      1165,
      1018,
      1514,
      1260,
      1411,
      1478,
      1316,
      1635
    ),
    c(
      -0.711792027907377,
      -0.578893855926013,
      -0.457292623663963,
      -0.399226543658018,
      -0.348478625063976,
      -0.280816726479214,
      -0.282833083733207,
      -0.233570127082276,
      -0.194029441167255,
      -0.186002451166171,
      -0.166916033629341,
      -0.126390301747902,
      -0.107125375206074,
      -0.116176514124248,
      -0.116613897618552,
      -0.0811485113692771,
      -0.0552501596504867,
      -0.000105576812395406,
      -0.0174117431748217,
      0.0149188721038136,
      0.0226999538837744,
      0.0464374061365272,
      0.0787090341193259,
      0.0960085904552361,
      0.108593745123183,
      0.139602860388728,
      0.157519600587737,
      0.195818940728144,
      0.247597642322602,
      0.26801990729495,
      0.265330940334051,
      0.330343664507751,
      0.374366270720777,
      0.471050387171474,
      0.580475986263097,
      0.728608735706007
    ),
    c(
      0.610603227939361,
      0.570357587947288,
      0.487584309883171,
      0.454222982105294,
      0.434112576353058,
      0.462122118497751,
      0.447777662447761,
      0.459491083363563,
      0.441661429902016,
      0.436077200320733,
      0.430416082094664,
      0.442961109519092,
      0.427420907718221,
      0.459441463605084,
      0.46443573120339,
      0.446264010052274,
      0.441710927275144,
      0.456965260306721,
      0.47146650652992,
      0.4533516196045,
      0.430383749145416,
      0.450423201525377,
      0.445725753265559,
      0.462956924562142,
      0.419116376287293,
      0.42152879917028,
      0.453165153945328,
      0.453180504788825,
      0.44596536054357,
      0.43298435214188,
      0.450116704189876,
      0.453250276679715,
      0.472483172300497,
      0.507939155843418,
      0.555679771600724,
      0.645774173625061
    )
  )
  expect_equivalent(ihs.cgu$frequency.class[, 1:3], expectedFreqClasses)
  expectedIhs <- data.frame(
    CHR = rep("1", 6),
    POSITION = c(113642,
                 244699, 369419, 447278, 487654, 524507),
    IHS = c(
      -0.558299239642611,
      0.272333661628506,
      0.481073631119111,
      1.06187102387775,
      0.818405971161885,
      -0.389702403570878
    ),
    LOGPVALUE = c(
      0.239095186782283,
      0.104928194005625,
      0.200339591023662,
      0.540164032852241,
      0.383918088706178,
      0.156918896931616
    ),
    row.names = c(
      "F0100190",
      "F0100220",
      "F0100250",
      "F0100270",
      "F0100280",
      "F0100290"
    ),
    check.names = FALSE,
    stringsAsFactors = FALSE
  )
  expect_equivalent(head(ihs.cgu$ihs), expectedIhs)
  
  # check Bonferoni correction
  ihs.cgu_adj <-
    ihh2ihs(wgscan.cgu, p.adjust.method = "bonferroni", verbose = FALSE)
  
  # multiply raw p-values by number of rows, cut-off at 1
  p.adj <-
    mapply(max, ihs.cgu$ihs$LOGPVALUE - log10(nrow(ihs.cgu$ihs)), 0)
  
  # identical only up to 1E-15
  expect_equivalent(ihs.cgu_adj$ihs$LOGPVALUE, p.adj)
})