In jtlovell/qtlTools: Data Processing and Plotting in Association with R/qtl

knitr::opts_chunk$set(echo = TRUE, fig.width = 7, fig.height = 7, warning = FALSE)
library(knitr)
library(qtl)
library(qtlTools)
library(plyr)

library(devtools)
install_github("jtlovell/qtlTools")

library(qtlTools)
library(qtl)
library(ggplot2)
library(lsmeans)
library(qtlpvl)

Simulate a few phenotypes in a five-chromosome map

For this tutorial, we will use the dataset multitrait from R/qtl.

data(multitrait)
cross<-multitrait
cross<-calc.genoprob(cross)
phes<-phenames(cross)

In R/qtl there is a pipeline to analyze multiple traits in a single model. Here, we will discuss analyzing them independently, then concatenating the results

Scanones

In the simplest case, we can run a bunch of scanones

s1<-scanone(cross, pheno.col = phes, method = "hk")
cols<-rainbow(length(phes))
plot(s1, type="n", ylim = c(0,max(as.matrix(s1[,-c(1:2)]))), 
     ylab = "LOD Score", main = "all phenotypes plotted together")
for(i in 1:length(phes)) plot(s1, add=T, lodcolumn = i, col = cols[i])

QTL peaks from scanones

This simplest way to look at these data is to pull out only the significant QTL peaks

perms<-scanone(cross, pheno.col = phes, method = "hk", n.perm = 100, verbose=F)
print(pullSigQTL(cross, pheno.col=phes,
                 s1.output = s1,  perm.output = perms, 
                 returnQTLModel=FALSE, alpha = 0.05,
                 controlAcrossCol = TRUE))

However, if we want to look at model statistics etc. it is best to create multiple qtl models for each phenotype

mods<-pullSigQTL(cross, pheno.col=phes,
                 s1.output = s1,  perm.output = perms, 
                 returnQTLModel=TRUE, alpha = 0.05,
                 controlAcrossCol = TRUE)

This returns a list of fitted QTL models, which we could use to calculate confidence intervals etc.

print(mods[1:5])

We can also plot confidence intervals directly from scanone data, however, it is important to note that multiple peaks on a single chromosome can violate assumptions of such intervals.

cis<-calcCis(cross,s1.output=s1, perm.output=perms)
head(cis)

We might be interested in where these QTL are - we can plot these on the genetic map

segmentsOnMap(cross, calcCisResults = cis, legendCex = .35)

Multiple QTL models

The approach above is nice because it is super fast, but it does not necessarily produce the best models (e.g. there can only be one QTL/chromosome). One could use stepwiseqtl to generate the models instead, however, scantwo permutations are needed for this appraoch, which is slow for many traits.

For the purposes of this tutorial, we will use the models generated through scanones to make statistical inference. The function qtlStats takes a qtl model and returns fitted statistics. For a single phenotype, we'd run it as follows:

mod<-mods[["X3.Methylsulfinylpropyl"]]
form<-paste("y ~",paste(mod$altname, collapse = " + "))
mod<-refineqtl(cross, qtl = mod, formula = form, verbose=F, method="hk")
qtlStats(cross, 
         pheno.col = "X3.Methylsulfinylpropyl", 
         form = form, 
         mod = mod)

Since mods is a list, we can loop through it using lapply First toss out an NULL QTL models

mods<-mods[sapply(mods,length)!=1] # a NULL QTL model has length = 1, toss these

Next do the loop, where x represents each phenotype.

stats<-lapply(names(mods), function(x){
  mod<-mods[[x]]
  form<-paste("y ~",paste(mod$altname, collapse = " + "))
  mod<-refineqtl(cross, qtl = mod, formula = form, 
                verbose=F, method="hk", pheno.col = x)
  qtlStats(cross, 
         pheno.col = x, 
         form = form, 
         mod = mod)
})
stats<-do.call(rbind, stats)
print(head(stats))

Plotting multiple QTL models

The most basic way to present these data is to show the position of QTL intervals on a genetic map. qtlTools::segmentsOnMap provides one method to do this.

segmentsOnMap(cross=cross, phe=stats$phenotype, 
              chr=stats$chr, 
              l = stats$lowposition, lwd=2,
              h =stats$highposition,legendCex=.35, leg.inset=.01)