\documentclass{article}

\usepackage{natbib}
\usepackage{graphics}
\usepackage{amsmath}
%\usepackage{indentfirst}
%\usepackage[utf8]{inputenc}

% \VignetteIndexEntry{LDheatmap adding tracks example}

\begin{document}

% Copied/modified next bit from Charlie Geyer's MCMC example vignette
<<foo,include=FALSE,echo=FALSE>>=
options(keep.source = TRUE, width = 60)
foo <- packageDescription("LDheatmap")
@

\title{LDheatmap (Version \Sexpr{foo$Version}): Example of Adding Tracks}
\author{Jinko Graham and Brad McNeney}
\maketitle

\section{Introduction}
As of version 0.9, \texttt{LDheatmap} allows users to ``flip'' the 
heatmap below a horizontal line in the style of Haploview. This 
extension was made so that ``tracks'' of
genomic annotation information could be placed above the genetic 
map. 
Facilities for adding tracks are illustrated in this vignette.
Limitations of the current implementation are discussed briefly in
Section~\ref{sec:limitations}

\section{Getting Started}

First load the package, the snpStats package
(to use its \texttt{SnpMatrix} data structure),
the example dataset \texttt{CIMAP5.CEU},
and an R data file that contains two
LDheatmap objects and one picture that will be used in the vignette.
<<load>>=
library(LDheatmap)
library(snpStats)
library(grid)
data(GIMAP5.CEU)
load(system.file("extdata/addTracks.RData",package="LDheatmap"))
@
The object \texttt{GIMAP5.CEU} is 
is a list with two elements: \texttt{snp.data} and \texttt{snp.support}.
\texttt{snp.data} is a \text{SnpMatrix} object from the
\texttt{snpStats} package, containing the SNP genotypes. 
Rows correspond to subjects and columns correspond to SNPs. 
\texttt{snp.support} is a data frame whose rows correspond to SNPs and
whose columns give information on the SNPs, such as their alleles and
genomic location. The help file \texttt{help("GIMAP5.CEU")} gives full
details.
In addition to \texttt{GIMAP5.CEU}, you should have the LDheatmap objects 
\texttt{llGenes} and
\texttt{llGenesRecomb}
in your workspace. These objects are the heatmap with tracks for
genes and recombination rates added incrementally. 
Adding these tracks requires
fetching information from the UCSC genome browser, which 
requires an internet connection and can be time-consuming.
These objects are provided so that users can see the track-annotated heatmaps
without having to wait.
The object \texttt{GIMAP5ideo} should also have been loaded. This 
is an array that represents a bitmap image of a chromosome ideogram
used in Section~\ref{sec:addGrob} of the vignette.

The basic LDheatmap for the \texttt{GIMAP5.CEU}
data to which tracks will be added can be created with the following 
R commands:
<<label=fig1com,include=FALSE>>= 
ll<-LDheatmap(GIMAP5.CEU$snp.data,GIMAP5.CEU$snp.support$Position,flip=TRUE)
@
% Above Sweave block shows the R commands that make the plot, but not the
% plot. Plot is shown by the Sweave block below. Doing the R commands and
% figure it outputs in one Sweave block makes the commands part of the Figure,
% which is hard for the reader to make sense of.
Figure~\ref{fig:fig1} displays the results.
\begin{figure}
\begin{center}
<<label=fig1,fig=TRUE,echo=FALSE>>= % do not echo commands -- already did above
<<fig1com>>
@
\end{center}
\caption{LDheatmap of GIMAP5.CEU data.}
\label{fig:fig1}
\end{figure}

\section{Adding tracks from the UCSC genome browser}

We will add ``UCSC Gene'' and ``Recomb Rate'' tracks 
to Figure~\ref{fig:fig1}. Data for these tracks is downloaded from the
UCSC genome browser using the \texttt{rtracklayer} package.
The genes track is added by the \texttt{LDheatmap.addGenes} function:
<<eval=FALSE>>=
llGenes <- LDheatmap.addGenes(ll, chr="chr7", genome="hg18")
@
It takes several minutes for this function to fetch the necessary data 
from the UCSC genome browser. 
To avoid this download time, users can draw the graphical object (grob) 
stored in \texttt{llGenes} directly (Figure~\ref{fig:fig2}) with 
the commands:
<<label=addGenes,eval=FALSE>>=
grid.newpage()
grid.draw(llGenes$LDheatmapGrob)
@
\begin{figure}
\begin{center}
<<label=fig2,fig=TRUE,echo=FALSE>>=
<<addGenes>>
@
\end{center}
\caption{LDheatmap of GIMAP5.CEU data with UCSC Genes track.}
\label{fig:fig2}
\end{figure}

We now add a track that illustrates recombination rates:
<<label=addRecomb,eval=FALSE>>=
llGenesRecomb <- LDheatmap.addRecombRate(llGenes, chr="chr7", genome="hg18")
@
The result is shown in Figure~\ref{fig:fig3}. Users who wish to 
avoid the download time may draw the
grob from the \texttt{llGenesRecomb} object with the commands:
<<label=addRecomb,eval=FALSE>>=
grid.newpage()
grid.draw(llGenesRecomb$LDheatmapGrob)
@
\begin{figure}
\begin{center}
<<label=fig3,fig=TRUE,echo=FALSE>>=
<<addRecomb>>
@
\end{center}
\caption{LDheatmap of GIMAP5.CEU data with UCSC Genes and Recomb Rate tracks.}
\label{fig:fig3}
\end{figure}

\section{Adding a scatterplot above the heatmap}

Above the two tracks fetched from the UCSC genome browser
we will superpose a scatterplot depicting the results of
association tests between each SNP in the {\it GIMAP5} gene
and a hypothetical disease.
We create a set of hypothetical p-values that includes
one very small value.
<<>>=
set.seed(1)
atests<-runif(nrow(GIMAP5.CEU$snp.support))
names(atests)<-rownames(GIMAP5.CEU$snp.support)
atests["rs6598"]<-1e-5
@
A Manhattan plot (scatterplot of the $-\log_{10}$ p-values versus genomic
location of the SNPs) 
can be added to the heatmap as follows:
<<label=addScatter,include=FALSE>>=
llGenesRecombScatter<-LDheatmap.addScatterplot(llGenesRecomb,-log10(atests),
ylab="-log10(p-values)")
@
\begin{figure}
\begin{center}
<<label=fig4,fig=TRUE,echo=FALSE>>=
grid.newpage()
grid.draw(llGenesRecombScatter$LDheatmapGrob)
@
\end{center}
\caption{LDheatmap of GIMAP5.CEU data with gene, recombination and Manhattan plot tracks.}
\label{fig:fig4}
\end{figure}
The result is shown in Figure~\ref{fig:fig4}.

\section{Adding other graphical objects}
\label{sec:addGrob}

Users may add an arbitrary grob above the heatmap with the 
\texttt{LDheatmap.addGrob} function. This function aligns the 
grob so that its horizontal dimension matches the gene map of the
heatmap. 
The idea is to allows users to create tracks
from grobs returned by 
other R packages that use grid graphics (e.g., \texttt{ggplot2}), 
or to add images using functions such as \texttt{rasterGrob}
or \texttt{pictureGrob}. 
Development of \texttt{LDheatmap.addGrob} is in its early stages. 

The first illustration of \texttt{LDheatmap.addGrob} is to add 
a Manhattan plot created with the \texttt{qplot} (quick plot) function 
of \texttt{ggplot2}.
We convert the output of \texttt{qplot} to a grob that can be added to our
heatmap.
<<label=addQplot,include=FALSE>>=
require("ggplot2")
posn<-GIMAP5.CEU$snp.support$Position
manhattan2<-ggplotGrob(
{
qplot(posn,-log10(atests),xlab="", xlim=range(posn),asp=1/10)
last_plot() + theme(axis.text.x=element_blank(),
              axis.title.y = element_text(size = rel(0.75)))
}
)
llQplot<-LDheatmap.addGrob(ll,manhattan2,height=.7)
@ 
The resulting heatmap is given in Figure~\ref{fig:fig5}
\begin{figure}
\begin{center}
<<label=fig5,fig=TRUE,echo=FALSE>>=
grid.newpage()
grid.draw(llQplot$LDheatmapGrob)
@
\end{center}
\caption{LDheatmap of GIMAP5.CEU data with Manhattan plot created with
the \texttt{qplot} function.}
\label{fig:fig5}
\end{figure}

In Figure~\ref{fig:fig5} the entire width of the
\texttt{manhattan2} grob has been aligned to the genetic map, but
the x-axis within this grob does not align. 
As a kludge, 
the Manhattan plot can be aligned manually by creating space on 
the display and using grid graphics functionality to place the 
grob: 
\begin{enumerate}
\item
Use \texttt{LDheatmap.addGrob} to make space for
\texttt{manhattan2} grob, by adding a track comprised of 
a rectangle with white lines.
\item Use grid graphics to create a viewport on the display to hold 
the \texttt{manhattan2} grob. 
This step takes some trial and error
with the placement and size of the viewport.
\item Draw the \texttt{manhattan2} grob. 
\item (optional) Clean up by popping the viewport created to hold 
the \texttt{manhattan2} grob off of the viewport tree.
\end{enumerate}
Alignment of the x-axis of the \texttt{manhattan2} 
is device-dependent. The following commands align the grob when displayed on
the \texttt{pdf()} device. The resulting plot does not render well
when the vignette is generated from its Sweave source, 
and so the plot is excluded from the vignette.
<<label=addQplot2,include=FALSE>>=
manhattan2<-ggplotGrob(
{
qplot(posn,-log10(atests),xlab="", xlim=range(posn),asp=1/10)
last_plot() + theme(axis.text.x=element_blank(),
              axis.title.y = element_text(size = rel(0.75)))
}
)
llQplot2<-LDheatmap.addGrob(ll,rectGrob(gp=gpar(col="white")),height=.2)
pushViewport(viewport(x=.48,y=.76,width=.99,height=.2))
grid.draw(manhattan2)
popViewport(1)
dev.off()
@

As another illustration of \texttt{LDheatmap.addGrob},
we add to a heatmap an ideogram of human chromosome 7 with the GIMAP5
region highlighted.
The ideogram was manually saved as a PDF file
from the UCSC genome browser 
and was then exported from PDF to a PNG file using 
Adobe Acrobat. The PNG file was read into R with the 
\texttt{readPNG} function from the \texttt{png} package and 
saved as the object \texttt{GIMAP5ideo}, which 
was loaded at the beginning of the vignette.
We package the image as a grob with the 
\texttt{rasterGrob} function and 
use \texttt{LDheatmap.addGrob} to add it to the heatmap:
<<label=addGrob,include=FALSE>>=
llImage<-LDheatmap.addGrob(ll,rasterGrob(GIMAP5ideo))
@ 
The results are shown in Figure~\ref{fig:fig6}.
The resolution of the ideogram image is rather poor, and the ideogram
is only readable when one zooms in on the plot.
\begin{figure}
\begin{center}
<<label=fig6,fig=TRUE,echo=FALSE>>=
grid.newpage()
grid.draw(llImage$LDheatmapGrob)
@
\end{center}
\caption{LDheatmap of GIMAP5.CEU data with chromosome 7 ideogram showing 
the GIMAP5 region as a small vertical red line.}
\label{fig:fig6}
\end{figure}

\section{Details and Limitations}

\subsection{Details}
Grobs representing the 
components of an \texttt{LDheatmap} object are stored as children of 
the main \texttt{LDheatmapGrob} grob:
<<pl>>=
names(ll$LDheatmapGrob$children)
@
Grobs that represent tracks are 
added to the list of children:
<<>>=
names(llGenesRecombScatter$LDheatmapGrob$children)
@
The \texttt{transcripts} grob is the genes track added 
by \texttt{LDheatmap.addGenes}, \texttt{recombRate}
is the graphic depicting recombination rates added by
\texttt{LDheatmap.addRecombRate}, and 
\texttt{association} is the Manhattan plot added by
\texttt{LDheatmap.addScatter}.
Tracks added by \texttt{LDheatmap.addGrob} have names
assigned by the grid graphics system:
<<>>=
names(llQplot$LDheatmapGrob$children)
names(llImage$LDheatmapGrob$children)
@

\subsection{Limitations}
\label{sec:limitations}
Adding tracks requires creating space for the track and
moving the gene map title (the one that 
reports the map length) and the main title of the graph. 
Currently the calculation of space needed to insert a track only
accounts for tracks previously added 
by \texttt{LDheatmap.addGenes},
\texttt{LDheatmap.addRecombRate},  or
\texttt{LDheatmap.addScatter}. Thus, a track added by
\texttt{LDheatmap.addGrob} will be over-written by subsequent tracks. 
Another limitation is that SNP names added with the \texttt{SNP.name}
argument of \texttt{LDheatmap} are over-written by any tracks added
to the heatmap.
We intend to fix these limitations in a future release of the package.


%\citep*{grg}
%\citet{geyer-temp}
%
%\begin{thebibliography}{}
%
%\bibitem[Gelman et al.(1996)Gelman, Roberts, and Gilks]{grg}
%Gelman, A., G.~O. Roberts, and W.~R. Gilks (1996).
%\newblock Efficient Metropolis jumping rules.
%\newblock In \emph{Bayesian Statistics, 5 (Alicante, 1994)}, pp.~599--607.
%  Oxford University Press.
%
%\bibitem[Geyer(1992)]{practical}
%Geyer, C.~J. (1992).
%\newblock Practical Markov chain Monte Carlo (with discussion).
%\newblock \emph{Statistical Science}, 7, 473--511.
%
%\bibitem[Geyer and Thompson(1995)]{geyer-temp}
%Geyer, C.~J. and E.~A. Thompson (1995).
%\newblock Annealing Markov chain Monte Carlo with applications to
%    ancestral inference.
%\newblock \emph{Journal of the American Statistical Association}, 90,
%909--920.
%
%\end{thebibliography}

\end{document}


