Mercurial > repos > brad-chapman > bam_to_bigwig

comparison bam_to_bigwig/bam_to_wiggle.py @ 1:0ff100a057ef

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Migrated tool version 0.0.2 from old tool shed archive to new tool shed repository
author brad-chapman
date Tue, 07 Jun 2011 16:26:46 -0400
parents d2c1af657010
children
comparison
equal deleted inserted replaced
0:d2c1af657010 1:0ff100a057ef
19 If not specified, these will be assumed to be present in the system path. 19 If not specified, these will be assumed to be present in the system path.
20 20
21 The script requires: 21 The script requires:
22 pysam (http://code.google.com/p/pysam/) 22 pysam (http://code.google.com/p/pysam/)
23 wigToBigWig from UCSC (http://hgdownload.cse.ucsc.edu/admin/exe/) 23 wigToBigWig from UCSC (http://hgdownload.cse.ucsc.edu/admin/exe/)
24 If a configuration file is used, then PyYAML is also required (http://pyyaml.org/)
24 """ 25 """
25 import os 26 import os
26 import sys 27 import sys
27 import subprocess 28 import subprocess
29 import tempfile
28 from optparse import OptionParser 30 from optparse import OptionParser
29 from contextlib import contextmanager 31 from contextlib import contextmanager, closing
30 32
31 import yaml
32 import pysam 33 import pysam
33 34
34 def main(bam_file, config_file=None, chrom='all', start=0, end=None, 35 def main(bam_file, config_file=None, chrom='all', start=0, end=None,
35 outfile=None): 36 outfile=None):
36 if config_file: 37 if config_file:
38 import yaml
37 with open(config_file) as in_handle: 39 with open(config_file) as in_handle:
38 config = yaml.load(in_handle) 40 config = yaml.load(in_handle)
39 else: 41 else:
40 config = {"program": {"ucsc_bigwig" : "wigToBigWig"}} 42 config = {"program": {"ucsc_bigwig" : "wigToBigWig"}}
41 if outfile is None: 43 if outfile is None:
43 if start > 0: 45 if start > 0:
44 start = int(start) - 1 46 start = int(start) - 1
45 if end is not None: 47 if end is not None:
46 end = int(end) 48 end = int(end)
47 regions = [(chrom, start, end)] 49 regions = [(chrom, start, end)]
48 if not os.path.exists(outfile): 50 if os.path.abspath(bam_file) == os.path.abspath(outfile):
49 wig_file = "%s.wig" % os.path.splitext(bam_file)[0] 51 sys.stderr.write("Bad arguments, input and output files are the same.\n")
50 with open(wig_file, "w") as out_handle: 52 sys.exit(1)
53 if not (os.path.exists(outfile) and os.path.getsize(outfile) > 0):
54 #Use a temp file to avoid any possiblity of not having write permission
55 out_handle = tempfile.NamedTemporaryFile(delete=False)
56 wig_file = out_handle.name
57 with closing(out_handle):
51 chr_sizes, wig_valid = write_bam_track(bam_file, regions, config, out_handle) 58 chr_sizes, wig_valid = write_bam_track(bam_file, regions, config, out_handle)
52 try: 59 try:
53 if wig_valid: 60 if wig_valid:
54 convert_to_bigwig(wig_file, chr_sizes, config, outfile) 61 convert_to_bigwig(wig_file, chr_sizes, config, outfile)
55 finally: 62 finally:
66 def write_bam_track(bam_file, regions, config, out_handle): 73 def write_bam_track(bam_file, regions, config, out_handle):
67 out_handle.write("track %s\n" % " ".join(["type=wiggle_0", 74 out_handle.write("track %s\n" % " ".join(["type=wiggle_0",
68 "name=%s" % os.path.splitext(os.path.split(bam_file)[-1])[0], 75 "name=%s" % os.path.splitext(os.path.split(bam_file)[-1])[0],
69 "visibility=full", 76 "visibility=full",
70 ])) 77 ]))
71 sizes = []
72 is_valid = False 78 is_valid = False
73 with indexed_bam(bam_file, config) as work_bam: 79 with indexed_bam(bam_file, config) as work_bam:
74 for ref_info in work_bam.header.get("SQ", []): 80 sizes = zip(work_bam.references, work_bam.lengths)
75 sizes.append((ref_info["SN"], ref_info["LN"]))
76 if len(regions) == 1 and regions[0][0] == "all": 81 if len(regions) == 1 and regions[0][0] == "all":
77 regions = [] 82 regions = [(name, 0, length) for name, length in sizes]
78 for ref_info in work_bam.header.get("SQ", []):
79 regions.append((ref_info["SN"], 0, None))
80 for chrom, start, end in regions: 83 for chrom, start, end in regions:
81 if end is None: 84 if end is None and chrom in work_bam.references:
82 for ref_info in work_bam.header.get("SQ", []): 85 end = work_bam.lengths[work_bam.references.index(chrom)]
83 if ref_info["SN"] == chrom:
84 end = int(ref_info["LN"])
85 break
86 assert end is not None, "Could not find %s in header" % chrom 86 assert end is not None, "Could not find %s in header" % chrom
87 out_handle.write("variableStep chrom=%s\n" % chrom) 87 out_handle.write("variableStep chrom=%s\n" % chrom)
88 for col in work_bam.pileup(chrom, start, end): 88 for col in work_bam.pileup(chrom, start, end):
89 out_handle.write("%s %s\n" % (col.pos+1, col.n)) 89 out_handle.write("%s %s\n" % (col.pos+1, col.n))
90 is_valid = True 90 is_valid = True