%\VignetteIndexEntry{GGtools overview}
%\VignetteDepends{}
%\VignetteKeywords{Genetical genomics,SNP,expression}
%\VignettePackage{GGtools}

%
% NOTE -- ONLY EDIT THE .Rnw FILE!!!  The .tex file is
% likely to be overwritten.
%
\documentclass[12pt]{article}

\usepackage{amsmath,pstricks}
\usepackage[authoryear,round]{natbib}
\usepackage{hyperref}


\textwidth=6.2in
\textheight=8.5in
%\parskip=.3cm
\oddsidemargin=.1in
\evensidemargin=.1in
\headheight=-.3in

\newcommand{\scscst}{\scriptscriptstyle}
\newcommand{\scst}{\scriptstyle}


\newcommand{\Rfunction}[1]{{\texttt{#1}}}
\newcommand{\Robject}[1]{{\texttt{#1}}}
\newcommand{\Rpackage}[1]{{\textit{#1}}}
\newcommand{\Rmethod}[1]{{\texttt{#1}}}
\newcommand{\Rfunarg}[1]{{\texttt{#1}}}
\newcommand{\Rclass}[1]{{\textit{#1}}}

\textwidth=6.2in

\bibliographystyle{plainnat} 
 
\begin{document}
%\setkeys{Gin}{width=0.55\textwidth}

\title{Overview of \Rpackage{GGtools} for genetical genomics}
\author{VJ Carey \texttt{<stvjc@channing.harvard.edu>}}
\maketitle

\section{Introduction}

We use the term \textit{genetical genomics} to refer to
data analysis activities that link genotypic information
such as SNP configurations to gene expression phenotype.
The \Rpackage{GGtools} package includes various
demonstration resources and analysis tools for these
activities.
@
We will attach the library and have a look at a basic
demonstration resource.
<<init>>=
library(GGtools)
data(chr20GGdem)
class(chr20GGdem)
chr20GGdem
@
The \Sexpr{class(chr20GGdem)} class is an extension of the
\Rclass{eSet} class.  It represents expression
data from the hgfocus chip on 48 individuals in the CEU
CEPH group, and SNP data obtained from their HapMap genotyping
results.

The data are organized into an \Sexpr{nrow(exprs(chr20GGdem))}
by \Sexpr{ncol(exprs(chr20GGdem))} matrix of expression values
accessible with the \texttt{exprs} method, and
an \Sexpr{nrow(snps(chr20GGdem))} by \Sexpr{ncol(snps(chr20GGdem))}
of rare allele counts:
<<lkc20>>=
dim(exprs(chr20GGdem))
dim(snps(chr20GGdem))
snps(chr20GGdem)[1:5,1:5]
@

We need some genetic metadata about SNPs; these are
culled from SNP genotyping panels released on a 
chromosome-by-chromosome basis for CEPH participants
by HapMap project:
<<lksmeta>>=
data(chr20meta)
chr20meta[1:4,]
@

A basic task is to compute a screen (over the genome, or,
more practically, over a chromosome) of genotypic determination
of expression.  The \Rfunction{snpScreen} method helps
with this; we illustrate an example related to
results in Cheung and Spielman 2005:

<<do20>>=
chr20GGdem = exclMono(chr20GGdem)
S100 = snpScreen(chr20GGdem, chr20meta, genesym("CPNE1"), ~., lm, gran=30) 
S100
@
A primitive display is obtained as follows.  We know that \texttt{lm}
was used, so the relevant p-values are in the coefficient component
of the summarized fit objects.
<<dofig,fig=TRUE>>=
ps = as.numeric(sapply(S100, function(x) try(summary(x)$coef[2,4])))
plot(S100@locs, -log10(ps), xlab="chromosomal location")
@

%\begin{center}
%<<dosc,fig=TRUE,results=hide>>=
%ss1 = snpScreen(chr20GGdem, 
@

%The \Rfunction{detScreen} function returns a list of fit results:
%<<lkds>>=
%names(ds1)
%ds1[[1]]
%length(ds1[[2]])
%names(ds1[[2]])
%ds1[[2]][[1]][[1]]
@

\section{Performance-oriented specialization}

The \texttt{snpScreen} method illustrated above is very
general (can accommodate and retain results of any R modeling function)
but fairly slow.  We have added an R function \texttt{fastAGM}
for fast fitting of an additive genetic model (equivalent to
but much faster than using \texttt{lm}).  
%We apply this
%for IRF5 on chromosome 7, and attempt to fit on all SNPs that
%have no missing values:
<<dofast>>=
ut = unix.time(sCPNE1 <- snpScreen(chr20GGdem, chr20meta, genesym("CPNE1"), ~.,
  fastAGM, 75))
ut
sCPNE1
wm = which.min(pp <- extract_p(sCPNE1))
pp[wm]
<<dofig5,fig=TRUE>>=
plot_mlp(sCPNE1, chr20meta)
@

%With the conservative approach to missing data taken thus far,
%the result does not agree with that of Cheung and Spielman,
%who found rs2280714 to be most significant.  The two SNPs
%are not completely concordant, but the basic patterns are similar:
<<lkevg,fig=TRUE>>=
# par(mfrow=c(1,2))
 plot_EvG(chr20GGdem, genesym("CPNE1"), "rs6060535")
# plot_EvG(chr7GGdem, genesym("IRF5"), "rs6969930")
@

 

\section{Appendix: Package documentation for \Rpackage{GGtools}}
<<getplop,echo=FALSE>>=
plop = function(x)
{
    if (!inherits(x, "packageInfo")) 
        stop("wrong class")
    #outConn <- file(outFile, open = "w")
    vignetteMsg <- gettextf("Further information is available in the following vignettes in directory '%s':", 
        file.path(x$path, "doc"))
    headers <- c(gettext("Description:\n\n"), gettext("Index:\n\n"), 
        paste(paste(strwrap(vignetteMsg), collapse = "\n"), "\n\n", 
            sep = ""))
    footers <- c("\n", "\n", "")
    formatDocEntry <- function(entry) {
        if (is.list(entry) || is.matrix(entry)) 
            formatDL(entry, style = "list")
        else entry
    }
    writeLines(gettextf("\n\t\tInformation on package '%s'\n", 
        x$name))#, outConn)
    for (i in which(!sapply(x$info, is.null))) {
        writeLines(headers[i])#, outConn, sep = "")
        writeLines(formatDocEntry(x$info[[i]]))#, outConn)
        writeLines(footers[i])#, outConn, sep = "")
    }
#    close(outConn)
}
<<gag,echo=FALSE>>=
plop(help(package="GGtools"))
@

Session information for this vignette build:
<<dos>>=
sessionInfo()
@

\end{document}
