Introductory comments

Description

A brief introduction to the R/qtl package.

New to R?

• Be sure to read the R documention and/or tutorials available from the R archive (http://cran.r-project.org).
• In order to use the R/qtl package, you must type (within R) library(qtl). You may wish to include this in a .First function or .Rprofile file.
• Use the help.start function to start the html version of the R help.
• Type library(help=qtl) to get a list of the functions in R/qtl.
• Use the example function to run examples of the various functions in R/qtl.

Walk-through of an analysis

You can use the library function to list available packages, list functions and data within a package, and load a package.

library() # List available packages
library(help=qtl) # List functions and data within "qtl"
library(qtl) # Load the qtl package

Use ? or help to get help about a function or dataset, or use help.start to start the html version of R help.

?calc.genoprob
help(calc.genoprob)
help.start()

Use the data function to load one of the example data sets provided with R/qtl, hyper.

data(hyper)

The hyper data set has class cross. The function plot.cross calls plot.missing (to plot a matrix displaying missing genotypes) and plot.map (to plot the genetic maps) and also displays histograms of all phenotypes. The plot.missing function can plot individuals ordered by their phenotypes; you can see that for most markers, only individuals with extreme phenotypes were genotyped.

plot(hyper)
plot.missing(hyper)
plot.missing(hyper,reorder=TRUE)
plot.map(hyper)

To re-estimate the genetic map for an experimental cross, use the function est.map. The function plot.map, in addition to plotting a single map, can plot the comparison of two genetic maps (as long as they are composed of the same numbers of chromosomes and markers per chromosome). The function replace.map map be used to replace the genetic map in a cross with a new one.

newmap <- est.map(hyper,error.prob=0.01,print.rf=TRUE)
plot.map(hyper,newmap)
hyper <- replace.map(hyper,newmap)

Before doing QTL analyses, a number of intermediate calculations may need to be performed. The function calc.genoprob calculates conditional genotype probabilities given the multipoint marker data. argmax.geno calculates the most likely sequence of underlying genotypes, given the observed marker data. sim.geno simulates sequences of genotypes from their joint distribution, given the observed marker data.

These three functions return a modified version of the input cross, with the intermediate calculations included.

hyper <- calc.genoprob(hyper,step=1,off.end=5,error.prob=0.01)
hyper <- argmax.geno(hyper,error.prob=0.01)
hyper <- sim.geno(hyper,step=1,off.end=5,n.draws=16,error.prob=0.01)

The function calc.errorlod may be used to assist in identifying possibly genotyping errors; it calculates the Lincoln and Lander error lod scores.

hyper <- calc.errorlod(hyper)
plot.errorlod(hyper)
top.errorlod(hyper,cutoff=3)

The function est.rf estimates the recombination fraction between each pair of markers, and calculates a LOD score testing r = 0.5. This is useful for identifying markers that are placed on the wrong chromosome. Note that since, for these data, many markers were typed only on recombinant individuals, the pairwise recombination fractions show rather odd patterns.

hyper <- est.rf(hyper)
plot.rf(hyper)
plot.rf(hyper,c(1,5,13))

The function scanone performs a genome scan with a single QTL model. With method="anova", analysis of variance at each marker is performed; any individuals with missing genotypes are discarded. With method="im", Lander and Botstein's interval mapping is perfomed, while with method="hk", the Haley-Knott approximation to interval mapping is performed. These two methods require the results from calc.genoprob.

out.anova <- scanone(hyper,method="anova")
out.im <- scanone(hyper,method="im") out.hk <- scanone(hyper,method="hk")

The output of scanone is a matrix with class scanone. The function plot.scanone may be used to plot the results, and may plot three sets of results against each other, as long as they conform appropriately.

plot(out.im)
plot(out.im,out.hk,c(1,4),lty=1,col=c("black","blue","red"))

The function summary.scanone may be used to list information on the peak LOD for each chromosome

summary(out.im,1.7)

Note

Also see the package BQTL (Bayesian QTL mapping toolkit) by Charles C. Berry. His package gave me the idea to have this "A starting point" help page.

Author(s)

Karl W Broman, kbroman@jhsph.edu
http://biosun01.biostat.jhsph.edu/~kbroman/software/qtl.html