Merge pull request #115 from kbroman/master

Add create_snpinfo() for creating snp information table from cross

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: qtl2
-Version: 0.19-11
+Version: 0.19-12
 Date: 2019-05-21
 Title: Quantitative Trait Locus Mapping in Experimental Crosses
 Description: R/qtl2 provides a set of tools to perform quantitative

NAMESPACE

@@ -66,6 +66,7 @@ export(convert2cross2)
 export(count_xo)
 export(covar_names)
 export(create_gene_query_func)
+export(create_snpinfo)
 export(create_variant_query_func)
 export(decomp_kinship)
 export(drop_markers)

NEWS.md

@@ -1,4 +1,4 @@
-## qtl2 0.19-11 (2019-05-21)
+## qtl2 0.19-12 (2019-05-21)
 
 ### Major changes
 
@@ -8,6 +8,9 @@
   `scan1blup()`). The previous behavior can be obtained with the
   argument `zerosum=FALSE`.
 
+- Add function `create_snpinfo()` for creating a SNP information table
+  from a cross2 object, for use with `scan1snps()`.
+
 ### Minor changes
 
 - Updated `extdata/mouse_genes_small.sqlite` using updated MGI

R/create_snpinfo.R

@@ -0,0 +1,76 @@
+#' Create snp information table for a cross
+#'
+#' Create a table of snp information from a cross, for use with [scan1snps()].
+#'
+#' @param cross Object of class `"cross2"`. For details, see the
+#' [R/qtl2 developer guide](https://kbroman.org/qtl2/assets/vignettes/developer_guide.html).
+#'
+#' @return A data frame of SNP information with the following columns:
+#' * `chr` - Character string or factor with chromosome
+#' * `pos` - Position (in same units as in the `"map"`
+#'     attribute in `genoprobs`.
+#' * `snp` - Character string with SNP identifier (if
+#'     missing, the rownames are used).
+#' * `sdp` - Strain distribution pattern: an integer, between
+#'     1 and \eqn{2^n - 2} where \eqn{n} is the number of strains, whose
+#'     binary encoding indicates the founder genotypes
+#' SNPs with missing founder genotypes are omitted.
+#'
+#' @examples
+#' \dontrun{
+#' # load example data and calculate genotype probabilities
+#' file <- paste0("https://raw.githubusercontent.com/rqtl/",
+#'                "qtl2data/master/DO_Recla/recla.zip")
+#' recla <- read_cross2(file)
+#' snpinfo <- create_snpinfo(recla)
+#'
+#' # calculate genotype probabilities
+#' pr <- calc_genoprob(recla, error_prob=0.002, map_function="c-f")
+#'
+#' # index the snp information
+#' snpinfo <- index_snps(recla$pmap, snpinfo)
+#'
+#' # sex covariate
+#' sex <- setNames((recla$covar$Sex=="female")*1, rownames(recla$covar))
+#'
+#' # perform a SNP scan
+#' out <- scan1snps(pr, recla$pmap, recla$pheno[,"bw"], addcovar=sex, snpinfo=snpinfo)
+#'
+#' # plot the LOD scores
+#' plot(out$lod, snpinfo, altcol="green3")
+#' }
+#'
+#' @seealso [index_snps()], [scan1snps()], [genoprob_to_snpprob()]
+#' @export
+create_snpinfo <-
+    function(cross)
+{
+    if(!("cross2" %in% class(cross))) {
+        stop("Input should be a cross2 object")
+    }
+
+    if("pmap" %in% names(cross)) map <- cross$pmap
+    else if("gmap" %in% names(cross)) map <- cross$gmap
+    else stop("cross contains neither pmap nor gmap")
+
+    # convert map to a data frame
+    nmar <- vapply(map, length, 1) # no. markers per chromosome
+    markers <- unlist(lapply(map, names))
+    map <- data.frame(chr=rep(names(map), nmar),
+                      pos=unlist(map),
+                      snp=markers,
+                      stringsAsFactors=FALSE)
+    rownames(map) <- map$snp
+
+    # founder genotypes -> SDP
+    if(!("founder_geno" %in% names(cross))) {
+        stop("cross does not contain founder_geno")
+    }
+    fg <- do.call("cbind", cross$founder_geno)
+
+    mar2drop <- (colSums(is.na(fg) | fg==0) > 0)
+
+    # drop markers with missing data; calculate founder strain distribution patterns (SDPs)
+    cbind(map[!mar2drop,,drop=FALSE],
+          sdp=calc_sdp(t(fg[,!mar2drop])))
+}

R/plot_snpasso.R

@@ -56,8 +56,9 @@
 #'
 #' @section Hidden graphics parameters:
 #' A number of graphics parameters can be passed via `...`. For
-#' example, `bgcolor` to control the background color and `altbgcolor`
-#' to control the background color on alternate chromosomes.
+#' example, `bgcolor` to control the background color,`altbgcolor`
+#' to control the background color on alternate chromosomes,
+#' `altcol` to control the point color on alternate chromosomes,
 #' `cex` for character expansion for the points (default 0.5),
 #' `pch` for the plotting character for the points (default 16),
 #' and `ylim` for y-axis limits.

tests/testthat/test-create_snpinfo.R

@@ -0,0 +1,23 @@
+context("create snpinfo from cross2 object")
+
+test_that("create_snpinfo works", {
+
+    if(isnt_karl()) skip("this test only run locally")
+
+    file <- paste0("https://raw.githubusercontent.com/rqtl/",
+                   "qtl2data/master/DOex/DOex.zip")
+
+    DOex <- read_cross2(file)
+
+    snpinfo <- create_snpinfo(DOex)
+
+    fg <- do.call("cbind", DOex$founder_geno)
+
+    expect_equal(sum(colSums(fg==0)==0), nrow(snpinfo))
+
+    n <- sapply(DOex$founder_geno, function(a) sum(colSums(a==0)==0))
+    expect_equivalent(unclass(table(snpinfo$chr)), n)
+
+    expect_equivalent(snpinfo$sdp, calc_sdp(t(fg[,colSums(fg==0)==0])))
+
+})