%\VignetteDepends{GO.db, org.Hs.eg.db, biomaRt}
%\VignetteIndexEntry{Annotation Lecture Presentation}
%\VignettePackage{SeattleIntro2010}

\documentclass{beamer}

\usepackage{BioconductorSlides}
\SweaveOpts{keep.source=TRUE}

\title{Annotation}
\author{Marc Carson}
\institute[FHCRC]{Fred Hutchinson Cancer Research Center}
\date{December 9, 2010}


\begin{document}

<<setup, echo=FALSE>>=
options(width = 64)
@ 

\section{Bioconductor Annotation Packages} 
  
\begin{frame}[fragile]
  \frametitle{Bioconductor annotation packages}
  Major types of annotation in Bioconductor.
  
  \structure{AnnotationDbi} packages:
  \begin{itemize}
    \item Organism level: \structure{org.Mm.eg.db}.
    \item Platform level: \structure{hgu133plus2.db}. 
    \item System-biology level: \structure{GO.db} or \structure{KEGG.db}.
  \end{itemize}
  
  \structure{biomaRt}: 
    \begin{itemize}
      \item Query web-based `biomart' resource for genes, sequence,
           SNPs, and etc. 
    \end{itemize}
  
  Other packages: 
    \begin{itemize}
      \item \structure{rtracklayer} -- export to UCSC web browsers.
      \item \structure{GenomicFeatures} -- for transcripts and soon for other tracks.
    \end{itemize}
\end{frame}




\section{AnnotationDbi}

\begin{frame}[fragile]
  \frametitle{AnnotationDbi}
  AnnotationDbi is a software package that enables the package annotations:
\begin{itemize}
  \item Each supported package contains a database.
  \item AnnotationDbi allows access to that data via Bimap objects.
  \item Some databases depend on the databases in other packages.
\end{itemize}
\end{frame}



\begin{frame}[fragile]
  \frametitle{Organism-level annotation}
  There are a number of organism annotation packages with names starting with
  org, e.g., \structure{org.Hs.eg.db} -- genome-wide annotation for human.
{\small
<<organnotation, results=hide>>=
library(org.Hs.eg.db)
org.Hs.eg()
org.Hs.eg_dbInfo()
org.Hs.egGENENAME
org.Hs.eg_dbschema()
@   
}
\end{frame}



\begin{frame}[fragile]
  \frametitle{platform based packages (chip packages)}
  There are a number of platform or chip specific annotation packages
  named after their respective platforms, e.g. \structure{hgu95av2.db}
  annotations for the hgu95av2 Affymetrix platform.
\begin{itemize}
  \item These packages appear to contain a lot of data but it's an illusion.
\end{itemize}
{\small
<<platannotation, results=hide>>=
library(hgu95av2.db)
hgu95av2()
hgu95av2_dbInfo()
hgu95av2GENENAME
hgu95av2_dbschema()
@   
}
\end{frame}




\begin{frame}[fragile]
  \frametitle{Kinds of annotation}
  What can you hope to extract from an annotation package? 
  
  \begin{itemize}
    \item GO IDs: \structure{GO}
    \item KEGG pathway IDs: \structure{PATH}
    \item Gene Symbols: \structure{SYMBOL}
    \item Chromosome start and stop locs: \structure{CHRLOC} and \structure{CHRLOCEND}
    \item Alternate Gene Symbols: \structure{ALIAS}
    \item Associated Pubmed IDs: \structure{PMID}
    \item RefSeq IDs: \structure{REFSEQ}
    \item Unigene IDs: \structure{UNIGENE}
    \item PFAM IDs: \structure{PFAM}
    \item Prosite IDs: \structure{PROSITE}
    \item ENSEMBL IDs: \structure{ENSEMBL}
  \end{itemize}  
\end{frame}



\section{AnnotationDbi Basics}

\begin{frame}[fragile]
  \frametitle{Basic Bimap structure and getters}
  Bimaps create a mapping from one set of keys to another.  And they can easily be searched.
\begin{itemize}
  \item \Rfunction{toTable}: converts a Bimap to a \structure{data.frame}
  \item \Rfunction{get}: pulls data from a Bimap
  \item \Rfunction{mget}: pulls data from a Bimap for multiple things at once
\end{itemize}
{\small
<<basic bimap1, results=hide>>=
head(toTable(hgu95av2SYMBOL))
get("38187_at",hgu95av2SYMBOL)
mget(c("38912_at","38187_at"),hgu95av2SYMBOL,ifnotfound=NA)
@   
}
\end{frame}


\begin{frame}[fragile]
  \frametitle{Reversing and subsetting Bimaps}
  Bimaps can also be reversed and subsetted:
\begin{itemize}
  \item \Rfunction{revmap}: reverses a Bimap
  \item \Rfunction{[[},\Rfunction{[}: Bimaps are subsettable.
\end{itemize}
{\small
<<basic bimap2, results=hide>>=
##revmap
mget(c("NAT1","NAT2"),revmap(hgu95av2SYMBOL),ifnotfound=NA)
##subsetting
head(toTable(hgu95av2SYMBOL[1:3]))
hgu95av2SYMBOL[["1000_at"]]
revmap(hgu95av2SYMBOL)[["MAPK3"]]
##Or you can combine things
toTable(hgu95av2SYMBOL[c("38912_at","38187_at")])
@   
}
\end{frame}


\begin{frame}[fragile]
  \frametitle{using merge, cbind}
  sometimes you will want to combine data
\begin{itemize}
  \item \Rfunction{cbind}: appends multiple columns (blindly by order)
  \item \Rfunction{merge}: "joins" a pair of data.frames based on a key 
\end{itemize}
{\small
<<merge and cbind, results=hide>>=
## 1st lets get some data
symbols = head(toTable(hgu95av2SYMBOL),n=3)
chrlocs = head(toTable(hgu95av2CHRLOC),n=3)
pmids = head(toTable(hgu95av2PMID),n=3)
##cbind
cbind(symbols, pmids, chrlocs)
##merge
merge(symbols, pmids, by.x="probe_id", by.y="probe_id") 
@   
}
\end{frame}



%EXERCISE 1
\begin{frame}[fragile]
  \frametitle{Annotation exercise 1}
  Find the gene symbol, chromosome position and KEGG pathway ID for "1003\_s\_at".
\end{frame}

\begin{frame}[fragile]
\end{frame}

\begin{frame}[fragile]
  \frametitle{Annotation exercise 1 solution}
{\small
<<exer1, results=hide>>=
get("1003_s_at",hgu95av2SYMBOL) 
get("1003_s_at",hgu95av2CHRLOC) 
get("1003_s_at",hgu95av2PATH) 
@
}
\end{frame}





\begin{frame}[fragile]
  \frametitle{Bimap keys}
  Bimaps create a mapping from one set of keys to another.  Some important methods include:
\begin{itemize}
  \item \Rfunction{keys}: centralID for the package (directional)
  \item \Rfunction{Lkeys}: centralID for the package (probe ID or gene ID)
  \item \Rfunction{Rkeys}: centralID for the package (attached data)
\end{itemize}
{\small
<<keys1 bimap, results=hide>>=
keys(hgu95av2SYMBOL[1:4])
Lkeys(hgu95av2SYMBOL[1:4])
Rkeys(hgu95av2SYMBOL)[1:4]
@   
}
\end{frame}


\begin{frame}[fragile]
  \frametitle{More Bimap structure}
  Not all keys have a partner (or are mapped)
\begin{itemize}
  \item \Rfunction{mappedkeys}: which of the key are mapped (directional)
  \item \Rfunction{mappedLkeys} \Rfunction{mappedRkeys}: which keys are mapped (absolute reference)
  \item \Rfunction{count.mappedkeys}: Number of mapped keys (directional)
  \item \Rfunction{count.mappedLkeys},\Rfunction{count.mappedRkeys}: Number of mapped keys (absolute)
\end{itemize}
{\small
<<keys2 Bimap, results=hide>>=
mappedkeys(hgu95av2SYMBOL[1:10])
mappedLkeys(hgu95av2SYMBOL[1:10])
mappedRkeys(hgu95av2SYMBOL[1:10])
count.mappedkeys(hgu95av2SYMBOL[1:100])
count.mappedLkeys(hgu95av2SYMBOL[1:100])
count.mappedRkeys(hgu95av2SYMBOL[1:100])
@   
}
\end{frame}



\begin{frame}[fragile]
  \frametitle{Bimap Conversions}
  How to handle conversions from Bimaps to lists
\begin{itemize}
  \item \Rfunction{as.list}: converts a Bimap to a \structure{list}
  \item \Rfunction{unlist2}: unlists a \structure{list} minus the name-mangling.
\end{itemize}
{\small
<<Bimap conversion, results=hide>>=
as.list(hgu95av2SYMBOL[c("38912_at","38187_at")])
unlist(as.list(hgu95av2SYMBOL[c("38912_at","38187_at")]))
unlist2(as.list(hgu95av2SYMBOL[c("38912_at","38187_at")]))

##but what happens when there are 
##repeating values for the left key?
unlist(as.list(revmap(hgu95av2SYMBOL)[c("STAT1","PTGER3")]))   

##unlist2 can help with this
unlist2(as.list(revmap(hgu95av2SYMBOL)[c("STAT1","PTGER3")]))   
@   
}
\end{frame}



%EXERCISE 2
\begin{frame}[fragile]
  \frametitle{Annotation exercise 2}
  Gene symbols are often recycled by other genes making them a poor
  choice for identifiers.  Using what you have learned, the SYMBOL
  Bimap, along with the \Rfunction{lapply}, \Rfunction{length} and
  \Rfunction{sort} functions, determine which gene symbols in hgu95av2
  are the worst offenders.
\end{frame}

\begin{frame}[fragile]
\end{frame}

\begin{frame}[fragile]
  \frametitle{Annotation exercise 2 solution}
{\small
<<exer2, results=hide>>=
badRank <- lapply(as.list(revmap(hgu95av2SYMBOL)), length)
tail(sort(unlist(badRank)))
@
}
\end{frame}



\begin{frame}[fragile]
  \frametitle{toggleProbes}
  How to hide/unhide ambiguous probes.
\begin{itemize}
  \item \Rfunction{toggleProbes}: hides or displays the probes that have multiple mappings to genes.
\end{itemize}
{\small
<<toggleProbes, results=hide>>=
## How many probes?
dim(hgu95av2ENTREZID)
## Make a mapping with multiple probes exposed 
multi <- toggleProbes(hgu95av2ENTREZID, "all")
## How many probes?
dim(multi)

## Make a mapping with ONLY multiple probes exposed 
multiOnly <- toggleProbes(multi, "multiple")
## How many probes?
dim(multiOnly)

## Then make a mapping with ONLY single mapping probes
singleOnly <- toggleProbes(multiOnly, "single")
## How many probes?
dim(singleOnly)  
@   
}
\end{frame}



%EXERCISE 3
\begin{frame}[fragile]
  \frametitle{Annotation exercise 3}
  Using the knowledge that Entrez IDs are good IDs that can be used to
  define genes uniquely, find the probe that maps to the largest
  number of different genes on hgu95av2.
\end{frame}

\begin{frame}[fragile]
\end{frame}

\begin{frame}[fragile]
  \frametitle{Annotation exercise 3 solution}
{\small
<<exer3, results=hide>>=
mult <- toggleProbes(hgu95av2ENTREZID, "multi")
dim(mult)
multRank <- lapply(as.list(mult), length)
tail(sort(unlist(multRank)))
@
}
\end{frame}





\section{Working with GO.db}


\begin{frame}[fragile]
  \frametitle{GO}
  Some important considerations about the Gene Ontology
\begin{itemize}
  \item GO is actually 3 ontologies (CC, BP and MF)
  \item Each ontology is a directed acyclic graph.
  \item The structure of GO is maintained separarately from the genes
    that these GO IDs are usually used to annotate.
\end{itemize}

\end{frame}




\begin{frame}[fragile]
  \frametitle{GO to gene mappings are stored in other packages}
  Mapping Entrez IDs to GO
  \begin{itemize}
    \item Each ENTREZ ID is associated with up to three GO categories.
    \item The objects returned from an ordinary GO mapping are complex.
  \end{itemize}
{\small
<<entrez to GO, results=hide>>=
go <- org.Hs.egGO[["1000"]]
length(go)
go[[2]]$GOID
go[[2]]$Ontology
@ 
}
\end{frame}


%TODO: give people a heads up on sapply(GOids,names)
%EXERCISE 4
\begin{frame}[fragile]
  \frametitle{Annotation exercise 4}
  Use what they have learned to write a function that gets the
  GOIDs for a particular entrez gene ID, and then returns only their
  GOID as a named vector.  Use \Rfunction{lapply} and \Rfunction(names).
\end{frame}

\begin{frame}[fragile]
\end{frame}

\begin{frame}[fragile]
  \frametitle{Annotation exercise 4 solution}
{\small
<<exer4, results=hide>>=
##get GOIDs from Hs package.
getGOIDs <- function(ids){
    require(org.Hs.eg.db)
    GOs = mget(ids, org.Hs.egGO, ifnotfound=NA)
    unlist2(lapply(GOs,names))
}
##usage example:
getGOIDs(c("1","10"))
@
}
\end{frame}




\begin{frame}[fragile]
  \frametitle{Working with GO.db}
  \begin{itemize}
    \item Encodes the hierarchical structure of GO terms.
    \item The mapping between \texttt{GO} terms and individual genes
      is maintained in the GO mappings from the other packages.
    \item the difference between children and offspring is how many
      generations are represented.  Children only nets you one step down 
      the graph.
  \end{itemize}
{\small  
<<working with GO, results=hide>>=
library(GO.db)
ls("package:GO.db")
## find children
as.list(GOMFCHILDREN["GO:0008094"])
## all the descendants (children, grandchildren, and so on)
as.list(GOMFOFFSPRING["GO:0008094"])
@ 
}  
\end{frame}



\begin{frame}[fragile]
  \frametitle{GO helper methods}
  Using the GO helper methods
  \begin{itemize}
    \item The GO terms are described in detail in the \structure{GOTERM} 
      mapping.
    \item The objects returned by GO.db are GOTerms objects, which can
      make use of helper methods like \Rfunction{GOID}, \Rfunction{Term},
      \Rfunction{Ontology} and \Rfunction{Definition} to retrieve various
      details.
    \item You can also pass GOIDs to these helper methods.
  \end{itemize}
{\small
<<GOTERMS, results=hide>>=
##Mapping a GOTerms object
go <- GOTERM[1]
GOID(go)
Term(go)

##OR you can supply GO IDs
id = c("GO:0007155","GO:0007156")
GOID(id)
Term(id)
Ontology(id)
Definition(id)
@ 
}
\end{frame}


%EXERCISE 5
\begin{frame}[fragile]
  \frametitle{Annotation exercise 5}
  Use what you have learned to write a function that calls your
  previous function so that it can get GOIDs and then returns the GO
  definitions.
\end{frame}

\begin{frame}[fragile]
\end{frame}

\begin{frame}[fragile]
  \frametitle{Annotation exercise 5 solution}
{\small
<<exer5, results=hide>>=
##get GOIDs from Hs package.
getGODefs <- function(ids){
    GOids <- getGOIDs(ids)
    defs <- Definition(GOids)
    names(defs) <- names(GOids)
    defs
}
##usage example:
getGODefs(c("1","10"))
@
}
\end{frame}


\section{Using biomaRt}

\begin{frame}[fragile]
  \frametitle{Using biomaRt}
  Setting up a biomaRt object
  \begin{itemize}
    \item biomaRt offers several "marts" to get data from
    \item each "mart" can have several datasets
    \item the mart object has to be configured with your choices
  \end{itemize}
{\small
<<biomaRt1, results=hide>>=
library(biomaRt)
##list the marts
head(listMarts())
## list the Datasets for a mart
head(listDatasets(useMart("ensembl")))
## now set up the fully qualified mart object
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
@
}
\end{frame}



\begin{frame}[fragile]
  \frametitle{Using biomaRt}
  Choosing biomaRt options
  \begin{itemize}
    \item filters are used to limit the query
    \item values are the values available for a specified filter
    \item attributes are information we want to retrieve
  \end{itemize}  
{\small
<<biomaRt2, results=hide>>=
## need to be able to list filters
head(listFilters(ensembl))
myFilter <- "chromosome_name"
## and list values that you expect back
head(filterOptions(myFilter, ensembl))
myValues <- c("21", "22")
## and list attributes
head(listAttributes(ensembl))
myAttributes <- c("ensembl_gene_id","chromosome_name")
@
}
\end{frame}


\begin{frame}[fragile]
  \frametitle{Using biomaRt}
  Calling \Rfunction{getBM} will extract the information
  \begin{itemize}
  \item \Rfunction{getBM} takes the information we have just shown you
    how to obtain as its parameters.
  \item With the exception of the mart object all these parameters are vectors 
    so you can request multiple values back if they are available etc.   
  \item If you should need to specify multiple filters, then you will need to 
    pass the values parameter in as a list of vectors instead of just a vector.
  \end{itemize}
{\small
<<getBM, results=hide, eval=FALSE>>=
## then you can assemble a query
res <- getBM(attributes =  myAttributes, 
             filters =  myFilter,
             values =  myValues,
             mart = ensembl)
head(res)
@
}
\end{frame}




%EXERCISE 6
\begin{frame}[fragile]
  \frametitle{Annotation exercise 6}
  Use what you have learned about biomaRt to find the gene symbol and name for 
  the entrez gene IDs 1, 10 and 100.
\end{frame}

\begin{frame}[fragile]
\end{frame}

\begin{frame}[fragile]
  \frametitle{Annotation exercise 6 solution}
{\small
<<exer6, results=hide>>=
res <- getBM(attributes =  c("entrezgene","hgnc_symbol"), 
             filters =  "entrezgene",
             values =  c("1","10","100"),
             mart = ensembl)
head(res)
@
}
\end{frame}




%EXERCISE 7
\begin{frame}[fragile]
  \frametitle{Annotation exercise 7}
  Use what you have learned to add annotations to the dataset that you
  used in the earlier session with gene symbols.  The code below will
  re-load the dataset into your R session.  The annotations for this
  package are found in the \structure{hgu95av2.db} chip package.
{\small
<<exer7question>>=
load(system.file("data", "result.rda",
                 package = "SeattleIntro2010"))
@
}
\end{frame}

\begin{frame}[fragile]
\end{frame}

\begin{frame}[fragile]
  \frametitle{Annotation exercise 7 solution}
{\small
<<exer7, results=hide>>=
ids <- result[,1]
library(hgu95av2.db)
merge(result, toTable(hgu95av2SYMBOL[ids]), by.x="ID", by.y="probe_id")
@
}
\end{frame}






%EXERCISE 8
\begin{frame}[fragile]
  \frametitle{Annotation exercise 8 (bonus round)}
  Use what you have learned about Annotations to come up with a more useful exercise to YOU, and a solution.  We will be here to help.
\end{frame}




%% Below is an example to create a slide as well as other features of
%% beamer-style vignette such as block and columns. If your slides
%% contain R chunck code, alway use \begain{frame}[fragile] to begin a slide
%% frame. Otherwise, you can use us \begin{frame}.
%% One more thing, \section{} allows the beamer-style vignette to create a
%% table of content, which will be placed at the first slide. 
% \section{AnnotationDbi}
% \begin{frame}[fragile]
%   \frametitle{haha}
%     An example of creating block 
%     \begin{block}{I am a block}
%       content of the block
%     \end{block}    
    
%     An example of creating columns
%     \begin{columns}
%       \begin{column}{0.5\textwidth}
%         column 1
%       \end{column}
%       \begin{column}{0.5\textwidth}
%         column 2
%       \end{column}
%     \end{columns}
    
% \end{frame}


\end{document}


%% TODO: configure emacs

