%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{An Introduction to Biofiles}
%\VignetteKeyword{DNA, Sequence, annotation, NCBI}
%\VignettePackage{biofiles}

\documentclass[10pt]{article}
\usepackage[utf8]{inputenc}

<<style-knitr, eval=TRUE, echo=FALSE, results="asis">>=
BiocStyle::latex(width = 80)
@

\title{Using the biofiles package}
\author{Gerhard Schöfl \\
        Dept. of Bioinformatics \\
        DKMS Life Science Lab,  \\
        Dresden, Germany}

\begin{document}
\maketitle
\date{\today}

\tableofcontents


\section{The \Rpackage{biofiles} package}

The \Rpackage{biofiles} package is an \R{} package for interfacing with \href{http://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html}{GenBank or GenPept} flat file records. The GenBank sequence format is a rich data format for storing sequences and associated annotations. \Rpackage{biofiles} includes utilities for parsing and writing GenBank files, and useful methods for interacting with the annotation and sequence data contained in these files.


\subsection{Classes in the \Rpackage{biofiles} package}

The classes in \Rpackage{biofiles} attempt to encapsulate the information contained in GenBank files as faithfully as possible. The design is based heavily on the GenBank/Embl feature tables. The central class for storing feature annotation information is the S4 class \Rclass{gbFeature}. Essentially, a \Robject{gbFeature} object contains a feature \emph{key}, a short description of the type of feature (e.g. ``CDS'' or ``gene''); the \emph{location} of the feature on the sequence string, and \emph{qualifiers}, additional information about the feature as key-value-pairs (e.g.: \emph{locus\_tag}, \emph{gene}, \emph{product}, or \emph{db\_xref} are frequently used qualifiers). A collection of \Rclass{gbFeature}s is organised in a \Rclass{gbFeatureTable}.

A \Rclass{gbFeatureTable} together with metadata and sequence information stored as \Rclass{XStringSet} objects (\Biocpkg{Biostrings}) make up a \Rclass{gbRecord}. A collection of \Rclass{gbRecord}s can be placed into a \Rclass{gbRecordList}.


\subsection{Importing sequence and annotation data}

\Rclass{gbRecord}s are generated by importing ``gb'' or ``gbk'' files. These files can be either downloaded and imported, or they can be directly retrieved and imported using the NCBI's public ENTREZ database server using the functionality provided the \CRANpkg{reutils} package. %Alternatively whole genomes can be imported directly from NCBI using the \Rfunction{genomeRecordFromNCBI} function.

As an example \Rpackage{biofiles} includes the annotation file for the mitochondrial genome of \emph{Saccharomyces cerevisiae}:

<<mito-path,>>=
mito.path <- system.file("extdata", "S_cerevisiae_mito.gb", package="biofiles")
@

<<read-rds, eval=TRUE, echo=FALSE>>=
mito.rds <- system.file("extdata", "NC_001224.rds", package="biofiles")
mito <- biofiles::loadRecord(mito.rds)
@

To read a GenBank file into \R{} use the command
<<read-mito, eval=FALSE, echo=TRUE>>=
mito <- biofiles::gbRecord(mito.path)
mito
@

<<show.rds, eval=TRUE, echo=FALSE>>=
mito
@

\subsection{Saving and loading \Rclass{gbRecord}s}

\Rclass{gbRecord} objects can be stored efficiently by writing the content to a compressed ``.rds'' file using \Rfunction{saveRecord} and re-imported using \Rclass{loadRecord}. Alternatively, \Rclass{gbRecord} objects can be exported as a GenBank flat file using \Rfunction{write.GenBank} or in GenBank feature table format using \Rfunction{write.FeatureTable}.
<<store-mito, eval=FALSE, echo=TRUE>>=
biofiles::saveRecord(mito)
rm(mito)
mito <- biofiles::loadRecord("NC_001224.rds")
biofiles::summary(mito, n = 3)
@

<<display-store-mito, eval=TRUE, echo=FALSE>>=
biofiles::summary(mito, n = 3)
@

\section{Summarising a Genbank record}

A quick overview over the content of a \Rclass{gbRecord} is provided by \Rfunction{summary}:
<<summarise>>=
biofiles::summary(mito)
@

The frequency of features/qualifiers can be quickly tablulated:
<<tabulate>>=
biofiles::qualifTable(mito)
biofiles::featureTable(mito)
@

\section{Exctracting subsets of the data}

There are a number of accessor methods for retrieving data from GenBank records.

We can access each component of the header, e.g.:
<<extract-header>>=
biofiles::getAccession(mito)
biofiles::getDefinition(mito)
biofiles::getGeneID(mito)
biofiles::getOrganism(mito)
biofiles::getLength(mito)
biofiles::getComment(mito)
@

We can extract the DNA sequence as a whole as a \Rclass{DNAStringSet}:
<<extract-sequence>>=
biofiles::getSequence(mito)
@

We can extract the Feature table:
<<extract-ft>>=
biofiles::ft(mito)
@

\subsection{Filtering the feature table}

The bulk of the data we are likely interested in is contained in the feature table. Features of interest can be filtered from a \Rclass{gbRecord} or a \Rclass{gbFeatureTable} using the function \Rfunction{filter}. For instence, to extract all ``CDS'' write:
<<extract-key>>=
cds <- biofiles::filter(mito, key = "CDS")
biofiles::summary(cds[1:2])
@

A handy shortcut for filtering by feature key is provided by using the \Rfunction{[} operator:
<<extact-key2>>=
cds <- mito["CDS"]
biofiles::summary(cds[3:4])
@

We can also filter features based on range:
<<extract-range>>=
f10000 <- biofiles::filter(mito, range = "..10000")
biofiles::summary(f10000)
@

Or, perhaps more usefully, based on arbitary qualifiers matching a regular expression. To find all CDSs for containing \emph{cytochrome b} as a product, we use:
<<extract-product>>=
cytb <- biofiles::filter(mito, key = "CDS", product = "^cytochrome b$")
cytb
@

\subsection{Accessing annotation data}

After having filtered the features of interest we might want to access the various annotation data in a form more emenable to further analysis. Accessor functions include \Rfunction{start}, \Rfunction{end}, \Rfunction{span}, \Rfunction{strand}, \Rfunction{key}, \Rfunction{locusTag}, \Rfunction{geneID}, \Rfunction{product}, \Rfunction{dbxref}, \Rfunction{proteinID}, or \Rfunction{translation}:
<<access-feature>>=
biofiles::start(cds[1:3])
biofiles::end(cds[1:3])
biofiles::span(cds[1:3])
biofiles::strand(cds[1:3])
biofiles::locusTag(cds[1:3])
biofiles::dbxref(cds[1:3])
biofiles::product(cds[1:3])
biofiles::translation(cds[1:3])
@


We can easily extract the DNA sequences for a bunch of CDSs by simply passing a feature table to the function \Rfunction{getSequence}:
<<access-sequence>>=
biofiles::getSequence(cds[1:6])
@


A general function for accessing GenBank qualifiers is \Rfunction{qualif}. If no further argument is provided \Rfunction{qualif} returns a \Rclass{data.frame} containing all qualifiers in the feature table. The argument \Rcode{which} can be used to targed specific qualifiers
<<access-qualif>>=
biofiles::qualif(cds[1:3])
biofiles::qualif(cds[1:3], which = c("gene", "locus_tag", "EC_number", "product", "db_xref.GeneID"))
@

To transform all annotation (``key'', ``location'', and ``qualifiers'') into a \Rclass{data.frame} we can use the function \Rfunction{select}:
<<access-select>>=
cols <- c("key", "gene", "locus_tag", "product")
biofiles::select(cds[1:4], .cols = cols)
@

Another useful format for annotation data are \Rclass{GRanges} provided by the \texttt{bioconductor} package \Biocpkg{GenomicRangs}. The function \Rfunction{ranges} allows to easily convert a \Rclass{gbFeatureTable} into a \Rclass{GRanges} object.
<<granges>>=
biofiles::ranges(cds)
@

A a default \Rfunction{ranges} will only include the feature key as metadata and will expand compound locations (such as typcal eukaryotic CDS). Compound locations can be merged using the argument \Rcode{join = TRUE}, and more qualifiers can be included using the argument \Rcode{include}:
<<granges-more>>=
biofiles::ranges(cds, join = TRUE, include = c("gene", "product", "db_xref"))
@


\section{Session Information}
All of the output in this vignette was produced under the following
conditions:
<<sessinfo>>=
utils::sessionInfo()
@

\end{document}
