comparison tools/protein_analysis/rxlr_motifs.py @ 20:a19b3ded8f33 draft

v0.2.11 Job splitting fast-fail; RXLR tools supports HMMER2 from BioConda; Capture more version information; misc internal changes

19:f3ecd80850e2

the predicted cleavage site, as this is expected to be more accurate.
Also note that the HMM score values have changed from v2.0 to v3.0.
Whisson et al. (2007) used SignalP v3.0 anyway.

Whisson et al. (2007) used HMMER 2.3.2, and althought their HMM model
can still be used with hmmsearch from HMMER 3 this this does give
slightly different results. We expect the hmmsearch from HMMER 2.3.2
(the last stable release of HMMER 2) to be present on the path under
the name hmmsearch2 (allowing it to co-exist with HMMER 3).
"""
import os
import sys
import re
import subprocess
from seq_analysis_utils import fasta_iterator

if "-v" in sys.argv:
    print("RXLR Motifs v0.0.10")
    sys.exit(0)

if len(sys.argv) != 5:
    sys.exit("Requires four arguments: protein FASTA filename, threads, model, and output filename")

    sys.exit("Did not recognise the model name %r\n"
             "Use Bhattacharjee2006, Win2007, or Whisson2007" % model)


def get_hmmer_version(exe, required=None):
    cmd = "%s -h" % exe
    try:
        child = subprocess.Popen([exe, "-h"], stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    except OSError:
        raise ValueError("Could not run %s" % exe)
    stdout, stderr = child.communicate()
    if required:
        return required in stdout

# Run hmmsearch for Whisson et al. (2007)
if model == "Whisson2007":
    hmm_file = os.path.join(os.path.split(sys.argv[0])[0],
                       "whisson_et_al_rxlr_eer_cropped.hmm")
    if not os.path.isfile(hmm_file):
        sys.exit("Missing HMM file for Whisson et al. (2007)")
    if not get_hmmer_version(hmmer_search, "HMMER 2.3.2 (Oct 2003)"):
        sys.exit("Missing HMMER 2.3.2 (Oct 2003) binary, %s" % hmmer_search)

# Cleanup
os.remove(signalp_input_file)
os.remove(signalp_output_file)

# Short summary to stdout for Galaxy's info display
print "%s for %i sequences:" % (model, total)
print ", ".join("%s = %i" % kv for kv in sorted(tally.iteritems()))

20:a19b3ded8f33

the predicted cleavage site, as this is expected to be more accurate.
Also note that the HMM score values have changed from v2.0 to v3.0.
Whisson et al. (2007) used SignalP v3.0 anyway.

Whisson et al. (2007) used HMMER 2.3.2, and althought their HMM model
can still be used with hmmsearch from HMMER 3, sadly this does give
slightly different results. We expect the hmmsearch from HMMER 2.3.2
(the last stable release of HMMER 2) to be present on the path under
the name hmmsearch2 (allowing it to co-exist with HMMER 3).

If using Conda, you should therefore install the special "hmmer2"
package from BioConda which provides "hmmsearch2" etc::

    conda install -c bioconda hmmer2

See https://bioconda.github.io/recipes/hmmer2/README.html and
https://anaconda.org/bioconda/hmmer2
"""

from __future__ import print_function

import os
import re
import subprocess
import sys

from seq_analysis_utils import fasta_iterator

if "-v" in sys.argv:
    print("RXLR Motifs v0.0.14")
    sys.exit(0)

if len(sys.argv) != 5:
    sys.exit("Requires four arguments: protein FASTA filename, threads, model, and output filename")

    sys.exit("Did not recognise the model name %r\n"
             "Use Bhattacharjee2006, Win2007, or Whisson2007" % model)


def get_hmmer_version(exe, required=None):
    try:
        child = subprocess.Popen([exe, "-h"],
                                 universal_newlines=True,
                                 stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    except OSError:
        raise ValueError("Could not run %s" % exe)
    stdout, stderr = child.communicate()
    if required:
        return required in stdout

# Run hmmsearch for Whisson et al. (2007)
if model == "Whisson2007":
    hmm_file = os.path.join(os.path.split(sys.argv[0])[0],
                            "whisson_et_al_rxlr_eer_cropped.hmm")
    if not os.path.isfile(hmm_file):
        sys.exit("Missing HMM file for Whisson et al. (2007)")
    if not get_hmmer_version(hmmer_search, "HMMER 2.3.2 (Oct 2003)"):
        sys.exit("Missing HMMER 2.3.2 (Oct 2003) binary, %s" % hmmer_search)

# Cleanup
os.remove(signalp_input_file)
os.remove(signalp_output_file)

# Short summary to stdout for Galaxy's info display
print("%s for %i sequences:" % (model, total))
print(", ".join("%s = %i" % kv for kv in sorted(tally.iteritems())))