QTL Analysis in R

ex qtl.r

#   ------------------------------------------------------------------
#  |PROGRAM NAME: EX_QTL_R
#  |DATE: 8-13-11
#  |CREATED BY: MATT BOGARD 
#  |PROJECT FILE: http://econometricsense.blogspot.com/2011/08/qtl-analysis-in-r.html             
#  |----------------------------------------------------------------
#  | PURPOSE: DEMONSTRATE STATISTICAL METHODS FOR QTL ANALYSIS IN R        
#  | 
#  |------------------------------------------------------------------
#  | REFERENCE:  R code for "A brief tour of R/qtl"
#  |  Karl W Broman, kbroman.wisc.edu http://www.rqtl.org
#  |-----------------------------------------------------------------
 
 
library(qtl) # call package qtl
 
ls() 
 
############################################################
#  exploring backcross data
############################################################
 
 
data(fake.bc) # load simulated backcross data (default data in qtl package)
ls()
 
summary(fake.bc) # summary of info in fake.bc which is and object of class 'cross'
 
# you can also get this information with specific function calls:
nind(fake.bc)   # number of individuals
nphe(fake.bc)   # number of phenotypes
nchr(fake.bc)   # number of chromosomes
totmar(fake.bc) # number of total  markers
nmar(fake.bc)   # list markers? 
 
############################################################
#   plotting maps, genotypes, phenotypes, and markers
###########################################################
 
 
plot(fake.bc) # gives plots data
 
# you can also call for the plots individually:
plot.missing(fake.bc) # just plot missing genotypes
plot.map(fake.bc)     # just plot the genetic map 
plot.pheno(fake.bc, pheno.col=1) # just plot phenotype 1 'phe 1' 
plot.map(fake.bc, chr=c(1, 4, 6, 7, 15), show.marker.names=TRUE) # specific chromosomes and marker names
plot.missing(fake.bc, reorder=TRUE) # order NA genotypes based on value of phenotype
 
fake.bc <- drop.nullmarkers(fake.bc)  # drop obs with missing genotypes
totmar(fake.bc) # total # markers left
 
 
##################################################################
#  specifying and estimating the likelihood used for QTL mapping
#################################################################
 
# From Broman:
# The function calc.genoprob calculates QTL genotype probabilities, conditional
# on the available marker data. These are needed for most of the QTL mapping 
# functions. The argument step indicates the step size (in cM) at which the
# probabilities are calculated, and determines the step size at which later 
# LOD scores are calculated.
 
 
fake.bc <- calc.genoprob(fake.bc, step=1, error.prob=0.01)
 
# function scanone performs single-QTL genome scan with a normal model. 
#  methods include: maximum likelihood via the EM algorithm 
#  and Haley-Knott regression 
 
out.em <- scanone(fake.bc)
out.hk <- scanone(fake.bc, method="hk")
 
# multiple imputation method using sim.geno utilizing the joint 
# genotype distribution, given the observed marker data.
 
fake.bc <- sim.geno(fake.bc, step=2, n.draws=16, error.prob=0.01)
out.imp <- scanone(fake.bc, method="imp")
 
# get the maximum LOD score on each chromosome 
# can also specify a threshold for LOD
 
summary(out.em) 
summary(out.em, threshold=3)
summary(out.hk, threshold=3)
summary(out.imp, threshold=3)
 
# function max.scanone returns just the highest peak from output of scanone.
max(out.em) # based on expectation maximization
max(out.hk) # based on Haley-Knott regression
max(out.imp) # based on multiple imputation
 
##################################################################
#  plot LOD scores by chromosome for QTL mapping
#################################################################
 
 
plot(out.em, chr=c(2,5)) # just based on em method
plot(out.em, out.hk, out.imp, chr=c(2,5)) # all methods
plot(out.em, chr=c(2)) # zoom in on specified chromosome