Mercurial > repos > peterjc > mira4_assembler

comparison tools/mira4_0/mira4_make_bam.py @ 32:56b421d59805 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
author peterjc
date Thu, 18 May 2017 13:36:08 -0400
parents fd95aaef8818
children 0785a6537f3e
comparison
equal deleted inserted replaced
31:fd95aaef8818 32:56b421d59805
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
10
9 11
10 def run(cmd, log_handle): 12 def run(cmd, log_handle):
11 try: 13 try:
12 child = subprocess.Popen(cmd, shell=True, 14 child = subprocess.Popen(cmd, shell=True,
13 stdout=subprocess.PIPE, 15 stdout=subprocess.PIPE,
14 stderr=subprocess.STDOUT) 16 stderr=subprocess.STDOUT)
15 except Exception, err: 17 except Exception as err:
16 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err)) 18 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err))
17 #TODO - call clean up? 19 # TODO - call clean up?
18 log_handle.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err)) 20 log_handle.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err))
19 sys.exit(1) 21 sys.exit(1)
20 #Use .communicate as can get deadlocks with .wait(), 22 # Use .communicate as can get deadlocks with .wait(),
21 stdout, stderr = child.communicate() 23 stdout, stderr = child.communicate()
22 assert not stderr #Should be empty as sent to stdout 24 assert not stderr # Should be empty as sent to stdout
23 if len(stdout) > 10000: 25 if len(stdout) > 10000:
24 #miraconvert can be very verbose (is holding stdout in RAM a problem?) 26 # miraconvert can be very verbose (is holding stdout in RAM a problem?)
25 stdout = stdout.split("\n") 27 stdout = stdout.split("\n")
26 stdout = stdout[:10] + ["...", "<snip>", "..."] + stdout[-10:] 28 stdout = stdout[:10] + ["...", "<snip>", "..."] + stdout[-10:]
27 stdout = "\n".join(stdout) 29 stdout = "\n".join(stdout)
28 log_handle.write(stdout) 30 log_handle.write(stdout)
29 return child.returncode 31 return child.returncode
30 32
33
31 def depad(fasta_file, sam_file, bam_file, log_handle): 34 def depad(fasta_file, sam_file, bam_file, log_handle):
32 log_handle.write("\n================= Converting MIRA assembly from SAM to BAM ===================\n") 35 log_handle.write("\n================= Converting MIRA assembly from SAM to BAM ===================\n")
33 #Also doing SAM to (uncompressed) BAM during depad 36 # Also doing SAM to (uncompressed) BAM during depad
34 bam_stem = bam_file + ".tmp" # Have write permissions and want final file in this folder 37 bam_stem = bam_file + ".tmp" # Have write permissions and want final file in this folder
35 cmd = 'samtools depad -S -u -T "%s" "%s" | samtools sort - "%s"' % (fasta_file, sam_file, bam_stem) 38 cmd = 'samtools depad -S -u -T "%s" "%s" | samtools sort - "%s"' % (fasta_file, sam_file, bam_stem)
36 return_code = run(cmd, log_handle) 39 return_code = run(cmd, log_handle)
37 if return_code: 40 if return_code:
38 return "Error %i from command:\n%s" % (return_code, cmd) 41 return "Error %i from command:\n%s" % (return_code, cmd)
39 if not os.path.isfile(bam_stem + ".bam"): 42 if not os.path.isfile(bam_stem + ".bam"):
46 return "Error %i from command:\n%s" % (return_code, cmd) 49 return "Error %i from command:\n%s" % (return_code, cmd)
47 if not os.path.isfile(bam_stem + ".bam.bai"): 50 if not os.path.isfile(bam_stem + ".bam.bai"):
48 return "samtools indexing of BAM file failed to produce BAI file" 51 return "samtools indexing of BAM file failed to produce BAI file"
49 52
50 shutil.move(bam_stem + ".bam", bam_file) 53 shutil.move(bam_stem + ".bam", bam_file)
51 os.remove(bam_stem + ".bam.bai") #Let Galaxy handle that... 54 os.remove(bam_stem + ".bam.bai") # Let Galaxy handle that...
52 55
53 56
54 def make_bam(mira_convert, maf_file, fasta_file, bam_file, log_handle): 57 def make_bam(mira_convert, maf_file, fasta_file, bam_file, log_handle):
55 if not os.path.isfile(mira_convert): 58 if not os.path.isfile(mira_convert):
56 return "Missing binary %r" % mira_convert 59 return "Missing binary %r" % mira_convert
69 if return_code: 72 if return_code:
70 return "Error %i from command:\n%s" % (return_code, cmd) 73 return "Error %i from command:\n%s" % (return_code, cmd)
71 if not os.path.isfile(sam_file): 74 if not os.path.isfile(sam_file):
72 return "Conversion from MIRA to SAM failed" 75 return "Conversion from MIRA to SAM failed"
73 76
74 #Also doing SAM to (uncompressed) BAM during depad 77 # Also doing SAM to (uncompressed) BAM during depad
75 msg = depad(fasta_file, sam_file, bam_file, log_handle) 78 msg = depad(fasta_file, sam_file, bam_file, log_handle)
76 if msg: 79 if msg:
77 return msg 80 return msg
78 81
79 os.remove(sam_file) 82 os.remove(sam_file)
80 os.rmdir(tmp_dir) 83 os.rmdir(tmp_dir)
81 84
82 return None #Good :) 85 return None # Good :)
86
83 87
84 if __name__ == "__main__": 88 if __name__ == "__main__":
85 mira_convert, maf_file, fasta_file, bam_file = sys.argv[1:] 89 mira_convert, maf_file, fasta_file, bam_file = sys.argv[1:]
86 msg = make_bam(mira_convert, maf_file, fasta_file, bam_file, sys.stdout) 90 msg = make_bam(mira_convert, maf_file, fasta_file, bam_file, sys.stdout)
87 if msg: 91 if msg: