Mercurial > repos > peterjc > mira4_assembler

changeset 21:4abe8d59a438 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded v0.0.4 preview 1; fix getting BAM without MAF
author peterjc
date Tue, 28 Oct 2014 08:29:59 -0400
parents aeb3e35f8236
children efeff76741b8
files test-data/header.mira tools/mira4/README.rst tools/mira4/mira4.py tools/mira4/mira4_convert.py tools/mira4/mira4_convert.xml tools/mira4/mira4_de_novo.xml tools/mira4/mira4_make_bam.py tools/mira4/mira4_mapping.xml tools/mira4/tool_dependencies.xml
diffstat 9 files changed, 397 insertions(+), 35 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/header.mira	Tue Oct 28 08:29:59 2014 -0400
@@ -0,0 +1,2 @@
+@Version	2	0
+@Program	MIRALIB
--- a/tools/mira4/README.rst	Tue Jun 10 10:11:58 2014 -0400
+++ b/tools/mira4/README.rst	Tue Oct 28 08:29:59 2014 -0400
@@ -44,9 +44,11 @@
 
 * ``mira4_de_novo.xml`` (the Galaxy tool definition for de novo usage)
 * ``mira4_mapping.xml`` (the Galaxy tool definition for mapping usage)
+* ``mira4_convert.xml`` (the Galaxy tool definition for converting MIRA files)
 * ``mira4_bait.xml`` (the Galaxy tool definition for mirabait)
 * ``mira4.py`` (the Python wrapper script)
-* ``mira_bait.py`` (the Python wrapper script for mirabait)
+* ``mira4_convert.py`` (the Python wrapper script for miraconvert)
+* ``mira4_bait.py`` (the Python wrapper script for mirabait)
 * ``mira4_validator.py`` (the XML parameter validation script)
 
 The suggested location is a new ``tools/mira4`` folder. You will also need to
@@ -75,6 +77,7 @@
     $ ./run_functional_tests.sh -id mira_4_0_bait
     $ ./run_functional_tests.sh -id mira_4_0_de_novo
     $ ./run_functional_tests.sh -id mira_4_0_mapping
+    $ ./run_functional_tests.sh -id mira_4_0_convert
 
 
 History
@@ -93,6 +96,7 @@
 v0.0.3  - Updated to target MIRA 4.0.2
 v0.0.4  - Using optparse for the Python wrapper script API
         - Made MAF and BAM outputs optional
+        - Include wrapper for ``miraconvert``
 ======= ======================================================================
 
 
@@ -105,7 +109,7 @@
 For making the "Galaxy Tool Shed" http://toolshed.g2.bx.psu.edu/ tarball use
 the following command from the Galaxy root folder::
 
-    $ tar -czf mira4_wrapper.tar.gz tools/mira4/README.rst tools/mira4/mira4_de_novo.xml tools/mira4/mira4_mapping.xml tools/mira4/mira4_bait.xml tools/mira4/mira4.py tools/mira4/mira4_make_bam.py tools/mira4/mira4_validator.py tools/mira4/mira4_bait.py tools/mira4/tool_dependencies.xml tools/mira4/repository_dependencies.xml test-data/U13small_m.fastq test-data/U13small_m.mira4_de_novo.fasta test-data/tvc_mini.fastq test-data/tvc_contigs.fasta test-data/tvc_map_ref_strain.fasta test-data/tvc_map_same_strain.fasta test-data/tvc_bait.fasta test-data/tvc_mini_bait_pos.fastq test-data/tvc_mini_bait_strict.fastq test-data/tvc_mini_bait_neg.fastq test-data/ecoli.fastq test-data/ecoli.mira4_de_novo.fasta test-data/empty_file.dat
+    $ tar -czf mira4_wrapper.tar.gz tools/mira4/README.rst tools/mira4/mira4_de_novo.xml tools/mira4/mira4_mapping.xml tools/mira4/mira4_bait.xml tools/mira4/mira4_convert.xml tools/mira4/mira4.py tools/mira4/mira4_make_bam.py tools/mira4/mira4_validator.py tools/mira4/mira4_convert.py tools/mira4/mira4_bait.py tools/mira4/tool_dependencies.xml tools/mira4/repository_dependencies.xml test-data/U13small_m.fastq test-data/U13small_m.mira4_de_novo.fasta test-data/tvc_mini.fastq test-data/tvc_contigs.fasta test-data/tvc_map_ref_strain.fasta test-data/tvc_map_same_strain.fasta test-data/tvc_bait.fasta test-data/tvc_mini_bait_pos.fastq test-data/tvc_mini_bait_strict.fastq test-data/tvc_mini_bait_neg.fastq test-data/ecoli.fastq test-data/ecoli.mira4_de_novo.fasta test-data/header.mira test-data/empty_file.dat
 
 Check this worked::
 
@@ -114,9 +118,11 @@
     tools/mira4/mira4_de_novo.xml
     tools/mira4/mira4_mapping.xml
     tools/mira4/mira4_bait.xml
+    tools/mira4/mira4_convert.xml
     tools/mira4/mira4.py
     tools/mira4/mira4_make_bam.py
     tools/mira4/mira4_validator.py
+    tools/mira4/mira4_convert.py
     tools/mira4/mira4_bait.py
     tools/mira4/tool_dependencies.xml
     tools/mira4/repository_dependencies.xml
@@ -132,6 +138,7 @@
     test-data/tvc_mini_bait_neg.fastq
     test-data/ecoli.fastq
     test-data/ecoli.mira4_de_novo.fasta
+    test-data/header.mira
     test-data/empty_file.dat
 
 
--- a/tools/mira4/mira4.py	Tue Jun 10 10:11:58 2014 -0400
+++ b/tools/mira4/mira4.py	Tue Oct 28 08:29:59 2014 -0400
@@ -58,6 +58,9 @@
 parser.add_option("--log", dest="log",
                   default="-", metavar="FILE",
                   help="MIRA logging output filename")
+parser.add_option("-v", "--version", dest="version",
+                  default=False, action="store_true",
+                  help="Show version and quit")
 options, args = parser.parse_args()
 manifest = options.manifest
 out_maf = options.maf
@@ -84,7 +87,7 @@
 mira_convert_ver = get_version(mira_convert)
 if not mira_convert_ver.strip().startswith("4.0"):
     stop_err("This wrapper is for MIRA V4.0, not:\n%s\n%s" % (mira_ver, mira_convert))
-if "-v" in sys.argv or "--version" in sys.argv:
+if options.version:
     print "%s, MIRA wrapper version %s" % (mira_ver, WRAPPER_VER)
     if mira_ver != mira_convert_ver:
         print "WARNING: miraconvert %s" % mira_convert_ver
@@ -192,7 +195,11 @@
     #For mapping mode, probably most people would expect a BAM file
     #using the reference FASTA file...
     if out_bam and out_bam != "-":
-        msg = make_bam(mira_convert, out_maf, ref_fasta, out_bam, handle)
+        if out_maf and out_maf != "-":
+            msg = make_bam(mira_convert, out_maf, ref_fasta, out_bam, handle)
+        else:
+            #Not collecting the MAF file, use original location        
+            msg = make_bam(mira_convert, old_maf, ref_fasta, out_bam, handle)
         if msg:
             stop_err(msg)
 
@@ -216,7 +223,9 @@
 
 assert os.path.isdir(temp)
 d = "%s_assembly" % name
-assert not os.path.isdir(d), "Path %s already exists" % d
+#This can fail on my development machine if stale folders exist
+#under Galaxy's .../database/job_working_directory/ tree:
+assert not os.path.isdir(d), "Path %r already exists:\n%s" % (d, os.path.abspath(d))
 try:
     #Check path access
     os.mkdir(d)
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/mira4/mira4_convert.py	Tue Oct 28 08:29:59 2014 -0400
@@ -0,0 +1,225 @@
+#!/usr/bin/env python
+"""A simple wrapper script to call MIRA and collect its output.
+
+This focuses on the miraconvert binary.
+"""
+import os
+import sys
+import subprocess
+import shutil
+import time
+import tempfile
+from optparse import OptionParser
+try:
+    from io import BytesIO
+except ImportError:
+    #Should we worry about Python 2.5 or older?
+    from StringIO import StringIO as BytesIO
+
+#Do we need any PYTHONPATH magic?
+from mira4_make_bam import depad
+
+WRAPPER_VER = "0.0.5" #Keep in sync with the XML file
+
+def stop_err(msg, err=1):
+    sys.stderr.write(msg+"\n")
+    sys.exit(err)
+
+def run(cmd):
+    #Avoid using shell=True when we call subprocess to ensure if the Python
+    #script is killed, so too is the child process.
+    try:
+        child = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
+    except Exception, err:
+        stop_err("Error invoking command:\n%s\n\n%s\n" % (" ".join(cmd), err))
+    #Use .communicate as can get deadlocks with .wait(),
+    stdout, stderr = child.communicate()
+    return_code = child.returncode
+    if return_code:
+        if stderr and stdout:
+            stop_err("Return code %i from command:\n%s\n\n%s\n\n%s" % (return_code, err, stdout, stderr))
+        else:
+            stop_err("Return code %i from command:\n%s\n%s" % (return_code, err, stderr))
+
+def get_version(mira_binary):
+    """Run MIRA to find its version number"""
+    # At the commend line I would use: mira -v | head -n 1
+    # however there is some pipe error when doing that here.
+    cmd = [mira_binary, "-v"]
+    try:
+        child = subprocess.Popen(cmd,
+                                 stdout=subprocess.PIPE,
+                                 stderr=subprocess.STDOUT)
+    except Exception, err:
+        sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (" ".join(cmd), err))
+        sys.exit(1)
+    ver, tmp = child.communicate()
+    del child
+    return ver.split("\n", 1)[0].strip()
+
+#Parse Command Line
+usage = """Galaxy MIRA4 wrapper script v%s - use as follows:
+
+$ python mira4_convert.py ...
+
+This will run the MIRA miraconvert binary and collect its output files as directed.
+""" % WRAPPER_VER
+parser = OptionParser(usage=usage)
+parser.add_option("--input", dest="input",
+                  default=None, metavar="FILE",
+                  help="MIRA input filename")
+parser.add_option("-x", "--min_length", dest="min_length",
+                  default="0",
+                  help="Minimum contig length")
+parser.add_option("-y", "--min_cover", dest="min_cover",
+                  default="0",
+                  help="Minimum average contig coverage")
+parser.add_option("-z", "--min_reads", dest="min_reads",
+                  default="0",
+                  help="Minimum reads per contig")
+parser.add_option("--maf", dest="maf",
+                  default="", metavar="FILE",
+                  help="MIRA MAF output filename")
+parser.add_option("--ace", dest="ace",
+                  default="", metavar="FILE",
+                  help="ACE output filename")
+parser.add_option("--bam", dest="bam",
+                  default="", metavar="FILE",
+                  help="Unpadded BAM output filename")
+parser.add_option("--fasta", dest="fasta",
+                  default="", metavar="FILE",
+                  help="Unpadded FASTA output filename")
+parser.add_option("--cstats", dest="cstats",
+                  default="", metavar="FILE",
+                  help="Contig statistics filename")
+parser.add_option("-v", "--version", dest="version",
+                  default=False, action="store_true",
+                  help="Show version and quit")
+options, args = parser.parse_args()
+if args:
+    stop_err("Expected options (e.g. --input example.maf), not arguments")
+
+input_maf = options.input
+out_maf = options.maf
+out_bam = options.bam
+out_fasta = options.fasta
+out_ace = options.ace
+out_cstats = options.cstats
+
+try:
+    mira_path = os.environ["MIRA4"]
+except KeyError:
+    stop_err("Environment variable $MIRA4 not set")
+mira_convert = os.path.join(mira_path, "miraconvert")
+if not os.path.isfile(mira_convert):
+    stop_err("Missing miraconvert under $MIRA4, %r\nFolder contained: %s"
+             % (mira_convert, ", ".join(os.listdir(mira_path))))
+
+mira_convert_ver = get_version(mira_convert)
+if not mira_convert_ver.strip().startswith("4.0"):
+    stop_err("This wrapper is for MIRA V4.0, not:\n%s\n%s" % (mira_ver, mira_convert))
+if options.version:
+    print "%s, MIRA wrapper version %s" % (mira_convert_ver, WRAPPER_VER)
+    sys.exit(0)
+
+if not input_maf:
+    stop_err("Input MIRA file is required")
+elif not os.path.isfile(input_maf):
+    stop_err("Missing input MIRA file: %r" % input_maf)
+
+if not (out_maf or out_bam or out_fasta or out_ace or out_cstats):
+    stop_err("No output requested")
+
+
+def check_min_int(value, name):
+    try:
+        i = int(value)
+    except:
+        stop_err("Bad %s setting, %r" % (name, value))
+    if i < 0:
+        stop_err("Negative %s setting, %r" % (name, value))
+    return i
+
+min_length = check_min_int(options.min_length, "minimum length")
+min_cover = check_min_int(options.min_cover, "minimum cover")
+min_reads = check_min_int(options.min_reads, "minimum reads")
+
+#TODO - Run MIRA in /tmp or a configurable directory?
+#Currently Galaxy puts us somewhere safe like:
+#/opt/galaxy-dist/database/job_working_directory/846/
+temp = "."
+
+
+cmd_list = [mira_convert]
+if min_length:
+    cmd_list.extend(["-x", str(min_length)])
+if min_cover:
+    cmd_list.extend(["-y", str(min_cover)])
+if min_reads:
+    cmd_list.extend(["-z", str(min_reads)])
+cmd_list.extend(["-f", "maf", input_maf, os.path.join(temp, "converted")])
+if out_maf:
+    cmd_list.append("maf")
+if out_bam:
+    cmd_list.append("samnbb")
+    if not out_fasta:
+        #Need this for samtools depad
+        out_fasta = os.path.join(temp, "depadded.fasta")
+if out_fasta:
+    cmd_list.append("fasta")
+if out_ace:
+    cmd_list.append("ace")
+if out_cstats:
+    cmd_list.append("cstats")
+run(cmd_list)
+
+def collect(old, new):
+    if not os.path.isfile(old):
+        stop_err("Missing expected output file %s" % old)
+    shutil.move(old, new)
+
+if out_maf:
+    collect(os.path.join(temp, "converted.maf"), out_maf)
+if out_fasta:
+    #Can we look at the MAF file to see if there are multiple strains?
+    old = os.path.join(temp, "converted_AllStrains.unpadded.fasta")
+    if os.path.isfile(old):
+        collect(old, out_fasta)
+    else:
+        #Might the output be filtered down to zero contigs?
+        old = os.path.join(temp, "converted.fasta")
+        if not os.path.isfile(old):
+            stop_err("Missing expected output FASTA file")
+        elif os.path.getsize(old) == 0:
+            print("Warning - no contigs (harsh filters?)")
+            collect(old, out_fasta)
+        else:
+            stop_err("Missing expected output FASTA file (only generic file present)")
+if out_ace:
+    collect(os.path.join(temp, "converted.maf"), out_ace)
+if out_cstats:
+    collect(os.path.join(temp, "converted_info_contigstats.txt"), out_cstats)
+
+if out_bam:
+    assert os.path.isfile(out_fasta)
+    old = os.path.join(temp, "converted.samnbb")
+    if not os.path.isfile(old):
+        old = os.path.join(temp, "converted.sam")
+    if not os.path.isfile(old):
+        stop_err("Missing expected intermediate file %s" % old)
+    h = BytesIO()
+    msg = depad(out_fasta, old, out_bam, h)
+    if msg:
+        print(msg)
+        print(h.getvalue())
+        h.close()
+        sys.exit(1)
+    h.close()
+    if out_fasta == os.path.join(temp, "depadded.fasta"):
+        #Not asked for by Galaxy, no longer needed
+        os.remove(out_fasta)
+
+if min_length or min_cover or min_reads:
+    print("Filtered.")
+else:
+    print("Converted.")
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/mira4/mira4_convert.xml	Tue Oct 28 08:29:59 2014 -0400
@@ -0,0 +1,114 @@
+<tool id="mira_4_0_convert" name="MIRA v4.0 miraconvert" version="0.0.5">
+    <description>Convert MIRA assembly to FASTA/SAM/BAM</description>
+    <requirements>
+        <requirement type="binary">miraconvert</requirement>
+        <requirement type="package" version="4.0">MIRA</requirement>
+    </requirements>
+    <version_command interpreter="python">mira4_convert.py --version</version_command>
+    <command interpreter="python">
+mira4_convert.py
+--input "$mira_file"
+--min_length $min_length
+--min_cover $min_cover
+--min_reads $min_reads
+#if str($maf_wanted)=="true":
+--maf "$out_maf"
+#end if
+#if str($fasta_wanted)=="true":
+--fasta "$out_fasta"
+#end if
+#if str($bam_wanted)=="true":
+--bam "$out_bam"
+#end if
+##Don't yet have a Galaxy datatype defined for ace:
+## #if str($ace_wanted)=="true":
+## --ace "$out_ace"
+## #end if
+#if str($cstats_wanted)=="true":
+--cstats "$out_cstats"
+#end if
+    </command>
+    <stdio>
+        <!-- Assume anything other than zero is an error -->
+        <exit_code range="1:" />
+        <exit_code range=":-1" />
+    </stdio>
+    <inputs>
+        <param name="mira_file" type="data" format="mira" required="true" label="MIRA Assembly Format input" />
+        <!-- TODO - top level select for contig versus read output? Or two Galaxy tools in different XML files? -->
+        <param name="min_length" type="integer" required="false" value="0" min="0"
+               label="Minimum contig length"
+               help="e.g. Set to 1000 to exclude small contigs. Default is to keep all contigs (minimum zero)" />
+        <param name="min_cover" type="integer" required="false" value="0" min="0"
+               label="Minimum average contig coverage"
+               help="e.g. Set to 10 to exclude low coverage contigs. Default is to keep all contigs (minimum zero)" />
+        <param name="min_reads" type="integer" required="false" value="0" min="0"
+               label="Minimum reads per contig"
+               help="e.g. Set to 5 to exclude low coverage contigs with only a few reads. Default is to keep all contigs (minimum zero)." />
+        <param name="maf_wanted" type="boolean" label="Output assembly in MIRA's own format? (useful if filtering)" checked="False" />
+        <param name="fasta_wanted" type="boolean" label="Convert assembly into (unpadded) FASTA?" checked="True" />
+        <param name="bam_wanted" type="boolean" label="Convert assembly into (upadded) BAM format?" checked="False" />
+        <!-- Don't yet have a Galaxy datatype defined for ace:
+        <param name="ace_wanted" type="boolean" label="Convert assembly in ACE format?" checked="False" />
+        -->
+        <param name="cstats_wanted" type="boolean" label="Assembly statistics file?" checked="False" />
+    </inputs>
+    <outputs>
+        <data name="out_maf" format="mira" label="$mira_file.name (filtered)">
+              <filter>maf_wanted is True</filter>
+        </data>
+        <data name="out_fasta" format="fasta" label="$mira_file.name (as FASTA)">
+              <filter>fasta_wanted is True</filter>
+        </data>
+        <data name="out_bam" format="bam" label="$mira_file.name (as BAM)">
+              <filter>bam_wanted is True</filter>
+        </data>
+        <!--
+        <data name="out_ace" format="ace" label="$mira_file.name (as ACE)">
+            <filter>ace_wanted is True</filter>
+        </data>
+        -->
+        <data name="out_cstats" format="tabular" label="$mira_file.name (filtered stats)">
+              <filter>cstats_wanted is True</filter>
+        </data>
+    </outputs>
+    <tests>
+        <!-- TODO -->
+    </tests>
+    <help>
+**What it does**
+
+Runs the ``miraconvert`` utility from MIRA v4.0 to filter and/or convert
+a MIRA Assembly Format file produced by a *mapping* or *de novo* assembly.
+
+**Example Usage**
+
+You want to remove all the low coverage contigs from a transcriptome
+assembly to focus on those with higher coverage.
+
+You want to convert your MIRA assembly into SAM/BAM to run a standard
+SNP finding tool.
+
+You've lost the FASTA consensus from your MIRA assembly and need to
+regenerate it.
+
+
+**Citation**
+
+If you use this Galaxy tool in work leading to a scientific publication please
+cite the following papers:
+
+Peter J.A. Cock, Björn A. Grüning, Konrad Paszkiewicz and Leighton Pritchard (2013).
+Galaxy tools and workflows for sequence analysis with applications
+in molecular plant pathology. PeerJ 1:e167
+http://dx.doi.org/10.7717/peerj.167
+
+Bastien Chevreux, Thomas Wetter and Sándor Suhai (1999).
+Genome Sequence Assembly Using Trace Signals and Additional Sequence Information.
+Computer Science and Biology: Proceedings of the German Conference on Bioinformatics (GCB) 99, pp. 45-56.
+http://www.bioinfo.de/isb/gcb99/talks/chevreux/main.html
+
+This wrapper is available to install into other Galaxy Instances via the Galaxy
+Tool Shed at http://toolshed.g2.bx.psu.edu/view/peterjc/mira4_assembler
+    </help>
+</tool>
--- a/tools/mira4/mira4_de_novo.xml	Tue Jun 10 10:11:58 2014 -0400
+++ b/tools/mira4/mira4_de_novo.xml	Tue Oct 28 08:29:59 2014 -0400
@@ -180,7 +180,8 @@
             <param name="bam_wanted" value="true"/>
             <output name="out_fasta" file="U13small_m.mira4_de_novo.fasta" ftype="fasta" />
             <output name="out_bam" file="empty_file.dat" compare="contains" />
-            <output name="out_maf" file="empty_file.dat" compare="contains" />
+            <!-- TODO: Suggest startswith as a compare method? -->
+            <output name="out_maf" file="header.mira" compare="contains" />
             <output name="out_log" file="empty_file.dat" compare="contains" />
         </test>
         <!-- Simple assembly based on MIRA's minidemo/solexa1 example
@@ -192,11 +193,9 @@
             <param name="job_quality" value="accurate" />
             <param name="type" value="none" />
             <param name="filenames" value="ecoli.fastq" ftype="fastqsanger" />
-            <param name="maf_wanted" value="true"/>
-            <param name="bam_wanted" value="true"/>
+            <param name="maf_wanted" value="false"/>
+            <param name="bam_wanted" value="false"/>
             <output name="out_fasta" file="ecoli.mira4_de_novo.fasta" ftype="fasta" />
-            <output name="out_bam" file="empty_file.dat" compare="contains" />
-            <output name="out_maf" file="empty_file.dat" compare="contains" />
             <output name="out_log" file="empty_file.dat" compare="contains" />
         </test>
     </tests>
--- a/tools/mira4/mira4_make_bam.py	Tue Jun 10 10:11:58 2014 -0400
+++ b/tools/mira4/mira4_make_bam.py	Tue Oct 28 08:29:59 2014 -0400
@@ -32,6 +32,29 @@
     log_handle.write(stdout)
     return child.returncode
 
+def depad(fasta_file, sam_file, bam_file, log_handle):
+    log_handle.write("\n================= Converting MIRA assembly from SAM to BAM ===================\n")
+    #Also doing SAM to (uncompressed) BAM during depad
+    bam_stem = bam_file + ".tmp" # Have write permissions and want final file in this folder
+    cmd = 'samtools depad -S -u -T "%s" "%s" | samtools sort - "%s"' % (fasta_file, sam_file, bam_stem)
+    return_code = run(cmd, log_handle)
+    if return_code:
+        return "Error %i from command:\n%s" % (return_code, cmd)
+    if not os.path.isfile(bam_stem + ".bam"):
+        return "samtools depad or sort failed to produce BAM file"
+
+    log_handle.write("\n====================== Indexing MIRA assembly BAM file =======================\n")
+    cmd = 'samtools index "%s.bam"' % bam_stem
+    return_code = run(cmd, log_handle)
+    if return_code:
+        return "Error %i from command:\n%s" % (return_code, cmd)
+    if not os.path.isfile(bam_stem + ".bam.bai"):
+        return "samtools indexing of BAM file failed to produce BAI file"
+
+    shutil.move(bam_stem + ".bam", bam_file)
+    os.remove(bam_stem + ".bam.bai") #Let Galaxy handle that...
+
+
 def make_bam(mira_convert, maf_file, fasta_file, bam_file, log_handle):
     if not os.path.isfile(mira_convert):
         return "Missing binary %r" % mira_convert
@@ -52,30 +75,14 @@
     if not os.path.isfile(sam_file):
         return "Conversion from MIRA to SAM failed"
 
-    log_handle.write("\n================= Converting MIRA assembly from SAM to BAM ===================\n")
     #Also doing SAM to (uncompressed) BAM during depad
-    bam_stem = bam_file + ".tmp" # Have write permissions and want final file in this folder
-    cmd = 'samtools depad -S -u -T "%s" "%s" | samtools sort - "%s"' % (fasta_file, sam_file, bam_stem)
-    return_code = run(cmd, log_handle)
-    if return_code:
-        return "Error %i from command:\n%s" % (return_code, cmd)
-    if not os.path.isfile(bam_stem + ".bam"):
-        return "samtools depad or sort failed to produce BAM file"
+    msg = depad(fasta_file, sam_file, bam_file, log_handle)
+    if msg:
+        return msg
 
     os.remove(sam_file)
     os.rmdir(tmp_dir)
 
-    log_handle.write("\n====================== Indexing MIRA assembly BAM file =======================\n")
-    cmd = 'samtools index "%s.bam"' % bam_stem
-    return_code = run(cmd, log_handle)
-    if return_code:
-        return "Error %i from command:\n%s" % (return_code, cmd)
-    if not os.path.isfile(bam_stem + ".bam.bai"):
-        return "samtools indexing of BAM file failed to produce BAI file"
-
-    shutil.move(bam_stem + ".bam", bam_file)
-    os.remove(bam_stem + ".bam.bai") #Let Galaxy handle that...
-
     return None #Good :)
 
 if __name__ == "__main__":
--- a/tools/mira4/mira4_mapping.xml	Tue Jun 10 10:11:58 2014 -0400
+++ b/tools/mira4/mira4_mapping.xml	Tue Oct 28 08:29:59 2014 -0400
@@ -186,7 +186,8 @@
             <param name="bam_wanted" value="true"/>
             <output name="out_fasta" file="tvc_map_ref_strain.fasta" ftype="fasta" />
             <output name="out_bam" file="empty_file.dat" compare="contains" />
-            <output name="out_maf" file="empty_file.dat" compare="contains" />
+            <!-- TODO: Suggest startswith as a compare method? -->
+            <output name="out_maf" file="header.mira" compare="contains" />
             <output name="out_log" file="empty_file.dat" compare="contains" />
         </test>
         <test>
@@ -196,11 +197,9 @@
             <param name="strain_setup" value="same" />
             <param name="type" value="none" />
             <param name="filenames" value="tvc_mini.fastq" ftype="fastqsanger" />
-            <param name="maf_wanted" value="true"/>
-            <param name="bam_wanted" value="true"/>
+            <param name="maf_wanted" value="false"/>
+            <param name="bam_wanted" value="false"/>
             <output name="out_fasta" file="tvc_map_same_strain.fasta" ftype="fasta" />
-            <output name="out_bam" file="empty_file.dat" compare="contains" />
-            <output name="out_maf" file="empty_file.dat" compare="contains" />
             <output name="out_log" file="empty_file.dat" compare="contains" />
         </test>
     </tests>
--- a/tools/mira4/tool_dependencies.xml	Tue Jun 10 10:11:58 2014 -0400
+++ b/tools/mira4/tool_dependencies.xml	Tue Oct 28 08:29:59 2014 -0400
@@ -1,7 +1,7 @@
 <?xml version="1.0"?>
 <tool_dependency>
     <package name="samtools" version="0.1.19">
-        <repository changeset_revision="40250a414486" name="package_samtools_0_1_19" owner="iuc" toolshed="https://testtoolshed.g2.bx.psu.edu" />
+        <repository changeset_revision="632f1a03db92" name="package_samtools_0_1_19" owner="iuc" toolshed="https://testtoolshed.g2.bx.psu.edu" />
     </package>
     <package name="MIRA" version="4.0">
         <install version="1.0">