Mercurial > repos > peterjc > mira4_assembler

annotate tools/mira4/mira4_make_bam.py @ 21:4abe8d59a438 draft

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