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),]) })