Mercurial > repos > greg > gregs

comparison bam_to_bigwig/bam_to_wiggle.py @ 0:bc71417f212c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author g2cmnty@test-web1.g2.bx.psu.edu
date Fri, 17 Jun 2011 15:25:25 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:bc71417f212c
1 #!/usr/bin/env python
2 """Convert BAM files to BigWig file format in a specified region.
3
4 Usage:
5 bam_to_wiggle.py <BAM file> [<YAML config>]
6 [--outfile=<output file name>
7 --chrom=<chrom>
8 --start=<start>
9 --end=<end>]
10
11 chrom start and end are optional, in which case they default to everything.
12
13 The config file is in YAML format and specifies the location of the wigToBigWig
14 program from UCSC:
15
16 program:
17 ucsc_bigwig: wigToBigWig
18
19 If not specified, these will be assumed to be present in the system path.
20
21 The script requires:
22 pysam (http://code.google.com/p/pysam/)
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/)
25 """
26 import os
27 import sys
28 import subprocess
29 import tempfile
30 from optparse import OptionParser
31 from contextlib import contextmanager, closing
32
33 import pysam
34
35 def main(bam_file, config_file=None, chrom='all', start=0, end=None,
36 outfile=None):
37 if config_file:
38 import yaml
39 with open(config_file) as in_handle:
40 config = yaml.load(in_handle)
41 else:
42 config = {"program": {"ucsc_bigwig" : "wigToBigWig"}}
43 if outfile is None:
44 outfile = "%s.bigwig" % os.path.splitext(bam_file)[0]
45 if start > 0:
46 start = int(start) - 1
47 if end is not None:
48 end = int(end)
49 regions = [(chrom, start, end)]
50 if os.path.abspath(bam_file) == os.path.abspath(outfile):
51 sys.stderr.write("Bad arguments, input and output files are the same.\n")
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):
58 chr_sizes, wig_valid = write_bam_track(bam_file, regions, config, out_handle)
59 try:
60 if wig_valid:
61 convert_to_bigwig(wig_file, chr_sizes, config, outfile)
62 finally:
63 os.remove(wig_file)
64
65 @contextmanager
66 def indexed_bam(bam_file, config):
67 if not os.path.exists(bam_file + ".bai"):
68 pysam.index(bam_file)
69 sam_reader = pysam.Samfile(bam_file, "rb")
70 yield sam_reader
71 sam_reader.close()
72
73 def write_bam_track(bam_file, regions, config, out_handle):
74 out_handle.write("track %s\n" % " ".join(["type=wiggle_0",
75 "name=%s" % os.path.splitext(os.path.split(bam_file)[-1])[0],
76 "visibility=full",
77 ]))
78 is_valid = False
79 with indexed_bam(bam_file, config) as work_bam:
80 sizes = zip(work_bam.references, work_bam.lengths)
81 if len(regions) == 1 and regions[0][0] == "all":
82 regions = [(name, 0, length) for name, length in sizes]
83 for chrom, start, end in regions:
84 if end is None and chrom in work_bam.references:
85 end = work_bam.lengths[work_bam.references.index(chrom)]
86 assert end is not None, "Could not find %s in header" % chrom
87 out_handle.write("variableStep chrom=%s\n" % chrom)
88 for col in work_bam.pileup(chrom, start, end):
89 out_handle.write("%s %s\n" % (col.pos+1, col.n))
90 is_valid = True
91 return sizes, is_valid
92
93 def convert_to_bigwig(wig_file, chr_sizes, config, bw_file=None):
94 if not bw_file:
95 bw_file = "%s.bigwig" % (os.path.splitext(wig_file)[0])
96 size_file = "%s-sizes.txt" % (os.path.splitext(wig_file)[0])
97 with open(size_file, "w") as out_handle:
98 for chrom, size in chr_sizes:
99 out_handle.write("%s\t%s\n" % (chrom, size))
100 try:
101 cl = [config["program"]["ucsc_bigwig"], wig_file, size_file, bw_file]
102 subprocess.check_call(cl)
103 finally:
104 os.remove(size_file)
105 return bw_file
106
107 if __name__ == "__main__":
108 parser = OptionParser()
109 parser.add_option("-o", "--outfile", dest="outfile")
110 parser.add_option("-c", "--chrom", dest="chrom")
111 parser.add_option("-s", "--start", dest="start")
112 parser.add_option("-e", "--end", dest="end")
113 (options, args) = parser.parse_args()
114 if len(args) not in [1, 2]:
115 print "Incorrect arguments"
116 print __doc__
117 sys.exit()
118 kwargs = dict(
119 outfile=options.outfile,
120 chrom=options.chrom or 'all',
121 start=options.start or 0,
122 end=options.end)
123 main(*args, **kwargs)