%\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 investigating
genetics of gene expression}
\author{VJ Carey \texttt{<stvjc@channing.harvard.edu>}}
\maketitle

\section{Introduction}

The \Rpackage{GGtools} package supports
data analysis activities that link genotypic information
such as SNP configurations to gene expression phenotype.
@
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
derived from 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, if seeking interactive performance with
commodity hardware, 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 is 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>>=
data(geneLocs_hsa)
plot_mlp(sCPNE1, chr20meta, geneLocDF=geneLocs_hsa)
@

%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{Creating new \texttt{racExSet} instances and supporting
objects}

Genotypes are available for a wide number of
reference strains of mice through the Wellcome Trust;
expression studies of such mice can be obtained using GEO.
The data structure \texttt{gse2031GG} is an example:
<<lkgse>>=
data(gse2031GG)
gse2031GG
@
The expression data are obtained by applying RMA to a
collection of CEL files readily obtainable from GEO.
The strains of the mice giving rise to the samples are
obtained from the WebQTL "contacts" link.

The genotype data for a large number of strains can be obtained
from the INBREDS distribution at Wellcome Trust.  The URL
is 
\url{www.well.ox.ac.uk/mouse/INBREDS}.  An example of the genotype
data 
can be seen using the following code:
<<lkgt>>=
ff = readLines(system.file("fileDemos/StrainInit.txt", package="GGtools"))
t(sapply(strsplit(ff[1:5]," "), function(x)x[c(1,2,3,55,56,57,101)]))
@
There are various representations provided in various places.
Tightening the path to accurate genotyping data for reference strains
would be useful.

Given an expression matrix and a genotyping file as provided at the INBREDS
site, the INBREDSworkflow function will help generate a \texttt{racExSet}
instance:
<<lkiw>>=
args(INBREDSworkflow)
@
The \texttt{estrains} argument is a vector of strings defining
the reference strains from which columns of \texttt{emat} were
obtained through microarray hybridization.  This vector will typically
be available from a \texttt{phenoData} variable on the
\texttt{ExpressionSet} from which the expression matrix was obtained.
For the \texttt{gse2031GG}, this is the variable \texttt{strain}:
<<lks>>=
as.character(gse2031GG$strain[1:5])
@
Now it appears that this nomenclature for strains is not always
adopted directly.  If we grep for a variation of BXD[nn] in
the INBREDS strains file, we find
<<lktag>>=
grep("BXD-", strsplit(ff[1]," ")[[1]], value=TRUE)[1:5]
@
and we see there are lexical variations introduced -- we want BXD6,
but the labels in the INBREDS file have hyphens.  A \texttt{fixup}
parameter allows inline reformatting of deviant strain identifiers
as column names of the INBREDS file.  Additionally, chromosome names
are sometimes postpended with ``-random'', and some steps may be needed
to simplify the chromosome nomenclature.  The user must create/find the R
code to obtain the desired results, or manually massage the source
data files so that the desired regularities are present.

Once genotype codes are available with a strain naming convention
that matches that of the expression samples, INBREDSworkflow will
compute rare allele counts, bind the genotyping,
expression, and phenotype data together, and generate a \texttt{racExSet} instance.

In summary, \texttt{make\_racExSet} is a specification of materials
that must be made mutually compatible for creation of a \texttt{GGtools}
resource with which investigations of genetics of gene expression
can be conveniently made.  \texttt{HMworkflow} and \texttt{INBREDSworkflow}
are utility functions that support this task for files coming from
HapMap and Wellcome repositories.




 

\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}
