shmlast

An improved implementation of Conditional Reciprocal Best Hits with LAST and Python

shmlast is a reimplementation of the Conditional Reciprocal Best Hits algorithm for finding potential orthologs between a transcriptome and a species-specific protein database. It uses the LAST aligner and the pydata stack to achieve much better performance while staying in the Python ecosystem.

About

Conditional Reciprocal Best Hits (CRBH) was originally described by Aubry et al. 2014 and implemented in the crb-blast package. CRBH builds on the traditional Reciprocal Best Hits (RBH) method for orthology assignment by training a simple model of the e-value cutoff for a particular length of sequence on an initial set of RBH's. From its github repository:

"Reciprocal best BLAST is a very conservative way to assign orthologs. The main innovation in
CRB-BLAST is to learn an appropriate e-value cutoff to apply to each pairwise alignment by taking
into account the overall relatedness of the two datasets being compared. This is done by fitting a
function to the distribution of alignment e-values over sequence lengths. The function provides the
e-value cutoff for a sequence of given length."

Unfortunately, the original implementation uses NCBI BLAST+ (which is incredibly slow), and is implemented in Ruby, which requires users to leave the Python-dominated bioinformatics software system. shmlast makes this algorithm available to users in Python-land, while also greatly improving performance by using LAST for initial homology searches. Additionally, shmlast outputs both the raw parameters and a plot of its model for inspection.

shmlast is designed for finding orthologs between transcriptomes and protein databases. As such, it currently does not support nucleotide-nucleotide or protein-protein alignments. This may be changed in a future version, but for now, it remains focused on that task.

Also note that RBH, and by extension CRBH, is meant for comparing between two species. Neither of these methods should be used for annotating a transcriptome with a mixed protein database (like, for example, uniref90).

Usage

For some transcriptome transcripts.fa and some protein database pep.faa, the basic usage is:

bash shmlast crbl -q transcripts.fa -d pep.faa

shmlast can be distributed across multiple cores using the --n_threads option.

bash shmlast crbl -q transcripts.fa -d pep.faa --n_threads 8

Another use case is to perform simple Reciprocal Best Hits; this can be done with the rbl subcommand. The maximum expectation-value can also be specified with -e.

bash shmlast rbl -q transcripts.fa -d pep.faa --e 0.000001

Output

shmlast outputs a plain CSV file with the CRBH's, which by default will be named $QUERY.x.$DATABASE.crbl.csv. This CSV file can be easily parsed with Pandas like so:

```Python import pandas as pd

crbldf = pd.readcsv('query.x.database.crbl.csv') ```

The columns are:

1. E: The e-value.
2. EG2: Expected alignments per square gigabase.
3. E_scaled: E-value rescaled for the model (see below for details).
4. ID: A unique ID for the alignment.
5. bitscore: The bitscore, calculated as (lambda * score - ln[K]) / ln[2].
6. qalnlen: Query alignment length.
7. q_frame: Frame in the query translation.
8. q_len: Length of the query sequence.
9. q_name: Name of the query sequence.
10. q_start: Start of query alignment. 11.q_strand: Strand of query alignment.
11. salnlen: Length of subject alignment.
12. s_len: Length of subject sequence.
13. s_name: Name of subject sequence.
14. s_start: Start of subject alignment.
15. s_strand: Strand of subject alignment.
16. score: The alignment score.

See http://last.cbrc.jp/doc/last-evalues.html for more information on e-values and scores.

Model Output

shmlast also outputs its model, both in CSV format and as a plot. The CSV file is named $QUERY.x.$DATABASE.crbl.model.csv, and has the following columns:

1. center: The center of the length bin.
2. size: The size of the bin.
3. left: The left of the bin.
4. right: The right of the bin.
5. fit: The scaled e-value cutoff for the bin.

To fit the model, the e-values are first scaled to a more suitable range using the equation Es = -log10(E), where Es is the scaled e-value. e-values of 0 are set to an arbitrarily small value to allow for log-scaling. The fit column of the model is this scaled value.

The model plot is named $QUERY.x.$DATABASE.crbl.model.plot.pdf by default.

Installation

via Conda

conda is the preferred installation method. shmlast is hosted on bioconda and it can be installed along with its dependencies using:

bash conda install shmlast -c bioconda

PyPI

If you really want to avoid conda, you can install via PyPI with:

bash pip install shmlast

After which you'll beed to install the third-party dependencies manually.

Third-party Dependencies

shmlast requires the LAST aligner and gnu-parallel. These will be installed automatically via conda if you choose that route; some other ways to install them follow.

Manually

LAST can be installed manually into your home directory like so:

bash cd curl -LO http://last.cbrc.jp/last-658.zip unzip last-658.zip pushd last-658 && make && make install prefix=~ && popd

And a recent version of gnu-parallel can be installed like so:

bash (wget -O - pi.dk/3 || curl pi.dk/3/ || fetch -o - http://pi.dk/3) | bash

Through a Package Manager

For Ubuntu 16.04 or newer, sufficiently new versions of both are available through the package manager:

bash sudo apt-get install last-align parallel

For OSX, you can get LAST through the homebrew-science channel:

bash brew tap homebrew/science brew install last

Library

shmlast is also a Python library. Each component of the pipeline is implemented as a pydoit task and can be used in doit workflows, and the implementations for calculating best hits, reciprocal best hits, and conditional reciprocal best hits are usable as Python classes. For example, the lastal task could be incorporated into a doit file like so:

```Python from shmlast.last import lastal_task

def tasklastal(): return lastaltask('query.fna', 'db.faa', translate=True) ```

Known Issues

There is currently an issue with IUPAC codes in RNA. This will be fixed soon.

Contributing

See CONTRIBUTING.md for guidelines.

References

1. Aubry S, Kelly S, Kümpers BMC, Smith-Unna RD, Hibberd JM (2014) Deep Evolutionary Comparison of Gene Expression Identifies Parallel Recruitment of Trans-Factors in Two Independent Origins of C4 Photosynthesis. PLoS Genet 10(6): e1004365. doi:10.1371/journal.pgen.1004365

2. O. Tange (2011): GNU Parallel - The Command-Line Power Tool, ;login: The USENIX Magazine, February 2011:42-47.

3. Kiełbasa, S. M., Wan, R., Sato, K., Horton, P., & Frith, M. C. (2011). Adaptive seeds tame genomic sequence comparison. Genome research, 21(3), 487-493.