Mercurial > repos > peterjc > mira4_assembler

comparison tools/mira4/mira4_make_bam.py @ 0:6a88b42ce6b9 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded v0.0.4, previously only on the TestToolShed
author peterjc
date Fri, 21 Nov 2014 06:42:56 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:6a88b42ce6b9
1 #!/usr/bin/env python
2 """Wrapper script using miraconvert & samtools to get BAM from MIRA.
3 """
4 import os
5 import sys
6 import shutil
7 import subprocess
8 import tempfile
9
10 def stop_err(msg, err=1):
11 sys.stderr.write(msg+"\n")
12 sys.exit(err)
13
14 def run(cmd, log_handle):
15 try:
16 child = subprocess.Popen(cmd, shell=True,
17 stdout=subprocess.PIPE,
18 stderr=subprocess.STDOUT)
19 except Exception, err:
20 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err))
21 #TODO - call clean up?
22 log_handle.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err))
23 sys.exit(1)
24 #Use .communicate as can get deadlocks with .wait(),
25 stdout, stderr = child.communicate()
26 assert not stderr #Should be empty as sent to stdout
27 if len(stdout) > 10000:
28 #miraconvert can be very verbose (is holding stdout in RAM a problem?)
29 stdout = stdout.split("\n")
30 stdout = stdout[:10] + ["...", "<snip>", "..."] + stdout[-10:]
31 stdout = "\n".join(stdout)
32 log_handle.write(stdout)
33 return child.returncode
34
35 def depad(fasta_file, sam_file, bam_file, log_handle):
36 log_handle.write("\n================= Converting MIRA assembly from SAM to BAM ===================\n")
37 #Also doing SAM to (uncompressed) BAM during depad
38 bam_stem = bam_file + ".tmp" # Have write permissions and want final file in this folder
39 cmd = 'samtools depad -S -u -T "%s" "%s" | samtools sort - "%s"' % (fasta_file, sam_file, bam_stem)
40 return_code = run(cmd, log_handle)
41 if return_code:
42 return "Error %i from command:\n%s" % (return_code, cmd)
43 if not os.path.isfile(bam_stem + ".bam"):
44 return "samtools depad or sort failed to produce BAM file"
45
46 log_handle.write("\n====================== Indexing MIRA assembly BAM file =======================\n")
47 cmd = 'samtools index "%s.bam"' % bam_stem
48 return_code = run(cmd, log_handle)
49 if return_code:
50 return "Error %i from command:\n%s" % (return_code, cmd)
51 if not os.path.isfile(bam_stem + ".bam.bai"):
52 return "samtools indexing of BAM file failed to produce BAI file"
53
54 shutil.move(bam_stem + ".bam", bam_file)
55 os.remove(bam_stem + ".bam.bai") #Let Galaxy handle that...
56
57
58 def make_bam(mira_convert, maf_file, fasta_file, bam_file, log_handle):
59 if not os.path.isfile(mira_convert):
60 return "Missing binary %r" % mira_convert
61 if not os.path.isfile(maf_file):
62 return "Missing input MIRA file: %r" % maf_file
63 if not os.path.isfile(fasta_file):
64 return "Missing padded FASTA file: %r" % fasta_file
65
66 log_handle.write("\n====================== Converting MIRA assembly to SAM =======================\n")
67 tmp_dir = tempfile.mkdtemp()
68 sam_file = os.path.join(tmp_dir, "x.sam")
69
70 # Note add nbb to the template name, possible MIRA 4.0 RC4 bug
71 cmd = '"%s" -f maf -t samnbb "%s" "%snbb"' % (mira_convert, maf_file, sam_file)
72 return_code = run(cmd, log_handle)
73 if return_code:
74 return "Error %i from command:\n%s" % (return_code, cmd)
75 if not os.path.isfile(sam_file):
76 return "Conversion from MIRA to SAM failed"
77
78 #Also doing SAM to (uncompressed) BAM during depad
79 msg = depad(fasta_file, sam_file, bam_file, log_handle)
80 if msg:
81 return msg
82
83 os.remove(sam_file)
84 os.rmdir(tmp_dir)
85
86 return None #Good :)
87
88 if __name__ == "__main__":
89 mira_convert, maf_file, fasta_file, bam_file = sys.argv[1:]
90 msg = make_bam(mira_convert, maf_file, fasta_file, bam_file, sys.stdout)
91 if msg:
92 stop_err(msg)