Mercurial > repos > peterjc > mira4_assembler

comparison tools/mira4/mira4_make_bam.py @ 9:302d13490b23 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded v0.0.2 preview 1, BAM output
author peterjc
date Thu, 28 Nov 2013 05:07:59 -0500
parents
children 4abe8d59a438
comparison
equal deleted inserted replaced
8:2ab1d6f6a8ec 9:302d13490b23
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 make_bam(mira_convert, maf_file, fasta_file, bam_file, log_handle):
36 if not os.path.isfile(mira_convert):
37 return "Missing binary %r" % mira_convert
38 if not os.path.isfile(maf_file):
39 return "Missing input MIRA file: %r" % maf_file
40 if not os.path.isfile(fasta_file):
41 return "Missing padded FASTA file: %r" % fasta_file
42
43 log_handle.write("\n====================== Converting MIRA assembly to SAM =======================\n")
44 tmp_dir = tempfile.mkdtemp()
45 sam_file = os.path.join(tmp_dir, "x.sam")
46
47 # Note add nbb to the template name, possible MIRA 4.0 RC4 bug
48 cmd = '"%s" -f maf -t samnbb "%s" "%snbb"' % (mira_convert, maf_file, sam_file)
49 return_code = run(cmd, log_handle)
50 if return_code:
51 return "Error %i from command:\n%s" % (return_code, cmd)
52 if not os.path.isfile(sam_file):
53 return "Conversion from MIRA to SAM failed"
54
55 log_handle.write("\n================= Converting MIRA assembly from SAM to BAM ===================\n")
56 #Also doing SAM to (uncompressed) BAM during depad
57 bam_stem = bam_file + ".tmp" # Have write permissions and want final file in this folder
58 cmd = 'samtools depad -S -u -T "%s" "%s" | samtools sort - "%s"' % (fasta_file, sam_file, bam_stem)
59 return_code = run(cmd, log_handle)
60 if return_code:
61 return "Error %i from command:\n%s" % (return_code, cmd)
62 if not os.path.isfile(bam_stem + ".bam"):
63 return "samtools depad or sort failed to produce BAM file"
64
65 os.remove(sam_file)
66 os.rmdir(tmp_dir)
67
68 log_handle.write("\n====================== Indexing MIRA assembly BAM file =======================\n")
69 cmd = 'samtools index "%s.bam"' % bam_stem
70 return_code = run(cmd, log_handle)
71 if return_code:
72 return "Error %i from command:\n%s" % (return_code, cmd)
73 if not os.path.isfile(bam_stem + ".bam.bai"):
74 return "samtools indexing of BAM file failed to produce BAI file"
75
76 shutil.move(bam_stem + ".bam", bam_file)
77 os.remove(bam_stem + ".bam.bai") #Let Galaxy handle that...
78
79 return None #Good :)
80
81 if __name__ == "__main__":
82 mira_convert, maf_file, fasta_file, bam_file = sys.argv[1:]
83 msg = make_bam(mira_convert, maf_file, fasta_file, bam_file, sys.stdout)
84 if msg:
85 stop_err(msg)