BLAST in Galaxy
Part of my work for the week was to start using Galaxy more extensively at SANBI. I.e. to make it more usable. Last week I wrote an authentication plugin to allow a Galaxy server to authenticate using PAM. This got accepted into the 15.07 release of Galaxy, so I updated our Galaxy server to that release. I had neglected to include the example
auth_conf.xml in the code I committed, but working off the example on my laptop I got PAM authentication working as a replacement for the previous HTTP authentication (which also spoke to PAM on the backend). I also took the opportunity to switch our server to using HTTPS using the SANBI wildcard certificate.
My first attempt at a practical use for our server came when I needed to run the BLAST step of the OrthoMCL pipeline. OrthoMCL uses an all-against-all BLAST as its input dataset, and based on the data I had from our colleagues, I had a collection of about 300,000 proteins to BLAST against each other. I started this off as an array job at CHPC but thought I could try and work locally as well, as a proof of concept. (Actually there was a previous step, filtering out poor proteins, but I’ll get to that below.) My first attempt at using BLAST hit a bug: “NotFound: cannot find ‘files_path’ while searching for ‘db_opts.histdb.files_path'”. This exception was thrown from
__build_command_line in Galaxy’s
lib/galaxy/tools/evaluation.py because the BLAST wrappers use an attribute called
files_path instead of
extra_files_path. Peter Cock and John Chilton discuss the problem in this Github issue and Peter quickly committed a workaround to the BLAST tools.
Having fixed that, and having prepared the protein set (outside Galaxy), I decided to take a chance on the Galaxy “parallelisation” code. This is enabled through appropriate tags in the tool XML, and in the case of the
blastp wrapper splits the query dataset into chunks of 1000 sequences each before submitting jobs (in Galaxy terms, actually tasks, not fully fledged jobs) to the cluster. Unfortunately these are individual jobs, not an array job, because array jobs are only implemented in the still-only-on-the-horizon DRMAA version 2. In any event, our cluster can handle thousands of job submissions so I hit go, saw the history item turn from grey to yellow, and waited. Unfortunately, after a day or so it went red (failed), but by then I was too busy with other stuff to debug it. To be continued…
(As an aside, the BLAST wrapper wraps BLAST+, whereas OrthoMCL uses legacy BLAST. I still need to check that the BLAST wrapper exposes enough flags in order to guarantee equivalence. A useful guide for some of the corresponding flags can be found on this page about ortholog finding).
OrthoMCL in and out of Galaxy
As mentioned previously, I was running BLAST as part of the OrthoMCL pipeline. OrthoMCL uses BLAST, MCL and a database (in the version we use, SQLite3) to compute the orthologs in a set of proteins. The pipeline has two steps before the BLAST stage (orthomclAdjustFasta and orthomclFilterFasta), five between the BLAST and MCL stages and a final step to process the MCL output. Currently I use a Makefile to execute the pipeline but at the GCC2015 Hackathon AJ started work on some wrappers for the steps in the pipeline. There has been previous work on executing OrthoMCL within Galaxy but that ran the entire workflow as a single tool. We want to implement the pipeline as a Galaxy workflow because that way we can (in theory at least) benefit from improvements in how BLAST is executed (e.g. parallelism) or even replace the BLAST step with a similar (but apparently faster) tool such as Diamond. The OrthoMCL pipeline is pretty linear so even given the limited capabilities of workflows in current Galaxy (as discussed by John Chilton at BOSC 2015) creating a OrthoMCL workflow should be pretty easy.
To that end we’ve now got a Github repository for the tool wrappers. I’m trying to follow the structure that groups like IUC use. AJ’s working on
orthomclAdjustFasta, so I decided to tackle
orthomclFilterFasta, a tool that takes a directory full of FASTA files as input, does some simple filtering and outputs a combined FASTA file. I’m not 100% sure on the requirements for the command line (I need to go back into the code and see how it is executed) so I’ve got a tool that generates a single shell command in the form:
mkdir inputs && /bin/bash orthomcl_prepare_dataset_for_filter.sh dataset1.dat && /bin/bash orthomcl prepare_dataset_for_filter.sh dataset2.dat && orthomclFilterFasta inputs/ <p1> <p2>
prepare_dataset_for_filter.sh is just a simple script to take a FASTA file, extract the tag that OrthoMCL uses to identify sets (added by
orthomclAdjustFasta) and renames the file according to that tag. The
orthomclFilterFasta tool insists that input files end in
.fasta and are named according to their tag.
In any event the tool runs fine on a local Galaxy install. The next step is to get tool dependencies right, which is where the stuff in the
package directory comes in. Galaxy can install packages for you (in an admin-configurable folder). For a tool the dependencies it needs are specified in a file called
tool_dependencies.xml that is in the same folder as the tool XML.
The tool dependencies specify packages to install. For OrthoMCL two new packages have been written (see here, one for OrthoMCL and one for the Perl DBD::SQLite module that it depends on. OrthoMCL in turn depends on Perl and DBD::SQLite. This is done using a
repository_dependencies.xml file – I’m still not sure if this is the correct approach, but in any event it follows the guide for simple repository dependencies in Galaxy. One limitation to repository dependencies is that they apparently only work within a single toolshed, so what to do if the package you require is in another toolshed?
Thus far the tool dependency stuff has been tested on a local Galaxy installation. It doesn’t work. Directories are created, but they are empty. Further testing and a better testing procedure is needed. Eric Rasche mentioned that Marius van den Beek has some Jenkins based testing framework that uses Docker to create sandboxes, and it is here – so perhaps getting this up and running is a next step.
And then finally, FastOrtho seems like a possibly viable alternative to OrthoMCL. The output seems roughly similar to OrthoMCLs and it is much faster (and as a single tool with no Perl dependencies, easier to package), but as with all new tools in bioinformatics, we’ll have to prove that it works well enough to replace OrthoMCL (which is somewhat of a standard in this domain). Well check back in a few weeks for updates…