tests/testthat/test-max_scan1.R

In rqtl/qtl2: Quantitative Trait Locus Mapping in Experimental Crosses

context("max_scan1")

test_that("max_scan1 works for intercross with two phenotypes", {

    # read data
    iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))

    # calculate genotype probabilities
    map <- insert_pseudomarkers(iron$gmap, step=1)
    probs <- calc_genoprob(iron, map, error_prob=0.002)

    # grab phenotypes and covariates; ensure that covariates have names attribute
    pheno <- iron$pheno
    covar <- match(iron$covar$sex, c("f", "m")) # make numeric
    names(covar) <- rownames(iron$covar)
    Xcovar <- get_x_covar(iron)

    # perform genome scan
    out <- scan1(probs, pheno, addcovar=covar, Xcovar=Xcovar)

    # maximum of first column
    expected <- data.frame(chr="16",
                           pos=28.6,
                           liver=max(out[,1]),
                           stringsAsFactors=FALSE)
    rownames(expected) <- "c16.loc29"
    expect_equal(max(out,map), expected)

    # maximum of spleen column
    expected <- data.frame(chr="9",
                           pos=56.6,
                           spleen=max(out[,2]),
                           stringsAsFactors=FALSE)
    rownames(expected) <- "c9.loc57"
    expect_equal(max(out, map, lodcolumn="spleen"), expected)

    # maximum of first column on chr 2
    expected <- data.frame(chr="2",
                           pos=56.8,
                           liver=maxlod(subset(out, map, "2", 1)),
                           stringsAsFactors=FALSE)
    rownames(expected) <- "D2Mit17"
    expect_equal(max(out, map, chr="2"), expected)

    # test that it works if output or map are subsetted
    expect_equal( max(out, map, chr="2"), max(subset(out, map, chr="2"), map) )
    expect_equal( max(out, map, chr="2"), max(out, map["2"]) )

    # test that it works if output has rows shuffled
    out_shuffled <- out[sample(1:nrow(out)),,drop=FALSE]
    class(out_shuffled) <- c("scan1", "matrix")
    expect_equal( max(out_shuffled, map), max(out, map))

    expect_equal(max_scan1(out), c(liver= 6.35264382891418))
    expect_equal(max_scan1(out, lodcolumn=2), c(spleen=12.5986057120873))
    expect_equal(max_scan1(out, lodcolumn="spleen"), c(spleen=12.5986057120873))

    # warning if you give lodcolumn as a vector
    expect_warning( expect_equal(max(out, lodcolumn=1:2), c(liver= 6.35264382891418)) )
    expected <- data.frame(chr="16",
                           pos=28.6,
                           liver=max(out[,1]),
                           stringsAsFactors=FALSE)
    rownames(expected) <- "c16.loc29"
    expect_warning( expect_equal( max(out, map, lodcolumn=1:2), expected) )

    # warning if you give chr but not map
    expect_warning( expect_equal( max(out, lodcolumn=2, chr=16), c(spleen=12.5986057120873)) )

    # results for all LOD score columns if lodcolumn=NULL
    expect_equal( max(out, lodcolumn=NULL),
                 c(liver= 6.35264382891418, spleen=12.5986057120873) )


    expected <- data.frame(lodindex=1:ncol(out),
                           lodcolumn=colnames(out),
                           chr=c("16","9"),
                           pos=c(28.6, 56.6),
                           lod=apply(out, 2, max),
                           stringsAsFactors=FALSE)
    rownames(expected) <- NULL
    expect_equal( max(out, map, lodcolumn=NULL), expected )

    expected <- data.frame(lodindex=1:ncol(out),
                           lodcolumn=colnames(out),
                           chr=c("16","16"),
                           pos=c(28.6, 30.6),
                           lod=apply(subset(out, map, chr=16), 2, max),
                           stringsAsFactors=FALSE)
    rownames(expected) <- NULL
    expect_equal( max(out, map, lodcolumn=NULL, chr=16), expected )

    expected <- data.frame(lodindex=1:ncol(out),
                           lodcolumn=colnames(out),
                           chr=c("3","9"),
                           pos=c(25.1, 56.6),
                           lod=apply(subset(out, map, chr=c(3,9)), 2, max),
                           stringsAsFactors=FALSE)
    rownames(expected) <- NULL
    expect_equal( max(out, map, lodcolumn=NULL, chr=c(3,9)), expected )
    expect_equal( max(out, map, lodcolumn=NULL, chr=c(9,3)), expected )

    # max handles NA (here, just guessing that I'm unlikely to NA-out the max
    set.seed(1331427)
    out_na <- out
    out_na[5,2] <- NA
    out_na[sample(nrow(out_na), 10),1] <- NA
    expect_equal(max(out), max(out_na))
    expect_equal(max(out, map), max(out_na, map))
    expect_equal(max(out, lodcolumn=NULL), max(out_na, lodcolumn=NULL, na.rm=TRUE))
    expect_equal(max(out, map, lodcolumn=NULL), max(out_na, map, lodcolumn=NULL, na.rm=TRUE))

    # max gives error if lodcolumn out of range
    expect_error(max(out, lodcolumn="blah"), 'LOD column "blah" not found.')
    expect_error(max(out, lodcolumn=0),
                 "column [0] out of range (should be in 1 - 2)", fixed=TRUE)
    expect_error(max(out, lodcolumn=3),
                 "column [3] out of range (should be in 1 - 2)", fixed=TRUE)

    # no column names
    colnames(out) <- colnames(out_na) <- NULL
    expect_equal(max(out), max(out_na))
    expect_equal(max(out, map), max(out_na, map))
    expect_equal(max(out, lodcolumn=2), max(out_na, lodcolumn=2))
    expect_equal(max(out, map, lodcolumn=2), max(out_na, map, lodcolumn=2))
    expect_equal(max(out, lodcolumn=NULL), max(out_na, lodcolumn=NULL, na.rm=TRUE))
    expect_equal(max(out, map, lodcolumn=NULL), max(out_na, map, lodcolumn=NULL, na.rm=TRUE))

})

test_that("maxlod works for intercross with two phenotypes", {

    # read data
    iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))

    # calculate genotype probabilities
    map <- insert_pseudomarkers(iron$gmap, step=1)
    probs <- calc_genoprob(iron, map, error_prob=0.002)

    # grab phenotypes and covariates; ensure that covariates have names attribute
    pheno <- iron$pheno
    covar <- match(iron$covar$sex, c("f", "m")) # make numeric
    names(covar) <- rownames(iron$covar)
    Xcovar <- get_x_covar(iron)

    # perform genome scan
    out <- scan1(probs, pheno, addcovar=covar, Xcovar=Xcovar)

    # 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")))))

    expect_equal(maxlod(out, map, "2"), max(unclass(subset(out, map,"2")), na.rm=TRUE))
    expect_equal(maxlod(out, map, "9"), max(unclass(subset(out, map,"9")), na.rm=TRUE))

    expect_equal( maxlod(out, map, c("2","8", "11")),
                 max( maxlod(out, map, "2"), maxlod(out, map, "8"), maxlod(out, map, "11") ) )

    # really could be using integers
    expect_equal( maxlod(out, map, c(2,8, 11)),
                 max( maxlod(out, map, 2), maxlod(out, map, 8), maxlod(out, map, 11) ) )

    # test that it works if output or map are subsetted
    expect_equal( max(out, map, chr="2"), max(subset(out, map, chr="2"), map) )
    expect_equal( max(out, map, chr="2"), max(out, map["2"]) )

    # test that it works if output has rows shuffled
    out_shuffled <- out[sample(1:nrow(out)),,drop=FALSE]
    class(out_shuffled) <- c("scan1", "matrix")
    expect_equal( max(out_shuffled, map), max(out, map))

    # maxlod handles NA (here, just guessing that I'm unlikely to NA-out the max
    set.seed(1331427)
    out_na <- out
    out_na[5,2] <- NA
    out_na[sample(nrow(out_na), 10),1] <- NA
    expect_equal(maxlod(out), maxlod(out_na))
    expect_equal(maxlod(out, map, 11), maxlod(out_na, map, 11))

})