Merge pull request #188 from kbroman/master

Several small changes
rqtl · Dec 15, 2020 · 896ba76 · 896ba76
2 parents 6db3a65 + 518bef1
commit 896ba76
Show file tree

Hide file tree

Showing 28 changed files with 121 additions and 51 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: qtl2
-Version: 0.23-13
-Date: 2020-12-10
+Version: 0.23-14
+Date: 2020-12-15
 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
@@ -13,6 +13,9 @@
 
 - `fit1()` now returns both fitted values and residuals.
 
+- `fit1()` can be run with genotype probabilities omitted, in which
+  case an intercept column of 1's is used (Issue #151).
+
 - Updated mouse gene database with 2020-09-07 data from
   [MGI](http://www.informatics.jax.org/downloads/mgigff3/archive/monthly/).
 
@@ -21,6 +24,10 @@
 - Make the [vdiffr](https://vdiffr.r-lib.org) package optional: only
   test the plots locally, and only if vdiffr is installed.
 
+- `calc_sdp()` can now take a plain vector (Issue #142).
+
+- Added a `lodcolumn` argument to `maxlod()` (Issue #137).
+
 ### Bug fixes
 
 - Fixed [Issue #181](https://github.com/rqtl/qtl2/issues/181), where

diff --git a/R/calc_sdp.R b/R/calc_sdp.R
@@ -23,9 +23,14 @@ calc_sdp <-
     function(geno)
 {
     # tolerate data frames, but convert to matrix
-    if(!is.matrix(geno) && is.data.frame(geno))
-        geno <- as.matrix(geno)
-    if(!is.matrix(geno)) stop("geno should be a matrix")
+    if(!is.matrix(geno)) {
+        if(is.data.frame(geno)) {
+            geno <- as.matrix(geno)
+        } else {
+            geno <- rbind(geno)
+            dimnames(geno) <- NULL
+        }
+    }
 
     n_str <- ncol(geno)
 

diff --git a/R/fit1.R b/R/fit1.R
@@ -3,7 +3,8 @@
 #' Fit a single-QTL model at a single putative QTL position and get detailed results
 #' about estimated coefficients and individuals contributions to the LOD score.
 #'
-#' @param genoprobs A matrix of genotype probabilities, individuals x genotypes
+#' @param genoprobs A matrix of genotype probabilities, individuals x genotypes.
+#'     If NULL, we create a single intercept column, matching the individual IDs in `pheno`.
 #' @param pheno A numeric vector of phenotype values (just one phenotype, not a matrix of them)
 #' @param kinship Optional kinship matrix.
 #' @param addcovar An optional numeric matrix of additive covariates.
@@ -122,9 +123,14 @@ fit1 <-
              contrasts=NULL, model=c("normal", "binary"),
              zerosum=TRUE, se=TRUE, hsq=NULL, reml=TRUE, blup=FALSE, ...)
 {
-    if(is.null(genoprobs)) stop("genoprobs is NULL")
     if(is.null(pheno)) stop("pheno is NULL")
 
+    if(missing(genoprobs) || is.null(genoprobs)) { # create matrix of 1's
+        if(!is.matrix(pheno)) pheno <- cbind(pheno)
+        genoprobs <- matrix(1, ncol=1, nrow=nrow(pheno))
+        dimnames(genoprobs) <- list(rownames(pheno), "intercept")
+    }
+
     model <- match.arg(model)
 
     if(blup) {

diff --git a/R/max_scan1.R b/R/max_scan1.R
@@ -163,6 +163,9 @@ max.scan1 <-
 #' @param map A list of vectors of marker positions, as produced by
 #' [insert_pseudomarkers()].
 #' @param chr Optional vector of chromosomes to consider.
+#' @param lodcolumn An integer or character string indicating the LOD
+#' score column, either as a numeric index or column name.
+#' If `NULL`, return maximum for all columns.
 #'
 #' @export
 #' @return A single number: the maximum LOD score across all columns and positions for
@@ -193,7 +196,7 @@ max.scan1 <-
 #' # maximum on chromosome 2
 #' maxlod(out, map, "2")
 maxlod <-
-    function(scan1_output, map=NULL, chr=NULL)
+    function(scan1_output, map=NULL, chr=NULL, lodcolumn=NULL)
 {
     if(is.null(scan1_output)) stop("scan1_output is NULL")
 
@@ -206,6 +209,10 @@ maxlod <-
         scan1_output <- subset(scan1_output, map=map, chr=chr)
     }
 
+    if(!is.null(lodcolumn)) {
+        scan1_output <- subset(scan1_output, lodcolumn=lodcolumn)
+    }
+
     # to handle output of either scan1() or scan1coef()
     # for coef(), look at the sign
     if(inherits(scan1_output, "scan1coef")) {

diff --git a/README.md b/README.md
@@ -1,4 +1,4 @@
-### [R/qtl2](https://kbroman.org/qtl2)
+### [R/qtl2](https://kbroman.org/qtl2/)
 
 [![R build status](https://github.com/rqtl/qtl2/workflows/R-CMD-check/badge.svg)](https://github.com/rqtl/qtl2/actions)
 [![CRAN_Status_Badge](https://www.r-pkg.org/badges/version/qtl2)](https://cran.r-project.org/package=qtl2)

diff --git a/man/fit1.Rd b/man/fit1.Rd
diff --git a/man/maxlod.Rd b/man/maxlod.Rd
diff --git a/tests/testthat/test-basic_summaries.R b/tests/testthat/test-basic_summaries.R
@@ -126,7 +126,7 @@ test_that("basic summaries give correct numbers for grav2 data", {
 
 test_that("basic summaries give correct numbers for recla data", {
 
-    skip_if(isnt_karl(), "This test only run locally")
+    skip_if(isnt_karl(), "this test only run locally")
 
     file <- paste0("https://raw.githubusercontent.com/rqtl/",
                    "qtl2data/master/DO_Recla/recla.zip")

diff --git a/tests/testthat/test-calc_het.R b/tests/testthat/test-calc_het.R
@@ -28,7 +28,7 @@ test_that("calc_het works for an intercross", {
 
 test_that("calc_het works with multi-core", {
 
-    skip_if(isnt_karl(), "This test only run locally")
+    skip_if(isnt_karl(), "this test only run locally")
 
     iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
     pr <- calc_genoprob(iron, err=0.002)

diff --git a/tests/testthat/test-calc_raw_summaries.R b/tests/testthat/test-calc_raw_summaries.R
@@ -2,7 +2,7 @@ context("calc raw summaries")
 
 test_that("calculation of raw summaries work", {
 
-    skip_if(isnt_karl(), "This test only run locally")
+    skip_if(isnt_karl(), "this test only run locally")
 
     # load example data
     file <- paste0("https://raw.githubusercontent.com/rqtl/",

diff --git a/tests/testthat/test-calc_sdp.R b/tests/testthat/test-calc_sdp.R
@@ -42,6 +42,10 @@ test_that("calc_sdp works", {
     g <- g[n_AA > 0 & n_AA < 8,]
     expect_equal(calc_sdp(g),
                  apply(g, 1, function(a) sum(((a-1)/2)*2^(seq(along=a)-1))))
+
+    expect_equal( calc_sdp( c(1,1,1,3,1,1,1,1) ), 8)
+    expect_equal( calc_sdp( data.frame(1,1,1,3,1,1,1,1)),  8)
+
 })
 
 

diff --git a/tests/testthat/test-fit1.R b/tests/testthat/test-fit1.R
@@ -105,7 +105,7 @@ test_that("fit1 by H-K works in intercross", {
 
 test_that("fit1 by H-K works in intercross, with weights", {
 
-    skip_if(isnt_karl(), "This test only run locally")
+    skip_if(isnt_karl(), "this test only run locally")
 
     iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
     iron <- iron[,c(18:19,"X")]
@@ -260,7 +260,7 @@ test_that("fit1 by H-K works in riself", {
 
 test_that("fit1 by LMM works in intercross", {
 
-    skip_if(isnt_karl(), "This test only run locally")
+    skip_if(isnt_karl(), "this test only run locally")
 
     iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
     map <- insert_pseudomarkers(iron$gmap, step=1)
@@ -341,7 +341,7 @@ test_that("fit1 by LMM works in intercross", {
 
 test_that("fit1 by LMM works in intercross, with weights", {
 
-    skip_if(isnt_karl(), "This test only run locally")
+    skip_if(isnt_karl(), "this test only run locally")
 
     iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
     map <- insert_pseudomarkers(iron$gmap, step=1)
@@ -514,7 +514,7 @@ test_that("fit1 handles contrasts properly in an intercross", {
 
 test_that("fit1 works with blup=TRUE", {
 
-    skip_if(isnt_karl(), "This test only run locally")
+    skip_if(isnt_karl(), "this test only run locally")
 
     iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
     map <- insert_pseudomarkers(iron$gmap, step=1)
@@ -543,7 +543,7 @@ test_that("fit1 works with blup=TRUE", {
 
 test_that("fit1 fitted values don't depend on order of individuals", {
 
-    skip_if(isnt_karl(), "This test only run locally")
+    skip_if(isnt_karl(), "this test only run locally")
 
     iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
     map <- insert_pseudomarkers(iron$gmap, step=5)
@@ -573,3 +573,34 @@ test_that("fit1 fitted values don't depend on order of individuals", {
 
     testthat::expect_equal(out_fit1$fitted[ind] , out_fit1b$fitted[ind])
 })
+
+test_that("fit1 works without genoprobs", {
+
+    set.seed(20201215)
+    n <- 100
+    nam <- paste0("mouse", sample(10*n, n))
+
+    phe <- setNames(rnorm(n), nam)
+    cov <- setNames(sample(0:1, n, replace=TRUE), nam)
+
+    lm_out <- lm(phe ~ cov)
+    lm_sum <- summary(lm_out)
+
+    coef_names <- c("intercept", "ac1", "intercept")
+    expected <- list(lod=0,
+                     ind_lod=setNames(rep(0, n), nam),
+                     coef=setNames(c(0, lm_out$coef[2], lm_out$coef[1]), coef_names),
+                     SE=setNames(lm_sum$coef[c(1,2,1),"Std. Error"], coef_names),
+                     fitted=lm_out$fitted,
+                     resid=lm_out$resid)
+
+    expect_equal(fit1(pheno=phe, addcovar=cov), expected)
+
+    k <- matrix(0.5,ncol=n, nrow=n)
+    diag(k) <- 1
+    dimnames(k) <- list(nam, nam)
+
+    # just test that this works
+    should_work <- fit1(pheno=phe, addcovar=cov, kinship=k)
+
+})
diff --git a/tests/testthat/test-max_scan1.R b/tests/testthat/test-max_scan1.R
@@ -149,6 +149,11 @@ test_that("maxlod works for intercross with two phenotypes", {
     # overall max
     expect_equal(maxlod(out), max(unclass(out)))
 
+    expect_equal(maxlod(out, lodcolumn="liver"), max(out[,"liver"], na.rm=TRUE))
+    expect_equal(maxlod(out, lodcolumn="spleen"), max(out[,"spleen"], na.rm=TRUE))
+    expect_equal(maxlod(out, lodcolumn=1), max(out[,"liver"], na.rm=TRUE))
+    expect_equal(maxlod(out, lodcolumn=2), max(out[,"spleen"], na.rm=TRUE))
+
     expect_equal(maxlod(out, map, c("2", "9")),
                  max(unclass(subset(out, map, c("2", "9")))))
 

diff --git a/tests/testthat/test-plot_coef.R b/tests/testthat/test-plot_coef.R
@@ -2,7 +2,7 @@ context("plot_coef")
 
 test_that("plot_coef works", {
 
-    skip_if(isnt_karl(), "plot tests only done locally")
+    skip_if(isnt_karl(), "plot tests only run locally")
 
     iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
     iron <- iron[,2]

diff --git a/tests/testthat/test-plot_genes.R b/tests/testthat/test-plot_genes.R
@@ -2,7 +2,7 @@ context("plot_genes")
 
 test_that("plot_genes works", {
 
-    skip_if(isnt_karl(), "plot tests only done locally")
+    skip_if(isnt_karl(), "plot tests only run locally")
 
     genes <- data.frame(chr = c("6", "6", "6", "6", "6", "6", "6", "6"),
                         start = c(139988753, 140680185, 141708118, 142234227, 142587862,

diff --git a/tests/testthat/test-plot_genoprob.R b/tests/testthat/test-plot_genoprob.R
@@ -2,7 +2,7 @@ context("plot_genoprob")
 
 test_that("plot_genoprob works", {
 
-    skip_if(isnt_karl(), "plot tests only done locally")
+    skip_if(isnt_karl(), "plot tests only run locally")
 
     iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
     iron <- iron[c("116", "232"),2]

diff --git a/tests/testthat/test-plot_peaks.R b/tests/testthat/test-plot_peaks.R
@@ -2,7 +2,7 @@ context("plot_peaks")
 
 test_that("plot_peaks works", {
 
-    skip_if(isnt_karl(), "plot tests only done locally")
+    skip_if(isnt_karl(), "plot tests only run locally")
 
     iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
     map <- insert_pseudomarkers(iron$gmap, step=1)

diff --git a/tests/testthat/test-plot_pxg.R b/tests/testthat/test-plot_pxg.R
@@ -2,7 +2,7 @@ context("plot_pxg")
 
 test_that("plot_pxg works", {
 
-    skip_if(isnt_karl(), "plot tests only done locally")
+    skip_if(isnt_karl(), "plot tests only run locally")
 
     iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
     map <- insert_pseudomarkers(iron$gmap, step=1)

diff --git a/tests/testthat/test-plot_scan1.R b/tests/testthat/test-plot_scan1.R
@@ -2,7 +2,7 @@ context("plot_scan1")
 
 test_that("plot_scan1 works", {
 
-    skip_if(isnt_karl(), "plot tests only done locally")
+    skip_if(isnt_karl(), "plot tests only run locally")
 
     iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
     map <- insert_pseudomarkers(iron$gmap, step=1)

diff --git a/tests/testthat/test-recode_snps.R b/tests/testthat/test-recode_snps.R
@@ -2,7 +2,7 @@ context("recode_snps")
 
 test_that("recode_snps works", {
 
-    skip_if(isnt_karl(), "This test only run locally")
+    skip_if(isnt_karl(), "this test only run locally")
 
     # load example data
     file <- paste0("https://raw.githubusercontent.com/rqtl/",

diff --git a/tests/testthat/test-scan1.R b/tests/testthat/test-scan1.R
@@ -205,7 +205,7 @@ test_that("scan1 for backcross with one phenotype", {
 
 test_that("scan1 for backcross with multiple phenotypes with NAs", {
 
-    skip_if(isnt_karl(), "This test only run locally")
+    skip_if(isnt_karl(), "this test only run locally")
 
     set.seed(20151202)
     library(qtl)
@@ -325,7 +325,7 @@ test_that("scan1 for backcross with multiple phenotypes with NAs", {
 
 test_that("scan1 works with NAs in the covariates", {
 
-    skip_if(isnt_karl(), "This test only run locally")
+    skip_if(isnt_karl(), "this test only run locally")
 
     set.seed(20151202)
     library(qtl)
@@ -372,7 +372,7 @@ test_that("scan1 works with NAs in the covariates", {
 
 test_that("scan1 aligns the individuals", {
 
-    skip_if(isnt_karl(), "This test only run locally")
+    skip_if(isnt_karl(), "this test only run locally")
 
     set.seed(20151202)
     library(qtl)
@@ -537,7 +537,7 @@ test_that("multi-core scan1 works", {
 
 test_that("scan1 LOD results don't depend on scale of x and y", {
 
-    skip_if(isnt_karl(), "This test only run locally")
+    skip_if(isnt_karl(), "this test only run locally")
 
     set.seed(20151202)
     library(qtl)
@@ -610,7 +610,7 @@ test_that("scan1 LOD results don't depend on scale of x and y", {
 
 test_that("scan1 deals with mismatching individuals", {
 
-    skip_if(isnt_karl(), "This test only run locally")
+    skip_if(isnt_karl(), "this test only run locally")
 
     iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
     map <- insert_pseudomarkers(iron$gmap, step=2.5)
@@ -633,7 +633,7 @@ test_that("scan1 deals with mismatching individuals", {
 
 test_that("scan1 can handle decomposed kinship matrix", {
 
-    skip_if(isnt_karl(), "This test only run locally")
+    skip_if(isnt_karl(), "this test only run locally")
 
     iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))