Mercurial > repos > peterjc > mira4_assembler

annotate tools/mira4/mira4.py @ 19:8487d70e82aa draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded v0.0.3 preview 1, target MIRA v4.0.2
author peterjc
date Wed, 21 May 2014 06:56:06 -0400
parents 381aa262c8cb
children aeb3e35f8236
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
1 #!/usr/bin/env python
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
2 """A simple wrapper script to call MIRA and collect its output.
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
3 """
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
4 import os
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
5 import sys
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
6 import subprocess
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
7 import shutil
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
8 import time
18
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
9 import tempfile
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
10
9
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
11 #Do we need any PYTHONPATH magic?
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
12 from mira4_make_bam import make_bam
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
13
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
14 WRAPPER_VER = "0.0.1" #Keep in sync with the XML file
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
15
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
16 def stop_err(msg, err=1):
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
17 sys.stderr.write(msg+"\n")
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
18 sys.exit(err)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
19
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
20
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
21 def get_version(mira_binary):
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
22 """Run MIRA to find its version number"""
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
23 # At the commend line I would use: mira -v | head -n 1
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
24 # however there is some pipe error when doing that here.
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
25 cmd = [mira_binary, "-v"]
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
26 try:
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
27 child = subprocess.Popen(cmd,
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
28 stdout=subprocess.PIPE,
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
29 stderr=subprocess.STDOUT)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
30 except Exception, err:
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
31 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (" ".join(cmd), err))
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
32 sys.exit(1)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
33 ver, tmp = child.communicate()
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
34 del child
9
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
35 return ver.split("\n", 1)[0].strip()
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
36
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
37
7
902f01c1084b Uploaded v0.0.1 preview 7, with mirabait wrapper
peterjc
parents: 5
diff changeset
38 try:
902f01c1084b Uploaded v0.0.1 preview 7, with mirabait wrapper
peterjc
parents: 5
diff changeset
39 mira_path = os.environ["MIRA4"]
12
a6a56440567c Uploaded v0.0.2 preview 4, now using MIRA 4.0 (not RC), and sets MIRA4 environment variable during install.
peterjc
parents: 11
diff changeset
40 except KeyError:
7
902f01c1084b Uploaded v0.0.1 preview 7, with mirabait wrapper
peterjc
parents: 5
diff changeset
41 stop_err("Environment variable $MIRA4 not set")
902f01c1084b Uploaded v0.0.1 preview 7, with mirabait wrapper
peterjc
parents: 5
diff changeset
42 mira_binary = os.path.join(mira_path, "mira")
902f01c1084b Uploaded v0.0.1 preview 7, with mirabait wrapper
peterjc
parents: 5
diff changeset
43 if not os.path.isfile(mira_binary):
13
7fcabeeca5df Uploaded v0.0.2 preview 5, fixes for MIRA 4.0 (final), more verbose error if $MIRA4 path wrong
peterjc
parents: 12
diff changeset
44 stop_err("Missing mira under $MIRA4, %r\nFolder contained: %s"
7fcabeeca5df Uploaded v0.0.2 preview 5, fixes for MIRA 4.0 (final), more verbose error if $MIRA4 path wrong
peterjc
parents: 12
diff changeset
45 % (mira_binary, ", ".join(os.listdir(mira_path))))
9
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
46 mira_convert = os.path.join(mira_path, "miraconvert")
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
47 if not os.path.isfile(mira_convert):
13
7fcabeeca5df Uploaded v0.0.2 preview 5, fixes for MIRA 4.0 (final), more verbose error if $MIRA4 path wrong
peterjc
parents: 12
diff changeset
48 stop_err("Missing miraconvert under $MIRA4, %r\nFolder contained: %s"
7fcabeeca5df Uploaded v0.0.2 preview 5, fixes for MIRA 4.0 (final), more verbose error if $MIRA4 path wrong
peterjc
parents: 12
diff changeset
49 % (mira_convert, ", ".join(os.listdir(mira_path))))
7
902f01c1084b Uploaded v0.0.1 preview 7, with mirabait wrapper
peterjc
parents: 5
diff changeset
50
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
51 mira_ver = get_version(mira_binary)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
52 if not mira_ver.strip().startswith("4.0"):
9
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
53 stop_err("This wrapper is for MIRA V4.0, not:\n%s\n%s" % (mira_ver, mira_binary))
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
54 mira_convert_ver = get_version(mira_convert)
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
55 if not mira_convert_ver.strip().startswith("4.0"):
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
56 stop_err("This wrapper is for MIRA V4.0, not:\n%s\n%s" % (mira_ver, mira_convert))
5
ffefb87bd414 Uploaded v0.0.1 preview 5, using MIRA 4.0 RC4, supports segment_placement (pairing type)
peterjc
parents: 4
diff changeset
57 if "-v" in sys.argv or "--version" in sys.argv:
ffefb87bd414 Uploaded v0.0.1 preview 5, using MIRA 4.0 RC4, supports segment_placement (pairing type)
peterjc
parents: 4
diff changeset
58 print "%s, MIRA wrapper version %s" % (mira_ver, WRAPPER_VER)
9
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
59 if mira_ver != mira_convert_ver:
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
60 print "WARNING: miraconvert %s" % mira_convert_ver
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
61 sys.exit(0)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
62
7
902f01c1084b Uploaded v0.0.1 preview 7, with mirabait wrapper
peterjc
parents: 5
diff changeset
63
11
02350bef2e99 Uploaded v0.0.2 preview 3, correction regression in URL
peterjc
parents: 9
diff changeset
64 try:
02350bef2e99 Uploaded v0.0.2 preview 3, correction regression in URL
peterjc
parents: 9
diff changeset
65 threads = int(os.environ.get("GALAXY_SLOTS", "1"))
02350bef2e99 Uploaded v0.0.2 preview 3, correction regression in URL
peterjc
parents: 9
diff changeset
66 except ValueError:
02350bef2e99 Uploaded v0.0.2 preview 3, correction regression in URL
peterjc
parents: 9
diff changeset
67 threads = 1
02350bef2e99 Uploaded v0.0.2 preview 3, correction regression in URL
peterjc
parents: 9
diff changeset
68 assert 1 <= threads, threads
7
902f01c1084b Uploaded v0.0.1 preview 7, with mirabait wrapper
peterjc
parents: 5
diff changeset
69
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
70
18
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
71 def override_temp(manifest):
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
72 """Override ``-DI:trt=/tmp`` in manifest with environment variable.
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
73
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
74 Currently MIRA 4 does not allow envronment variables like ``$TMP``
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
75 inside the manifest, which is a problem if you need to override
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
76 the default at run time.
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
77
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
78 The tool XML will ``/tmp`` and we replace that here with
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
79 ``tempfile.gettempdir()`` which will respect $TMPDIR, $TEMP, $TMP
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
80 as explained in the Python standard library documentation:
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
81 http://docs.python.org/2/library/tempfile.html#tempfile.tempdir
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
82
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
83 By default MIRA 4 would write its temporary files within the output
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
84 folder, which is a problem if that is a network drive.
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
85 """
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
86 handle = open(manifest, "r")
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
87 text = handle.read()
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
88 handle.close()
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
89
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
90 #At time of writing, this is at the end of a file,
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
91 #but could be followed by a space in future...
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
92 text = text.replace("-DI:trt=/tmp", "-DI:trt=" + tempfile.gettempdir())
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
93
19
8487d70e82aa Uploaded v0.0.3 preview 1, target MIRA v4.0.2
peterjc
parents: 18
diff changeset
94 #Want to try to ensure this gets written to disk before MIRA attempts
8487d70e82aa Uploaded v0.0.3 preview 1, target MIRA v4.0.2
peterjc
parents: 18
diff changeset
95 #to open it - any networked file system may impose a delay...
18
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
96 handle = open(manifest, "w")
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
97 handle.write(text)
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
98 handle.flush()
19
8487d70e82aa Uploaded v0.0.3 preview 1, target MIRA v4.0.2
peterjc
parents: 18
diff changeset
99 os.fsync(handle.fileno())
18
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
100 handle.close()
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
101
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
102
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
103 def log_manifest(manifest):
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
104 """Write the manifest file to stderr."""
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
105 sys.stderr.write("\n%s\nManifest file\n%s\n" % ("="*60, "="*60))
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
106 with open(manifest) as h:
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
107 for line in h:
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
108 sys.stderr.write(line)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
109 sys.stderr.write("\n%s\nEnd of manifest\n%s\n" % ("="*60, "="*60))
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
110
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
111
4
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
112 def collect_output(temp, name, handle):
9
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
113 """Moves files to the output filenames (global variables)."""
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
114 n3 = (temp, name, name, name)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
115 f = "%s/%s_assembly/%s_d_results" % (temp, name, name)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
116 if not os.path.isdir(f):
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
117 log_manifest(manifest)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
118 stop_err("Missing output folder")
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
119 if not os.listdir(f):
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
120 log_manifest(manifest)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
121 stop_err("Empty output folder")
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
122 missing = []
4
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
123
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
124 old_maf = "%s/%s_out.maf" % (f, name)
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
125 if not os.path.isfile(old_maf):
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
126 #Triggered extractLargeContigs.sh?
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
127 old_maf = "%s/%s_LargeContigs_out.maf" % (f, name)
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
128
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
129 #De novo or single strain mapping,
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
130 old_fasta = "%s/%s_out.unpadded.fasta" % (f, name)
9
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
131 ref_fasta = "%s/%s_out.padded.fasta" % (f, name)
4
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
132 if not os.path.isfile(old_fasta):
9
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
133 #Mapping (StrainX versus reference) or de novo
4
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
134 old_fasta = "%s/%s_out_StrainX.unpadded.fasta" % (f, name)
9
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
135 ref_fasta = "%s/%s_out_StrainX.padded.fasta" % (f, name)
4
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
136 if not os.path.isfile(old_fasta):
9
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
137 old_fasta = "%s/%s_out_ReferenceStrain.unpadded.fasta" % (f, name)
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
138 ref_fasta = "%s/%s_out_ReferenceStrain.padded.fasta" % (f, name)
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
139
4
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
140
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
141 missing = False
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
142 for old, new in [(old_maf, out_maf),
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
143 (old_fasta, out_fasta)]:
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
144 if not os.path.isfile(old):
4
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
145 missing = True
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
146 else:
4
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
147 handle.write("Capturing %s\n" % old)
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
148 shutil.move(old, new)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
149 if missing:
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
150 log_manifest(manifest)
4
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
151 sys.stderr.write("Contents of %r:\n" % f)
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
152 for filename in sorted(os.listdir(f)):
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
153 sys.stderr.write("%s\n" % filename)
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
154
9
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
155 #For mapping mode, probably most people would expect a BAM file
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
156 #using the reference FASTA file...
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
157 msg = make_bam(mira_convert, out_maf, ref_fasta, out_bam, handle)
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
158 if msg:
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
159 stop_err(msg)
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
160
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
161 def clean_up(temp, name):
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
162 folder = "%s/%s_assembly" % (temp, name)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
163 if os.path.isdir(folder):
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
164 shutil.rmtree(folder)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
165
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
166 #TODO - Run MIRA in /tmp or a configurable directory?
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
167 #Currently Galaxy puts us somewhere safe like:
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
168 #/opt/galaxy-dist/database/job_working_directory/846/
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
169 temp = "."
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
170 #name, out_fasta, out_qual, out_ace, out_caf, out_wig, out_log = sys.argv[1:8]
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
171 name = "MIRA"
9
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
172 manifest, out_maf, out_bam, out_fasta, out_log = sys.argv[1:]
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
173
18
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
174 override_temp(manifest)
381aa262c8cb Uploaded v0.0.2 preview 10, override /tmp via environment variable
peterjc
parents: 13
diff changeset
175
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
176 start_time = time.time()
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
177 #cmd_list =sys.argv[8:]
11
02350bef2e99 Uploaded v0.0.2 preview 3, correction regression in URL
peterjc
parents: 9
diff changeset
178 cmd_list = [mira_binary, "-t", str(threads), manifest]
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
179 cmd = " ".join(cmd_list)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
180
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
181 assert os.path.isdir(temp)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
182 d = "%s_assembly" % name
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
183 assert not os.path.isdir(d), "Path %s already exists" % d
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
184 try:
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
185 #Check path access
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
186 os.mkdir(d)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
187 except Exception, err:
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
188 log_manifest(manifest)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
189 sys.stderr.write("Error making directory %s\n%s" % (d, err))
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
190 sys.exit(1)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
191
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
192 #print os.path.abspath(".")
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
193 #print cmd
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
194
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
195 handle = open(out_log, "w")
4
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
196 handle.write("======================== MIRA manifest (instructions) ========================\n")
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
197 m = open(manifest, "rU")
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
198 for line in m:
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
199 handle.write(line)
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
200 m.close()
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
201 del m
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
202 handle.write("\n")
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
203 handle.write("============================ Starting MIRA now ===============================\n")
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
204 handle.flush()
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
205 try:
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
206 #Run MIRA
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
207 child = subprocess.Popen(cmd_list,
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
208 stdout=handle,
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
209 stderr=subprocess.STDOUT)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
210 except Exception, err:
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
211 log_manifest(manifest)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
212 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err))
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
213 #TODO - call clean up?
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
214 handle.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err))
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
215 handle.close()
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
216 sys.exit(1)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
217 #Use .communicate as can get deadlocks with .wait(),
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
218 stdout, stderr = child.communicate()
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
219 assert not stdout and not stderr #Should be empty as sent to handle
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
220 run_time = time.time() - start_time
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
221 return_code = child.returncode
4
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
222 handle.write("\n")
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
223 handle.write("============================ MIRA has finished ===============================\n")
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
224 handle.write("MIRA took %0.2f hours\n" % (run_time / 3600.0))
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
225 if return_code:
9
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
226 print "MIRA took %0.2f hours" % (run_time / 3600.0)
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
227 handle.write("Return error code %i from command:\n" % return_code)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
228 handle.write(cmd + "\n")
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
229 handle.close()
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
230 clean_up(temp, name)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
231 log_manifest(manifest)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
232 stop_err("Return error code %i from command:\n%s" % (return_code, cmd),
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
233 return_code)
4
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
234 handle.flush()
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
235
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
236 if os.path.isfile("MIRA_assembly/MIRA_d_results/ec.log"):
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
237 handle.write("\n")
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
238 handle.write("====================== Extract Large Contigs failed ==========================\n")
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
239 e = open("MIRA_assembly/MIRA_d_results/ec.log", "rU")
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
240 for line in e:
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
241 handle.write(line)
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
242 e.close()
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
243 handle.write("============================ (end of ec.log) =================================\n")
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
244 handle.flush()
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
245
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
246 #print "Collecting output..."
9
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
247 start_time = time.time()
4
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
248 collect_output(temp, name, handle)
9
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
249 collect_time = time.time() - start_time
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
250 handle.write("MIRA took %0.2f hours; collecting output %0.2f minutes\n" % (run_time / 3600.0, collect_time / 60.0))
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
251 print("MIRA took %0.2f hours; collecting output %0.2f minutes\n" % (run_time / 3600.0, collect_time / 60.0))
4
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
252
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
253 if os.path.isfile("MIRA_assembly/MIRA_d_results/ec.log"):
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
254 #Treat as an error, but doing this AFTER collect_output
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
255 sys.stderr.write("Extract Large Contigs failed\n")
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
256 handle.write("Extract Large Contigs failed\n")
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
257 handle.close()
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
258 sys.exit(1)
0
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
259
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
260 #print "Cleaning up..."
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
261 clean_up(temp, name)
32f693f6e741 Uploaded v0.0.1 preview0, very much a work in progress, primarily checking mira_datatypes dependency
peterjc
parents:
diff changeset
262
9
302d13490b23 Uploaded v0.0.2 preview 1, BAM output
peterjc
parents: 7
diff changeset
263 handle.write("\nDone\n")
4
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
264 handle.close()
df86ed992a1b Uploaded preview 4, lots of work on mapping
peterjc
parents: 0
diff changeset
265 print("Done")