Mercurial > repos > peterjc > mira4_assembler
annotate tools/mira4_0/mira4_make_bam.py @ 32:56b421d59805 draft
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 |
| 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, | |
| 15 stdout=subprocess.PIPE, | |
| 16 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
|
17 except Exception as err: |
| 25 | 18 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
|
19 # TODO - call clean up? |
| 25 | 20 log_handle.write("Error invoking command:\n%s\n\n%s\n" % (cmd, err)) |
| 21 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
|
22 # Use .communicate as can get deadlocks with .wait(), |
| 25 | 23 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
|
24 assert not stderr # Should be empty as sent to stdout |
| 25 | 25 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
|
26 # miraconvert can be very verbose (is holding stdout in RAM a problem?) |
| 25 | 27 stdout = stdout.split("\n") |
| 28 stdout = stdout[:10] + ["...", "<snip>", "..."] + stdout[-10:] | |
| 29 stdout = "\n".join(stdout) | |
| 30 log_handle.write(stdout) | |
| 31 return child.returncode | |
| 32 | |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
33 |
| 25 | 34 def depad(fasta_file, sam_file, bam_file, log_handle): |
| 35 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
|
36 # 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
|
37 bam_stem = bam_file + ".tmp" # Have write permissions and want final file in this folder |
| 25 | 38 cmd = 'samtools depad -S -u -T "%s" "%s" | samtools sort - "%s"' % (fasta_file, sam_file, bam_stem) |
| 39 return_code = run(cmd, log_handle) | |
| 40 if return_code: | |
| 41 return "Error %i from command:\n%s" % (return_code, cmd) | |
| 42 if not os.path.isfile(bam_stem + ".bam"): | |
| 43 return "samtools depad or sort failed to produce BAM file" | |
| 44 | |
| 45 log_handle.write("\n====================== Indexing MIRA assembly BAM file =======================\n") | |
| 46 cmd = 'samtools index "%s.bam"' % bam_stem | |
| 47 return_code = run(cmd, log_handle) | |
| 48 if return_code: | |
| 49 return "Error %i from command:\n%s" % (return_code, cmd) | |
| 50 if not os.path.isfile(bam_stem + ".bam.bai"): | |
| 51 return "samtools indexing of BAM file failed to produce BAI file" | |
| 52 | |
| 53 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
|
54 os.remove(bam_stem + ".bam.bai") # Let Galaxy handle that... |
| 25 | 55 |
| 56 | |
| 57 def make_bam(mira_convert, maf_file, fasta_file, bam_file, log_handle): | |
| 58 if not os.path.isfile(mira_convert): | |
| 59 return "Missing binary %r" % mira_convert | |
| 60 if not os.path.isfile(maf_file): | |
| 61 return "Missing input MIRA file: %r" % maf_file | |
| 62 if not os.path.isfile(fasta_file): | |
| 63 return "Missing padded FASTA file: %r" % fasta_file | |
| 64 | |
| 65 log_handle.write("\n====================== Converting MIRA assembly to SAM =======================\n") | |
| 66 tmp_dir = tempfile.mkdtemp() | |
| 67 sam_file = os.path.join(tmp_dir, "x.sam") | |
| 68 | |
| 69 # Note add nbb to the template name, possible MIRA 4.0 RC4 bug | |
| 70 cmd = '"%s" -f maf -t samnbb "%s" "%snbb"' % (mira_convert, maf_file, sam_file) | |
| 71 return_code = run(cmd, log_handle) | |
| 72 if return_code: | |
| 73 return "Error %i from command:\n%s" % (return_code, cmd) | |
| 74 if not os.path.isfile(sam_file): | |
| 75 return "Conversion from MIRA to SAM failed" | |
| 76 | |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
77 # Also doing SAM to (uncompressed) BAM during depad |
| 25 | 78 msg = depad(fasta_file, sam_file, bam_file, log_handle) |
| 79 if msg: | |
| 80 return msg | |
| 81 | |
| 82 os.remove(sam_file) | |
| 83 os.rmdir(tmp_dir) | |
| 84 | |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
85 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
|
86 |
| 25 | 87 |
| 88 if __name__ == "__main__": | |
| 89 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) | |
| 91 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
|
92 sys.exit(msg) |