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