Mercurial > repos > peterjc > mira4_assembler

comparison 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
comparison
equal deleted inserted replaced
20:aeb3e35f8236 21:4abe8d59a438
30 stdout = stdout[:10] + ["...", "<snip>", "..."] + stdout[-10:] 30 stdout = stdout[:10] + ["...", "<snip>", "..."] + stdout[-10:]
31 stdout = "\n".join(stdout) 31 stdout = "\n".join(stdout)
32 log_handle.write(stdout) 32 log_handle.write(stdout)
33 return child.returncode 33 return child.returncode
34 34
35 def depad(fasta_file, sam_file, bam_file, log_handle):
36 log_handle.write("\n================= Converting MIRA assembly from SAM to BAM ===================\n")
37 #Also doing SAM to (uncompressed) BAM during depad
38 bam_stem = bam_file + ".tmp" # Have write permissions and want final file in this folder
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)
55 os.remove(bam_stem + ".bam.bai") #Let Galaxy handle that...
56
57
35 def make_bam(mira_convert, maf_file, fasta_file, bam_file, log_handle): 58 def make_bam(mira_convert, maf_file, fasta_file, bam_file, log_handle):
36 if not os.path.isfile(mira_convert): 59 if not os.path.isfile(mira_convert):
37 return "Missing binary %r" % mira_convert 60 return "Missing binary %r" % mira_convert
38 if not os.path.isfile(maf_file): 61 if not os.path.isfile(maf_file):
39 return "Missing input MIRA file: %r" % maf_file 62 return "Missing input MIRA file: %r" % maf_file
50 if return_code: 73 if return_code:
51 return "Error %i from command:\n%s" % (return_code, cmd) 74 return "Error %i from command:\n%s" % (return_code, cmd)
52 if not os.path.isfile(sam_file): 75 if not os.path.isfile(sam_file):
53 return "Conversion from MIRA to SAM failed" 76 return "Conversion from MIRA to SAM failed"
54 77
55 log_handle.write("\n================= Converting MIRA assembly from SAM to BAM ===================\n")
56 #Also doing SAM to (uncompressed) BAM during depad 78 #Also doing SAM to (uncompressed) BAM during depad
57 bam_stem = bam_file + ".tmp" # Have write permissions and want final file in this folder 79 msg = depad(fasta_file, sam_file, bam_file, log_handle)
58 cmd = 'samtools depad -S -u -T "%s" "%s" | samtools sort - "%s"' % (fasta_file, sam_file, bam_stem) 80 if msg:
59 return_code = run(cmd, log_handle) 81 return msg
60 if return_code:
61 return "Error %i from command:\n%s" % (return_code, cmd)
62 if not os.path.isfile(bam_stem + ".bam"):
63 return "samtools depad or sort failed to produce BAM file"
64 82
65 os.remove(sam_file) 83 os.remove(sam_file)
66 os.rmdir(tmp_dir) 84 os.rmdir(tmp_dir)
67
68 log_handle.write("\n====================== Indexing MIRA assembly BAM file =======================\n")
69 cmd = 'samtools index "%s.bam"' % bam_stem
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(bam_stem + ".bam.bai"):
74 return "samtools indexing of BAM file failed to produce BAI file"
75
76 shutil.move(bam_stem + ".bam", bam_file)
77 os.remove(bam_stem + ".bam.bai") #Let Galaxy handle that...
78 85
79 return None #Good :) 86 return None #Good :)
80 87
81 if __name__ == "__main__": 88 if __name__ == "__main__":
82 mira_convert, maf_file, fasta_file, bam_file = sys.argv[1:] 89 mira_convert, maf_file, fasta_file, bam_file = sys.argv[1:]