Mercurial > repos > peterjc > mira4_assembler

comparison tools/mira4_0/mira4_make_bam.py @ 4:1713289d9908 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
v0.0.11 tweak for use with bioconda dependencies
author peterjc
date Thu, 10 Aug 2017 11:09:10 -0400
parents 4eb32a3d67d1
children
comparison
equal deleted inserted replaced
3:a4f602cc3aa9 4:1713289d9908
1 #!/usr/bin/env python 1 #!/usr/bin/env python
2 """Wrapper script using miraconvert & samtools to get BAM from MIRA. 2 """Wrapper script using miraconvert & samtools to get BAM from MIRA.
3 """ 3 """
4
4 import os 5 import os
5 import sys
6 import shutil 6 import shutil
7 import subprocess 7 import subprocess
8 import sys
8 import tempfile 9 import tempfile
9 10
10 def sys_exit(msg, err=1):
11 sys.stderr.write(msg+"\n")
12 sys.exit(err)
13 11
14 def run(cmd, log_handle): 12 def run(cmd, log_handle):
15 try: 13 try:
16 child = subprocess.Popen(cmd, shell=True, 14 child = subprocess.Popen(cmd, shell=True,
15 universal_newlines=True,
17 stdout=subprocess.PIPE, 16 stdout=subprocess.PIPE,
18 stderr=subprocess.STDOUT) 17 stderr=subprocess.STDOUT)
19 except Exception, err: 18 except Exception as err:
20 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err)) 19 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err))
21 #TODO - call clean up? 20 # TODO - call clean up?
22 log_handle.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err)) 21 log_handle.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err))
23 sys.exit(1) 22 sys.exit(1)
24 #Use .communicate as can get deadlocks with .wait(), 23 # Use .communicate as can get deadlocks with .wait(),
25 stdout, stderr = child.communicate() 24 stdout, stderr = child.communicate()
26 assert not stderr #Should be empty as sent to stdout 25 assert not stderr # Should be empty as sent to stdout
27 if len(stdout) > 10000: 26 if len(stdout) > 10000:
28 #miraconvert can be very verbose (is holding stdout in RAM a problem?) 27 # miraconvert can be very verbose (is holding stdout in RAM a problem?)
29 stdout = stdout.split("\n") 28 stdout = stdout.split("\n")
30 stdout = stdout[:10] + ["...", "<snip>", "..."] + stdout[-10:] 29 stdout = stdout[:10] + ["...", "<snip>", "..."] + stdout[-10:]
31 stdout = "\n".join(stdout) 30 stdout = "\n".join(stdout)
32 log_handle.write(stdout) 31 log_handle.write(stdout)
33 return child.returncode 32 return child.returncode
34 33
34
35 def depad(fasta_file, sam_file, bam_file, log_handle): 35 def depad(fasta_file, sam_file, bam_file, log_handle):
36 log_handle.write("\n================= Converting MIRA assembly from SAM to BAM ===================\n") 36 log_handle.write("\n================= Converting MIRA assembly from SAM to BAM ===================\n")
37 #Also doing SAM to (uncompressed) BAM during depad 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 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) 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) 40 return_code = run(cmd, log_handle)
41 if return_code: 41 if return_code:
42 return "Error %i from command:\n%s" % (return_code, cmd) 42 return "Error %i from command:\n%s" % (return_code, cmd)
43 if not os.path.isfile(bam_stem + ".bam"): 43 if not os.path.isfile(bam_stem + ".bam"):
50 return "Error %i from command:\n%s" % (return_code, cmd) 50 return "Error %i from command:\n%s" % (return_code, cmd)
51 if not os.path.isfile(bam_stem + ".bam.bai"): 51 if not os.path.isfile(bam_stem + ".bam.bai"):
52 return "samtools indexing of BAM file failed to produce BAI file" 52 return "samtools indexing of BAM file failed to produce BAI file"
53 53
54 shutil.move(bam_stem + ".bam", bam_file) 54 shutil.move(bam_stem + ".bam", bam_file)
55 os.remove(bam_stem + ".bam.bai") #Let Galaxy handle that... 55 os.remove(bam_stem + ".bam.bai") # Let Galaxy handle that...
56 56
57 57
58 def make_bam(mira_convert, maf_file, fasta_file, bam_file, log_handle): 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): 59 if not os.path.isfile(maf_file):
62 return "Missing input MIRA file: %r" % maf_file 60 return "Missing input MIRA file: %r" % maf_file
63 if not os.path.isfile(fasta_file): 61 if not os.path.isfile(fasta_file):
64 return "Missing padded FASTA file: %r" % fasta_file 62 return "Missing padded FASTA file: %r" % fasta_file
65 63
73 if return_code: 71 if return_code:
74 return "Error %i from command:\n%s" % (return_code, cmd) 72 return "Error %i from command:\n%s" % (return_code, cmd)
75 if not os.path.isfile(sam_file): 73 if not os.path.isfile(sam_file):
76 return "Conversion from MIRA to SAM failed" 74 return "Conversion from MIRA to SAM failed"
77 75
78 #Also doing SAM to (uncompressed) BAM during depad 76 # Also doing SAM to (uncompressed) BAM during depad
79 msg = depad(fasta_file, sam_file, bam_file, log_handle) 77 msg = depad(fasta_file, sam_file, bam_file, log_handle)
80 if msg: 78 if msg:
81 return msg 79 return msg
82 80
83 os.remove(sam_file) 81 os.remove(sam_file)
84 os.rmdir(tmp_dir) 82 os.rmdir(tmp_dir)
85 83
86 return None #Good :) 84 return None # Good :)
85
87 86
88 if __name__ == "__main__": 87 if __name__ == "__main__":
89 mira_convert, maf_file, fasta_file, bam_file = sys.argv[1:] 88 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) 89 msg = make_bam(mira_convert, maf_file, fasta_file, bam_file, sys.stdout)
91 if msg: 90 if msg:
92 sys_exit(msg) 91 sys.exit(msg)