Econometric Sense

# ------------------------------------------------------------------
# |PROGRAM NAME: EX_QTL_R
# |DATE: 8-13-11
# |CREATED BY: MATT BOGARD 
# |PROJECT FILE: 
# |----------------------------------------------------------------
# | 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