Nothing
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)
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.