Mercurial > repos > peterjc > mira4_assembler
annotate tools/mira4_0/mira4_make_bam.py @ 38:cee8f9005e43 draft
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit 206259620376b322fc8ed99a6efdd3712f38764b
| author | peterjc |
|---|---|
| date | Wed, 11 Jul 2018 12:35:35 -0400 |
| parents | 259891fce7fd |
| children | bbf14bb9607b |
| rev | line source |
|---|---|
| 25 | 1 #!/usr/bin/env python |
|
38
cee8f9005e43
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit 206259620376b322fc8ed99a6efdd3712f38764b
peterjc
parents:
35
diff
changeset
|
2 """Wrapper script using miraconvert & samtools to get BAM from MIRA.""" |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
3 |
| 25 | 4 import os |
| 5 import shutil | |
| 6 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
|
7 import sys |
| 25 | 8 import tempfile |
| 9 | |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
10 |
| 25 | 11 def run(cmd, log_handle): |
| 12 try: | |
| 13 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
|
14 universal_newlines=True, |
| 25 | 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(maf_file): | |
| 59 return "Missing input MIRA file: %r" % maf_file | |
| 60 if not os.path.isfile(fasta_file): | |
| 61 return "Missing padded FASTA file: %r" % fasta_file | |
| 62 | |
| 63 log_handle.write("\n====================== Converting MIRA assembly to SAM =======================\n") | |
| 64 tmp_dir = tempfile.mkdtemp() | |
| 65 sam_file = os.path.join(tmp_dir, "x.sam") | |
| 66 | |
| 67 # Note add nbb to the template name, possible MIRA 4.0 RC4 bug | |
| 68 cmd = '"%s" -f maf -t samnbb "%s" "%snbb"' % (mira_convert, maf_file, sam_file) | |
| 69 return_code = run(cmd, log_handle) | |
| 70 if return_code: | |
| 71 return "Error %i from command:\n%s" % (return_code, cmd) | |
| 72 if not os.path.isfile(sam_file): | |
| 73 return "Conversion from MIRA to SAM failed" | |
| 74 | |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
75 # Also doing SAM to (uncompressed) BAM during depad |
| 25 | 76 msg = depad(fasta_file, sam_file, bam_file, log_handle) |
| 77 if msg: | |
| 78 return msg | |
| 79 | |
| 80 os.remove(sam_file) | |
| 81 os.rmdir(tmp_dir) | |
| 82 | |
|
32
56b421d59805
planemo upload for repository https://github.com/peterjc/galaxy_mira/tree/master/tools/mira4_0 commit fd979d17340cde155de176604744831d9597c6b6
peterjc
parents:
31
diff
changeset
|
83 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
|
84 |
| 25 | 85 |
| 86 if __name__ == "__main__": | |
| 87 mira_convert, maf_file, fasta_file, bam_file = sys.argv[1:] | |
| 88 msg = make_bam(mira_convert, maf_file, fasta_file, bam_file, sys.stdout) | |
| 89 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
|
90 sys.exit(msg) |