%\VignetteIndexEntry{QTLBIM Prototype Slides Analyzing Hyper Data}
%\VignetteDepends{qtlbim}
%\VignetteKeywords{QTL}
%\VignettePackage{qtlbim}
\documentclass{beamer}

\usepackage{beamerthemesplit}

\usepackage{Sweave}

\begin{document}

<<echo=FALSE,label=Initialization>>=
## Initialization if not called from qb.sweave.
.qb.Package <- !exists(".qb.name")
if(.qb.Package) { ## Called by Sweave directly or in R build.
  .qb.name <- "hyper"
  require(qtl, quietly = TRUE)
  data(hyper)
  hyper <- subset(clean(hyper),
    chr = ("X" != unlist(lapply(hyper$geno, class))))
  .qb.cross <- hyper
  .qb.pheno <- 1
  .qb.niter <- 3000
  .qb.draws <- 8
  .qb.scan.type <- "2logBF"
  .qb.hpd.level <- 0.5
  .qb.threshold <- c(upper=2)
  .qb.remove <- TRUE
  .qb.SweaveFile <- system.file("doc", "hyperslide.Rnw", package="qtlbim")
  .qb.SweaveExtra <- system.file("external", "hyperslideextra.Rnw", package="qtlbim")
  .qb.PDFDir <- paste(names(.qb.cross$pheno)[.qb.pheno], "PDF", sep = "")
}
.qb.pheno.name <- names(.qb.cross$pheno)[.qb.pheno]
if(!file.exists(.qb.PDFDir)) {
  dir.create(.qb.PDFDir)
  warning(paste("Creating PDF directory", .qb.PDFDir),
    call. = FALSE, immediate. = TRUE)
}
## Make sure Sweave.sty is locally available.
invisible(file.copy(system.file("../../share/texmf","Sweave.sty"),"."))
## Assign visible names for script.
cross <- .qb.cross
pheno.col <- .qb.pheno
hpd.level <- .qb.hpd.level
scan.type <- .qb.scan.type
threshold <- .qb.threshold
n.iter <- .qb.niter
n.draws <- .qb.draws
remove.qb <- .qb.remove
if(!is.null(.qb.SweaveExtra)) {
  .qb.ExtraTex <- basename(.qb.SweaveExtra)
  .qb.ExtraTex <- substring(.qb.ExtraTex,1,nchar(.qb.ExtraTex)-4)
}
@


\SweaveOpts{eps=FALSE}

\title{Prototype QTL Strategy: Phenotype
\Sexpr{.qb.pheno.name}
in Cross
\Sexpr{.qb.name}}
\author{Brian S. Yandell, W. Whipple Neely, Nengjun Yi}
\date{\today}

\frame{\titlepage}

\section[Outline]{}
\frame{\tableofcontents}

\section{Overview}

\frame
{
  \frametitle{Automated Strategy}

  \begin{itemize}
  \item Estimate positions and effects of main QTL.
  \item Find chromosomes with epistasis.
  \item Estimate epistatic pair positions and effects.
  \item Confirm genetic architecture with ANOVA.
  \end{itemize}
}

\begin{frame}[fragile]
  \frametitle{Running Sweave}


\tiny

<<echo=TRUE,label=LoadQtlbim:CallSweave>>=
library(qtlbim)
@
\begin{Schunk}
\begin{Sinput}
> qb.sweave(\Sexpr{.qb.name}, pheno.col = \Sexpr{.qb.pheno},
+  n.iter = \Sexpr{.qb.niter}, n.draws = \Sexpr{.qb.draws},
+  scan.type = "\Sexpr{.qb.scan.type}", hpd.level = \Sexpr{.qb.hpd.level},
+  threshold = c(\Sexpr{paste(names(.qb.threshold), "=", .qb.threshold, collapse = ", ")}),
+  SweaveFile = "\Sexpr{.qb.SweaveFile}",
+  SweaveExtra = "\Sexpr{.qb.SweaveExtra}",
+  PDFDir = "\Sexpr{.qb.PDFDir}",
+  remove.qb = \Sexpr{.qb.remove})
\end{Sinput}
\end{Schunk}

\end{frame}

\subsection{Initialization}

\begin{frame}[fragile]
  \frametitle{Cross Object}

\tiny

<<echo=TRUE,label=SummaryCrossObject>>=
summary(cross)
@

\end{frame}

\begin{frame}[fragile]
  \frametitle{Create MCMC runs}

\tiny

\begin{Schunk}
\begin{Sinput}
> cross <- qb.genoprob(cross,step=2)
> cross.qb <- qb.mcmc(cross, pheno.col = pheno.col,
+   genoupdate=TRUE, n.iter = \Sexpr{n.iter}, verbose=FALSE)
\end{Sinput}
\end{Schunk}
<<echo=FALSE,label=MCMCSamples>>=
if(.qb.Package) {
  qb.load(hyper, qbHyper)
  cross <- hyper
  cross.qb <- qbHyper
  cross.qb$cross.name <- "cross"
  rm(hyper, qbHyper)
} else {
  cross <- qb.genoprob(cross,step=2)
}

## Create cross.qb if it does not exist.
if(exists("cross.qb")) {
  remove.qb <- FALSE
} else {
  cross.qb <- qb.mcmc(cross, pheno.col = pheno.col, genoupdate=TRUE,
    n.iter = n.iter, verbose=FALSE)
}
@

\end{frame}

\section{1-D \& 2-D Scans}

\begin{frame}[fragile]
  \frametitle{1-D \Sexpr{.qb.scan.type} Scan}

\tiny

<<echo=TRUE,label=HPDSummary>>=
hpd.level
scan.type
cross.hpd <- qb.hpdone(cross.qb, hpd.level, scan.type)
sum.one <- summary(cross.hpd)
sum.one
chrs <- as.vector(sum.one[, "chr"])
pos <- sum.one[, "pos"]
@
\begin{Schunk}
\begin{Sinput}
> plot(cross.hpd)
\end{Sinput}
\end{Schunk}

\end{frame}

\begin{frame}[fragile]
  \frametitle{1-D Scan: \Sexpr{.qb.scan.type} Profile}

<<results=tex,echo=FALSE,label=HPDPlot>>=
file <- paste(.qb.PDFDir, "/slide1hpd.pdf", sep = "")
pdf(file = file, paper = "special", width = 9, height = 6)
plot(cross.hpd)
invisible(dev.off())
cat("\\includegraphics{", file, "}\n\n", sep = "")
@

\end{frame}

\begin{frame}[fragile]
  \frametitle{2-D: find epistatic pairs}

\tiny

<<echo=TRUE,label=ScanTwoSummary>>=
two <- qb.scantwo(cross.qb, chr = chrs, type = scan.type)
sum.two <- summary(two,sort="upper",threshold=threshold,
  refine = TRUE)
sum.two
@

\end{frame}

\begin{frame}[fragile]
  \frametitle{Initial Genetic Architecture}

\tiny

<<echo=TRUE,label=InitialArchitecture>>=
cross.arch <- qb.arch(sum.two, chrs, pos)
cross.arch
@

\end{frame}
 
\section{Anova Fit}

\begin{frame}[fragile]
  \frametitle{Construct QTL Object}

\tiny

use R/qtl tools to check model fit\\
first simulate missing markers\\
then construct QTL object

<<echo=TRUE,label=InitializeFitQTL>>=
cross.sub <- subset(cross, chr = unique(cross.arch$qtl$chr))
n.draws
cross.sub <- sim.geno(cross.sub, n.draws=n.draws, step=2, error=0.01)
qtl <- makeqtl(cross.sub, cross.arch$qtl$chr, cross.arch$qtl$pos)
cross.sub <- clean(cross.sub)
@
\end{frame}

\begin{frame}[fragile]
  \frametitle{Stepwise Reduction}

\tiny

<<echo=TRUE,label=StepwiseFitQTL>>=
cross.step <- step.fitqtl(cross.sub, qtl, pheno.col, cross.arch)
@
\begin{Schunk}
\begin{Sinput}
> summary(cross.step$fit)
\end{Sinput}
\end{Schunk}
<<echo=FALSE,label=SummaryFitQTL>>=
sum.fit <- summary(cross.step$fit)
print(sum.fit$result.full, quote = FALSE, na.print = "")
@
\end{frame}

<<results=tex,echo=FALSE>>=
if(!is.null(sum.fit$result.drop)) {
  cat("\\begin{frame}[fragile]\n")
  cat("\\frametitle{Stepwise Reduction}\n\n")
  cat("\\tiny\n\n")
  cat("\\begin{Schunk}\n")
  cat("\\begin{Soutput}\n")
  printCoefmat(sum.fit$result.drop[,-6], digits = 4, cs.ind = 1, P.values = TRUE, 
               has.Pvalue = TRUE, signif.legend = FALSE)
  cat("\\end{Soutput}\n\n")
  cat("\\end{Schunk}\n")
  cat("\\end{frame}\n\n")
}
@

\begin{frame}[fragile]
  \frametitle{Reduced Genetic architecture}

\tiny

<<echo=TRUE,label=FinalArchitecture>>=
cross.arch <- cross.step$arch
cross.arch
@
\end{frame}

\begin{frame}[fragile]
  \frametitle{2-D Plots}

2-D plots by cliques (if any epistasis)

\tiny

\begin{Schunk}
\begin{Sinput}
> for(i in names(cross.arch$chr.by.set))
+   plot(two, chr = cross.arch$chr.by.set[[i]], smooth = 3,
+     col = "gray", contour = 3)
\end{Sinput}
\end{Schunk}

\end{frame}

<<results=tex,echo=FALSE,label=ScanTwoPlotByGroup>>=
## Note extra R overhead to produce an arbitrary number of plots.
if(!is.null(cross.arch$chr.by.set)) {
  for(i in names(cross.arch$chr.by.set)) {
    file <- paste(.qb.PDFDir, "/slide2LOD-", i, ".pdf", sep = "")
    pdf(file = file, paper = "special", width = 8, height = 6)
    plot(two, chr = cross.arch$chr.by.set[[i]], smooth = 3,
      col = "gray", contour = 3)
    invisible(dev.off())
    cat("\\begin{frame}[fragile]\n")
    cat("\\frametitle{2-D Plots: clique", i, "}\n\n")
    cat("\\includegraphics{", file, "}\n\n", sep = "")
    cat("\n\\end{frame}\n\n")
    warning(paste("writing", file), call. = FALSE, immediate. = TRUE)
  }
  cat("\n\n")
}

@

\begin{frame}[fragile]
  \frametitle{Slice Each Epistatic Pair}

show detail plots for epistatic pairs (if any)

\tiny

\begin{Schunk}
\begin{Sinput}
> if(!is.null(cross.arch$pair.by.chr)) {
+  for(i in seq(nrow(cross.arch$pair.by.chr$chr))) {
+    chri <- cross.arch$pair.by.chr$chr[i,]
+    posi <- cross.arch$pair.by.chr$pos[i,]
+    if(chri[1] != chri[2])
+      plot(qb.slicetwo(cross.qb, chri, posi, scan.type))
+  }
+}
\end{Sinput}
\end{Schunk}

\end{frame}

<<echo=FALSE>>=
if(!is.null(cross.arch$pair.by.chr))
  warning(paste("creating", nrow(cross.arch$pair.by.chr$chr), "epistatic pair plots"),
    call. = FALSE, immediate. = TRUE)
@

<<results=tex,echo=FALSE,label=PairPlots>>=
## Note extra R overhead to produce an arbitrary number of plots.
if(!is.null(cross.arch$pair.by.chr)) {
  cross <- sim.geno(cross, step = qb.get(cross.qb, "step"))
  for(i in seq(nrow(cross.arch$pair.by.chr$chr))) {
    chri <- cross.arch$pair.by.chr$chr[i,]
    posi <- cross.arch$pair.by.chr$pos[i,]
    if(chri[1] != chri[2]) {
      file <- paste(.qb.PDFDir, "/slide-", chri[[1]], "-", chri[[2]], ".pdf",
        sep = "")
      pdf(file = file, paper = "special", width = 9, height = 6)
      tmp <- qb.slicetwo(cross.qb, chri, posi, scan.type)
      plot(tmp)
      invisible(dev.off())
      cat("\\begin{frame}[fragile]\n")
      cat("\\frametitle{Epistatic Pair", chri[[1]], "and", chri[[2]], "}\n\n")
      cat("\\includegraphics{", file, "}\n\n", sep = "")
      cat("\n\\end{frame}\n\n")
      warning(paste("writing", file), call. = FALSE, immediate. = TRUE)
    }
  }
  cat("\n\n")
}
@

<<results=tex,echo=FALSE,label=UserExtraSweave>>=
if(!is.null(.qb.SweaveExtra)) {
  warning(paste("Running Sweave on Extra to create ",
      .qb.ExtraTex, ".tex", sep = ""),
    call. = FALSE, immediate. = TRUE)
  Sweave(.qb.SweaveExtra, quiet = TRUE)
  cat("\n\\input{", .qb.ExtraTex, "}\n\n", sep = "")
}
@

\section{Conclusion}

\begin{frame}[fragile]

<<echo=FALSE>>=
.qb.sweave.tex <- basename(.qb.SweaveFile)
.qb.sweave.tex <- paste(substring(.qb.sweave.tex, 1,
   nchar(.qb.sweave.tex) - 4), "tex", sep = ".")
.qb.pheno.tex <- paste(.qb.pheno.name, "tex", sep = ".")
@
final tasks:\\
externally rename file \Sexpr{.qb.sweave.tex} to \Sexpr{.qb.pheno.tex}\\
and run pdflatex twice on it\\
remove objects created by {\tt R/qtlbim} if desired

\tiny

\begin{Schunk}
\begin{Sinput}
> file.rename("\Sexpr{.qb.sweave.tex}", "\Sexpr{.qb.pheno.tex}")
> invisible(system("pdflatex \Sexpr{.qb.pheno.tex}",intern=TRUE))
> invisible(system("pdflatex \Sexpr{.qb.pheno.tex}",intern=TRUE))
\end{Sinput}
\end{Schunk}
<<echo=TRUE,label=RemoveObjects>>=
remove.qb
if(remove.qb) {
  qb.remove(cross.qb)
  rm(cross, cross.sub, pheno.col, threshold, n.iter, n.draws, remove.qb)
}
@

<<echo=FALSE>>=
remove(list=objects(all.names = TRUE, pattern="^\\.qb\\..*"),
       pos=".GlobalEnv")
@

\end{frame}

\end{document}
