<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/restez-query-genbank-locally/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://ropensci.org/technotes/2018/12/03/restez/"> rOpenSci - open tools for open science</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>What is <code>restez</code>?</h2>
<p>R packages for interacting with the National Center for Biotechnology Information (NCBI) have, to-date, depended on <a rel="nofollow" href="https://en.wikipedia.org/wiki/Application_programming_interface" target="_blank">API</a> query calls via <a rel="nofollow" href="https://www.ncbi.nlm.nih.gov/search/" target="_blank">NCBI’s Entrez</a>.<br>
For computational analyses that require the automated look-up of reams of biological sequence data, piecemeal querying via bandwith-limited requests is evidently not ideal. These queries are not only slow, but they depend on network connections and the remote server’s consistent behaviour. Additionally, users who make very large requests over extended periods of time run the risk of being blocked.</p>
<p><code>restez</code> attempts to make large queries to <a rel="nofollow" href="https://www.ncbi.nlm.nih.gov/genbank/" target="_blank">NCBI GenBank</a> more efficient by allowing users to download whole sections of GenBank, create a local database from these downloaded files and then query this mini-GenBank version instead.<br>
This process is far more efficient as the downloaded files are compressed and users can limit the size of the database by only creating it<br>
from sequences of interest (limiting by taxonomic domain and/or sequence size).</p>
<p><code>restez</code> tries to be user-friendly: a database can be set up in just a few function calls (set path, download and create),<br>
a database can be queried with a consistent set of functions (the <code>gb_*_get()</code> functions),<br>
the number of arguments per function is limited, and the package is designed to integrate with pre-exisiting R packages<br>
that interact with NCBI (<a rel="nofollow" href="https://github.com/ropensci/rentrez" target="_blank"><code>rentrez</code></a> and<br>
<a rel="nofollow" href="https://github.com/ropensci/phylotaR" target="_blank"><code>phylotaR</code></a>)</p>
<p>For more a detailed description and for tutorials of the package, please visit the<br>
<a rel="nofollow" href="https://ropensci.github.io/restez/" target="_blank"><code>restez</code> website</a>.</p>
<p><img alt="restez_outline" src="https://i2.wp.com/ropensci.org/img/blog-images/2018-12-03-restez/outline.png?w=456&ssl=1" data-recalc-dims="1"><br>
<em>Figure 1. Diagrammatic outline of the <code>restez</code> functions and folder structure. Data is downloaded from NCBI into a file path<br>
set by the user. Raw downloads are stored in “downloads/” the generated database is stored in “sql_db/“. This database can<br>
then be queried with a series of <code>gb_*_get()</code> functions as well as some additional wrappers.</em></p>
<h2>Installation</h2>
<p><code>restez</code> (v1.0.0) is available from CRAN.</p>
<pre><code>install.packages('restez')
</code></pre>
<p>Alternatively, the latest development version can be downloaded from <code>restez</code>’s <a rel="nofollow" href="https://github.com/ropensci/restez" target="_blank">GitHub page</a>.</p>
<pre><code>devtools::install_github(repo = 'ropensci/restez')
</code></pre>
<h2>Usage</h2>
<p>To walk us through the basics of <code>restez</code> let’s pretend we’re a microbiologist interested in the sequence diversity among all the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Bacteriophage" target="_blank">bacteriophages</a> that infect <a rel="nofollow" href="https://en.wikipedia.org/wiki/Escherichia" target="_blank"><em>Escherichia</em> bacteria</a>. First, we will first need to create our sequence database. Second, we will need to identify the <em>Escherichia</em> phage sequences in this database. And then finally we will need to write out the sequences in a suitable format for a sequence diversity analysis. Thankfully, <code>restez</code> can perform all these steps.</p>
<h3>Set-up</h3>
<p>To get started with <code>restez</code> we first have to download and create a database. This set-up consists of three steps:</p>
<ul>
<li>set a filepath where we would like to store the database, using<br>
<a rel="nofollow" href="https://ropensci.github.io/restez/reference/restez_path_set.html" target="_blank"><code>restez_path_set()</code></a></li>
<li>download sections of GenBank,<br>
using <a rel="nofollow" href="https://ropensci.github.io/restez/reference/db_download.html" target="_blank"><code>db_download()</code></a></li>
<li>create the database from the downloaded files,<br>
using <a rel="nofollow" href="https://ropensci.github.io/restez/reference/db_create.html" target="_blank"><code>db_create()</code></a></li>
</ul>
<p>Depending on how many GenBank files you select to download, the above process can take up to several hours.<br>
In this example, however, we will only download and set up a database for ‘phage’ sequences<br>
which should take between 5-10 mintues depending on your machine and internet connection.</p>
<h4>Path</h4>
<pre><code>library(restez)
# create a new folder in your working directory to host a database
restez_path &lt;- file.path(getwd(), 'phages')
dir.create(restez_path)
# set the restez path
restez_path_set(restez_path)
</code></pre>
<p>The above code will set the <code>restez</code> path. All downloaded files and the created database will be stored in this path.<br>
(Keep a note of the <code>restez</code> path, you will need it later.)</p>
<h4>Download</h4>
<p>To download sequences, run the interactive function <code>db_download()</code>.</p>
<pre><code>db_download()
</code></pre>
<p>This will produce a list of options, like this:</p>
<pre><code>Which sequence file types would you like to download?
Choose from those listed below:
● 1  - 'Bacterial'
        542 files and 123 GB
● 2  - 'EST (expressed sequence tag)'
        486 files and 243 GB
● 3  - 'Constructed'
        372 files and 84.6 GB
● 4  - 'Patent'
        340 files and 77.9 GB
● 5  - 'GSS (genome survey sequence)'
        308 files and 117 GB
● 6  - 'TSA (transcriptome shotgun assembly)'
        234 files and 53.9 GB
● 7  - 'Plant sequence entries (including fungi and algae)'
        232 files and 62.8 GB
● 8  - 'HTGS (high throughput genomic sequencing)'
        155 files and 36.7 GB
● 9  - 'Invertebrate'
        112 files and 24.8 GB
● 10 - 'Environmental sampling'
        103 files and 24.1 GB
● 11 - 'Other vertebrate'
        94 files and 19.4 GB
● 12 - 'Primate'
        59 files and 13.7 GB
● 13 - 'Viral'
        58 files and 13.1 GB
● 14 - 'Other mammalian'
        55 files and 9.44 GB
● 15 - 'Rodent'
        31 files and 7.32 GB
● 16 - 'STS (sequence tagged site)'
        20 files and 4.45 GB
● 17 - 'HTC (high throughput cDNA sequencing)'
        15 files and 3.43 GB
● 18 - 'Synthetic and chimeric'
        10 files and 2.42 GB
● 19 - 'Phage'
        5 files and 1.12 GB
● 20 - 'Unannotated'
        1 files and 0.00111 GB

Provide one or more numbers separated by spaces.
e.g. to download all Mammalian sequences, type: &quot;12 14 15&quot; followed by Enter

Which files would you like to download?
(Press Esc to quit) 
</code></pre>
<p>We can download all phage sequences by typing <code>19</code>. After pressing Enter, we will be told of the likely total file size<br>
of the download. If you have enough free space, push any key to continue. This will then initiate a download process<br>
for all phage sequences files on GenBank.</p>
<h4>Create</h4>
<p>After the download process has completed, we can create the database with <code>db_create()</code>.</p>
<pre><code>db_create()
</code></pre>
<p>This will add all of the downloaded files to the database. It will take a while to complete. When it finishes, we can then identify our <em>Escherichia</em> phage sequences!</p>
<h3>Querying</h3>
<h4>Status</h4>
<p>After we have built the database, we can query it! For every new R session, we will always need to point <code>restez</code> to the<br>
database using <code>restez_path_set()</code> and then connect to the database with <code>restez_connect()</code>. To get started, let’s see<br>
the database status, is it ready for querying?</p>
<pre><code>library(restez)
restez_path &lt;- file.path(getwd(), 'phages')
restez_path_set(restez_path)
restez_connect()
restez_status()
</code></pre>
<pre><code>Checking setup status at  ...
────────────────────────────────────────────────────────────────────────────────
Restez path ...
... Path '/phages/restez'
... Does path exist? 'Yes'
────────────────────────────────────────────────────────────────────────────────
Download ...
... Path '/phages/restez/downloads'
... Does path exist? 'Yes'
... N. files 6
... N. GBs 0.34
... GenBank division selections 'Phage'
... GenBank Release 228
... Last updated '2018-11-16 14:46:56'
────────────────────────────────────────────────────────────────────────────────
Database ...
... Path '/phages/restez/sql_db'
... Does path exist? 'Yes'
... N. GBs 1.12
... Is database connected? 'Yes'
... Does the database have data? 'Yes'
... Number of sequences 14911
... Min. sequence length 0
... Max. sequence length Inf
... Last_updated '2018-11-16 14:54:32'
</code></pre>
<p>The above status report tells us the database, exists, has data and is connected – which means it’s ready for queries.<br>
(To get a simple TRUE or FALSE for whether the database is ready, use<br>
<a rel="nofollow" href="https://ropensci.github.io/restez/reference/restez_ready.html" target="_blank"><code>restez_ready()</code></a>.)</p>
<h4>Get-tools</h4>
<p><code>restez</code> comes with a series of <a rel="nofollow" href="https://ropensci.github.io/restez/reference/index.html#section-genbank" target="_blank"><code>gb_*_get()</code></a> functions for parsing the GenBank records to pull out specific elements (sequences, definition lines, whole records). We can find records in the database using <a rel="nofollow" href="https://www.ncbi.nlm.nih.gov/Sequin/acc.html" target="_blank">Accession IDs</a>. To list all Accession IDs in a database, we can use <a rel="nofollow" href="https://ropensci.github.io/restez/reference/list_db_ids.html" target="_blank"><code>list_db_ids()</code></a>.</p>
<pre><code># get a random accession ID from the database
id &lt;- sample(list_db_ids(), 1)
# definitions
def &lt;- gb_definition_get(id)[[1]]
print(def)
#> [1] &quot;Unidentified clone B15 DNA sequence from ocean beach sand&quot;
</code></pre>
<h3><em>Escherichia</em> phage sequences</h3>
<p>In our scenario, we’re interested in finding and writing out all the <em>Escherichia</em> phage sequences. We can do this by looking up the organism names of the sequence sources of all the sequences in the database. We can then parse these names for <code>&quot;escherichia&quot;</code> and write out the resulting list of sequences.</p>
<pre><code># get list of ALL IDs in database
ids &lt;- list_db_ids(n = NULL)
# get all organisms
organisms &lt;- gb_organism_get(ids)
# parse for Escherichia
is_escherichia &lt;- grepl('escherichia', organisms, ignore.case = TRUE)
# fetch fasta formatted sequences
fastas &lt;- gb_fasta_get(ids[is_escherichia])
# check ...
cat(fastas[1])
#> >AB000833.1 Bacteriophage Mu DNA for ORF1, sheath protein gpL, ORF2, ORF3, complete cds
#> ACGGTCAGACGTTTGGCCCGACCACCGGGATGAGGCTGACGCAGGTCAGAAATCTTTGTGACGACAACCG
#> TATCAATGCCGGTGTGGTGCTTTACGGCGTTCTGTTCAGTGGCACAACCCCGCTGCCGTCCGTAGTGGAC
#> CTGGATTCGCTGGATGATTACGAGCGTCACTGGCAGACCTGGAAATTCCCGGACGAAACCCCGGAATTTG
#> CCGCACATATCAATGTGAATCAGGAAAAGGATCATGATGCTGAAAATTAAACCCGCAGCGGGAAAAGCCA
#> ....
# write out
write(fastas, file = 'phage_seqs.fasta')
</code></pre>
<p>In this little example, we could identify our sequences of interest using <code>restez</code> itself. Ordinarily, however, because sequences can only be looked up via Accession IDs, users will probably not use <code>restez</code> for sequence discovery, only retrieval. For a more adaptable example of searching and fetching sequences, see <a rel="nofollow" href="https://ropensci.github.io/restez/articles/2_search_and_fetch.html" target="_blank">How to search and fetch sequences</a></p>
<h2>Integrations</h2>
<p>To minimise the coding effort on the part of a user, <code>restez</code> has been built to work with R packages that already connect to<br>
<a rel="nofollow" href="https://www.ncbi.nlm.nih.gov/books/NBK49540/" target="_blank">NCBI’s Entrez</a>. After setting up a <code>restez</code> database the same functions of these<br>
other packages can be used to query NCBI Entrez. Internally, <code>restez</code> will query its local database and if it cannot find all<br>
of the requested sequences it will pass these arguments on to these other packages.</p>
<p>For example, users can use the <code>entrez_fetch()</code> function of the <a rel="nofollow" href="https://github.com/ropensci/rentrez" target="_blank"><code>rentrez</code></a> package. Running<br>
this function through <code>restez</code> means a user can first check the local database rather than make lots of queries over the internet.<br>
The function arguments are exactly the same.<br>
Additionally, users can set up up a <code>restez</code> database before launching a <a rel="nofollow" href="https://github.com/ropensci/phylotaR" target="_blank"><code>phylotaR</code></a> run.<br>
<code>phylotaR</code> searches NCBI for orthologous sequence clusters for a given taxonomic ID. If a <code>restez</code> database is set-up, <code>phylotaR</code><br>
will first search the local database before downloading via Entrez.</p>
<p>For more information on these integrations see the additional documentation:<br>
<a rel="nofollow" href="https://ropensci.github.io/restez/articles/restez.html#entrez" target="_blank"><code>rentrez</code> and <code>restez</code></a> and<br>
<a rel="nofollow" href="https://ropensci.github.io/restez/articles/4_phylotar.html" target="_blank"><code>phylotaR</code> and <code>restez</code></a></p>
<h2>Future</h2>
<p>We have many ideas for improving <code>restez</code> and we welcome forks and pull requests! Our current list of ideas for improvement, include:</p>
<ul>
<li><strong>Protein database</strong> – the current code could be easily duplicated for working with protein databases, not just GenBank.</li>
<li><strong>Taxonomy</strong> – integration of existing taxonomic packages with <code>restez</code>.</li>
<li><strong>Retmodes</strong> – <code>restez</code> only supports text-based return modes, it could be expanded to include XML.</li>
</ul>
<p>Please see the <a rel="nofollow" href="https://github.com/ropensci/restez/blob/master/CONTRIBUTING.md" target="_blank">contributing page</a> for more details and any updates.</p>
<p>If you have any ideas of your own for new features than please open a <a rel="nofollow" href="https://github.com/ropensci/restez/issues" target="_blank">new issue</a>.</p>
<h2>Acknowledgements</h2>
<p>Big thanks to <a rel="nofollow" href="https://github.com/eveskew" target="_blank">Evan Eskew</a> and <a rel="nofollow" href="https://github.com/naupaka" target="_blank">Naupaka Zimmerman</a> for reviewing the package;<br>
and to <a rel="nofollow" href="https://github.com/cboettig" target="_blank">Carl Boettiger</a> and <a rel="nofollow" href="https://github.com/noamross" target="_blank">Noam Ross</a> for useful comments during the review;<br>
and, of course, to <a rel="nofollow" href="https://github.com/sckott" target="_blank">Scott Chamberlain</a> for editing!</p>
<h2>Useful Links</h2>
<ul>
<li><a rel="nofollow" href="https://ropensci.github.io/restez/" target="_blank">The <code>restez</code> website</a></li>
<li><a rel="nofollow" href="https://github.com/ropensci/rentrez#rentrez" target="_blank"><code>rentrez</code></a></li>
<li><a rel="nofollow" href="https://www.ncbi.nlm.nih.gov/genbank/" target="_blank">GenBank</a></li>
<li><a rel="nofollow" href="https://github.com/ropensci/restez" target="_blank"><code>restez</code> on GitHub</a></li>
<li><a rel="nofollow" href="https://github.com/ropensci/phylotaR" target="_blank"><code>phylotaR</code></a></li>
</ul>
<h2>Reference</h2>
<p>Bennett, D.J., Hettling, H., Silvestro, D., Vos, R. and Antonelli, A. <strong>2018</strong>. 2018. restez: Create and Query a Local Copy of GenBank in R. <em>Journal of Open Source Software</em>, 3(31), 1102, <a rel="nofollow" href="https://doi.org/10.21105/joss.01102" target="_blank">https://doi.org/10.21105/joss.01102</a></p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://ropensci.org/technotes/2018/12/03/restez/"> rOpenSci - open tools for open science</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/C9UHplkRbz0" height="1"></a></div>