From 3f33b7783fed71e071e981e2ec2b929decfbafd6 Mon Sep 17 00:00:00 2001
From: Karl Broman <kbroman@gmail.com>
Date: Wed, 18 Nov 2020 09:55:03 -0600
Subject: [PATCH] Fix Issue #181 (calc_het wasn't working with qtl2fst probs)

- Fixed a similar bug in calc_geno_freq()
- Added related tests in R/qtl2fst
---
 DESCRIPTION        | 4 ++--
 NEWS.md            | 7 ++++++-
 R/calc_geno_freq.R | 8 +++++---
 R/calc_het.R       | 8 ++++----
 4 files changed, 17 insertions(+), 10 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index fd4e9d73..4b9e2963 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: qtl2
-Version: 0.23-6
-Date: 2020-10-22
+Version: 0.23-7
+Date: 2020-11-18
 Title: Quantitative Trait Locus Mapping in Experimental Crosses
 Description: Provides a set of tools to perform quantitative
     trait locus (QTL) analysis in experimental crosses. It is a
diff --git a/NEWS.md b/NEWS.md
index 681e689c..6228c48d 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,4 +1,4 @@
-## qtl2 0.23-6 (2020-10-22)
+## qtl2 0.23-7 (2020-11-18)
 
 ### Major changes
 
@@ -18,6 +18,11 @@
 
 ### Bug fixes
 
+- Fixed [Issue #181](https://github.com/rqtl/qtl2/issues/181), where
+  `calc_het()` gave values > 1 when used with
+  [R/qtl2fst](https://github.com/rqtl/qtl2fst)-based probabilities.
+  Also fixed a similar bug in `calc_geno_freq()`.
+
 - Fixed [Issue #172](https://github.com/rqtl/qtl2/issues/172), where
   `fit1()` gave incorrect fitted values when `kinship` is provided,
   because they weren't "rotated back".
diff --git a/R/calc_geno_freq.R b/R/calc_geno_freq.R
index 12b6058b..a42b86fd 100644
--- a/R/calc_geno_freq.R
+++ b/R/calc_geno_freq.R
@@ -61,12 +61,14 @@ calc_geno_freq <-
 
     # for rest, can assume that they're all one group
 
+    n_chr <- length(probs)
+
     if(by=="individual") {
         # total markers
-        total_mar <- sum( vapply(probs, function(a) dim(a)[3], 1) )
+        total_mar <- sum(dim(probs)[3,])
 
         # summarize each chromosome
-        result <- lapply(probs, apply, 1:2, sum)
+        result <- lapply(seq_len(n_chr), function(chr) apply(probs[[chr]], 1:2, sum))
 
         if(length(result)>1) {
             for(i in seq_along(result)[-1])
@@ -76,5 +78,5 @@ calc_geno_freq <-
     }
 
     # else: by marker
-    t(do.call("cbind", lapply(probs, apply, 2:3, mean)))
+    t(do.call("cbind", lapply(seq_len(n_chr), function(chr) apply(probs[[chr]], 2:3, mean))))
 }
diff --git a/R/calc_het.R b/R/calc_het.R
index d8155839..797a9fc1 100644
--- a/R/calc_het.R
+++ b/R/calc_het.R
@@ -41,17 +41,17 @@ calc_het <-
 
     # determine which columns are het
     het_col <- vector("list", n_chr)
+    geno <- dimnames(probs)[[2]]
     for(chr in seq_len(n_chr)) {
-        geno <- colnames(probs[[chr]])
-        a1 <- substr(geno, 1, 1)
-        a2 <- substr(geno, 2, 2)
+        a1 <- substr(geno[[chr]], 1, 1)
+        a2 <- substr(geno[[chr]], 2, 2)
         if(is_x_chr[chr]) het_col[[chr]] <- (a1 != a2 & a2 != "Y")
         else het_col[[chr]] <- (a1 != a2)
     }
 
     if(by=="individual") {
         # total markers
-        total_mar <- sum( vapply(probs, function(a) dim(a)[3], 1) )
+        total_mar <- sum(dim(probs)[3,])
 
         # summarize each chromosome
         result <- lapply(seq_len(n_chr), function(chr) apply(probs[[chr]][,het_col[[chr]],,drop=FALSE], 1, sum))