%
% NOTE -- ONLY EDIT THE .Rnw FILE!!!  The .tex file is
% likely to be overwritten.
%

%\VignetteIndexEntry{gwascat -- exploring NHGRI GWAS catalog}
%\VignetteDepends{GenomicRanges, DO.db}
%\VignetteKeywords{genetics, GWAS}
%\VignettePackage{gwascat}

\documentclass[12pt]{article}

\usepackage{amsmath,pstricks}
\usepackage{auto-pst-pdf}
\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{\textit{gwascat}: structuring and querying the NHGRI GWAS catalog}
\author{VJ Carey\footnote{Generous support of Robert Gentleman and the
Computational Biology Group of
Genentech, Inc. is gratefully acknowledged}.}
\maketitle
\tableofcontents

\section{Introduction}

NHGRI maintains and routinely updates a database of selected genome-wide
association studies.  This document describes R/Bioconductor facilities for
working with contents of this database.

\subsection{Installation}
The package can be installed using Bioconductor's \textit{BiocInstaller}
package, with the sequence
\begin{verbatim}
library(BiocInstaller)
biocLite("gwascat")
\end{verbatim}

\subsection{Attachment and access to documentation}
Once the package has been installed, use \verb+library(gwascat)+ to 
obtain interactive access to all the facilities.  After executing
this command, use \verb+help(package="gwascat")+ to obtain an overview.
The current version of this vignette can always be accessed at
www.bioconductor.org, or by suitably navigating the web pages generated
with \texttt{help.start()}.

<<getlib,echo=FALSE,results=hide>>=
library(gwascat)
@

Some noteworthy limitations:  As of
\Sexpr{gwrngs@extractDate},
there are \Sexpr{sum(gwcat$SNPs == "NR")}
records in the database for which no SNP is identified.
There are \Sexpr{sum(gwcat$Chr_pos == "")} records for which no
chromosomal position of the associated locus is given.


\subsection{Illustrations}

<<lk1,echo=FALSE,results=hide>>=
if (length(grep("gwascat", search()))>0) detach("package:gwascat")
@
Available functions are:
<<lk2>>=
library(gwascat)
objects("package:gwascat")
@

The GRanges instance with all SNP-disease associations is:
<<lkgr>>=
gwrngs
@

To determine the most frequently occurring traits:
<<lktr>>=
topTraits(gwrngs)
@

For a given trait, obtain a GRanges with all recorded associations; here
only three associations are shown:
<<lklocs>>=
subsetByTraits(gwrngs, tr="LDL cholesterol")[1:3]
@

\section{SNP sets and trait sets}

\subsection{SNPs by name}

We can regard the content of a SNP chip as a set of SNP, referenced
by name.  The
pd.genomewidesnp.6 package describes the Affymetrix SNP 6.0 chip.
We can determine which traits are associated with loci interrogated
by the chip as follows.  We work with a subset of the 1 million loci
for illustration.

The \texttt{locon6} data frame has information on 10000 probes,
acquired through the following code (not executed here to reduce
dependence on the pd.genomewidesnp.6 package, which is very large.
<<dobig,eval=FALSE>>=
library(pd.genomewidesnp.6)
con = pd.genomewidesnp.6@getdb()
locon6 = dbGetQuery(con, 
   "select dbsnp_rs_id, chrom, physical_pos from featureSet limit 10000")
@
Instead use the serialized information:
<<doloc>>=
data(locon6)
rson6 = as.character(locon6[[1]])
rson6[1:5]
@

We subset the GWAS ranges structure with rsids that
are common to both the chip and the GWAS catalog.
We then tabulate the diseases associated with the
common loci.
<<lkdtab>>=
intr = gwrngs[ intersect(getRsids(gwrngs), rson6) ]
sort(table(getTraits(intr)), decreasing=TRUE)[1:10]
@

\subsection{Traits by genomic location}

We will assemble genomic coordinates for SNP on the Affymetrix 6.0 chip
and show the effects of identifying the trait-associated
loci with regions of width 1000bp instead
of 1bp.

The following code retrieves coordinates for SNP interrogated
on 10000 probes (to save time)
on the 6.0 chip, and stores the results in a GRanges instance.
<<lkexp,keep.source=TRUE>>=
gr6.0 = GRanges(seqnames=ifelse(is.na(locon6$chrom),0,locon6$chrom), 
       IRanges(ifelse(is.na(locon6$phys),1,locon6$phys), width=1))
elementMetadata(gr6.0)$rsid = as.character(locon6$dbsnp_rs_id)
seqlevels(gr6.0) = paste("chr", seqlevels(gr6.0), sep="")
@

Here we compute overlaps with both the raw disease-associated locus
addresses, and with the locus address $\pm$ 500bp.
<<dosub>>=
ag = function(x) as(x, "GRanges")
ovraw = subsetByOverlaps(ag(gwrngs), gr6.0)
length(ovraw)
ovaug = subsetByOverlaps(ag(gwrngs+500), gr6.0)
length(ovaug)
@

To acquire the subset of the catalog to which 6.0 probes are
within 500bp, use:
<<dosub>>=
rawrs = elementMetadata(ovraw)$SNPs
augrs = elementMetadata(ovaug)$SNPs
gwrngs[augrs]
@

Relaxing the intersection criterion in this
limited case leads to a larger set of traits.

<<lkrelax>>=
setdiff( getTraits(gwrngs[augrs]), getTraits(gwrngs[rawrs]) )
@

\section{Counting alleles associated with traits}

We can use \texttt{riskyAlleleCount} to count
risky alleles enumerated in the GWAS catalog.  This
particular function assumes that we have genotyped at
the catalogued loci.  Below we will discuss how to impute
from non-catalogued loci to those enumerated in the catalog.
<<lkcout>>=
data(gg17N) # translated from GGdata chr 17 calls using ABmat2nuc
gg17N[1:5,1:5]
@
This function can use genotype information in the A/B format,
assuming that B denotes the alphabetically later nucleotide.
Because we have direct nucleotide coding in our matrix, we set
the \texttt{matIsAB} parameter to false in this call.

<<dorun>>=
h17 = riskyAlleleCount(gg17N, matIsAB=FALSE, chr="ch17")
h17[1:5,1:5]
table(as.numeric(h17))
@

It is of interest to bind the counts back to the catalog data.
<<domo>>=
gwr = gwrngs
gwr = gwr[colnames(h17),]
elementMetadata(gwr) = cbind(elementMetadata(gwr), DataFrame(t(h17)))
sn = rownames(h17)
gwr[,c("Disease.Trait", sn[1:4])]
@

Now by programming on the elementMetadata, we can identify
individuals with particular risk profiles.

\section{Imputation to unobserved loci}

If we lack information on a specific locus $s$, but have
reasonably dense genotyping on a subject, population genetics
may allow a reasonable guess at the genotype at $s$ for
this subject.  Many algorithms for genotype imputation have been proposed.
Here we use a very simple approach due to
David Clayton in the \textit{snpStats} package.

We use the ``low coverage'' 1000 genomes genotypes for
the CEU (central European) HapMap cohort as a base for constructing
imputation rules.  We focus on chromosome 17 for illustration.

The base data are
<<getbase,cache=TRUE>>=
data(low17)
low17
@

A somewhat sparser set of genotypes (HapMap phase II,
genomewide 4 million loci) on chromosome 17 is 
archived as \texttt{g17SM}.
This has a compact SnpMatrix encoding of genotypes.
<<lkggd,cache=TRUE>>=
data(g17SM)
g17SM
@
For a realistic demonstration, we use the subset of these
loci that are present on the Affy 6.0 SNP array.
<<dog,cache=TRUE>>=
data(gw6.rs_17)
g17SM = g17SM[, intersect(colnames(g17SM), gw6.rs_17)]
dim(g17SM)
@

The base data were used to create a set of rules allowing
imputation from genotypes in the sparse set to the richer set.
Some rules involve only a single locus, some as many as 4.
The construction of rules involves tuning of modeling parameters.
See snp.imputation in snpStats for details.
<<lkrul,cache=TRUE>>=
data(rules_6.0_1kg_17)
rules_6.0_1kg_17[1:5,]
@
The summary of rules shows the degree of association between
the predictors and predictands in terms of $R^2$.  Many potential
targets are not imputed.
<<lksum>>=
summary(rules_6.0_1kg_17)
@
The overlap between the 6.0-resident g17SM loci and the catalog is
<<lkov>>=
length(intersect(colnames(g17SM), values(gwrngs)$SNPs))
@

The new expected B allele counts are
<<doimp>>=
exg17 = impute.snps(rules_6.0_1kg_17, g17SM)
@
The number of new loci that coincide with risk loci in the catalog is:
<<lkl>>=
length(intersect(colnames(exg17), values(gwrngs)$SNPs))
@

\section{Formalizing disease traits with Disease Ontology}

The Disease Ontology project \cite{osborne} formalizes a vocabulary
for human diseases.  Bioconductor's DO.db package is a curated
representation.

<<getdo>>=
library(DO.db)
DO()
@

All tokens of the ontology are acquired via:
<<getallt>>=
alltob = unlist(mget(mappedkeys(DOTERM), DOTERM))
allt = sapply(alltob, Term)
allt[1:5]
@

Direct mapping from disease trait tokens in the catalog
to this vocabulary succeeds for a modest proportion of
records.
<<dohit>>=
cattra = elementMetadata(gwrngs)$Disease.Trait
mat = match(tolower(cattra), tolower(allt))
catDO = names(allt)[mat]
catDO[1:50]
mean(is.na(catDO))
@

Approximate matching of unmatched tokens can proceed
by various routes.  Some traits are not diseases, and
will not be mappable.  However, consider
<<dogr>>=
unique(cattra[is.na(catDO)])[1:20]
nomatch = cattra[is.na(catDO)]
@
Manual searching shows that a number of these have
very close matches.

% NB the graph stuff can be used for other OBO without .db packages..
%The \textit{gwascat} package includes an OBO-formatted
%image of the ontology and a model for it as a \texttt{graphNEL} instance
%as defined in the Bioconductor \textit{graph} package \cite{carey2005}.

%The graph model is constructed as follows.
%<<dogrmod>>=
%doobo = dir(system.file("obo", package="gwascat"), full=TRUE, patt="HumanDO")
%DOgraph = obo2graphNEL(doobo)
%DOgraph
%@
%
%Nodes are the formal disease term identifiers; node data provides additional
%metadata.
%<<lkno>>=
%nodeData(DOgraph)[1:5]
%@

\section{Appendix: Adequacy of location annotation}

A basic question concerning the use of archived
SNP identifiers is durability
of the association between asserted location
and SNP identifier.  The \texttt{chklocs}
function uses a current Bioconductor
SNPlocs package to check this.

For example, to verify that locations asserted on chromosome
20 agree between the Bioconductor dbSNP image and the
gwas catalog,
<<getrs,keep.source=TRUE,cache=TRUE>>=
if ("SNPlocs.Hsapiens.dbSNP.20110815" %in% installed.packages()[,1])
suppressWarnings(chklocs("20"))
@

This is not a fast procedure but has succeeded for 
all chromosomes 1-22 when checked off line.


\end{document}