This is a lab tutorial for construction of a genetic map using r/qtl. Example scripts are from http://www.rqtl.org/tutorials/geneticmaps.R. The tutorial can be downloaded at http://www.rqtl.org/tutorials/geneticmaps.pdf

Data summary and filtering

library(qtl)
data(mapthis)
summary(mapthis) ### summarize data
##     F2 intercross
## 
##     No. individuals:    300 
## 
##     No. phenotypes:     1 
##     Percent phenotyped: 100 
## 
##     No. chromosomes:    1 
##         Autosomes:      1 
## 
##     Total markers:      100 
##     No. markers:        100 
##     Percent genotyped:  95.4 
##     Genotypes (%):      AA:26.2  AB:48.2  BB:25.6  not BB:0.0  not AA:0.0
### missing data
par(las=1, cex=0.8, mar=c(4.1, 4.1, 3, 1.1))
plotMissing(mapthis, main="Missing data")

### individuals / markers
plot(ntyped(mapthis), ylab="No. typed markers", main="No. genotypes by individual")

plot(ntyped(mapthis, "mar"), ylab="No. typed individuals", 
     main="No. genotypes by marker")

Segregation Distortion

The function geno.table to inspect the segregation patterns. It calculates the genotype frequencies and also a P-value for a test of departure from the 1:2:1 expected ratios. We will focus on those markers that show significant distortion at the 5% level, after a Bonferroni correction for the multiple tests.

gt <- geno.table(mapthis)
bon.pval.cutoff <- 0.05/totmar(mapthis) ### Bonferroni corrected p-value cutoff
gt[gt$P.value < bon.pval.cutoff, ]
##       chr missing  AA  AB  BB not.BB not.AA       P.value
## C4M2    1       5 100 144  51      0      0  2.686681e-04
## C1M4    1       8   8 210  74      0      0  2.173881e-19
## C2M9    1       5 290   3   2      0      0 2.581933e-184
## C1M21   1      10 200  10  80      0      0  7.058444e-77
## C2M15   1       9   0   1 290      0      0 1.455993e-188
## C2M27   1       7   2 218  73      0      0  2.360267e-23
todrop <- rownames(gt[gt$P.value < bon.pval.cutoff, ])
mapthis <- drop.markers(mapthis, todrop) ### drop markers showing significant segregation distortion

Plot of LOD scores versus estimated recombination fractions for all marker pairs.

mapthis <- est.rf(mapthis)
## Warning in est.rf(mapthis): Alleles potentially switched at markers 
##   C3M16 C2M16 C1M2 C3M9 C2M14 C1M24 C1M1 C1M36 C3M1 C2M25 C1M22 C5M5 C5M7 C1M17 C5M1 C3M5 C5M2 C1M15 C2M24 C2M17 C1M23 C5M6 C1M16 C3M2 C3M10 C3M6 C2M13
rf <- pull.rf(mapthis)
lod <- pull.rf(mapthis, what="lod")
par(mar=c(4.1,4.1,0.6,0.6), las=1, cex=0.8)
plot(as.numeric(rf), as.numeric(lod), xlab="Recombination fraction", ylab="LOD score")

Linkage groups

The author said “The function has two main arguments: max.rf and min.lod. Two markers will be placed in the same linkage groups if they have estimated recombination fraction ≤ max.rf and LOD score ≥ min.lod. The linkage groups are closed via the transitive property. That is, if markers a and b are linked and markers b and c are linked, then all three are placed in the same linkage group. I generally start with min.lod relatively large (say 6 or 8 or even 10). The appropriate value depends on the number of markers and chromosomes (and individuals). The aim is to get as many truly linked markers together, but avoid completely putting unlinked markers in the same linkage group. It is usually easier to combine linkage groups after the fact rather than to have to split linkage groups apart.”

mapthis <- formLinkageGroups(mapthis, max.rf=0.35, min.lod=6, reorgMarkers=TRUE)

Order markers

The function orderMarkers() may be used to establish an initial order. It picks an arbitrary pair of markers, and then adds an additional marker (chosen at random), one at a time, in the best supported position. With the argument use.ripple=TRUE, the function ripple() is used after the addition of each marker, to consider all possible orders in a sliding window of markers, to attempt to improve marker order. The argument window defines the length of the window. Larger values will explore more possible orders but will require considerably more computation time. The value window=7 is usually about the largest one would ever consider;

Here is an example for chr5 (group 5).

mapthis <- orderMarkers(mapthis, use.ripple=TRUE, window = 7, chr=5)
pull.map(mapthis, chr=5)
##    C3M16    C3M10     C3M7     C3M6     C3M4     C3M3     C3M2 
##  0.00000 36.33398 46.57328 53.45508 76.44823 78.73700 90.65764

switch alleles (specifically for heterozygous parents)

plotRF(mapthis, alternate.chrid=TRUE)

toswitch <- markernames(mapthis, chr=c(5, 7:11))
mapthis <- switchAlleles(mapthis, toswitch)
mapthis <- formLinkageGroups(mapthis, max.rf=0.35, min.lod=6, reorgMarkers=TRUE)
## Warning in formLinkageGroups(mapthis, max.rf = 0.35, min.lod = 6,
## reorgMarkers = TRUE): Running est.rf.
mapthis <- orderMarkers(mapthis, use.ripple=TRUE, window = 7, chr=5)

Look for problem individuals

plot(countXO(mapthis), ylab="Number of crossovers")

Estimate genotyping error rate, bugs contained

#loglik <- err <- c(0.0025, 0.005, 0.0075, 0.01, 0.0125, 0.015, 0.0175, 0.02)
#for(i in seq(along = err)) {
#  cat(i, "of", length(err), "\n")
#  tempmap <- est.map(mapthis, error.prob = err[i])
#  loglik[i] <- sum(sapply(tempmap, attr, "loglik"))
#}
#lod <- (loglik - max(loglik))/log(10)
### plot MLE
#plot(err, lod, xlab="Genotyping error rate", xlim=c(0,0.02), ylab=expression(paste(log[10], " likelihood")))

Look for genotyping errors and zero out these suspicious genotypes (that is, make them missing (NA))

mapthis <- calc.errorlod(mapthis, error.prob=0.005)  ### determine error LOD
toperr <- top.errorlod(mapthis, cutoff=4)
print(toperr)
##   chr    id marker errorlod
## 1   5 id113   C4M7 5.953700
## 2   5 id233   C4M9 4.949502
mapthis.clean <- mapthis
for(i in 1:nrow(toperr)) {
  chr <- toperr$chr[i]
  id <- toperr$id[i]
  mar <- toperr$marker[i]
  mapthis.clean$geno[[chr]]$data[mapthis$pheno$id==id, mar] <- NA
}

Summary and plot map

summaryMap(mapthis)
##         n.mar length ave.spacing max.spacing
## 1          34  330.0        10.0        10.0
## 2          25  240.0        10.0        10.0
## 3          16  150.0        10.0        10.0
## 4          10   90.0        10.0        10.0
## 5           9   60.2         7.5        20.8
## overall    94  870.2         9.8        20.8
plotRF(mapthis, alternate.chrid=TRUE)

par(las=1, mar=c(4.1,4.1,1.1,0.1), cex=0.3)
plotMap(mapthis, main="", show.marker.names=TRUE)