Merge pull request #187 from kbroman/master

Fix Issue #146 so predict_snpgeno() works with RIL

Author: kbroman

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: qtl2
-Version: 0.23-11
-Date: 2020-12-08
+Version: 0.23-12
+Date: 2020-12-09
 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

NEWS.md

@@ -1,4 +1,4 @@
-## qtl2 0.23-11 (2020-12-08)
+## qtl2 0.23-12 (2020-12-09)
 
 ### Major changes
 
@@ -40,19 +40,22 @@
   `genoprob_to_alleleprob()` for general AIL crosses. We had not
   implemented the `geno2allele_matrix()` function.
 
-- Fix Issue #164, so `plot_pxg()` can handle a phenotype that is a
+- Fixed Issue #164, so `plot_pxg()` can handle a phenotype that is a
   single-column data frame.
 
-- Fix Issue #135, so `plot_scan1()` can take vector input (which is
+- Fixed Issue #135, so `plot_scan1()` can take vector input (which is
   then converted to a single-column matrix).
 
-- Fix Issue #157, to have `calc_genoprob()` give a better error
+- Fixed Issue #157, to have `calc_genoprob()` give a better error
   message about missing genetic map.
 
-- Fix Issue #178, to have `read_cross2()` give a warning not an error
+- Fixed Issue #178, to have `read_cross2()` give a warning not an error
   if incorrect number of alleles.
 
-- Fix Issue #180 re `scan1()` error if phenotypes' rownames have rownames.
+- Fixed Issue #180 re `scan1()` error if phenotypes' rownames have rownames.
+
+- Fixed Issue #146, revising `predict_snpgeno()` so that it works for
+  homozygous populations, like MAGIC lines or the Collaborative Cross.
 
 
 ## qtl2 0.22-11 (2020-07-09)

R/predict_snpgeno.R

@@ -43,21 +43,27 @@ function(cross, geno, cores=1)
    # ensure same chromosomes
    if(n_chr(cross) != length(geno) ||
       any(chr_names(cross) != names(geno))) {
-       chr <- get_common_ids(setNames(chr_names(cross), NULL), setNames(names(geno), NULL))
+       chr <- get_common_ids(setNames(names(geno), NULL), setNames(chr_names(cross), NULL))
        cross <- cross[,chr]
        geno <- geno[,chr]
    }
 
    # ensure same individuals
    if(n_ind_geno(cross) != nrow(geno[[1]]) ||
       any(ind_ids_geno(cross) != rownames(geno[[1]]))) {
-       ind <- get_common_ids(setNames(ind_ids_geno(cross), NULL), setNames(rownames(geno[[1]]), NULL))
+       ind <- get_common_ids(setNames(rownames(geno[[1]]), NULL), setNames(ind_ids_geno(cross), NULL))
        cross <- cross[ind,]
        geno <- geno[ind,]
    }
 
+   phase_known <- is_phase_known(cross)
+
    # guess phase
-   ph <- qtl2::guess_phase(cross, geno, cores=cores)
+   if(!phase_known) {
+       ph <- guess_phase(cross, geno, cores=cores)
+   } else {
+       ph <- cross$geno
+   }
 
    # subset cross chromosomes
    cross <- cross[,names(ph)]
@@ -72,8 +78,12 @@ function(cross, geno, cores=1)
    by_chr_func <- function(chr) {
        # ensure the same markers
        fg <- cross$founder_geno[[chr]]
-       ph1 <- ph[[chr]][,,1]
-       ph2 <- ph[[chr]][,,2]
+       if(phase_known) {
+           ph1 <- ph2 <- geno[[chr]]
+       } else {
+           ph1 <- ph[[chr]][,,1]
+           ph2 <- ph[[chr]][,,2]
+       }
        if(ncol(fg) != ncol(ph1) ||
           any(colnames(fg) != colnames(ph1))) {
            mar <- get_common_ids(setNames(colnames(fg), NULL), setNames(colnames(ph1), NULL))

tests/testthat/test-predict_snpgeno.R

@@ -22,3 +22,36 @@ test_that("predict_snpgeno works", {
     expect_equivalent(infg[[1]][11,51:60], c(NA,NA,NA,NA,1,3,1,2,2,3))
 
 })
+
+
+test_that("predict_snpgeno works for magic lines", {
+
+    skip_if(isnt_karl(), "this test only run locally")
+
+    # load example data and calculate genotype probabilities
+    file <- paste0("https://raw.githubusercontent.com/rqtl/",
+                   "qtl2data/master/ArabMAGIC/arabmagic_tair9.zip")
+    magic <- read_cross2(file)
+    ind <- paste0("MAGIC", ".", 1:20)
+    chr <- "2"
+    probs <- calc_genoprob(magic[ind,chr], error_prob=0.002)
+    m <- maxmarg(probs)
+
+    infg <- predict_snpgeno(magic, m)
+
+    expect_equal(class(infg), "list")
+    expect_equal(length(infg), 1)
+    expect_equal(names(infg), chr)
+    expect_equal(dim(infg[[1]]), c(20,211))
+    expect_equivalent(infg[[1]][,1], c(3,1,3,1,NA,3,1,3,1,1,3,3,NA,1,3,NA,1,1,NA,1))
+    expect_equivalent(infg[[1]][9,1:10], c(1,1,1,3,3,1,1,3,3,1))
+    expect_equivalent(infg[[1]][10,101:110],  c(rep(1,9), 3))
+    expect_equivalent(infg[[1]][11,51:60], c(1,1,3,1,1,1,1,1,1,1))
+
+    fg <- magic$founder_geno[[chr]]
+    fg[fg==0] <- NA
+    infg_hard <- t(sapply(1:20, function(wh_ind) sapply(seq_along(m[[chr]][wh_ind,]), function(i) {
+        g <- m[[chr]][wh_ind,i]; ifelse(is.na(g), NA, fg[g,i]) })))
+    expect_equivalent(infg[[chr]], infg_hard)
+
+})