############################################################## # R code for "R/qtl tutorial" # # Karl W Broman, kbroman@biostat.wisc.edu # University of Wisconsin Madison # # http://www.rqtl.org # # 30 November 2012 ############################################################## ############################################################ # Preliminaries ############################################################ library(qtl) ls() help(read.cross) ?read.cross # url.show("http://www.rqtl.org/rqtltour3.R") ############################################################ # Data import ############################################################ sug <- read.cross("csv", "http://www.rqtl.org", "sug.csv", genotypes=c("CC", "CB", "BB"), alleles=c("C", "B")) ############################################################ # Summaries ############################################################ summary(sug) nind(sug) nchr(sug) totmar(sug) nmar(sug) nphe(sug) plot(sug) plotMissing(sug) plotMap(sug) plotPheno(sug, pheno.col=1) plotPheno(sug, pheno.col=2) plotPheno(sug, pheno.col=3) plotPheno(sug, pheno.col=4) plotPheno(sug, pheno.col=5) plotPheno(sug, pheno.col=6) plotPheno(sug, pheno.col="bp") plotPheno(sug, pheno.col="bw") ############################################################ # Single-QTL analysis ############################################################ sug <- calc.genoprob(sug, step=1) out.em <- scanone(sug) summary(out.em) summary(out.em, threshold=3) plot(out.em) out.hk <- scanone(sug, method="hk") plot(out.em, out.hk, col=c("blue", "red")) plot(out.em, col="blue") plot(out.hk, col="red", add=TRUE) plot(out.em, out.hk, col=c("blue", "red"), chr=c(7,15)) plot(out.em, col="blue", chr=c(7,15)) plot(out.hk, col="red", chr=c(7,15), add=TRUE) plot(out.hk - out.em, ylim=c(-0.3, 0.3), ylab="LOD(HK)-LOD(EM)") ############################################################ # Permutation tests ############################################################ load(url("http://www.rqtl.org/various.RData")) operm <- scanone(sug, method="hk", n.perm=1000) plot(operm) summary(operm) summary(operm, alpha=c(0.05, 0.2)) summary(out.hk, perms=operm, alpha=0.2, pvalues=TRUE) ############################################################ # Interval estimates of QTL location ############################################################ lodint(out.hk, chr=7) bayesint(out.hk, chr=7) lodint(out.hk, chr=7, expandtomarkers=TRUE) bayesint(out.hk, chr=7, expandtomarkers=TRUE) lodint(out.hk, chr=7, drop=2) bayesint(out.hk, chr=7, prob=0.99) lodint(out.hk, chr=15) bayesint(out.hk, chr=15) ############################################################ # QTL effects ############################################################ max(out.hk) mar <- find.marker(sug, chr=7, pos=47.7) plotPXG(sug, marker=mar) sug <- sim.geno(sug, n.draws=64, step=1) effectplot(sug, mname1=mar) effectplot(sug, mname1="7@47.7") max(out.hk, chr=15) mar2 <- find.marker(sug, chr=15, pos=12) plotPXG(sug, marker=mar2) effectplot(sug, mname1="15@12") plotPXG(sug, marker=c(mar, mar2)) plotPXG(sug, marker=c(mar2, mar)) effectplot(sug, mname1="7@47.7", mname2="15@12") effectplot(sug, mname2="7@47.7", mname1="15@12") ############################################################ # Other phenotypes ############################################################ out.hr <- scanone(sug, pheno.col=2, method="hk") out.bw <- scanone(sug, pheno.col="bw", method="hk") out.logbw <- scanone(sug, pheno.col=log(sug$pheno$bw), method="hk") out.all <- scanone(sug, pheno.col=1:4, method="hk") summary(out.all, threshold=3) summary(out.all, threshold=3, lodcolumn=4) summary(out.all, threshold=3, format="allpeaks") summary(out.all, threshold=3, format="allpheno") summary(out.all, threshold=3, format="tabByCol") summary(out.all, threshold=3, format="tabByChr") ############################################################ # Data diagnostics ############################################################ bad <- read.cross("csv", "http://www.rqtl.org", "bad.csv") summary(bad) par(mfrow=c(3,1)) for(i in 1:3) plotPheno(bad, pheno.col=i) par(mfrow=c(3,1)) for(i in 1:3) plot(bad$pheno[,i], ylab=names(bad$pheno)[i]) par(mfrow=c(1,1)) pairs(bad$pheno) par(mfrow=c(2,1)) plot(ntyped(bad, "ind"), main="No. genotypes, by ind'l", ylab="No. genotypes") plot(ntyped(bad, "mar"), main="No. genotypes, by marker", ylab="No. genotypes") bad <- subset(bad, ind=(ntyped(bad, "ind") > 100)) nt.mar <- ntyped(bad, "mar") bad <- drop.markers(bad, names(nt.mar)[nt.mar < 180]) cg <- comparegeno(bad) lowertri <- cg[lower.tri(cg)] par(mfrow=c(1,1)) hist(lowertri, breaks=seq(0, 1, by=0.01)) rug(lowertri) sort(lowertri, decreasing=TRUE)[1:10] wh <- which(!is.na(cg) & cg > 0.9, arr.ind=TRUE) bad <- subset(bad, ind = -wh[,2]) gt <- geno.table(bad) gt[gt$P.value < 0.05/nrow(gt),] gt2 <- geno.table(bad, scanone.output=TRUE) par(mfrow=c(2,1)) plot(gt2) plot(gt2, lod=3:5, ylab="Genotype frequency") abline(h=c(0.25, 0.5), lty=2, col="green") legend("bottomleft", colnames(gt2)[5:7], lwd=2, col=c("black","blue","red")) geno.table(bad, chr=3) geno.table(bad, chr=6) gt3 <- geno.table(bad, chr=-3) badmar <- rownames(gt3)[gt3$P.value < 0.05/nrow(gt3)] bad <- drop.markers(bad, badmar) bad <- est.rf(bad) par(mfrow=c(1,1)) plotRF(bad) plotRF(bad, chr=c(7,10)) markernames(bad, chr=7) markernames(bad, chr=7)[4] newmap <- est.map(bad) par(mfrow=c(1,1)) plotMap(bad, newmap) plotMap(bad, newmap, chr=7, show.marker.names=TRUE) out <- tryallpositions(bad, "D7M4", chr=10) summary(out) bad <- movemarker(bad, "D7M4", 10, 9.245) newmap <- est.map(bad) plotMap(bad, newmap) plotRF(bad, chr=18) rip <- ripple(bad, chr=18, window=7) summary(rip) compareorder(bad, chr=18, rip[2,]) mar <- markernames(bad, chr=18)[3] out <- tryallpositions(bad, mar, 18) summary(out) bad <- movemarker(bad, mar, 18, 53.7) newmap <- est.map(bad) plotMap(bad, newmap) xo <- countXO(bad) hist(xo, breaks=50) rug(xo) bad <- subset(bad, ind=(xo < 60)) bad <- calc.errorlod(bad, err=0.01) top.errorlod(bad) ############################################################ # Two-dimensional, two-QTL scans ############################################################ sug <- read.cross("csv", "http://www.rqtl.org", "sug.csv", genotypes=c("CC", "CB", "BB"), alleles=c("C", "B")) load(url("http://www.rqtl.org/various.RData")) sug <- calc.genoprob(sug, step=2) # out2 <- scantwo(sug, method="hk") plot(out2) plot(out2, lower="fv1") plot(out2, lower="fv1", upper="av1") # operm2 <- scantwo(sug, method="hk", n.perm=5) summary(out2, perms=operm2, alpha=0.2, pvalues=TRUE) ############################################################ # Multiple-QTL analyses ############################################################ sug <- calc.genoprob(sug, step=1) qtl <- makeqtl(sug, chr=c(7,15), pos=c(47.7, 12), what="prob") out.fq <- fitqtl(sug, qtl=qtl, method="hk") summary(out.fq) summary(fitqtl(sug, qtl=qtl, method="hk", get.ests=TRUE, dropone=FALSE)) out.fqi <- fitqtl(sug, qtl=qtl, method="hk", formula=y~Q1*Q2) out.fqi <- fitqtl(sug, qtl=qtl, method="hk", formula=y~Q1+Q2+Q1:Q2) summary(out.fqi) addint(sug, qtl=qtl, method="hk") rqtl <- refineqtl(sug, qtl=qtl, method="hk") rqtl summary(out.fqr <- fitqtl(sug, qtl=rqtl, method="hk")) plotLodProfile(rqtl) out.hk <- scanone(sug, method="hk") plot(out.hk, chr=c(7,15), col="red", add=TRUE) out.aq <- addqtl(sug, qtl=rqtl, method="hk") plot(out.aq) print(pen <- calc.penalties(operm2)) out.sq <- stepwiseqtl(sug, max.qtl=5, penalties=pen, method="hk", verbose=2) out.sq # end of rqtltour3.R