Reciprocal Smallest Distance
Introduction
Citation: Wall, D.P., Fraser, H.B. and Hirsh, A.E. (2003) Detecting putative orthologs, Bioinformatics, 19, 1710-1711.
The reciprocal smallest distance (RSD) algorithm (Wall, et al., 2003. http://bioinformatics.oxfordjournals.org/content/19/13/1710) accurately infers orthologs between pairs of genomes by considering global sequence alignment and maximum likelihood evolutionary distance between sequences. Orthologs inferred with RSD for many species are available at Roundup (http://roundup.hms.harvard.edu/), which provides multi-species clusters of orthologous genes, output in formats for other phylogenetics packages, and sequence metadata such as Gene Ontology terms and database cross-references.
This package contains source code, scripts for running RSD, and example input and output files.
- README.md: the file you are reading now
- bin/rsd_search: a script that runs the reciprocal smallest distance (RSD) algorithm to search for orthologs.
- bin/rsd_blast: a script that computes and saves BLAST hits for use in multiple runs of RSD.
- bin/rsd_format: a script that turns FASTA-formatted genomes into BLAST-formatted indexes.
- rsd/: python package implementing the RSD algorithm.
- rsd/jones.dat, rsd/codeml.ctl: used by codeml/paml to compute the evolutionary distance between two sequences.
- examples/: a directory containing examples of inputs and outputs to rsd, including fasta-formatted genome protein sequence files, a query sequence id file (for --ids), and an orthologs output file.
Installing RSD
Prerequisites
RSD depends on Python, NCBI BLAST, PAML, and Kalign. It has been tested to work with the versions below. It might work with other versions too.
Install:
- Python 2.7: http://www.python.org/download/
- NCBI BLAST 2.2.24: ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/
- PAML 4.4: http://abacus.gene.ucl.ac.uk/software/paml.html
- Kalign 2.04: http://msa.sbc.su.se/cgi-bin/msa.cgi
Add the executables for python (version 2.7), makeblastdb, blastp, codeml, and kalign, to your PATH.
Install Using Pip
Use pip ( http://www.pip-installer.org/ ) to easily install the reciprocal_smallest_distance package:
pip install reciprocal_smallest_distance
Install From Github
Installing from github is an easy way to make the examples/ directory
accessible for running the examples in this README file.
Clone the latest version from GitHub:
cd ~ # Optional
git clone https://github.com/todddeluca/reciprocal_smallest_distance
Install reciprocal_smallest_distance, making sure to use Python 2.7:
cd reciprocal_smallest_distance
python setup.py install
Using RSD to Find Othologs
The following example commands demonstrate the main ways to run rsd_search.
Every invocation of rsd_search requires specifying the location of a
FASTA-formatted sequence file for two genomes, called the query and subject
genomes. Their order is arbitrary, but if you use the --ids option, the ids
must come from the query genome. You must also specify a file to write the
results of the orthologs found by the RSD algorithm. The
formats of the output file are described in detail below.
You can optionally specify a file containing ids using the --ids option.
Then rsd will only search for orthologs for those ids. Using --de gives you
the option of using different divergence and evalue thresholds from the
defaults.
Get help on how to run rsd_search, rsd_blast, or rsd_format:
rsd_search -h
rsd_blast -h
rsd_format -h
Find orthologs between all the sequences in the query and subject genomes, using the default divergence and evalue thresholds:
rsd_search -q examples/genomes/Mycoplasma_genitalium.aa/Mycoplasma_genitalium.aa \
--subject-genome=examples/genomes/Mycobacterium_leprae.aa/Mycobacterium_leprae.aa \
-o Mycoplasma_genitalium.aa_Mycobacterium_leprae.aa_0.8_1e-5.orthologs.txt
Find orthologs using several divergence and evalue threshold combinations:
rsd_search -q examples/genomes/Mycoplasma_genitalium.aa/Mycoplasma_genitalium.aa \
--subject-genome=examples/genomes/Mycobacterium_leprae.aa/Mycobacterium_leprae.aa \
-o Mycoplasma_genitalium.aa_Mycobacterium_leprae.aa.several.orthologs.txt \
--de 0.2 1e-20 --de .5 0.00001 --de 0.8 0.1
It is not necessary to format a FASTA file for BLAST or compute BLAST hits
because rsd_search does it for you. However if you plan on running
rsd_search multiple times for the same genomes, especially for large genomes,
you can save time by using rsd_format to preformatting the FASTA files and
rsd_blast to precomputing the BLAST hits. When running rsd_blast, make
sure to use an --evalue as large as the largest evalue threshold you intend to
give to rsd_search.
Here is how to format a pair of FASTA files in place:
rsd_format -g examples/genomes/Mycoplasma_genitalium.aa/Mycoplasma_genitalium.aa
rsd_format -g examples/genomes/Mycobacterium_leprae.aa/Mycobacterium_leprae.aa
And here is how to format the FASTA files, putting the results in another directory (the current directory in this case):
rsd_format -g examples/genomes/Mycoplasma_genitalium.aa/Mycoplasma_genitalium.aa -d .
rsd_format -g examples/genomes/Mycobacterium_leprae.aa/Mycobacterium_leprae.aa -d .
Here is how to compute forward and reverse blast hits (using the default evalue):
rsd_blast -v -q examples/genomes/Mycoplasma_genitalium.aa/Mycoplasma_genitalium.aa \
--subject-genome=examples/genomes/Mycobacterium_leprae.aa/Mycobacterium_leprae.aa \
--forward-hits q_s.hits --reverse-hits s_q.hits
Here is how to compute forward and reverse blast hits for rsd_search, using
genomes that have already been formatted for blast and a non-default evalue:
rsd_blast -v -q Mycoplasma_genitalium.aa \
--subject-genome=Mycobacterium_leprae.aa \
--forward-hits q_s.hits --reverse-hits s_q.hits \
--no-format --evalue 0.1
Find orthologs between all the sequences in the query and subject genomes using genomes that have already been formatted for blast:
rsd_search -q Mycoplasma_genitalium.aa \
--subject-genome=Mycobacterium_leprae.aa \
-o Mycoplasma_genitalium.aa_Mycobacterium_leprae.aa_0.8_1e-5.orthologs.txt \
--no-format
Find orthologs between all the sequences in the query and subject genomes using hits that have already been computed. Notice that --no-format is included, because since the blast hits have already been computed the genomes do not need to be formatted for blast:
rsd_search -v --query-genome Mycoplasma_genitalium.aa \
--subject-genome=Mycobacterium_leprae.aa \
-o Mycoplasma_genitalium.aa_Mycobacterium_leprae.aa.default.orthologs.txt \
--forward-hits q_s.hits --reverse-hits s_q.hits --no-format
Find orthologs for specific sequences in the query genome. For finding
orthologs for only a few sequences, using --no-blast-cache can speed up
computation. YMMV.
rsd_search -q examples/genomes/Mycoplasma_genitalium.aa/Mycoplasma_genitalium.aa \
--subject-genome=examples/genomes/Mycobacterium_leprae.aa/Mycobacterium_leprae.aa \
-o examples/Mycoplasma_genitalium.aa_Mycobacterium_leprae.aa_0.8_1e-5.orthologs.txt \
--ids examples/Mycoplasma_genitalium.aa.ids.txt --no-blast-cache
Orthologs can be saved in several different formats using the --outfmt option
of rsd_search. The default format, --outfmt -1, refers to --outfmt 3.
Inspired by Uniprot dat files, a set of orthologs starts with a parameters
line, then has 0 or more ortholog lines, then has an end line. The parametes
are the query genome name, subject genome name, divergence threshold, and
evalue threshold. Each ortholog is on a single line listing the query sequence
id, the subject sequence id, and the maximum likelihood distance estimate.
This format can represent orthologs for multiple sets of parameters in a single
file as well as sets of parameters with no orthologs. Therefore it is suitable
for use with rsd_search when specifying multiple divergence and evalue
thresholds.
Here is an example containing 2 parameter combinations, one of which has no orthologs:
PA\tLACJO\tYEAS7\t0.2\t1e-15
OR\tQ74IU0\tA6ZM40\t1.7016
OR\tQ74K17\tA6ZKK5\t0.8215
//
PA\tMYCGE\tMYCHP\t0.2\t1e-30
//
The original format of RSD, --outfmt 1, is provided for backward
compatibility. Each line contains an ortholog, represented as subject sequence
id, query sequence id, and maximum likelihood distance estimate. It can only
represent a single set of orthologs in a file.
Example:
A6ZM40\tQ74IU0\t1.7016
A6ZKK5\tQ74K17\t0.8215
Also provided for backward compatibility is a format used internally by Roundup (http://roundup.hms.harvard.edu/) which is like the original RSD format, except the query sequence id column is before the subject sequence id.
Example:
Q74IU0\tA6ZM40\t1.7016
Q74K17\tA6ZKK5\t0.8215
Testing
This assumes you have cloned the repository.
For convenience, define some environment variables:
export TMP=$HOME/tmp
mkdir -p $TMP
export REPO_PATH=<path/to/reciprocal_smallest_distance> # e.g. export REPO_PATH=~/git/reciprocal_smallest_distance
Create a virtualenv in which to install RSD:
cd $TMP
virtualenv venv
Install RSD into the virtualenv:
venv/bin/pip install -e $REPO_PATH
Run an example command:
venv/bin/rsd_search -q $REPO_PATH/examples/genomes/Mycoplasma_genitalium.aa/Mycoplasma_genitalium.aa \
--subject-genome=$REPO_PATH/examples/genomes/Mycobacterium_leprae.aa/Mycobacterium_leprae.aa \
-o Mycoplasma_genitalium.aa_Mycobacterium_leprae.aa_0.8_1e-5.orthologs.txt
Examine the results:
less Mycoplasma_genitalium.aa_Mycobacterium_leprae.aa_0.8_1e-5.orthologs.txt
