Working with messy genetic maps in R/qtl

JT Lovell

2016-10-14

library(qtl)

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

Before begining QTL mapping, we must generate a genetic map. This is often the hardest part of the whole analysis. For example, we may start with a badly ordered map, with far more markers than are necessary - such genotype matrices are typical of RAD or other NGS data.

Below is a short workflow to go from a map that has bad marker order and several regions with more markers than are necessary to one that is smaller and correctly ordered.

To do this, we employ the tools of R/qtl and a few wrapper tools from the package qtlTools. The following tutorial gives an example of this, and complements the very helpful tutorial from R/qtl: http://www.rqtl.org/tutorials/geneticmaps.pdf

Part 1: Make a map with some marker order and marker density problems and

set.seed(42)
map<-sim.map(len = c(50,50,20,30), n.mar = c(25, 10, 10, 50), include.x=F)
plot(map)

cross0<-sim.cross(map, n.ind=50, type="f2",
          error.prob=0.001, missing.prob=0.001, map.function="kosambi")
plot.rf(cross0)

##########
jitterMarkerOrder<-function(cross, chr){
  mars<-1:nmar(cross)[chr]
  set.seed(42)
  badorder<-order(jitter(mars, factor = 10))
  cross<-switch.order(cross, chr = chr, order = badorder, 
                      error.prob=0.001, map.function="kosambi")
}
##########
cross1<-cross0
for(i in 1:nchr(cross1))   cross1<-jitterMarkerOrder(cross=cross1, chr = i)
plot.rf(cross1)

newmap<-est.map(cross1, error.prob=0.001, map.function="kosambi")
cross1 <- replace.map(cross1, newmap)
cross1<-est.rf(cross1)

Part 2: Toss out markers that are not informative

cross2<-dropSimilarMarkers(cross1, rf.threshold = 0.03)
## initial n markers:  95 
## final n markers:  32
plot.map(cross1, cross2, main = "comparison of full and culled maps")

Part 3: Reorder markers iteratively

cross3<-repRipple(cross2, error.prob=0.001, map.function="kosambi",window = 6)
## running ripple for chromosome:  1 
## ripple run # 1 ... n crossovers reduced by 12 
## ripple run # 2 ... n crossovers not reduced
## running ripple for chromosome:  2 
## ripple run # 1 ... n crossovers reduced by 17 
## ripple run # 2 ... n crossovers not reduced
## running ripple for chromosome:  3 
## ripple run # 1 ... n crossovers not reduced
## running ripple for chromosome:  4 
## ripple run # 1 ... n crossovers not reduced
## re-estimating genetic map
plot.rf(cross2, main = "recombination fractions before ripple")

plot.rf(cross3, main = "recombination fractions after ripple")