Sample R file (hyper.R)

hyper.R

# R script for analyzing QTL data
# Bob Gotwals
# hyper.qtl script
# October 31, 2012
#
# clean things up
rm(list=ls())
setwd("~/Desktop/hyper")
#
# load the QTL library
library(qtl)
#
# Load the data!
cross <- read.cross("csv", file="hyper.csv", genotypes = c("A", "H", "B"), na.strings="-", alleles = c("A", "B"))
jittermap(cross)
summary(cross)
names(cross$pheno)
# take a look at my data, make sure it's pretty clean
cross <- est.rf(cross)
plot.rf(cross)
plot.map(cross)
plot.missing(cross)
hist(cross$pheno$bp)
bp <- cross$pheno
qqnorm(cross$pheno$bp)
qqline(cross$pheno$bp, col="red")
#
cross <- calc.genoprob(cross, step=2.0, off.end=0.0, error.prob=1.0e-4, map.function= "haldane", stepwidth ="fixed")
#
# Perform the mainscan for the QTL
cross.scanBP <- scanone(cross, pheno.col=1, model = "normal", method="em")
cross.scanBP.perm <- scanone(cross, pheno.col=1, model="normal", method="em", n.perm=100)
#
# plot the mainscan
thresh <- summary(cross.scanBP.perm, alpha = c(0.63, 0.1, 0.05))

plot(cross.scanBP, main="Mainscan plot of BP")
abline(h=thresh[1], col="blue")
abline(h=thresh[2], col="red")
abline(h=thresh[3], col="green")
#
summary(cross.scanBP, perm=cross.scanBP.perm, lodcolumn=1, alpha=0.05)
#
# do an effect plot
mname1 <- find.marker(cross, chr=1, pos=49.2)
effectplot(cross, pheno.col=1, mname1=mname1)
#
mname2 <- find.marker(cross, chr=4, pos=29.5)
effectplot(cross, pheno.col=1, mname1=mname2)

CIchr1 <- bayesint(cross.scanBP, chr=1, prob=0.95)
plot(cross.scanBP, chr=1, lodcolumn=1, main="Confidence interval for Chr1")
lines(x=CIchr1[c(1,3),2], y= c(0,0), type = "l", col="green", lwd=4)