Mercurial > repos > peterjc > mira4_assembler
annotate tools/mira4_0/mira4_make_bam.py @ 35:259891fce7fd draft
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit e4c56df75150c82d3e9c4ac487c4209108e52412
| author | peterjc |
|---|---|
| date | Wed, 09 Aug 2017 13:32:11 -0400 |
| parents | 0785a6537f3e |
| children | cee8f9005e43 |
| rev | line source |
|---|---|
| 25 | 1 #!/usr/bin/env python |
| 2 """Wrapper script using miraconvert & samtools to get BAM from MIRA. | |
| 3 """ | |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
4 |
| 25 | 5 import os |
| 6 import shutil | |
| 7 import subprocess | |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
8 import sys |
| 25 | 9 import tempfile |
| 10 | |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
11 |
| 25 | 12 def run(cmd, log_handle): |
| 13 try: | |
| 14 child = subprocess.Popen(cmd, shell=True, | |
|
34
0785a6537f3e
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit 6405ba93fcec7ea93452bf54d559c7507ee7a57c
peterjc
parents:
32
diff
changeset
|
15 universal_newlines=True, |
| 25 | 16 stdout=subprocess.PIPE, |
| 17 stderr=subprocess.STDOUT) | |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
18 except Exception as err: |
| 25 | 19 sys.stderr.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err)) |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
20 # TODO - call clean up? |
| 25 | 21 log_handle.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err)) |
| 22 sys.exit(1) | |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
23 # Use .communicate as can get deadlocks with .wait(), |
| 25 | 24 stdout, stderr = child.communicate() |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
25 assert not stderr # Should be empty as sent to stdout |
| 25 | 26 if len(stdout) > 10000: |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
27 # miraconvert can be very verbose (is holding stdout in RAM a problem?) |
| 25 | 28 stdout = stdout.split("\n") |
| 29 stdout = stdout[:10] + ["...", "<snip>", "..."] + stdout[-10:] | |
| 30 stdout = "\n".join(stdout) | |
| 31 log_handle.write(stdout) | |
| 32 return child.returncode | |
| 33 | |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
34 |
| 25 | 35 def depad(fasta_file, sam_file, bam_file, log_handle): |
| 36 log_handle.write("\n================= Converting MIRA assembly from SAM to BAM ===================\n") | |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
37 # Also doing SAM to (uncompressed) BAM during depad |
|
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
38 bam_stem = bam_file + ".tmp" # Have write permissions and want final file in this folder |
| 25 | 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) | |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
55 os.remove(bam_stem + ".bam.bai") # Let Galaxy handle that... |
| 25 | 56 |
| 57 | |
| 58 def make_bam(mira_convert, maf_file, fasta_file, bam_file, log_handle): | |
| 59 if not os.path.isfile(maf_file): | |
| 60 return "Missing input MIRA file: %r" % maf_file | |
| 61 if not os.path.isfile(fasta_file): | |
| 62 return "Missing padded FASTA file: %r" % fasta_file | |
| 63 | |
| 64 log_handle.write("\n====================== Converting MIRA assembly to SAM =======================\n") | |
| 65 tmp_dir = tempfile.mkdtemp() | |
| 66 sam_file = os.path.join(tmp_dir, "x.sam") | |
| 67 | |
| 68 # Note add nbb to the template name, possible MIRA 4.0 RC4 bug | |
| 69 cmd = '"%s" -f maf -t samnbb "%s" "%snbb"' % (mira_convert, maf_file, sam_file) | |
| 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(sam_file): | |
| 74 return "Conversion from MIRA to SAM failed" | |
| 75 | |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
76 # Also doing SAM to (uncompressed) BAM during depad |
| 25 | 77 msg = depad(fasta_file, sam_file, bam_file, log_handle) |
| 78 if msg: | |
| 79 return msg | |
| 80 | |
| 81 os.remove(sam_file) | |
| 82 os.rmdir(tmp_dir) | |
| 83 | |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
84 return None # Good :) |
|
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
85 |
| 25 | 86 |
| 87 if __name__ == "__main__": | |
| 88 mira_convert, maf_file, fasta_file, bam_file = sys.argv[1:] | |
| 89 msg = make_bam(mira_convert, maf_file, fasta_file, bam_file, sys.stdout) | |
| 90 if msg: | |
|
31
fd95aaef8818
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit bc3d484c5cd68ddcf456db2fff489d584aa2034c
peterjc
parents:
25
diff
changeset
|
91 sys.exit(msg) |