% NOTE -- ONLY EDIT THE .Rnw FILE!!!  The .tex file is
% likely to be overwritten.
%
%\VignetteIndexEntry{EMBO03 Lab 2}
%\VignetteDepends{Biobase,annotate,tkWidgets,geneplotter,golubEsets,hu6800,GO}
%\VignetteKeywords{Annotation}
\documentclass[12pt]{article}

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


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

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

\bibliographystyle{plainnat}


\title{EMBO03 Lab 2: Introduction to Bioconductor Annotation Packages}

\author{Sandrine Dudoit and Robert Gentleman}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}

\maketitle

One of the largest challenges in analyzing genomic data is associating
the experimental data with the available biological {\em metadata}, e.g.,
sequence, gene annotation, chromosomal maps, literature. Bioconductor
provides two main packages for this purpose: \texttt{annotate}
(end-user) and \texttt{AnnBuilder} (developer). These packages can be
used for querying databases such as GenBank, GO, LocusLink, and
PubMed, from R, and for processing the query results in R.\\

In this lab, we focus on annotation resources for Affymetrix
chips. Since Affymetrix chips have a standard layout, Bioconductor
can provide annotation data packages for the main types of Affymetrix
chips, e.g., hu6800, hgu33, hgu95, mgu74, and rgu34 series. These data packages are built using the
\texttt{AnnBuilder} package. By contrast, there is no standard
array design for two-color spotted microarray experiments; each lab tends
to have its own custom design. In this case, specific annotation data packages
will have to be created for each facility using
\texttt{AnnBuilder}. Once annotation data packages are constructed to provide
mappings between different sets of gene identifiers, the tools in
\texttt{annotate} can be used in a similar manner for both platforms.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Mappings between different gene identifiers}

A main task in annotation is to associate manufacturer (e.g., for
Affymetrix chips) or in-house (e.g., for custom cDNA spotted arrays)
probe identifiers to other available identifiers
(e.g., PMID, GenBank accession number).\\
 
In this lab, we will use the Bioconductor annotation data packages to map
between Affymetrix identifiers and identifiers for biological metadata
available on the WWW.
Because these data sources are very large, are constantly evolving, and
are similar across species, we have adopted the strategy of distributing
data in regularly updated R packages. Each mapping is contained in an {\em environment}.
For our purposes, an environment is simply a {\em hash table}.
It provides a very fast way of looking up {\em values} for specific
text {\em keys}.
This implementation is not the best and we hope to develop real hash tables
for R but that will take some time.\\


We will again use the Golub et al. (1999) dataset as a case study and
will explore annotation resources for the Affymetrix HU6800 chip used
in this study. The Bioconductor
annotation data package for the HU6800 chip is \texttt{hu6800} and can
be downloaded from the "Data packages" section of the Bioconductor
website. To load the packages needed for this lab
<<loadAnnot1,eval=FALSE>>=
library(Biobase)
library(annotate)
library(tkWidgets)
library(geneplotter)
library(golubEsets)
library(hu6800)
library(GO)
@
Alternately, since all these packages are required by the
\texttt{EMBO03} course package, simply use
<<loadAnnot2,eval=TRUE>>=
library(EMBO03)
@

After loading the \texttt{hu6800} package, we have access to a number of
different datasets (environments). These provide the many different
mappings from Affymetrix identifiers to other probe identifiers, such
as chromosomal location (in the \texttt{hu6800CHR} environment) and PMID
(in the \texttt{hu6800PMID} environment). Quality control data, counts,
etc., for the mappings are available by calling the function
\texttt{hu6800()} and by examining the help page, \texttt{?hu6800}.
<<hu6800>>=
hu6800()
@

The main R functions we will use for dealing with environments are \texttt{ls} (to list the names of objects in a
specified environment) and  \texttt{get} (to search for an R object with
a given name and return its value, if found, in the specified environment).
We obtain the Affymetrix identifiers using \texttt{ls} and note that there are 7,129 identifiers for the
HU6800 chip (this matches the Golub \texttt{exprSets}).
<<lsID>>=
affyID <- ls(env=hu6800CHR)
length(affyID)
@

We now arbitrarily select the probe set with Affymetrix identifier
\texttt{"U18237\_at"} and obtain a number of other IDs corresponding to
this probe set. Note that there may be several PMIDs, i.e., PubMed abstracts, associated
with a given gene.

<<getID>>=
 mygene <- affyID[4001]
 mygene
 get(mygene, env = hu6800ACCNUM)
 get(mygene, env = hu6800LOCUSID)
 get(mygene, env = hu6800SYMBOL)
 get(mygene, env = hu6800GENENAME)
 get(mygene, env = hu6800SUMFUNC)
 get(mygene, env = hu6800UNIGENE)
 get(mygene, env = hu6800CHR)
 get(mygene, env = hu6800CHRLOC)
 get(mygene, env = hu6800CHRORI)
 get(mygene, env = hu6800MAP)
 get(mygene, env = hu6800PMID)
 get(mygene, env = hu6800GO)
@

In some cases, we need to obtain data on several genes at once. We
wrote a special function for this purpose: \texttt{multiget}.
<<multigetID>>=
fivegenes<-affyID[6:10]
fivegenes
multiget(fivegenes, env = hu6800PMID)
@

Instead of relying on the general R functions for environments (e.g., \texttt{get}), the
development version of the \texttt{annotate} package also provides new and more user-friendly
functions for accessing specific identifiers. 
<<easyGetID>>=
getSYMBOL(fivegenes, data="hu6800")
getLL(fivegenes, data="hu6800")
getPMID(fivegenes, data="hu6800")
gg<- getGO(fivegenes, data="hu6800")
getGOdesc(gg[[2]], "MF")
@

A variety of information about the probes can now easily be obtained
by, for example, querying PubMed, GenBank, and LocusLink, as described below.  In
addition, one can also use R to compute on these data and perform a
number of quality control and exploratory analyses.
<<chromTable>>=
 whChrom <- multiget(affyID, env=hu6800CHR)
 table(unlist(whChrom))
 vv<-sapply(whChrom, length)
 table(vv)
 whChrom[vv==2]
@
This yields the distribution of probe sets by chromosome.
Note that there are 10 genes that have been assigned two chromosome locations.
We might want to explore these further by examining resources at the NCBI.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Querying PubMed}

The National Library of Medicine (NLM) provides a great deal of information
on biological data. We have
only just started to explore these resources. We have already
developed tools for interacting with PubMed, but need to develop
further functionality on MeSH and other available resources.
A mapping has been done between LocusLink and PubMed.  This associates specific
articles with specific genes. The articles can be queried interactively using
tools available in R and other systems.
Full text abstracts are readily available; in some cases full text
articles are also available (we don't have the tools to make use of
the later).
We next demonstrate how to interact with PubMed using functions
in \texttt{annotate}.\\

Bioconductor functions for querying PubMed and other WWW databases
from R rely on the \texttt{browseURL} function and the R \texttt{XML}
package to parse query results. 
For more details, please consult articles in  R News (Gentleman \&
Gentry (2002), {\it R News} 2(2); Temple Lang (2001), {\it R News} 1(1)).

<<browseURL,eval=FALSE>>=
browseURL("www.r-project.org")
@

A \texttt{pubMedAbst} class structure was defined for handling 
PubMed abstracts in R. The slots are
<<pubMedAbst>>=
slotNames("pubMedAbst")
@

The basic engine for talking to PubMed is the function
\texttt{pubmed}. Given a vector of PMIDs, the function either has a
browser display a URL showing the results of the PubMed query for those
identifiers (\texttt{disp="browser"}) or creates an \texttt{XMLdoc}
object with the same data (\texttt{disp="data"}). 
<<pubmed1,eval=FALSE>>=
pmids<-getPMID(mygene,data="hu6800")
pubmed(pmids, disp="browser")
absts1<-pubmed(pmids, disp="data")
@

The function \texttt{pm.getabst} provides a simpler way to 
download the specified PubMed abstracts (stored in XML) and create a
list of  \texttt{pubMedAbst} objects. The following commands can be
used to store the PubMed abstracts for 5 genes in a list of objects of
class \texttt{pubMedAbst}, to compute the number of abstracts retrieved
for each gene, and to print the abstract(s) for the first gene. We can
see, for example, that one of the genes has 21 abstracts associated with it.
<<pubmed2>>=
absts2<-pm.getabst(fivegenes,"hu6800")
lapply(absts2,length)
absts2[[1]]
@

The functions 
\texttt{pm.titles} and  \texttt{pm.abstGrep} can then be used to
extract the titles from a set of PubMed abstracts and for regular
expression matching on these abstracts, respectively. For the genes with Affy IDs \texttt{fivegenes}, the following
commands extract the abstract titles and search the
abstracts for the word "protein", with lower case "p" or upper case "P".
<<pubmed3>>=
pm.titles(absts2)
sapply(absts2, function(x) pm.abstGrep("[Pp]rotein", x)) 
@

The new function \texttt{pmAbst2html} from the development version of
\texttt{annotate} takes a list of
\texttt{pubMedAbst} objects and generates an HTML report with
the titles of the abstracts and links to their full page on
PubMed. The report will be stored in the file \texttt{pm.html} in the
working directory.


%%%***<<pmAbst2html,eval=FALSE>>=
\begin{Sinput}
> pmAbst2html(absts2[[2]],filename="pm.html")
\end{Sinput}
%%%@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Querying GenBank}

Another source of data is GenBank. We can interact with GenBank in
much the same way as with PubMed; the main function for this purpose
is \texttt{genbank}.

<<genbank,eval=FALSE>>=
gbacc<-multiget(fivegenes,hu6800ACCNUM)
genbank(gbacc,disp="browser")
gb<-genbank(gbacc[1],disp="data")
@


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Querying LocusLink}

Finally, we consider interactions with LocusLink to obtain sequence
and other biological information on probe sets of interest.

<<locuslink,eval=FALSE>>=
llid<-getLL(fivegenes, data="hu6800")
locuslinkByID(llid)
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{HTML reports}

The function \texttt{ll.htlmpage} generates HTML reports, with one row per gene and a clickable entry which opens the LocusLink webpage for that gene. Additional information can be displayed using the 
\texttt{othernames} argument as shown below. Such HTML reports provide an
easy way to communicate results with other researchers working on the same analysis.

<<ll.htmlpage>>=
tengenes <- affyID[1:10]
llid<-getLL(tengenes, data="hu6800")
symb<-getSYMBOL(tengenes, data="hu6800")
map<-multiget(tengenes,hu6800MAP)
res<-data.frame(tengenes,cbind(unlist(symb),unlist(map)))
names(res)<-c("Affy ID","Gene symbol", "Chromosomal location")
ll.htmlpage(llid, filename="ll.html",
  title="HTML report for 10 genes",
  othernames=res, table.head=c("LocusID",names(res)),
  table.center = TRUE)
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Plotting genomic data}

The classes \texttt{chromLoc} and \texttt{chromLocation} are used to
keep track of location information for a single gene and a set of
genes (entire genome), respectively. The object
\texttt{hu6800ChrClass} from the \texttt{EMBO03} course package is an
instance of the class \texttt{chromLocation} for the Affymetrix HU6800
chip. It was constructed using the "byChroms" vignette in the \texttt{geneplotter} package.

<<chromLoc>>=
slotNames("chromLoc")
data(hu6800ChrClass)
species(hu6800ChrClass)
chromNames(hu6800ChrClass)
@

The \texttt{geneplotter} package provides two main functions for
plotting genomic data: \texttt{cPlot} and \texttt{alongChrom}. These
functions operate on instances of the class \texttt{chromLocation}.
With \texttt{cPlot}, we can render all chromosomes the same length or
we can scale them by their relative lengths. Note that the \texttt{mva} package
in the next release of R (1.7.0) will have new functions for heat maps and dendrograms.

\newpage

<<cPlot1,fig=TRUE>>=
cPlot(hu6800ChrClass,scale="relative")
@

<<cPlot2,fig=TRUE>>=
cPlot(hu6800ChrClass,scale="max")
@

<<alongChrom1,fig=TRUE>>=
data(golubTrain)
cols<-as.numeric(pData(golubTrain)$ALL.AML)+1
alongChrom(golubTrain,chrom=1,specChrom=hu6800ChrClass,colors=cols)
@

<<alongChrom2,fig=TRUE>>=
alongChrom(golubTrain,chrom=22,specChrom=hu6800ChrClass,plotFormat="image")
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}
