<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/phylotar-retrieve-orthologous-sequences-from-genbank/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://ropensci.org/technotes/2018/08/08/phylotar/"> 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>
<p>In this technote I will outline what phylotaR was developed for, how to install it and how to run it with some simple examples.</p>
<h2>What is phylotaR?</h2>
<p>In any phylogenetic analysis it is important to identify sequences that share the same <a rel="nofollow" href="https://en.wikipedia.org/wiki/Sequence_homology#Orthology" target="_blank">orthology</a> – homologous sequences separated by speciation events. This is often performed by simply searching an online sequence repository using sequence labels. Relying solely on sequence labels, however, can miss sequences that have either not been labelled, have unanticipated names or have been mislabelled.</p>
<p>The phylotaR package, like its earlier inspiration <a rel="nofollow" href="http://phylota.net/" target="_blank">PhyLoTa</a>, identifies sequences of shared orthology not by naming matching but instead through the use of an alignment search tool, <a rel="nofollow" href="https://en.wikipedia.org/wiki/BLAST" target="_blank">BLAST</a>. As a result, a user of phylotaR is able to download more relevant sequence data for their phylogenetic analysis than they otherwise would.</p>
<p>From the taxon of interest the phylotaR pipeline traverses down the taxonomy and downloads entire subsets of sequences representing different nodes in the taxonomy. For taxa with too many sequences (e.g. <em>Homo sapiens</em>, <em>Arabidopsis thaliana</em>), a random subset is instead downloaded. For each of these subsets of downloaded sequences all-vs-all BLAST is then performed. The BLAST results of these subsets are parsed to identify all the independent, homologous clusters of sequences as determined using E-values and coverages. A secondary clustering step performs BLAST again to identify higher-level taxonomic clusters, above the level at which the all-vs-all BLAST searches were performed. This pipeline is automated and all a user needs to supply is a taxonomic group of interest from which the pipeline will begin its traverse. The pipeline is broken down into four stages: retrieve taxonomic information (‘taxise’), download sequences, identify clusters and identify clusters among the clusters.</p>
<p><img alt="phylotaR pipeline stages" src="https://i1.wp.com/ropensci.org/img/blog-images/2018-08-08-phylotar/stages.png?w=456&ssl=1" data-recalc-dims="1"><br>
<em>Figure 1. The stages of phylotaR pipeline: taxise, download, cluster and cluster^2. Note, ‘taxise’ is the name of a stage and does not relate to the package <code>taxize</code>.</em></p>
<p>After a pipeline has completed, the package provides a series of tools for interacting with the results to identify potential ortholgous clusters suitable for phylogenetic analysis. These tools include methods for filtering based on sequence lengths, sequence annotations, taxonomy …etc. As well as for visualisations to, for example, check the relative abundance of sequences by taxon.</p>
<p>For more information on how the pipeline works, please see the open-access scientific article: <a rel="nofollow" href="https://doi.org/10.3390/life8020020" target="_blank">phylotaR: An Automated Pipeline for Retrieving Orthologous DNA Sequences from GenBank in R</a></p>
<h2>Installation</h2>
<p>phylotaR is available from CRAN.</p>
<pre><code>install.packages('phylotaR')
</code></pre>
<p>Alternatively, the latest development version can be downloaded from phylotaR’s <a rel="nofollow" href="https://github.com/ropensci/phylotaR" target="_blank">GitHub page</a>.</p>
<pre><code>devtools::install_github(repo = 'ropensci/phylotaR', build_vignettes = TRUE)
</code></pre>
<h3>BLAST+</h3>
<p>In addition to the R package you will also need to have installed BLAST+ – a local copy of the well-known <a rel="nofollow" href="https://blast.ncbi.nlm.nih.gov/Blast.cgi" target="_blank">BLAST software</a> on your machine. Unfortunately this can be a little complex as its installation depends on your operating system (Windows, Mac or Linux). Fortunately, NCBI provides detailed installation instructions for each flavour of operating system: <a rel="nofollow" href="https://www.ncbi.nlm.nih.gov/books/NBK279671/" target="_blank">NCBI installation instructions</a>.</p>
<p>Once BLAST+ is installed you will need to record the location of the BLAST+ file system path where the exectuable programs are located. This should be something like <code>C:\users\tao\desktop\blast-2.2.29+\bin\</code> on a Windows machine or <code>/usr/local/ncbi/blast/bin/</code> on a Unix.</p>
<h2>Quick guide</h2>
<h3>Taxonomic ID</h3>
<p>After installation, a phylotaR pipeline can be launched with the <code>setup()</code> and <code>run()</code> functions. Before launching a run, a user must first look up the NCBI taxonomic ID for their chosen taxonomic group. To find the ID of your taxon, search <a rel="nofollow" href="https://www.ncbi.nlm.nih.gov/taxonomy" target="_blank">NCBI taxonomy</a> using the scientific name or an English common name. For example, upon searching for the primate genus ‘Aotus’ we should come across this <a rel="nofollow" href="https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=9504" target="_blank">page</a> which states the taxonomic ID as <code>9504</code>.</p>
<p>Because names are not unique phylotaR does not allow names to be used instead of taxonomic IDs. Aotus is a good example as there is also a plant genus of the same name. If we were to search NCBI with just the name then we would be downloading a mixture of primate and plant DNA. Such is the importance of IDs. (Alternatively, you could automate this step yourself using the <a rel="nofollow" href="https://github.com/ropensci/taxize" target="_blank">taxize</a> package.)</p>
<h3>Running</h3>
<p>With the package installed, BLAST+ available and ID to hand, we can run our pipeline.</p>
<pre><code># load phylotaR
library(phylotaR)
# create a folder to hold the cached data and results
dir.create('aotus')
# setup the folder and check that BLAST is available
setup(wd = wd, txid = 9504, ncbi_dr = '/usr/local/ncbi/blast/bin/')
# run the pipeline
run(wd = wd)
# ^^ running should take about 5 minutes to complete
</code></pre>
<p>For more details on running the pipeline, such as changing the <a rel="nofollow" href="https://ropensci.github.io/phylotaR/reference/parameters.html" target="_blank">parameters</a> or understanding the results, see the vignette, <code>vignette('phylotaR')</code> or visit the <a rel="nofollow" href="https://ropensci.github.io/phylotaR/" target="_blank">website</a>.</p>
<h3>Timings</h3>
<p>The time it takes for a pipeline to complete depends on the number of taxonomic groups a taxon contains and the number of sequences it represents. Below are some examples of the time taken for various taxonomic groups.</p>
<table>
<thead>
<tr>
<th align="left">Taxon</th>
<th align="right">N. taxa</th>
<th align="right">N. sequences</th>
<th align="right">Time taken (mins.)</th>
</tr>
</thead>
<tbody>
<tr>
<td align="left">Anisoptera (“dragonflies”)</td>
<td align="right">1175</td>
<td align="right">11432</td>
<td align="right">72</td>
</tr>
<tr>
<td align="left">Acipenseridae (“sturgeons”)</td>
<td align="right">51</td>
<td align="right">2407</td>
<td align="right">13</td>
</tr>
<tr>
<td align="left">Tinamiformes (“tinamous”)</td>
<td align="right">25</td>
<td align="right">251</td>
<td align="right">2.7</td>
</tr>
<tr>
<td align="left">Aotus (“night/owl monkeys”)</td>
<td align="right">13</td>
<td align="right">1499</td>
<td align="right">3.9</td>
</tr>
<tr>
<td align="left">Bromeliaceae (“bromeliads”)</td>
<td align="right">1171</td>
<td align="right">9833</td>
<td align="right">66</td>
</tr>
<tr>
<td align="left">Cycadidae (“cycads”)</td>
<td align="right">353</td>
<td align="right">8331</td>
<td align="right">37</td>
</tr>
<tr>
<td align="left">Eutardigrada (“water bears”)</td>
<td align="right">261</td>
<td align="right">960</td>
<td align="right">14</td>
</tr>
<tr>
<td align="left">Kazachstania (a group of yeast)</td>
<td align="right">40</td>
<td align="right">623</td>
<td align="right">23</td>
</tr>
<tr>
<td align="left">Platyrrhini (“new world monkeys”)</td>
<td align="right">212</td>
<td align="right">12731</td>
<td align="right">60</td>
</tr>
</tbody>
</table>
<h3>Reviewing the results</h3>
<p>After a pipeline has completed, the <code>wd</code> will contain all the results files. The information contained within these files can be read into the R environment using the <code>read_phylota()</code> function. It will generate a <code>Phylota</code> object that contains information on all the identified ‘clusters’ – i.e. groups of ortholgous sequences. A user can interact with the object to filter out unwanted sequences, clusters or taxonomic groups. The package comes with completed results for playing with. In the example below, we demonstrate how to generate a presence/absence matrix of all the genera in the top 10 clusters identifed for cycads, <a rel="nofollow" href="https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=1445963" target="_blank">1445963</a>.</p>
<pre><code>library(phylotaR)
data(cycads)
# drop all but first ten
cycads &lt;- drop_clstrs(cycads, cycads@cids[1:10])
# get genus-level taxonomic names
genus_txids &lt;- get_txids(cycads, txids = cycads@txids, rnk = 'genus')
genus_txids &lt;- unique(genus_txids)
# dropping missing
genus_txids &lt;- genus_txids[genus_txids !=  '']
genus_nms &lt;- get_tx_slot(cycads, genus_txids, slt_nm = 'scnm')
# make alphabetical for plotting
genus_nms &lt;- sort(genus_nms, decreasing = TRUE)
# generate geom_object
p &lt;- plot_phylota_pa(phylota = cycads, cids = cycads@cids, txids = genus_txids,
                     txnms = genus_nms)
# plot
print(p)
</code></pre>
<p><img alt="presenceabsence of cycad genera" src="https://i2.wp.com/ropensci.org/img/blog-images/2018-08-08-phylotar/pa.png?w=456&ssl=1" data-recalc-dims="1"><br>
<em>Figure 2. The presence (dark block) or absence (light block) for different cycad genera across the top ten clusters in the example dataset.</em></p>
<p>In this next example, we create a treemap showing the differences in the number of sequences and clusters identifed between genera in tinamous, <a rel="nofollow" href="https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=8802" target="_blank">8802</a>. (For the unintiated, <a rel="nofollow" href="https://en.wikipedia.org/wiki/Tinamou" target="_blank">tinamous</a> are semi-flightless birds found in South America and members of the ratities – the same group comprising ostrichs and kiwis.)</p>
<pre><code>library(phylotaR)
data(tinamous)
# Plot taxa, size by n. sq, fill by ncl
txids &lt;- get_txids(tinamous, txids = tinamous@txids, rnk = 'genus')
txids &lt;- txids[txids !=  '']
txids &lt;- unique(txids)
txnms &lt;- get_tx_slot(tinamous, txids, slt_nm = 'scnm')
p &lt;- plot_phylota_treemap(phylota = tinamous, txids = txids, txnms = txnms,
                          area = 'nsq', fill = 'ncl')
print(p)
</code></pre>
<p><img alt="treemap of tinmous genera" src="https://i0.wp.com/ropensci.org/img/blog-images/2018-08-08-phylotar/treemap.png?w=456&ssl=1" data-recalc-dims="1"><br>
<em>Figure 3. Relative number of sequences and clusters per tinamous genus. The larger the size of the box, the more sequences are represented for the genus. The lighter the blue colour, the more clusters are represented for the genus.</em></p>
<p>Through interacting with the Phylota object and using the various functions for manipulating it, a user can extract the specific ortholgous sequences of interest and write out the sequences in <a rel="nofollow" href="https://en.wikipedia.org/wiki/FASTA" target="_blank">fasta format</a> with the <a rel="nofollow" href="https://ropensci.github.io/phylotaR/reference/write_sqs.html" target="_blank"><code>write_sqs()</code></a> function.</p>
<pre><code># get sequences for a cluster and write out
data('aotus')
# take the first cluser ID
cid &lt;- aotus@cids[[1]]
# extract its sequence IDs from Phylota object
sids &lt;- aotus[[cid]]@sids
# design sequence names for definition line
sid_txids &lt;- get_sq_slot(phylota = aotus, sid = sids, slt_nm = 'txid')
sid_spnms &lt;- get_tx_slot(phylota = aotus, txid = sid_txids, slt_nm = 'scnm')
sq_nms &lt;- paste0(sid_spnms, ' | ', sids)
# write out
write_sqs(phylota = aotus, outfile = file.path(tempdir(), 'test.fasta'),
          sq_nm = sq_nms, sid = sids)
</code></pre>
<p>For more information and examples on manipulating the Phylota object, see the <a rel="nofollow" href="https://ropensci.github.io/phylotaR/" target="_blank">phylotaR website</a>.</p>
<h2>Future</h2>
<p>We have many ideas for improving phylotaR, such as making use of the BLAST API – instead of relying on users installing BLAST+ on their machine – and allowing users to introduce their own non-GenBank sequences. Please see the <a rel="nofollow" href="https://github.com/ropensci/phylotaR/blob/master/CONTRIBUTING.md" target="_blank">contributing page</a> for a complete and current list of development options. Fork requests are welcome!</p>
<p>Alternatively, if you have any ideas of your own for new features than please open a <a rel="nofollow" href="https://github.com/ropensci/phylotaR/issues" target="_blank">new issue</a>.</p>
<h2>Acknowledgements</h2>
<p>Big thanks to <a rel="nofollow" href="https://github.com/arendsee" target="_blank">Zebulun Arendsee</a>, <a rel="nofollow" href="https://github.com/naupaka" target="_blank">Naupaka Zimmerman</a> and <a rel="nofollow" href="https://github.com/sckott" target="_blank">Scott Chamberlain</a> for reviewing/editing the package for ROpenSci.</p>
<h2>Useful Links</h2>
<ul>
<li><a rel="nofollow" href="https://ropensci.github.io/phylotaR/" target="_blank">The phylotaR wesbite</a></li>
<li><a rel="nofollow" href="http://phylota.net/" target="_blank">The original PhyLoTa website</a></li>
<li><a rel="nofollow" href="https://www.ncbi.nlm.nih.gov/books/NBK279671/" target="_blank">BLAST+ installation and usage instructions</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/phylotaR" target="_blank">phylotaR on GitHub</a></li>
</ul>
<h2>References</h2>
<p>Bennett, D., Hettling, H., Silvestro, D., Zizka, A., Bacon, C., Faurby, S., … Antonelli, A. (2018). phylotaR: An Automated Pipeline for Retrieving Orthologous DNA Sequences from GenBank in R. <em>Life</em>, <strong>8</strong>(2), 20. <a rel="nofollow" href="https://doi.org/10.3390/life8020020" target="_blank">DOI:10.3390/life8020020</a></p>
<p>Sanderson, M. J., Boss, D., Chen, D., Cranston, K. A., &amp; Wehe, A. (2008). The PhyLoTA Browser: Processing GenBank for molecular phylogenetics research. <em>Systematic Biology</em>, <strong>57</strong>(3), 335–346. <a rel="nofollow" href="https://doi.org/10.1080/10635150802158688" target="_blank">DOI:10.<sup>1080</sup>⁄<sub>10635150802158688</sub></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/08/08/phylotar/"> 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/sEMcEKZSNfg" height="1"></a></div>