%\VignetteIndexEntry{03 Annotation -- Exercises}
%\VignetteEngine{knitr::knitr}

\documentclass{article}

<<style, eval=TRUE, echo=FALSE, results='asis'>>=
options(max.print=1000)
stopifnot(BiocInstaller::biocVersion() == "2.13")
BiocStyle::latex()
library(knitr)
opts_chunk$set(cache=TRUE, tidy=FALSE)
@ 

<<packages, eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE>>=
suppressPackageStartupMessages({
    library(org.Hs.eg.db)
    library(TxDb.Hsapiens.UCSC.hg19.knownGene)
    library(BSgenome.Hsapiens.UCSC.hg19)
    library(rtracklayer)
    library(biomaRt)
})
@ 

\title{Practical: Annotations}
\author{Martin Morgan (\url{mtmorgan@fhcrc.org})}
\date{3 February 2014}

\newcommand{\Hsap}{\emph{H.~sapiens}}
\newcommand{\Dmel}{\emph{D.~melanogaster}}

\usepackage{theorem}
\newtheorem{Ext}{Exercise}
\newenvironment{Exercise}{
  \renewcommand{\labelenumi}{\alph{enumi}.}\begin{Ext}%
}{\end{Ext}}
\newenvironment{Solution}{%
  \noindent\textbf{Solution:}\renewcommand{\labelenumi}{\alph{enumi}.}%
}{\bigskip}

\setlength{\abovecaptionskip}{6pt}
\setlength{\belowcaptionskip}{6pt}

\begin{document}

\maketitle
\tableofcontents

\section{Gene annotation}

\subsection{Data packages}

Organism-level (`org') packages contain mappings between a central
identifier (e.g., Entrez gene ids) and other identifiers (e.g. GenBank
or Uniprot accession number, RefSeq id, etc.).  The name of an org
package is always of the form \texttt{org.<Sp>.<id>.db}
(e.g. \Biocannopkg{org.Hs.eg.db}) where \texttt{<Sp>} is a 2-letter
abbreviation of the organism (e.g. \texttt{Hs} for \emph{Homo
  spaiens}) and \texttt{<id>} is an abbreviation (in lower-case)
describing the type of central identifier (e.g. \texttt{eg} for ENTREZ
gene identifiers).  The ``How to use the `.db' annotation packages''
vignette in the \Biocpkg{AnnotationDbi} package (org packages are only
one type of ``.db'' annotation packages) is a key reference.  The
`.db' and most other \Bioconductor{} annotation packages are updated
every 6 months.

Annotation packages usually contain an object named after the package
itself.  These objects are collectively called \Rclass{AnnotationDb}
objects, with more specific classes named \Rclass{OrgDb},
\Rclass{ChipDb} or \Rclass{TranscriptDb} objects.  Methods that can be
applied to these objects include \Rfunction{cols}, \Rfunction{keys},
\Rfunction{keytypes} and \Rfunction{select}.  Common operations for
retrieving annotations are summarized in Table~\ref{tab:select-ops}.

\begin{table}
  \centering
  \caption{Common operations for retrieving and manipulating annotations.}
  \label{tab:select-ops}
  \begin{tabular}{lll}
    Category & Function & Description \\
    \hline\noalign{\smallskip}

    Discover & \Rfunction{columns} & List the kinds of columns that can be returned \\
    & \Rfunction{keytypes} & List columns that can be used as keys \\
    & \Rfunction{keys} & List values that can be expected for a given keytype \\
    & \Rfunction{select} & Retrieve annotations matching
    \Rcode{keys}, \Rcode{keytype} and \Rcode{columns} \\
    
    Manipulate & \Rfunction{setdiff}, \Rfunction{union}, \Rfunction{intersect} & Operations on sets \\
    & \Rfunction{duplicated}, \Rfunction{unique} & Mark or remove duplicates \\
    & \Rfunction{\%in\%},  \Rfunction{match} &  Find matches  \\
    & \Rfunction{any}, \Rfunction{all} &  Are any \Rcode{TRUE}?  Are all? \\
    & \Rfunction{merge} & Combine two different \Robject{data.frames} based on shared keys \\

    \Rclass{GRanges*} & \Rfunction{transcripts}, \Rfunction{exons}, \Rfunction{cds} & 
        Features (transcripts, exons, coding sequence) as \Rclass{GRanges}. \\

    & \Rfunction{transcriptsBy} , \Rfunction{exonsBy} & 
    Features group by  gene, transcript, etc., as \Rclass{GRangesList}.\\
    & \Rfunction{cdsBy}\\
    
    \hline
  \end{tabular}
\end{table}

\begin{Exercise}
  This exercise illustrates basic use of the `select' interface to
  annotation packages.

  \begin{enumerate}
  \item What is the name of the org package for \emph{Homo sapiens}?
    Load it.  Display the \Rclass{OrgDb} object for the
    \Biocpkg{org.Hs.eg.db} package.  Use the \Rfunction{keytypes} and
    \Rfunction{columns} methods to discover which sorts of annotations
    can be queried and extracted.
  \item Here are some ENTREZID values.
<<select-setup>>=
egids <- c("3183", "91828", "81537", "4776", "283624", "4053", "85446", 
           "10484", "55701", "1112")
@ 
%% 
These are the most strongly differentially expressed genes from a
subset of an RNA-seq differential expression analysis that you will
encounter later in the course. The biological background is provided
in \cite{pmid23374342}; see the ArrayExpress entry for
\href{https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-1147/}{E-MTAB-1147}. The
data are from chromosome 14 only

Use the ENTREZIDs in the \Rfunction{select} method in such a way that
you extract the SYMBOL (gene symbol) and GENENAME information for
each. To what extent do the differentially expressed genes make
biological sense?
  \end{enumerate}
\end{Exercise}
\begin{Solution}
  The `org' package for humans (\emph{Homo sapiens}) is
  \Biocannopkg{org.Hs.eg.db}.  Load the \Biocannopkg{org.Hs.eg.db}
  package. 
<<load-org>>=
library(org.Hs.eg.db)
@ 
Discover the key types and columns  in the annotation package.
<<select>>=
keytypes(org.Hs.eg.db)
columns(org.Hs.eg.db)
cols <- c("SYMBOL", "GENENAME")
select(org.Hs.eg.db, keys=egids, columns=cols, keytype="ENTREZID")
@
\end{Solution}

\begin{Exercise}
  This exercise annotates a larger selection of the differential
  expression results, merging the statistics of differential
  expression with annotation.
  \begin{enumerate}
  \item Read the comma-separate value file
    \texttt{E-MTAB-1147-toptable.csv} in to \R{} using
    \Rfunction{read.csv}; include the argument \Rcode{row.names=1} to
    name the rows of the data (these are ENTREZIDs). Perform basic
    \R{} operations to discover the object dimensions, and to view the
    head (first 6 rows) of the data frame. The columns of this data
    frame are statistics of differential representation, to be
    discussed later in the course.
  \item Use \Rfunction{select} to map \emph{all} ENTREZIDs
    (\Rfunction{rownames}) to their CHR (chromosome), SYMBOL and
    GENENAME. Verify that in fact all genes are located on chromosome 14.
  \item Use \Rfunction{merge} to add the SYMBOL and GENENAME
    annotations to the differential expression statistics.
  \item Print out the 10 rows of the merged data frame with highest
    absolute log fold change. To do this, you'll need to (a) take the
    absolute value of the log fold change column; (b) determine the
    order, decreasing from largest to smallest, of the absolute
    values; (c) select the rows of the merged data frame in decreasing
    order of fold change; and (d) select, using \Rfunction{head}, just
    the first 10 rows of the ordered data frame.
  \end{enumerate}
\end{Exercise}

\begin{Solution}
  \begin{enumerate}
  \item Find the file \texttt{E-MTAB-1147-toptable.csv} of differentially expressed genes.
<<merge-setup, echo=FALSE>>=
fl <- system.file("extdata", "E-MTAB-1147-toptable.csv", package="SummerX")
@ 
<<merge-user-setup, eval=FALSE>>=
fl <- file.choose()
@ 
%%
Verify that the file exists!
<<file-exists>>=
file.exists(fl)
@ 
%% 
Input the file, specifying that the the first column should be used
for row names; look at basic properties of the data.
<<merge>>=
csv <- read.csv(fl, row.names=1)
class(csv)                              # data.frame
dim(csv)                                # 528 genes x 6 columns
head(csv)                               # log (base 2) fold change, # adjusted P-values
@
%% 
\item Select CHR, SYMBOL, and GENENAME annotations for each ENTREZID,
  and verify that all genes are from CHR 14.
<<anno>>=
cols <- c("CHR", "SYMBOL", "GENENAME")
anno <- select(org.Hs.eg.db, rownames(csv), cols)
class(anno)
dim(anno)
head(anno)
all(anno$CHR %in% "14")                 # all on CHR 14?
@ 
%% 
\item Merge the differential expression and annotation data frames,
  specifying that rows are to be matched based on row names in the
  first data frame (\Rcode{by.x=0}) and by the ENTREZID in the second
  data frame (\Rcode{by.y="ENTREZID"}). Verify that the result is as
  expected using standard \R{} commands.
<<merge-annotated>>=
annotated <- merge(csv, anno, by.x=0, by.y="ENTREZID")
class(annotated)
dim(annotated)
head(annotated)
@ 
\item Re-organize the annotated data frame by (a) taking the absolute
  value of the log fold change column and (b) determine the order,
  decreasing from largest to smallest, of the absolute values of the
  fold change. Here is a 'one-liner'; explain it to your
  neighbor. What is \Rcode{o}?
<<order-annotated>>=
o <- order(abs(annotated$log2FoldChange), decreasing=TRUE)
@ 
%% 
Finally, order the annotated data frame and view the 10 genes with
largest differential expression.
<<head-annotated>>=
annotated[head(o),]
@ 
%% 
Verify that the \emph{most} differentially expressed genes are at the
top of your table!
\end{enumerate}
\end{Solution}

\subsection{Internet resources}

A short summary of select \Bioconductor{} packages enabling web-based
queries is in Table~\ref{tab:webannoservices}.

\begin{table}
  \centering
  \caption{Selected packages querying web-based annotation services.}
  \label{tab:webannoservices}
  \begin{tabular}{ll}
    Package & Description\\
    \hline\noalign{\smallskip}
    \Biocpkg{AnnotationHub} & Ensembl, Encode, dbSNP, UCSC data objects \\
    \Biocpkg{biomaRt} & \url{http://biomart.org}, Ensembl and other annotations\\
    \Biocpkg{PSICQUIC} & \url{https://code.google.com/p/psicquic.org}, protein interactions \\
    \Biocpkg{uniprot.ws} & \url{http://uniprot.org}, protein annotations\\
    \Biocpkg{KEGGREST} & \url{http://www.genome.jp/kegg}, KEGG pathways\\
    \Biocpkg{SRAdb} & \url{http://www.ncbi.nlm.nih.gov/sra}, sequencing experiments.\\
    \Biocpkg{rtracklayer} & \url{http://genome.ucsc.edu}, genome tracks.\\
    \Biocpkg{GEOquery} & \url{http://www.ncbi.nlm.nih.gov/geo/}, array and other data\\
    \Biocpkg{ArrayExpress} & \url{http://www.ebi.ac.uk/arrayexpress/}, array and other data\\
    \hline
  \end{tabular}
\end{table}

\paragraph{Using biomaRt}

The \Biocpkg{biomaRt} package offers access to the online
\href{http://www.biomart.org}{biomart} resource. this consists of
several data base resources, referred to as `marts'.  Each mart allows
access to multiple data sets; the \Biocpkg{biomaRt} package provides
methods for mart and data set discovery, and a standard method
\Rfunction{getBM} to retrieve data.

\begin{Exercise}
  \warning{This exercise requires INTERNET ACCESS}

  \begin{enumerate}
  \item Load the \Biocpkg{biomaRt} package and list the available
    marts.  Choose the \emph{ensembl} mart and list the datasets for
    that mart.  Set up a mart to use the \emph{ensembl} mart and the
    \emph{hsapiens\_gene\_ensembl} dataset.
  \item A \Biocpkg{biomaRt} dataset can be accessed via
    \Rfunction{getBM}. In addition to the mart to be accessed, this
    function takes filters and attributes as arguments.  Use
    \Rfunction{filterOptions} and \Rfunction{listAttributes} to
    discover values for these arguments.  Call \Rfunction{getBM} using
    filters and attributes of your choosing.
  \end{enumerate}
\end{Exercise}

\begin{Solution}
<<biomaRt1, eval=FALSE, results="hide">>=
## NEEDS INTERNET ACCESS !!
library(biomaRt)
head(listMarts(), 3)                      ## list the marts
head(listDatasets(useMart("ensembl")), 3) ## mart datasets
ensembl <-                                ## fully specified mart
    useMart("ensembl", dataset = "hsapiens_gene_ensembl")

head(listFilters(ensembl), 3)             ## filters
myFilter <- "chromosome_name"
head(filterOptions(myFilter, ensembl), 3) ## return values
myValues <- c("21", "22")
head(listAttributes(ensembl), 3)          ## attributes
myAttributes <- c("ensembl_gene_id","chromosome_name")

## assemble and query the mart
res <- getBM(attributes =  myAttributes, filters =  myFilter,
             values =  myValues, mart = ensembl)
@
Use \Rcode{head(res)} to see the results.
\end{Solution}

\begin{Exercise}
  As an optional exercise, annotate the genes that are differentially
  expressed in the DESeq2 laboratory, e.g., find the \texttt{GENENAME}
  associated with the five most differentially expressed genes. Do
  these make biological sense? Can you \Rfunction{merge} the
  annotation results with the `top table' results to provide a
  statistically and biologically informative summary?
\end{Exercise}

\paragraph{Using PSICQUIC}
\href{https://code.google.com/p/psicquic/}{PSICQUIC} is a really
useful effort to provide programmatic access to molecular interaction
data bases.  The \Biocpkg{PSICQUIC} package provides an \R{} /
\Bioconductor{} interace to PSICQUIC. 
\begin{enumerate}
\item Follow instructions on the PSICQUIC package landing
  page\footnote{\url{http://bioconductor.org/packages/release/bioc/html/PSICQUIC.html}}
  to install the package.
\item Work through section 2 `Quick Start' of the PSICQUIC
  vignette\footnote{\url{http://bioconductor.org/packages/release/bioc/vignettes/PSICQUIC/inst/doc/PSICQUIC.pdf}}, discovering documented interactions between Myc and TP53.
\item If interested, explore more of the PSICQUIC vignette. Save
  yourself typing by using the R script from the package landing page.
\end{enumerate}



\section{Genome annotation}

There are a diversity of packages and classes available for
representing large genomes. Several include:
\begin{description}
\item [\Rpackage{TxDb.*}] For transcript and other genome / coordinate
  annotation.
\item [\Biocpkg{BSgenome}] For whole-genome representation. See
  \Rfunction{available.packages} for pre-packaged genomes, and the
  vignette `How to forge a BSgenome data package' in the
\item [\Biocannopkg{Homo.sapiens}] For integrating \Rpackage{TxDb*} and
  \Rpackage{org.*} packages.
\item [\Rpackage{SNPlocs.*}] For model organism SNP locations derived
  from dbSNP.
\item [\Rfunction{FaFile}] (\Biocpkg{Rsamtools}) for accessing indexed
  FASTA files.
\item [\Rpackage{SIFT.*}, \Rpackage{PolyPhen}, \Rpackage{ensemblVEP}]
  Variant effect scores.
\end{description}

\subsection{Transcript annotation packages}

Genome-centric packages are very useful for annotations involving
genomic coordinates. It is straight-forward, for instance, to discover
the coordinates of coding sequences in regions of interest, and from
these retrieve corresponding DNA or protein coding sequences. Other
examples of the types of operations that are easy to perform with
genome-centric annotations include defining regions of interest for
counting aligned reads in RNA-seq experiments and retrieving DNA
sequences underlying regions of interest in ChIP-seq analysis, e.g.,
for motif characterization.

\begin{Exercise}
  This exercise uses annotation packages to go from gene identifiers
  to coding sequences.
  \begin{enumerate}
  \item Map from an informal gene SYMBOL, e.g., BRCA1, to ENTREZID
    gene identifiers using the \Biocannopkg{org.Hs.eg.db} package and
    the \Rfunction{select} function, use the
    \Biocannopkg{TxDb.Hsapiens.UCSC.hg19.knownGene} package and a
    second map to go from ENTREZID to TXNAME.
  \item Extract the coding sequence grouped by transcript using the
    \Biocannopkg{TxDb.Hsapiens.UCSC.hg19.knownGene} package and
    \Rfunction{cdsBy} function; select just those transcripts we are
    interested in.
  \item Retrieve the nucleotide sequence from the
    \Biocannopkg{BSgenome.Hsapiens.UCSC.hg19} package using the
    function \Rfunction{extractTranscriptsFromGenome}.
  \item Verify that the coding sequences are all multiples of 3, and
    translate from nucleotide to amino acid sequence.
  \end{enumerate}
\end{Exercise}

\begin{Solution}
  Map from gene SYMBOL to ENTREZID, and from ENTREZID to TXNAME
<<SYMBOL-to-ENTREZID>>=
library(org.Hs.eg.db)
egid <- select(org.Hs.eg.db, "BRCA1", "ENTREZID", "SYMBOL")$ENTREZID
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
egToTx <- select(txdb, egid, "TXNAME", "GENEID")
@ 
\noindent Extract the releveant coding sequence, grouped by transcript
<<cdsBy>>=
cdsByTx <- cdsBy(txdb, "tx", use.names=TRUE)[egToTx$TXNAME]
@ 
\noindent Retrieve the sequence
<<getsequence>>=
library(BSgenome.Hsapiens.UCSC.hg19)
txx <- extractTranscriptsFromGenome(Hsapiens, cdsByTx)
@ 
\noindent Translate to amino acid sequence
<<translate>>=
all(width(txx) %% 3 == 0)  # sanity check
translate(txx)             # amino acid sequence
@ 
\end{Solution}

\subsection{\Rpackage{rtracklayer}}
\label{subsec:rtracklayer}

The \Biocpkg{rtracklayer} package allows us to query the UCSC genome
browser, as well as providing \Rfunction{import} and
\Rfunction{export} functions for common annotation file formats like
GFF, GTF, and BED.

\begin{Exercise}
  \warning{This exercise requires INTERNET ACCESS}

  Here we use \Biocpkg{rtracklayer} to retrieve estrogen receptor
  binding sites identified across cell lines in the ENCODE project. We
  focus on binding sites in the vicinity of a particularly interesting
  region of interest.
  \begin{enumerate}
  \item Define our region of interest by creating a \Rclass{GRanges}
    instance with appropriate genomic coordinates. Our region
    corresponds to 10Mb up- and down-stream of a particular gene.
  \item Create a session for the UCSC genome browser
  \item Query the UCSC genome browser for ENCODE estrogen receptor ERalpha$_a$
    transcription marks; identifying the appropriate track, table, and
    transcription factor requires biological knowledge and detective
    work.
  \item Visualize the location of the binding sites and their scores;
    annotate the mid-point of the region of interest.
  \end{enumerate}
\end{Exercise}

\begin{Solution}
Define the region of interest
<<rtracklayer-roi>>=
roi <- GRanges("chr10", IRanges(92106877, 112106876, names="ENSG00000099194"))
@ 
\noindent Create a session
<<rtracklayer-session>>=
library(rtracklayer) 
session <- browserSession()
@ 
\noindent Query the UCSC for a particular track, table, and
transcription factor, in our region of interest
<<rtracklayer-marks>>=
trackName <- "wgEncodeRegTfbsClusteredV2"
tableName <- "wgEncodeRegTfbsClusteredV2"
trFactor <- "ERalpha_a"
ucscTable <- getTable(ucscTableQuery(session, track=trackName,
    range=roi, table=tableName, name=trFactor))
@ 
\noindent Visualize the result
<<rtracklayer-plot, fig.height=3>>=
plot(score ~ chromStart, ucscTable, pch="+")
abline(v=start(roi) + (end(roi) - start(roi) + 1) / 2, col="blue")
@ 
\end{Solution}

\appendix

\bibliography{SummerX}

\end{document}