The internal data structures for R/qtl2 version 0.5 have undergone major changes. In addition, in several places, QTL analysis now requires additional steps and the need to keep track of additional objects.

The main differences are that the genotype probabilities from qtl2geno::calc_genoprob and the genome scan results from qtl2scan::scan1 (as well as estimated QTL effects from qtl2scan::scan1coef and qtl2scan::scan1blup) no longer include the genetic map of marker/pseudomarker locations.

To convert data structures from the old format to the new format, we provide functions convert_probs() and convert_scan1() in the file convert_04_to_05.R.

Genotype probabilities

First, the genotype probabilities produced by qtl2geno::calc_genoprob no longer include the genetic map of marker/pseudomarker positions, and calculating these probabilities generally requires one extra step.

Let’s load some sample data.

library(qtl2geno)
iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2geno"))

The calc_genoprob function no longer takes arguments step, stepwidth, and off_end. Instead, you need to first call insert_pseudomarkers() to insert pseudomarker positions.

map <- insert_pseudomarkers(iron$gmap, step=1)

You then call calc_genoprob using that map, with the pseudomarker positions inserted. Alternatively, you could use iron$gmap directly, and do the calculations solely at the markers.

pr <- calc_genoprob(iron, map, error_prob=0.002, map_function="c-f")

Genome scan

A genome scan with qtl2scan::scan1 is performed as before, but the result no longer includes the genetic map of marker/pseudomarker positions, and so later use of these results generally require that you provide both the results and the map.

We first get the special covariates for the X chromosome.

Xcovar <- get_x_covar(iron)

Then we perform the genome scan.

library(qtl2scan)
out <- scan1(pr, iron$pheno, Xcovar=Xcovar)

We can use qtl2scan::find_peaks to get a summary of inferred QTL meeting some LOD threshold, but this function needs both the scan1 output and the map of marker/pseudomarker positions.

find_peaks(out, map, threshold=4)
##   lodindex lodcolumn chr  pos       lod
## 1        1     liver   2 56.8  4.956559
## 2        1     liver   7 50.1  4.050449
## 3        1     liver  16 28.6  7.682563
## 4        2    spleen   8 13.6  4.302929
## 5        2    spleen   9 56.6 12.062549

Similarly, to plot the results, we need to provide both the scan1 output and the map.

library(qtl2plot)
plot(out, map, ylim=c(0, maxlod(out)*1.05))
plot(out, map, lodcolumn=2, col="violetred", add=TRUE)

Moreover, to subset scan1 results, you need to provide both the scan1 object and the map.

out_selectedchr <- subset(out, map, chr=c(2, 7, 8, 9, 16))

Previously, when subsetting scan1 objects with brackets ([ ]), it was done with scan1_output[chr, lodcolumns], but we have dropped this feature for now. To subset by chromosome, you need to explicitly use subset(), while the bracket notation just subsets the object as a matrix scan1_output[positions, lodcolumns].

SNP association analyses

In analysis of SNP association, with a table of SNP information and the function qtl2scan::genoprob_to_snpprob(), there’s now an additional step, and the SNP information table is no longer included in the results, but needs to be provided separately.

Let’s first load the DOex example, which is a subset of the data from subset of the data from Recla et al. (2014).

file <- paste0("https://raw.githubusercontent.com/rqtl/",
               "qtl2data/master/DOex/DOex.zip")
DOex <- read_cross2(file)

Let’s first do a single-QTL genome scan with the haplotype dosages.

pr <- calc_genoprob(DOex, error_prob=0.002)
apr <- genoprob_to_alleleprob(pr)
k <- calc_kinship(apr, "loco")
sex <- (DOex$covar$Sex == "male")*1
names(sex) <- rownames(DOex$covar)
out <- scan1(apr, DOex$pheno, k, sex)

Consider the peak marker on chromosome 2. We find it’s physical location:

marker <- rownames(max(out, DOex$gmap, chr="2"))
peak_Mbp <- DOex$pmap[["2"]][marker]

We now grab SNP information for the region.

tmpfile <- tempfile()
file <- "https://raw.githubusercontent.com/rqtl/qtl2data/master/DOex/c2_snpinfo.rds"
download.file(file, tmpfile, quiet=TRUE)
snpinfo <- readRDS(tmpfile)
unlink(tmpfile)

We convert the founder genotypes to a “strain distribution pattern” (SDP): an integer whose binary encoding corresponds to the 8 founders’ genotypes.

snpinfo$sdp <- calc_sdp(snpinfo[,-(1:4)])

Previously, qtl2scan::genoprob_to_snpprob would first partition the SNPs into groups of equivalent SNPs (within the same marker interval and having the same strain distribution pattern) and then collapse the genotype or haplotype probabilities to SNP genotype or allele probabilities.

The partitioning/indexing of SNPs is now a separate step:

snpinfo <- index_snps(DOex$pmap, snpinfo)

This adds a few columns, indicating which SNPs are the non-equivalent ones, and which marker intervals they belong to.

We then can call qtl2scan::genoprob_to_snpprob(), giving both the probabilities and the augmented SNP info table.

snp_apr <- genoprob_to_snpprob(apr, snpinfo)

We can then perform a genome scan.

out_snps <- scan1(snp_apr, DOex$pheno, k[["2"]], sex)

The function qtl2plot::plot_snpasso() can be used to plot the results. You need to provide both the scan1 output and the SNP information table (augmented by index_snps, with the indexing information).

plot_snpasso(out_snps, snpinfo)

To get a table of the SNPs with the largest LOD scores, again use the function top_snps(), but now you need to provide the SNP information table.

top_snps(out_snps, snpinfo)
##            snp_id chr  pos_Mbp alleles AJ B6 129 NOD NZO CAST PWK WSB sdp index interval on_map      lod
## 3264  rs212414861   2 96.75489     C|G  1  1   1   1   3    3   1   1  48  3264       64  FALSE 6.931185
## 3265  rs229578122   2 96.75494     T|A  1  1   1   1   3    3   1   1  48  3264       64  FALSE 6.931185
## 3266  rs254318131   2 96.75494     C|T  1  1   1   1   3    3   1   1  48  3264       64  FALSE 6.931185
## 3267  rs217679969   2 96.75495     G|T  1  1   1   1   3    3   1   1  48  3264       64  FALSE 6.931185
## 3268  rs238404461   2 96.75496     T|G  1  1   1   1   3    3   1   1  48  3264       64  FALSE 6.931185
## 3269  rs262749925   2 96.75497     C|G  1  1   1   1   3    3   1   1  48  3264       64  FALSE 6.931185
## 3271  rs231282689   2 96.75498   C|G/T  1  1   1   1   3    3   1   1  48  3264       64  FALSE 6.931185
## 3274  rs260286709   2 96.75518     G|A  1  1   1   1   3    3   1   1  48  3264       64  FALSE 6.931185
## 3275   rs27396282   2 96.75522     T|C  1  1   1   1   3    3   1   1  48  3264       64  FALSE 6.931185
## 3283  rs263841533   2 96.75541     T|C  1  1   1   1   3    3   1   1  48  3264       64  FALSE 6.931185
## 3288  rs231205920   2 96.75557     G|A  1  1   1   1   3    3   1   1  48  3264       64  FALSE 6.931185
## 3289  rs242885221   2 96.75656     T|C  1  1   1   1   3    3   1   1  48  3264       64  FALSE 6.931185
## 5242  rs586746690   2 96.88443     C|T  1  1   3   1   3    3   1   1  52  5242       64  FALSE 5.904999
## 5426   rs49002164   2 96.89204     G|A  1  1   3   1   3    3   1   1  52  5242       64  FALSE 5.904999
## 7170  rs244595995   2 96.96972     A|T  1  1   3   1   3    3   1   1  52  5242       64  FALSE 5.904999
## 16182 rs220351620   2 97.57962     A|G  1  1   1   1   3    1   1   1  16 16182       65  FALSE 5.994374
## 22474  rs52579091   2 98.07299     C|T  1  1   1   1   3    1   1   1  16 16182       65  FALSE 5.994374
## 23017 rs243489710   2 98.11061     C|T  1  1   1   1   3    1   1   1  16 16182       65  FALSE 5.994374
## 23184 rs244316263   2 98.12529     C|A  1  1   1   1   3    1   1   1  16 16182       65  FALSE 5.994374
## 23186 rs219729956   2 98.12534     C|A  1  1   1   1   3    1   1   1  16 16182       65  FALSE 5.994374
## 23360 rs235315566   2 98.20433     G|T  1  1   1   1   3    1   1   1  16 16182       65  FALSE 5.994374
## 23494 rs250167663   2 98.21481     C|T  1  1   1   1   3    1   1   1  16 16182       65  FALSE 5.994374
## 23518 rs234831418   2 98.21704     T|G  1  1   1   1   3    1   1   1  16 16182       65  FALSE 5.994374
## 23601 rs240832432   2 98.22253     C|A  1  1   1   1   3    1   1   1  16 16182       65  FALSE 5.994374
## 23649 rs220815439   2 98.22511     G|A  1  1   1   1   3    1   1   1  16 16182       65  FALSE 5.994374
## 23830 rs579950897   2 98.24240     C|T  1  1   1   1   3    1   1   1  16 16182       65  FALSE 5.994374
## 26015 rs224994851   2 98.39510     A|T  1  1   1   1   3    1   1   1  16 16182       65  FALSE 5.994374
## 26213 rs248208898   2 98.42249     C|T  1  1   1   1   3    1   1   1  16 16182       65  FALSE 5.994374
## 26554 rs245525925   2 98.45411     C|A  1  1   1   1   3    1   1   1  16 16182       65  FALSE 5.994374
## 26891 rs229631954   2 98.47723     C|T  1  1   1   1   3    1   1   1  16 16182       65  FALSE 5.994374