tests/testthat/test-top_snps.R

context("top snps from snp association analysis")

test_that("top_snps() works", {
    skip_if(isnt_karl(), "this test only run locally")

    # load example DO data from web
    file <- paste0("https://raw.githubusercontent.com/rqtl/",
                   "qtl2data/main/DOex/DOex.zip")
    DOex <- read_cross2(file)

    # subset to chr 2
    DOex <- DOex[,"2"]

    # calculate genotype probabilities and convert to allele probabilities
    pr <- calc_genoprob(DOex, error_prob=0.002, cores=2) # multi-core
    apr <- genoprob_to_alleleprob(pr)

    # download snp info from web
    tmpfile <- tempfile()
    file <- paste0("https://raw.githubusercontent.com/rqtl/",
                   "qtl2data/main/DOex/c2_snpinfo.rds")
    download.file(file, tmpfile, quiet=TRUE)
    snpinfo <- readRDS(tmpfile)
    unlink(tmpfile)

    # identify equivalent SNPs
    snpinfo <- index_snps(DOex$pmap, snpinfo)

    # convert to snp probabilities
    snp_pr <- genoprob_to_snpprob(apr, snpinfo)

    # perform SNP association analysis (here, ignoring residual kinship)
    out_snps <- scan1(snp_pr, DOex$pheno)

    # table with top SNPs
    result <- top_snps(out_snps, snpinfo)

    expected <- structure(list(snp_id = c("rs264175039", "rs227493750", "rs218001398", "rs212907587",
                                          "rs245211479", "rs580781679", "rs213082844", "rs224780822",
                                          "rs247650681", "rs241598130", "rs265940269", "rs232479852",
                                          "rs212414861", "rs229578122", "rs254318131", "rs217679969",
                                          "rs238404461", "rs262749925", "rs231282689", "rs260286709",
                                          "rs27396282",  "rs245362008", "rs265321723", "rs263841533",
                                          "rs231205920", "rs242885221", "rs240030573", "2:96756871_T/A",
                                          "rs230778665", "rs216514125", "rs578994109", "rs224938158",
                                          "rs216200700", "rs257147584", "rs263558150", "rs241319958",
                                          "rs245545514", "rs586746690", "rs49002164", "rs229778455",
                                          "rs252775618", "2:96945646_C/A", "2:96945653_C/A", "rs232602612",
                                          "rs244595995", "rs221707344", "rs251046484", "rs258630238",
                                          "rs235073020", "rs237252713", "rs230537549", "2:97362128_T/C",
                                          "rs220351620", "rs52579091", "rs243489710", "rs244316263",
                                          "rs219729956", "rs235315566", "rs250167663", "rs234831418",
                                          "rs240832432", "rs220815439", "rs579950897", "rs224994851",
                                          "rs248208898", "rs245525925", "rs229631954"),
                               chr = rep("2", 67),
                               pos_Mbp = c(96.500224, 96.500276, 96.506074, 96.651629, 96.653345,
                                           96.685074, 96.752351, 96.752682, 96.753044, 96.754757,
                                           96.754758, 96.754759, 96.754889, 96.754937, 96.75494,
                                           96.754955, 96.754963, 96.75497,  96.754979, 96.75518,
                                           96.755219, 96.755357, 96.755358, 96.755411, 96.755571,
                                           96.756558, 96.75669,  96.756871, 96.757316, 96.757973,
                                           96.758174, 96.758487, 96.77603,  96.777354, 96.808596,
                                           96.879951, 96.883644, 96.884435, 96.892038, 96.892047,
                                           96.936957, 96.945646, 96.945653, 96.967673, 96.969724,
                                           96.981743, 96.982601, 96.984073, 96.991667, 96.995521,
                                           97.010408, 97.362128, 97.579618, 98.072994, 98.110607,
                                           98.125288, 98.125343, 98.204326, 98.214808, 98.217039,
                                           98.222529, 98.225114, 98.242395, 98.395104, 98.422488,
                                           98.45411, 98.477226),
                               alleles = c("A|C", "C|T", "C|T", "G|A", "C|A", "T|C", "C|T",
                                           "T|A", "G|A", "A|G", "A|G", "G|A", "C|G", "T|A", "C|T", "G|T",
                                           "T|G", "C|G", "C|G/T", "G|A", "T|C", "T|C", "G|A", "T|C", "G|A",
                                           "T|C", "G|A", "T|A", "A|C", "A|T", "G|A", "T|C", "C|A", "C|T",
                                           "T|C", "A|G", "C|T", "C|T", "G|A", "T|G", "C|A", "C|A", "C|A",
                                           "T|G", "A|T", "C|T", "G|A/T", "T|A", "G|C/A", "A|G", "C|T/A",
                                           "T|C", "A|G", "C|T", "C|T", "C|A", "C|A", "G|T", "C|T", "T|G",
                                           "C|A", "G|A", "C|T", "A|T", "C|T", "C|A", "C|T"),
                               AJ = rep(1, 67),
                               B6 = rep(1, 67),
                               `129` = c(1, 1, 1, 1, 1, 1, 1, 1,
                                         1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                                         1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 1, 1, 1, 1, 1, 3, 1, 1, 1, 1, 1,
                                         1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
                               NOD = rep(1, 67),
                               NZO = rep(3, 67),
                               CAST = c(1, 1, 1, 1, 1, 1, 1,
                                        1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 1, 1, 3, 3, 3, 1, 1,
                                        1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 1, 1, 1, 1, 1, 3, 1, 1, 1, 1,
                                        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
                               PWK = rep(1, 67),
                               WSB = rep(1, 67),
                               sdp = c(16L, 16L, 16L, 16L, 16L,
                                       16L, 16L, 16L, 16L, 16L, 16L, 16L, 48L, 48L, 48L, 48L, 48L, 48L,
                                       48L, 48L, 48L, 16L, 16L, 48L, 48L, 48L, 16L, 16L, 16L, 16L, 16L,
                                       16L, 16L, 16L, 16L, 16L, 16L, 52L, 52L, 16L, 16L, 16L, 16L, 16L,
                                       52L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L,
                                       16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L),
                               index = c(2L,
                                         2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3264L, 3264L, 3264L,
                                         3264L, 3264L, 3264L, 3264L, 3264L, 3264L, 2L, 2L, 3264L, 3264L,
                                         3264L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 5152L, 5152L, 5242L,
                                         5242L, 5152L, 5152L, 5152L, 5152L, 5152L, 5242L, 5152L, 5152L,
                                         5152L, 5152L, 5152L, 5152L, 5152L, 5152L, 5152L, 5152L, 5152L,
                                         5152L, 5152L, 5152L, 5152L, 5152L, 5152L, 5152L, 26015L, 26015L,
                                         26015L, 26015L),
                               interval = c(64L, 64L, 64L, 64L, 64L, 64L, 64L,
                                            64L, 64L, 64L, 64L, 64L, 64L, 64L, 64L, 64L, 64L, 64L, 64L, 64L,
                                            64L, 64L, 64L, 64L, 64L, 64L, 64L, 64L, 64L, 64L, 64L, 64L, 64L,
                                            64L, 64L, 65L, 65L, 65L, 65L, 65L, 65L, 65L, 65L, 65L, 65L, 65L,
                                            65L, 65L, 65L, 65L, 65L, 65L, 65L, 65L, 65L, 65L, 65L, 65L, 65L,
                                            65L, 65L, 65L, 65L, 66L, 66L, 66L, 66L),
                               on_map = rep(FALSE, 67),
                               lod = c(5.29789577677169, 5.29789577677169, 5.29789577677169,
                                       5.29789577677169, 5.29789577677169, 5.29789577677169, 5.29789577677169,
                                       5.29789577677169, 5.29789577677169, 5.29789577677169, 5.29789577677169,
                                       5.29789577677169, 6.5663685577182,  6.5663685577182,  6.5663685577182,
                                       6.5663685577182,  6.5663685577182,  6.5663685577182,  6.5663685577182,
                                       6.5663685577182,  6.5663685577182,  5.29789577677169, 5.29789577677169,
                                       6.5663685577182,  6.5663685577182,  6.5663685577182,  5.29789577677169,
                                       5.29789577677169, 5.29789577677169, 5.29789577677169, 5.29789577677169,
                                       5.29789577677169, 5.29789577677169, 5.29789577677169, 5.29789577677169,
                                       5.59194100993414, 5.59194100993414, 5.88250557969092, 5.88250557969092,
                                       5.59194100993414, 5.59194100993414, 5.59194100993414, 5.59194100993414,
                                       5.59194100993414, 5.88250557969092, 5.59194100993414, 5.59194100993414,
                                       5.59194100993414, 5.59194100993414, 5.59194100993414, 5.59194100993414,
                                       5.59194100993414, 5.59194100993414, 5.59194100993414, 5.59194100993414,
                                       5.59194100993414, 5.59194100993414, 5.59194100993414, 5.59194100993414,
                                       5.59194100993414, 5.59194100993414, 5.59194100993414, 5.59194100993414,
                                       5.18076131028991, 5.18076131028991, 5.18076131028991, 5.18076131028991)),
                          row.names = c(2L, 3L, 40L, 1593L, 1652L, 2128L, 3197L, 3205L,
                                        3215L, 3258L, 3259L, 3260L, 3264L, 3265L, 3266L, 3267L, 3268L,
                                        3269L, 3271L, 3274L, 3275L, 3276L, 3277L, 3283L, 3288L, 3289L,
                                        3291L, 3296L, 3302L, 3311L, 3316L, 3328L, 3672L, 3691L, 4271L,
                                        5152L, 5229L, 5242L, 5426L, 5427L, 6610L, 6813L, 6814L, 7135L,
                                        7170L, 7253L, 7272L, 7317L, 7457L, 7579L, 7769L, 12874L, 16182L,
                                        22474L, 23017L, 23184L, 23186L, 23360L, 23494L, 23518L, 23601L,
                                        23649L, 23830L, 26015L, 26213L, 26554L, 26891L), class = "data.frame")

    expect_equal(result, expected)

    # top SNPs among the distinct subset at which calculations were performed
    result <- top_snps(out_snps, snpinfo, show_all_snps=FALSE)
    expect_equal(result, expected[c(1,13,36,38,64),])

    # top SNPs within 0.5 LOD
    result <- top_snps(out_snps, snpinfo, drop=0.5)
    expect_equal(result, expected[c(13:21,24:26),])

})
rqtl/qtl2 documentation built on March 20, 2024, 6:35 p.m.