Skip to content

Commit

Permalink
Fix Issue rqtl#180 re scan1 error if phenotype rownames have names
Browse files Browse the repository at this point in the history
  • Loading branch information
kbroman committed Dec 8, 2020
1 parent ae53a5c commit a888e91
Show file tree
Hide file tree
Showing 11 changed files with 47 additions and 13 deletions.
2 changes: 2 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,8 @@
- Fix 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.


## qtl2 0.22-11 (2020-07-09)

Expand Down
2 changes: 1 addition & 1 deletion R/est_herit.R
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ est_herit <-
# find individuals in common across all arguments
# and drop individuals with missing covariates or missing *all* phenotypes
ind2keep <- get_common_ids(kinshipIDs, addcovar, weights, complete.cases=TRUE)
ind2keep <- get_common_ids(ind2keep, rownames(pheno)[rowSums(is.finite(pheno)) > 0])
ind2keep <- get_common_ids(ind2keep, pheno[rowSums(is.finite(pheno)) > 0,,drop=FALSE])
if(length(ind2keep)<=2) {
if(length(ind2keep)==0)
stop("No individuals in common.")
Expand Down
3 changes: 2 additions & 1 deletion R/kinship_util.R
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,7 @@ check_kinship_onechr <-

# multiply kinship by weights (from and back)
# assuming weights are really square-root weights
#' @importFrom stats setNames
weight_kinship <-
function(kinship, weights=NULL, tol=1e-8)
{
Expand All @@ -132,7 +133,7 @@ weight_kinship <-
}

# line them up
ind2keep <- get_common_ids(rownames(kinship), names(weights))
ind2keep <- get_common_ids(setNames(rownames(kinship), NULL), setNames(names(weights), NULL))
weights <- weights[ind2keep]
kinship <- kinship[ind2keep, ind2keep, drop=FALSE]

Expand Down
7 changes: 4 additions & 3 deletions R/predict_snpgeno.R
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
#' @keywords utilities
#' @seealso [maxmarg()], [viterbi()], [calc_errorlod()]
#' @export
#' @importFrom stats setNames
#'
#' @examples
#' \dontrun{
Expand All @@ -42,15 +43,15 @@ function(cross, geno, cores=1)
# ensure same chromosomes
if(n_chr(cross) != length(geno) ||
any(chr_names(cross) != names(geno))) {
chr <- get_common_ids(chr_names(cross), names(geno))
chr <- get_common_ids(setNames(chr_names(cross), NULL), setNames(names(geno), 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(ind_ids_geno(cross), rownames(geno[[1]]))
ind <- get_common_ids(setNames(ind_ids_geno(cross), NULL), setNames(rownames(geno[[1]]), NULL))
cross <- cross[ind,]
geno <- geno[ind,]
}
Expand All @@ -75,7 +76,7 @@ function(cross, geno, cores=1)
ph2 <- ph[[chr]][,,2]
if(ncol(fg) != ncol(ph1) ||
any(colnames(fg) != colnames(ph1))) {
mar <- get_common_ids(colnames(fg), colnames(ph1))
mar <- get_common_ids(setNames(colnames(fg), NULL), setNames(colnames(ph1), NULL))
fg <- fg[,mar,drop=FALSE]
ph1 <- ph1[,mar,drop=FALSE]
ph2 <- ph2[,mar,drop=FALSE]
Expand Down
2 changes: 1 addition & 1 deletion R/scan1.R
Original file line number Diff line number Diff line change
Expand Up @@ -210,7 +210,7 @@ scan1 <-
# and drop individuals with missing covariates or missing *all* phenotypes
ind2keep <- get_common_ids(genoprobs, addcovar, Xcovar, intcovar,
weights, complete.cases=TRUE)
ind2keep <- get_common_ids(ind2keep, rownames(pheno)[rowSums(is.finite(pheno)) > 0])
ind2keep <- get_common_ids(ind2keep, pheno[rowSums(is.finite(pheno)) > 0,,drop=FALSE])
if(length(ind2keep)<=2) {
if(length(ind2keep)==0)
stop("No individuals in common.")
Expand Down
2 changes: 1 addition & 1 deletion R/scan1_pg.R
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ scan1_pg <-
# and drop individuals with missing covariates or missing *all* phenotypes
ind2keep <- get_common_ids(genoprobs, addcovar, Xcovar, intcovar,
kinshipIDs, weights, complete.cases=TRUE)
ind2keep <- get_common_ids(ind2keep, rownames(pheno)[rowSums(is.finite(pheno)) > 0])
ind2keep <- get_common_ids(ind2keep, pheno[rowSums(is.finite(pheno)) > 0,,drop=FALSE])
if(length(ind2keep)<=2) {
if(length(ind2keep)==0)
stop("No individuals in common.")
Expand Down
2 changes: 1 addition & 1 deletion R/scan1max.R
Original file line number Diff line number Diff line change
Expand Up @@ -147,7 +147,7 @@ scan1max <-
# and drop individuals with missing covariates or missing *all* phenotypes
ind2keep <- get_common_ids(genoprobs, addcovar, Xcovar, intcovar,
weights, complete.cases=TRUE)
ind2keep <- get_common_ids(ind2keep, rownames(pheno)[rowSums(is.finite(pheno)) > 0])
ind2keep <- get_common_ids(ind2keep, pheno[rowSums(is.finite(pheno)) > 0,,drop=FALSE])
if(length(ind2keep)<=2) {
if(length(ind2keep)==0)
stop("No individuals in common.")
Expand Down
2 changes: 1 addition & 1 deletion R/scan1max_pg.R
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ scan1max_pg <-
# and drop individuals with missing covariates or missing *all* phenotypes
ind2keep <- get_common_ids(genoprobs, addcovar, Xcovar, intcovar,
kinshipIDs, weights, complete.cases=TRUE)
ind2keep <- get_common_ids(ind2keep, rownames(pheno)[rowSums(is.finite(pheno)) > 0])
ind2keep <- get_common_ids(ind2keep, pheno[rowSums(is.finite(pheno)) > 0,,drop=FALSE])
if(length(ind2keep)<=2) {
if(length(ind2keep)==0)
stop("No individuals in common.")
Expand Down
2 changes: 1 addition & 1 deletion R/scan1perm.R
Original file line number Diff line number Diff line change
Expand Up @@ -223,7 +223,7 @@ scan1perm <-
# and drop individuals with missing covariates or missing *all* phenotypes
ind2keep <- get_common_ids(genoprobs, addcovar, Xcovar, intcovar,
kinshipIDs, weights, perm_strata, complete.cases=TRUE)
ind2keep <- get_common_ids(ind2keep, rownames(pheno)[rowSums(is.finite(pheno)) > 0])
ind2keep <- get_common_ids(ind2keep, pheno[rowSums(is.finite(pheno)) > 0,,drop=FALSE])
if(length(ind2keep)<=2) {
if(length(ind2keep)==0)
stop("No individuals in common.")
Expand Down
9 changes: 6 additions & 3 deletions R/weights_util.R
Original file line number Diff line number Diff line change
Expand Up @@ -36,18 +36,20 @@ sqrt_weights <-


# multiply a vector by a set of weights
#' @importFrom stats setNames
weight_vector <-
function(vec, weights, tol=1e-12)
{
if(is_null_weights(weights, tol) || is.null(vec)) return(vec)

# align and multiply
id <- get_common_ids(names(vec), names(weights))
id <- get_common_ids(setNames(names(vec), NULL), setNames(names(weights), NULL))
vec[id] * weights[id]
}


# multiply a matrix by a set of weights
#' @importFrom stats setNames
weight_matrix <-
function(mat, weights, tol=1e-12)
{
Expand All @@ -57,11 +59,12 @@ weight_matrix <-
if(!is.matrix(mat)) mat <- as.matrix(mat)

# align and multiply
id <- get_common_ids(rownames(mat), names(weights))
id <- get_common_ids(setNames(rownames(mat), NULL), setNames(names(weights), NULL))
mat[id,,drop=FALSE] * weights[id]
}

# multiply an array by a set of weights
#' @importFrom stats setNames
weight_array <-
function(arr, weights, tol=1e-12)
{
Expand All @@ -72,6 +75,6 @@ weight_array <-
stop("arr should be a 3-dimensional array")

# align and multiply
id <- get_common_ids(rownames(arr), names(weights))
id <- get_common_ids(setNames(rownames(arr), NULL), setNames(names(weights), NULL))
arr[id,,,drop=FALSE] * weights[id]
}
27 changes: 27 additions & 0 deletions tests/testthat/test-scan1.R
Original file line number Diff line number Diff line change
Expand Up @@ -664,3 +664,30 @@ test_that("weights derived from table() are okay", {
expect_equal(out1k, out2k)

})

test_that("scan1 can handle names in the phenotype rownames", {

iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
pr <- calc_genoprob(iron, error_prob=0.002)
pheno <- iron$pheno
out <- scan1(pr, pheno)

names(rownames(pheno)) <- paste0("mouse", rownames(pheno))
expect_equal(scan1(pr, pheno), out)

names(rownames(pheno)) <- rep(NA, nrow(pheno))
expect_equal(scan1(pr, pheno), out)


# the same with kinship matrix
pheno <- iron$pheno
k <- calc_kinship(pr)
out <- scan1(pr, pheno, k) # no error

names(rownames(pheno)) <- paste0("mouse", rownames(pheno))
expect_equal(scan1(pr, pheno, k), out)

names(rownames(pheno)) <- rep(NA, nrow(pheno))
expect_equal(scan1(pr, pheno, k), out)

})

0 comments on commit a888e91

Please sign in to comment.