###################################################################### # Code for "A guide to QTL mapping with R/qtl" # Karl W Broman and Saunak Sen # # Chapter 10: Case study I ####################################################################### ###################################################################### # 10.1 Diagnostics ###################################################################### # load the data library(qtl) library(qtlbook) data(ovar) # summary of the data summary(ovar) # summary plot of the data plot(ovar) # scatterplot of on1 and on2 plot(jitter(on2) ~ jitter(on1), data=ovar$pheno, xlab="on1", ylab="on2", cex=0.6, xlim=c(6.66, 17.53), ylim=c(6.66, 17.53)) # correlation between on1 and on2 cor(pull.pheno(ovar,"on1"), pull.pheno(ovar,"on2"), use="complete") # check that onm = (on1+on2)/2 max(abs(pull.pheno(ovar,"onm") - (pull.pheno(ovar,"on1")+pull.pheno(ovar,"on2"))/2), na.rm=TRUE) # check whether onm is missing when on1 or on2 are table(apply(is.na(pull.pheno(ovar,1:3)), 1, paste, collapse=":")) # individuals with onm but not on1 ovar$pheno[!is.na(ovar$pheno$onm) & is.na(ovar$pheno$on1),] # Make onm phenotype missing when on1 is missing ovar$pheno$onm[is.na(ovar$pheno$on1)] <- NA # box plot of ovariole counts in the two crosses boxplot(onm ~ cross, data=ovar$pheno, horizontal=TRUE, xlab="Ovariole count", ylab="Cross") # t-test for difference in average phenotype between the two crosses t.test(onm ~ cross, data=ovar$pheno) # pull out first cross ovar1 <- subset(ovar, ind=(pull.pheno(ovar, "cross")==1)) # plot phenotype against number of typed genotypes plot(jitter(ntyped(ovar1)),jitter(pull.pheno(ovar1,"onm")), xlab="No. genotypes", ylab="Ovariole count") # plot pairwise recombination fractions ovar1 <- est.rf(ovar1) plot.rf(ovar1) # reestimate map and plot it against that in data newmap <- est.map(ovar1, error.prob=0.001, verbose=FALSE) plot.map(ovar1, newmap) # estimate map with Kosambi map function newmap.k <- est.map(ovar1, err=0.001, map.function="kosambi") # summary of the two maps summary(newmap) summary(newmap.k) # two-locus genotypes for markers st and e in 2nd cross ovar2 <- subset(ovar, ind=(pull.pheno(ovar,"cross")==2)) geno.crosstab(ovar2, "st", "e") # two-locus genotypes for markers st and e in 1st cross geno.crosstab(ovar1, "st", "e") # plot chr 3 genotypes for a random set of individuals from the 2nd cross toplot <- sort(sample(nind(ovar2), 15)) plot.geno(ovar2, chr=3, ind=toplot) ###################################################################### # 10.2 Initial cross ###################################################################### # Omit individuals with missing genotype data ovar1 <- subset(ovar1, ind=!is.na(pull.pheno(ovar1,"onm"))) # Multiple imputations ovar1 <- sim.geno(ovar1, n.draws=512, step=2, err=0.001, map.function="kosambi") # genome scan by multiple imputation out1 <- scanone(ovar1, method="imp") # summary of genome scan summary(out1) # plot of genome scan plot(out1, ylab="LOD score") # effect plot for chr 3 effectplot(ovar1, mname1="cpo") # fitqtl to estimate QTL effect qtl <- makeqtl(ovar1, chr=3, pos=82.4) summary(fitqtl(ovar1, qtl=qtl, dropone=FALSE, get.ests=TRUE)) # scan adjusting for chr 3 locus out1.c3 <- addqtl(ovar1, qtl=qtl) # summary of genome scan, adjusting for chr 3 locus summary(out1.c3) # plot of genome scan, adjusting for chr 3 locus plot(out1.c3, ylab="LOD score") # 2d scan on chromosome 3 qtl2 <- makeqtl(ovar1, 2, 90) out1.ap <- addpair(ovar1, qtl=qtl2, chr=3, verbose=FALSE) # summary of the 2d scan summary(out1.ap) # refine QTL locations in the 3-QTL model qtl <- makeqtl(ovar1, c(2, 3, 3), c(90, 82.8, 121)) rqtl <- refineqtl(ovar1, qtl=qtl, verbose=FALSE) # print the results of refineqtl rqtl # plot LOD profiles for the 3-QTL model plotLodProfile(rqtl, col=c("black","red","blue"), ylab="Profile LOD score") # approximate confidence intervals for the location of the proximal c3 QTL lodint(rqtl, qtl.index=2) bayesint(rqtl, qtl.index=2) # permutation test with 2-dim, 2-QTL genome scan strat <- (nmissing(ovar1) < 15) operm <- scantwo(ovar1, method="imp", n.perm=1000, perm.strata=strat) # Significance thresholds summary(operm) # Calculate penalties print(pen <- calc.penalties(operm)) # stepwise selection starting at a 3-qtl model stepout1 <- stepwiseqtl(ovar1, penalties=pen, qtl=rqtl, max.qtl=8, verbose=FALSE) # Print stepwiseqtl results stepout1 # Plots to study the inferred pair of QTL on chr 2 mar <- find.marker(ovar1, 2, c(92.4, 94)) par(mfrow=c(1,2)) effectplot(ovar1, mname1="2@92.4", mname2="2@94") plot.pxg(ovar1, marker=mar) # drop a QTL, refine locations, do dropone analysis qtl <- dropfromqtl(stepout1, 2) qtl <- refineqtl(ovar1, qtl=qtl, verbose=FALSE) # Print the result qtl ###################################################################### # 10.3 Combined data ###################################################################### # clean up memory object.size(ovar1) object.size(ovar1)/1024^2 ovar1 <- clean(ovar1) object.size(ovar1) object.size(ovar1)/1024 # remove individuals with missing phenotype ovar <- subset(ovar, ind=!is.na(pull.pheno(ovar,"onm"))) # Multiple imputations ovar <- sim.geno(ovar, n.draws=512, step=2, err=0.001, map.function="kosambi") # genome scan with combined data cross <- pull.pheno(ovar, "cross") outc <- scanone(ovar, method="imp", addcovar=cross) # genome scan with just the second cross out2 <- scanone(subset(ovar, ind=(cross==2)), method="imp") # plot genome scan for combined data and the two individual crosses plot(outc, out1, out2, ylab="LOD score") # location of maximum LOD in combined data max(outc) cross <- data.frame(cross=cross) qtlc <- makeqtl(ovar, chr=3, pos=72.9) outc.c3 <- addqtl(ovar, qtl=qtlc, covar=cross) qtl2 <- makeqtl(subset(ovar, ind=(cross==2)), chr=3, pos=72.9) out2.c3 <- addqtl(subset(ovar, ind=(cross==2)), qtl=qtl2) # plot genome scan for combined data and the two individual crosses plot(outc.c3, out1.c3, out2.c3, ylab="LOD score") # create vector indicating the 3 strata strat <- pull.pheno(ovar, "cross") strat[strat==1 & nmissing(ovar)<15] <- 3 table(strat) # the permutation test opermc <- scantwo(ovar, method="imp", addcovar=cross, n.perm=n.perm, perm.strata=strat) # significance thresholds summary(opermc) print(penc <- calc.penalties(opermc)) # 2d scan on chr 3, controlling for chr 2 locus qtl.c2 <- makeqtl(ovar, 2, 115) out.ap.c3 <- addpair(ovar, qtl=qtl.c2, chr=3, covar=cross, verbose=FALSE) # summary of the results summary(out.ap.c3) # plot of the results plot(out.ap.c3, lower="cond-int") # 2d scan on chr 2, controlling for the chr 3 locus out.ap.c2 <- addpair(ovar, qtl=qtlc, chr=2, covar=cross, verbose=FALSE) # summary of the results summary(out.ap.c2, perms=opermc, pval=TRUE) # plot of the results plot(out.ap.c2, lower="cond-int", upper="cond-add") # refine qtl locations qtl <- makeqtl(ovar, c(2, 2, 3, 3), c(4, 114, 62.8, 74.8)) rqtl <- refineqtl(ovar, qtl=qtl, covar=cross, verbose=FALSE, formula=y~cross+Q1+Q2+Q3+Q4+Q3:Q4) # print result rqtl # drop-one analysis with fitqtl summary(fitqtl(ovar, qtl=rqtl, covar=cross, formula=y~cross+Q1+Q2+Q3+Q4+Q3:Q4), pvalues=FALSE) # plot lod profiles plotLodProfile(rqtl, col=c("red","blue","red","blue"), ylab="Profile LOD score") addint(ovar, qtl=rqtl, covar=cross, formula=y~cross+Q1+Q2+Q3+Q4+Q3:Q4, pvalues=FALSE) # scan for one more QTL onemore <- addqtl(ovar, qtl=rqtl, covar=cross, formula=y~cross+Q1+Q2+Q3+Q4+Q3:Q4) # maximum LOD on each chromosome summary(onemore) qtl2 <- addtoqtl(ovar, rqtl, 2, 89.1) qtl2 <- reorderqtl(qtl2) rqtl2 <- refineqtl(ovar, qtl=qtl2, covar=cross, verbose=FALSE, formula=y~cross+Q1+Q2+Q3+Q4+Q5+Q4:Q5) # print result rqtl2 summary(fitqtl(ovar, qtl=rqtl2, covar=cross, formula=y~cross+Q1+Q2+Q3+Q4+Q5+Q4:Q5), pvalues=FALSE) # stepwise QTL analysis stepout1 <- stepwiseqtl(ovar, covar=cross, pen=pen, max.qtl=8, verbose=FALSE) # print results stepout1 # stepwise QTL analysis with exclusively heavy penalties stepout2 <- stepwiseqtl(ovar, covar=cross, pen=pen[1:2], max.qtl=8, verbose=FALSE) # print results stepout2 ###################################################################### # 10.4 Discussion ###################################################################### # end of chap10.R