tests/testthat/test-calcerrorlod.R

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

context("calc_errorlod")

test_that("calc_errorlod works for an intercross", {

    # this is not much more than a regression test, really
    iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))

    pr <- calc_genoprob(iron, err=0.002, map_function="c-f")
    err <- calc_errorlod(iron, pr)

    expected_1to5_chr1 <- structure(c(-2.82278691244427, 0, -2.82278691244427, -3.24248113822165,
                                      0, -2.76132216098372, 0, -2.76132216098372, -3.23746097636371,
                                      0, -3.00415918819379, 0, -3.00415918819379, -2.99576369311289,
                                      0), .Dim = c(5L, 3L),
                                    .Dimnames = list(c("1", "2", "3", "4", "5"),
                                                     c("D1Mit18", "D1Mit80", "D1Mit17")))
    expect_equal(err[[1]][1:5,], expected_1to5_chr1)

    ind <- c(1:5, 244:249)
    expected_chrX <- structure(c(-3.10992625978842, 0, -3.10992625978842, -3.10992625978842,
                                 0, -3.10992625978842, -3.10992625978842, 0, 0, -3.10992625978842,
                                 -3.10992625978842, -3.10992625978842, 0, -3.10992625978842, -3.10992625978842,
                                 0, -3.10992625978842, -3.10992625978842, 0, 0, -3.10992625978842,
                                 -3.10992625978842), .Dim = c(11L, 2L),
                               .Dimnames = list(c("1", "2", "3", "4", "5", "244", "245", "246", "247", "248", "249"),
                                                c("DXMit16", "DXMit186")))
    expect_equal(err$X[ind,], expected_chrX)

    # missing values -> error lod == 0
    for(i in seq(along=err)) {
        if(any(iron$geno[[i]]==0)) {
            expect_true(max(abs(err[[i]][iron$geno[[i]]==0])) < 1e-15)
        }
    }

    # check a locus on chr 3
    p <- pr[[7]][,,4]
    g <- iron$geno[[7]][,4]
    expected <- g
    expected[g==1] <- log10((1-p[,1])/p[,1]*0.25/0.75)[g==1]
    expected[g==2] <- log10((1-p[,2])/p[,2])[g==2]
    expected[g==3] <- log10((1-p[,3])/p[,3]*0.25/0.75)[g==3]

    expect_equal(err[[7]][,4], expected)

    # check a locus on the X chr
    p <- pr$X[,,2]
    g <- iron$geno$X[,2]
    male <- (iron$covar$sex=="m") # genotypes 1/3
    female_f <- (iron$covar$sex=="f" & iron$covar$cross_direction=="(SxB)x(SxB)") # genotypes 1/2
    female_r <- (iron$covar$sex=="f" & iron$covar$cross_direction=="(BxS)x(BxS)") # genotypes 2/3
    expected <- g
    expected[g==1 & male] <- log10((1-p[,5])/p[,5])[g==1 & male]
    expected[g==3 & male] <- log10((1-p[,6])/p[,6])[g==3 & male]
    expected[g==1 & female_f] <- log10((1-p[,1])/p[,1])[g==1 & female_f]
    expected[g==2 & female_f] <- log10((1-p[,2])/p[,2])[g==2 & female_f]
    expected[g==2 & female_r] <- log10((1-p[,3])/p[,3])[g==2 & female_r]
    expected[g==3 & female_r] <- log10((1-p[,4])/p[,4])[g==3 & female_r]

    expect_equal(err$X[,2], expected)

})


test_that("calc_errorlod works for RIL", {

    # this is not much more than a regression test, really
    grav2 <- read_cross2(system.file("extdata", "grav2.zip", package="qtl2"))

    pr <- calc_genoprob(grav2, err=0.002, map_function="c-f")
    err <- calc_errorlod(grav2, pr)

    ind <- c(9L, 23L, 25L, 42L, 44L, 49L, 54L, 66L, 134L, 158L)
    mar <- c(32L, 35L, 36L, 48L, 56L)

    expected <- structure(c(-5.82983567140969, -5.82983567140969, -5.82983567140969,
                            -5.82983567140969, -5.82983567140969, -5.82983288817781, -5.82983567140969,
                            -5.82983567140969, -5.82983567140043, -5.82983567140969, -5.60155654532719,
                            -5.60155654532719, -5.60155654532719, -5.60155654532719, -5.60155654532719,
                            -2.61600519342757, -5.60155654532719, -5.60155654532719, -5.60016580984213,
                            -5.60155654266195, -8.29839428989781, -8.29839428989781, -8.29839428989781,
                            -8.29839428989781, -8.29839428989781, -5.31493687487081, -8.29839428989475,
                            -8.29839428989775, -5.47724186085683, -8.29784500171819, -14.3504722627229,
                            -14.350472262709, -14.3504722627229, -14.3504722627229, -14.3504722627229,
                            -14.3504722627229, -14.3504722578884, -14.3504499660042, -14.3504722627229,
                            -14.3504722627119, -5.86578939924429, -4.92654601370327, -5.86578939924429,
                            -7.62403763042528, -7.62403763042528, -5.86578939924429, -7.62403763042528,
                            -7.62392376500139, -7.62403763042527, -7.62392376500139),
                          .Dim = c(10L, 5L),
                          .Dimnames = list(c("9", "23", "25", "42", "44", "49", "54", "66", "134", "158"),
                                           c("DF.300C", "CD.179L", "CH.88L", "BH.81L-Col", "FD.345C")))

    expect_equal(err[[5]][ind,mar], expected)

    # missing values -> error lod == 0
    for(i in seq(along=err)) {
        if(any(grav2$geno[[i]]==0)) {
            expect_true(max(abs(err[[i]][grav2$geno[[i]]==0])) < 1e-15)
        }
    }

    # check a locus on chr 4
    p <- pr[[4]][,,35]
    g <- grav2$geno[[4]][,35]
    expected <- g
    for(i in 1:2) expected[g==i] <- log10((1-p[,i])/p[,i])[g==i]
    expect_equal(err[[4]][,35], expected)

})


test_that("calc_errorlod works for RIL", {

    skip_if(isnt_karl(), "this test only run locally")

    # this is not much more than a regression test, really
    file <- paste0("https://raw.githubusercontent.com/rqtl/",
                   "qtl2data/main/DO_Recla/recla.zip")
    recla <- read_cross2(file)
    recla <- recla[c(1:2,53:54), c("19","X")] # subset to 4 mice and 2 chromosomes
    probs <- calc_genoprob(recla, err=0.002)

    err <- calc_errorlod(recla, probs)

    expected_c19 <- structure(c(-9.32903376530038, -8.47970241856038, 2.32305910951409,
                                -5.89845329201031, -9.47388004360523, -8.33322922086628, -3.99257600876601,
                                -4.20445567758459, -7.50493914753867, -9.14775695748182, -9.09524930864587,
                                -6.79076628467775, -5.15284031339008, -6.99865791362416, -9.76892294431087,
                                -5.40229996503646), .Dim = c(4L, 4L),
                              .Dimnames = list(c("1", "4", "58", "59"),
                                               c("UNC190139327", "UNC190279660", "backupUNC190139992", "UNC190140427")))

    expect_equal(err[[1]][,7:10], expected_c19)

    expected_cX <- structure(c(-6.08215766199779, -18.6736643714498, -6.5183772856749,
                               -7.07293927915729, 0, -16.5645778105845, -0.442200926998983,
                               -1.14542766590646, -17.9466342307662, -12.3112920811831, -16.8345767687497,
                               7.53542542841645, -16.9579570002804, -16.9636066624958, -18.3869460431974,
                               -9.65729043655822, -18.2175770053576, -9.87948829125263, -17.2496166243865,
                               -7.59183451077769), .Dim = 4:5, .Dimnames = list(c("1", "4",
                               "58", "59"), c("UNC200034449", "UNC200261375", "UNC200062736",
                               "backupUNC200189557", "UNC200192385")))

    expect_equal(err$X[,266:270], expected_cX)

})