only calc ihs for focal snp if given

statgenlmu · Jun 3, 2015 · 19ef074 · 19ef074
1 parent 43de9ff
commit 19ef074
Show file tree

Hide file tree

Showing 4 changed files with 24 additions and 10 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: coala
-Version: 0.0.9104
-Date: 2015-05-20
+Version: 0.0.9105
+Date: 2015-06-03
 License: GPL (>= 3)
 Title: A Framework for Coalescent Simulation
 Authors@R: c(

diff --git a/R/sumstat_ihh.R b/R/sumstat_ihh.R
@@ -35,15 +35,26 @@ SumstatIhh <- R6Class('sumstat_ihh', inherit = Sumstat, #nolint
           return(private$empty_matrix)
         }
         snps <- private$get_snp(pos[[locus]], locus, model)
-        assert_that(length(snps) > 0)
         haplohh <- self$segsites_to_rehh_data(seg_sites[[locus]],
                                               pos[[locus]],
                                               ind)
-        rehh::scan_hh(haplohh)[snps , -(1:3), drop=FALSE] #nolint
+        if (!is.na(private$position)) {
+          assert_that(length(snps) == 1)
+          ihh <- matrix(0, 1, 3)
+          colnames(ihh) <- c("IHHa", "IHHd", "IES")
+          ihh[1, 1:2] <- rehh::calc_ehh(haplohh, mrk = snps, plotehh = FALSE)$ihh
+          ihh[1, 3] <- rehh::calc_ehhs(haplohh, mrk = snps, plotehh = FALSE)$ies
+          return(ihh)
+        }
+        rehh::scan_hh(haplohh)[snps , -(1:3), drop = FALSE] #nolint
       })
     },
     segsites_to_snp_map = function(seg_sites, pos) {
-      map <- data.frame(name=seq(along = pos), chr=1, pos = pos, anc=0, der=1)
+      map <- data.frame(name = seq(along = pos),
+                        chr = 1,
+                        pos = pos,
+                        anc = 0,
+                        der = 1)
       file <- tempfile('snp_map')
       write.table(map, file, row.names = FALSE, col.names = FALSE)
       file
@@ -57,9 +68,11 @@ SumstatIhh <- R6Class('sumstat_ihh', inherit = Sumstat, #nolint
     segsites_to_rehh_data = function(seg_sites, pos, ind) {
       haplo <- self$segsites_to_haplo(seg_sites, ind)
       snp_map <- self$segsites_to_snp_map(seg_sites, pos)
-      capture.output(rehh <- rehh::data2haplohh(haplo,
-                                                snp_map,
-                                                recode.allele = TRUE))
+      suppressWarnings(
+        capture.output(rehh <- rehh::data2haplohh(haplo,
+                                                  snp_map,
+                                                  recode.allele = TRUE))
+      )
       unlink(c(snp_map, haplo))
       rehh
     }

diff --git a/tests/testthat/test-sumstat-iHH.R b/tests/testthat/test-sumstat-iHH.R
@@ -44,7 +44,8 @@ test_that('calculation of ihh works', {
   expect_equal(length(ihh2), 1)
   expect_that(ihh2[[1]], is_a('matrix'))
   expect_equal(dim(ihh2[[1]]), c(1, 3))
-  expect_equal(ihh[[1]][3, , drop=FALSE], ihh2[[1]])
+  expect_equivalent(ihh[[1]][3, , drop = FALSE], ihh2[[1]])
+  expect_equal(rownames(ihh), rownames(ihh2))
 
   model <- coal_model(4, 3, 337)
   ihh2 <- stat_ihh$calculate(list(seg_sites, seg_sites, seg_sites), NULL, model)

diff --git a/tests/testthat/test-sumstat-nSL.R b/tests/testthat/test-sumstat-nSL.R
@@ -23,7 +23,7 @@ test_that('calculation of nSL works', {
   expect_equal(length(nsl2), 1)
   expect_that(nsl2[[1]], is_a('numeric'))
   expect_equal(length(nsl2[[1]]), 1)
-  expect_equal(nsl[[1]][3], nsl2[[1]])
+  expect_equivalent(nsl[[1]][3], nsl2[[1]])
 
   model <- coal_model(4, 3, 337)
   nsl <- stat_nsl$calculate(list(seg_sites, seg_sites, seg_sites), NULL, model)