\documentclass[a4paper]{article} \usepackage{Rd} % \VignetteIndexEntry{Annotation mapping functions} % \VignetteDepends{AnnotationDbi} \title{Annotation mapping functions} \author{Stefan McKinnon Edwards \\ Faculty of Agricultural Sciences, Aarhus University} \date{\today} \begin{document} \maketitle \tableofcontents \section{Introduction} The Bioconductor\cite{BioC} contains over 500 packages for annotational data (as of release 2.7), containing all gene data, chip data for microarrays or homology data. These packages can be used to map between different types of annotations, such as RefSeq\cite{RefSeq}, Entrez or Ensembl\cite{Hubbard01012009}, but also pathways such as KEGG\cite{Kanehisa}\cite{Kanehisa01012000}\cite{Kanehisa01012010} og Gene Ontology\cite{GO}. This package contains functions for mapping between different annotations using the Bioconductors annotations packages and can in most cases be done with a single line of code. \subsection{Types of Annotation Packages} Bioconductors website on the annotation packages\footnote{http://www.bioconductor.org/help/workflows/annotation-data/} lists five different types of annotation packages. Here, we list some of them: \begin{itemize} \item org.``XX''.``yy''.db - organism annotation packages containing all the gene data for en entire organism, where ``XX'' is the abbreviation for Genus and species. ``yy'' is the \emph{central ID} that binds the data together, such as Entrez gene (eg). \item KEGG.db and GO.db containing systems biology data / pathways. \item hom.``XX''.inp.db - homology packages mapping gene between an organism and 35 other. \item Chip annotation packages for accessing data for e.g. Affymetrix chips. \end{itemize} \subsection{Why this package?} Querying the annotation packages for data is easy, as demonstrated below: <<>>= library(org.Bt.eg.db) genes <- c('280705', '280706', '100327208') symbols <- org.Bt.egSYMBOL[genes] toTable(symbols) @ However, if we were to query for a non-existing gene, \begin{Schunk} \begin{Sinput} > org.Bt.egSYMBOL[c("280705", "280706", "100327208", "123")] \end{Sinput} \begin{Soutput} Error in .checkKeys(value, Lkeys(x), x@ifnotfound) : value for "123" not found \end{Soutput} \end{Schunk} or try to map from a symbol to e.g. RefSeq, things become a bit more difficult: <<>>= symbols <- c('BCSE','TF','TG') entrez <- org.Bt.egSYMBOL2EG[symbols] entrez <- toTable(entrez)[,'gene_id'] entrez refseq <- org.Bt.egREFSEQ[entrez] toTable(refseq) @ Note here, that the central ID in \code{org.Bt.eg.db} is \emph{Entrez Gene ID} as indicated by the ''eg'' in the package name. This means that all annotational data is linked together by an Entrez Gene ID. Therefore to map from e.g. symbols to RefSeq, it is neccesary to map to Entrez first and then to RefSeq. The consequence of this is that elements can be lost in the mapping. With this package, it can be done with one-liners: <<>>= library(AnnotationFuncs) translate(c(280705, 280706, 100327208, 123), org.Bt.egSYMBOL) translate(c('BCSE','TF','TG'), from=org.Bt.egSYMBOL2EG, to=org.Bt.egREFSEQ) @ Note that the elements that could not be mapped to the end product are removed (by option). \section{Usage guide} First things first, we will start by describing the most usefull of functions: \code{translate}. After this we will cover the other functions that are intended for Gene Ontology pathways (GO) and RefSeq. Please note, that throughout this chapter, we use \code{org.Bt.eg.db}, but the annotation package is more or less interchangable with the other annotation packages from Bioconductor. However, when we refer to the annotation types (e.g. \code{org.Bt.egCHR} or \code{org.Bt.egREFSEQ}), the package cannot be directly replaced by another, so you must check which types are available for the package. \subsection{translate} To map between different types of annotation data, you must first choose an annotation package for the appropiate organism. For cows there is the \code{org.Bt.eg.db} that maps between different identifers or \code{bovine.db} that maps the probes on an Affymatrix bovine chip. For humans there is also \code{org.Hs.eg.db}, \code{hgu95av2.db} for Affymetrix Human Genome U95 Set annotation data or \code{hom.Hs.inp.db} that maps homologues genes from human to 35 other organims. For \code{org.Bt.eg.db} there is a data object for each set of data, such as mapping from Entrez Gene IDs to chromosome, \code{org.Bt.egCHR}. Some of the objects comes in pairs, such as Entrez and RefSeq there is \code{org.Bt.egREFSEQ} and \code{org.Bt.egREFSEQ2EG}. The latter can also be obtained with \code{revmap(org.Bt.egREFSEQ)}. Note that the mapping is generally from the central ID (Entrez Gene ID in this case) to the end product. We start of by getting some basic information about some Entrez Gene IDs: <<>>= library(AnnotationFuncs) library(org.Bt.eg.db) genes <- c(280705, 280706, 100327208) translate(genes, org.Bt.egGENENAME) translate(genes, org.Bt.egCHR) translate(genes, org.Bt.egENSEMBL) translate(genes, org.Bt.egREFSEQ, remove.missing=FALSE) @ This is however trivial, but not the last example how the non-mapped was included in the result. The first argument is coerced into a character vector of unique entries, so providing a list will do no good (unless you use the function \code{unlist}). Now we would like to map from something else than the central ID to the central ID and to a third annotation: <<>>= symbols <- c("SERPINA1","KERA","CD5") translate(symbols, org.Bt.egSYMBOL2EG) translate(symbols, revmap(org.Bt.egSYMBOL)) @ The two results should be exactly the same. To map to a third annotation, we specify a object for the mapping to the central ID and another object for the mapping from the central ID to the end product: <<>>= translate(symbols, from=org.Bt.egSYMBOL2EG, to=org.Bt.egGENENAME) # As a curiosity, if you specify another organism, you are warned: library(org.Hs.eg.db) translate(symbols, from=org.Bt.egSYMBOL2EG, to=org.Hs.egGENENAME) warnings() @ If you specify an annotation type that can map to multiple entries, you can use the argument \code{reduce} to specify how the result is reduced. The options are \code{all} (default), that does not reduce, \code{first} that only selects the arbitrarily first element and \code{last} that does the same on the last element in the set. <<>>= symbols <- c("SERPINA1","KERA","CD5") translate(symbols, org.Bt.egSYMBOL2EG, org.Bt.egREFSEQ) translate(symbols, org.Bt.egSYMBOL2EG, org.Bt.egREFSEQ, reduce='first') translate(symbols, org.Bt.egSYMBOL2EG, org.Bt.egREFSEQ, reduce='last') @ And finally, if you for some reason needed the result on matrix form, there is the argument \code{return.list}: <<>>= translate(symbols, org.Bt.egSYMBOL2EG, org.Bt.egREFSEQ, return.list=FALSE) @ \subsubsection{Combining with other lists and using lapply and sapply} Sometimes your input might not be a vector of values, but a list object, mapping a grouping of genes: <<>>= groups <- list('a'=c('ACR','ASM','S','KERA'), 'IL'=c('IL1','IL2','IL3','IL10'), 'bwahh'=c('ACR','SERPINA1','IL1','IL10','CD5')) @ For these three groups, you want the entrez gene id for each gene. For this we have two options: \code{lapply/sapply} and \code{mapLists}. When using \code{lapply} or \code{sapply}, we have added an argument, \code{simplify} which reduces the result from translate to a simple character vector: <<>>= lapply(groups, translate, from=org.Bt.egSYMBOL2EG, simplify=TRUE) @ Of cause, if there are many recurring IDs, we do not want redundant look ups for the genes, so we do it a bit smarter and just start by asking for a translation of all the genes. This just leaves us with the task of mapping the groups to the result. <<>>= symbols <- unlist(groups, use.names=FALSE) ent <- translate(symbols, org.Bt.egSYMBOL2EG) ent mapLists(groups, ent) @ Tip: When unlisting a list (\code{unlist}), set the argument \code{use.names} to \code{FALSE}, as this will greatly speed up the process. But, naturally, only if you do not need the names. Tip 2: \code{unlist} concatenate the names with a number for each entry, with might make a mess, if the names already end with a number. There is an alternative, \code{unlist\textbf{2}} in the package \code{AnnotationDbi} which preserves the names. \subsection{RefSeq} As noted in the previous section, a single gene can map to several RefSeq identifiers. The format of the RefSeq identifier differs depending on the type, but they can be recognized on their prefix. For instance, RefSeq identifiers starting with `NM\_' are mature mRNA transcripts and `XP\_' are model proteins. A thorough explanation can be found by executing \code{?org.Bt.egREFSEQ} or on NCBIs website for RefSeq\footnote{http://www.ncbi.nlm.nih.gov/projects/RefSeq/key.html}. To be able to sort the result from \code{translate} (or a character vector of RefSeq identifiers in general), the function \code{pickRefSeq}. It uses string comparisons to select the correct identifiers, so options are not limited to those shown here. <<>>= symbols <- c("SERPINA1","KERA","CD5") refseq <- translate(symbols, from=org.Bt.egSYMBOL2EG, to=org.Bt.egREFSEQ) mRNA <- pickRefSeq(refseq, priorities=c('NM','XM')) # Equals pickRefSeq.mRNA mRNA proteins <- pickRefSeq(refseq, priorities=c('NP','XP')) # Equals pickRefSeq.Protein proteins @ If used in a context where each element must be mapped to exactly one, the argument \code{reduce} can be supplied here as in \code{translate}. \subsection{Gene Ontology} The annotation packages can map to pathways; either Gene Ontology (GO) or KEGG. Using org.Bt.egGO makes things a bit more complicated, as each GO identifier is accompanied by an evidence and category code. We therefore supply the function \code{pickGO} to pick the correct pathway based on category (biological process, molecular function and/or cellular component) and evidence (typically how they were inferred). <>= symbols <- c("SERPINA1","KERA","CD5") GOs <- translate(symbols, from=org.Bt.egSYMBOL2EG, to=org.Bt.egGO) # Pick biological process: pickGO(GOs, category='BP') # Pick only those biological processes for Entrez Gene ID 280730 # which have been inferred from sequence similarity or electronic annotation: pickGO(translate(280730, org.Bt.egGO), category='BP', evidence=c('ISS','IEA')) @ If a GO object is used as the \code{from} argument in \code{translate} for further mapping, the arguments \code{evidence} and \code{category} can be applied to \code{translate} to filter the intermediate product before mapping to the end product. The evidence codes can be found by executing \code{?org.Bt.egGO} or \code{getEvidenceCodes()}. \subsection{Finding orthologs / using the INPARANOID data packages} The Stockholm Bioinformatics Centre has compiled some data packages containing orthologs and paralogs (collectively called homologs) \cite{inparanoid6,inparanoid,Reem2001}. These are very useful if you need to map ids from one species to another. At the current time, there are about 38 species in the data packages. These data packages can be downloaded from BioConductor and are named \code{hom.Xx.inp.db}, where Xx is the central species, such as At, Ce, Dm, Dr, Hs, Mm, Rn and Sc (guess what the abbreviation stands for). The central species are of the same concept as for the organism annotation packages mentioned earlier. For the human package, \code{hom.Hs.inp.db}, this means that the ortholog mappings are between Humans and cattle, humans and mice and so on. The builtin method of performing a ortholog mapping is as such: <<>>=@ library(hom.Hs.inp.db) submap <- hom.Hs.inpBOSTA[c('ENSP00000356224','ENSP00000384582','ENSP00000364403')] toTable(submap) @ This one was quick. But if you query for lots of ids (e.g. the entire set), it takes a long time, and even longer if you are querying on the reversed map. Furthermore, if you included values that could not be mapped, an error is returned: \begin{Schunk} \begin{Sinput} > hom.Hs.inpBOSTA['bwahh'] \end{Sinput} \begin{Soutput} Error in .checkKeys(value, Lkeys(x), x@ifnotfound) : value for "bwahh" not found \end{Soutput} \end{Schunk} Note: This can be resolved by using \code{AnnotationDbi::mget(values, hom.Hs.inpBOSTA, \textbf{ifnotfound=NA})} which also returns a list object. We therefore present the function \code{getOrthologs}. It is faster and more powerful. Faster, as retrieving all mapped ids from bovine to human took approx. 11 minutes with the method shown above, while \code{getOrthologs} did it in just 1.17 seconds. More powerful, as it can translate the input values \textit{before} mapping \textit{and} the mapped ids \textit{after} mapping (but not automatically - you have to tell it which translations to use). To reproduce the above (in list form): <<>>=@ getOrthologs(c('ENSP00000356224','ENSP00000384582','ENSP00000364403'), hom.Hs.inpBOSTA, 'BOSTA') @ Unfortunately we require the INPARANOID style genus names ('BOSTA') to be provided as a character string. These are the same as the five last characters in the mapping object, if you use the original names from the package. However, the INPARANOID data packages works primarily with Ensembl protein ids. If for instance your ids are entrez, refseq, etc. you are required to translate them into Ensembl before the mapping. This can be done in the same call to the function: <<>>=@ symbols <- c("SERPINA1","KERA","CD5") getOrthologs(symbols, hom.Hs.inpBOSTA, 'BOSTA', pre.from=org.Hs.egSYMBOL2EG, pre.to=org.Hs.egENSEMBLPROT, post.from=org.Bt.egENSEMBLPROT2EG, post.to=org.Bt.egSYMBOL) @ Are you surprised about the result? The four arguments \code{pre.from}, \code{pre.to}, \code{post.from} and \code{post.to} corresponds directly to the \code{from} and \code{to} arguments in \code{translate}, with "pre" for the translation \textit{before} mapping and "post" for the translation \textit{after} mapping. In this example we used the \code{org.Hs.eg.db} to map the input from the symbol to Ensembl. As \code{org.Hs.eg.db} is centred around Entrez gene ids, it requires the two-step translation from symbol to Entrez to Ensembl, and likewise for the post mapping translation. If your input was Entrez gene ids, you would only require the \code{pre.from} argument. \textbf{Note!} Some of the species in the INPARANOID homology data packages are represented with other ids than Ensembl. For instance, \textit{Arabidopsis thaliana} (\code{ARATH}) has ids such as \code{At1g80070.1}, while \textit{Trichoplax adhaerens} (\code{TRIAD}) as ids such as \code{24890} which should not be confused with Entrez gene IDs. \section{What was used to create this document} The version number of R and the packages and their versions that were used to generate this document are listed below. <<>>= sessionInfo() @ \section{References} \begin{thebibliography}{10} \bibitem{inparanoid6} AC~Berglund, E~Sjolund, G~Ostlund, and Sonnhammer ELL. \newblock {InParanoid 6: eukaryotic ortholog clusters with inparalogs}. \newblock {\em Nucleic Acids Research}, 36:D263--266, 2008. \bibitem{GO} The Gene~Ontology Consortium. \newblock Gene ontology: tool for the unification of biology. \newblock {\em Nat Genet}, 25:25--29, 2000. \bibitem{BioC} Robert~C Gentleman, Vincent~J. Carey, Douglas~M. Bates, et~al. \newblock Bioconductor: Open software development for computational biology and bioinformatics. \newblock {\em Genome Biology}, 5:R80, 2004. \bibitem{Hubbard01012009} T.~J.~P. Hubbard, B.~L. Aken, S.~Ayling, B.~Ballester, K.~Beal, E.~Bragin, S.~Brent, Y.~Chen, P.~Clapham, L.~Clarke, G.~Coates, S.~Fairley, S.~Fitzgerald, J.~Fernandez-Banet, L.~Gordon, S.~Graf, S.~Haider, M.~Hammond, R.~Holland, K.~Howe, A.~Jenkinson, N.~Johnson, A.~Kahari, D.~Keefe, S.~Keenan, R.~Kinsella, F.~Kokocinski, E.~Kulesha, D.~Lawson, I.~Longden, K.~Megy, P.~Meidl, B.~Overduin, A.~Parker, B.~Pritchard, D.~Rios, M.~Schuster, G.~Slater, D.~Smedley, W.~Spooner, G.~Spudich, S.~Trevanion, A.~Vilella, J.~Vogel, S.~White, S.~Wilder, A.~Zadissa, E.~Birney, F.~Cunningham, V.~Curwen, R.~Durbin, X.~M. Fernandez-Suarez, J.~Herrero, A.~Kasprzyk, G.~Proctor, J.~Smith, S.~Searle, and P.~Flicek. \newblock {Ensembl 2009}. \newblock {\em Nucleic Acids Research}, 37(suppl 1):D690--D697, 2009. \bibitem{Kanehisa01012000} Minoru Kanehisa and Susumu Goto. \newblock {KEGG: Kyoto Encyclopedia of Genes and Genomes}. \newblock {\em Nucleic Acids Research}, 28(1):27--30, 2000. \bibitem{Kanehisa01012010} Minoru Kanehisa, Susumu Goto, Miho Furumichi, Mao Tanabe, and Mika Hirakawa. \newblock {KEGG for representation and analysis of molecular networks involving diseases and drugs}. \newblock {\em Nucleic Acids Research}, 38(suppl 1):D355--D360, 2010. \bibitem{Kanehisa} Minoru Kanehisa, Susumu Goto, Masahiro Hattori, Kiyoko~F. Aoki-Kinoshita, Masumi Itoh, Shuichi Kawashima, Toshiaki Katayama, Michihiro Araki, and Mika Hirakawa. \newblock {From genomics to chemical genomics: new developments in KEGG}. \newblock {\em Nucleic Acids Research}, 34(suppl 1):D354--D357. \bibitem{inparanoid} Kevin~P O'Brien, Maido Remm, and Erik~L.L Sonnhammer. \newblock {Inparanoid: A Comprehensive Database of Eukaryotic Orthologs}. \newblock {\em NAR}, 33:D476--D480, 2005. \bibitem{RefSeq} Kim~D Pruitt, Tatiana Tatusova, and Donna~R Maglott. \newblock Ncbi reference sequences (refseq): a curated non-redundant sequence database of genomes, transcripts and proteins. \newblock {\em Nucleic Acids Res}, 35 (Database issue(:D61--65, 2007. \bibitem{Reem2001} Maido Remm, Christian E.~V. Storm, and Erik L.~L. Sonnhammer. \newblock {Automatic clustering of orthologs and in-paralogs from pairwise species comparisons}. \newblock {\em J. Mol. Biol.}, 314:1041--1052, 2001. \end{thebibliography} \end{document}