# HG changeset patch # User brad-chapman # Date 1307476201 14400 # Node ID 20373941a69186d6d864cb4d4baaf2edb7db09ef Migrated tool version 0.0.1 from old tool shed archive to new tool shed repository diff -r 000000000000 -r 20373941a691 bam_to_bigwig/README.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/bam_to_bigwig/README.txt Tue Jun 07 15:50:01 2011 -0400 @@ -0,0 +1,15 @@ +Convert a BAM file into a BigWig coverage file. This can be used directly from +Galaxy for display at UCSC. The advantage over standard Wiggle format is that +the data is stored in a compressed format and can be retrieved by genome +region. This allows you to view regions of arbitrarily large Wiggle file data +at UCSC while avoiding the upload costs. + +The latest version of the bam_to_wiggle.py script is available at: + +https://github.com/chapmanb/bcbb/blob/master/nextgen/scripts/bam_to_wiggle.py + +This requires: + +Python2, version 2.6 or better +pysam (http://code.google.com/p/pysam/) +wigToBigWig from UCSC (http://hgdownload.cse.ucsc.edu/admin/exe/) diff -r 000000000000 -r 20373941a691 bam_to_bigwig/bam_to_bigwig.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/bam_to_bigwig/bam_to_bigwig.xml Tue Jun 07 15:50:01 2011 -0400 @@ -0,0 +1,25 @@ + + Calculates coverage from a BAM alignment file + bam_to_wiggle.py $align --outfile=$out + + + + + + + + +**What it does** + +Creates a coverage file in BigWig format, given a BAM alignment file. + +**Input** + +A BAM alignment file. This needs to have the genome database build used in alignment annotated. If your file has '?' for the database build, click on the pencil icon to edit the alignment attributes, and specify the organism used to align against. + +**Output** + +BigWig files can be loaded directly from Galaxy into the UCSC browser. They can be loaded incrementally by UCSC, so a single file can be used to represent the entire genome without having to upload the entire thing as a custom track. + + + diff -r 000000000000 -r 20373941a691 bam_to_bigwig/bam_to_wiggle.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/bam_to_bigwig/bam_to_wiggle.py Tue Jun 07 15:50:01 2011 -0400 @@ -0,0 +1,123 @@ +#!/usr/bin/env python +"""Convert BAM files to BigWig file format in a specified region. + +Usage: + bam_to_wiggle.py [] + [--outfile= + --chrom= + --start= + --end=] + +chrom start and end are optional, in which case they default to everything. + +The config file is in YAML format and specifies the location of the wigToBigWig +program from UCSC: + +program: + ucsc_bigwig: wigToBigWig + +If not specified, these will be assumed to be present in the system path. + +The script requires: + pysam (http://code.google.com/p/pysam/) + wigToBigWig from UCSC (http://hgdownload.cse.ucsc.edu/admin/exe/) +""" +import os +import sys +import subprocess +from optparse import OptionParser +from contextlib import contextmanager + +import yaml +import pysam + +def main(bam_file, config_file=None, chrom='all', start=0, end=None, + outfile=None): + if config_file: + with open(config_file) as in_handle: + config = yaml.load(in_handle) + else: + config = {"program": {"ucsc_bigwig" : "wigToBigWig"}} + if outfile is None: + outfile = "%s.bigwig" % os.path.splitext(bam_file)[0] + if start > 0: + start = int(start) - 1 + if end is not None: + end = int(end) + regions = [(chrom, start, end)] + if not os.path.exists(outfile): + wig_file = "%s.wig" % os.path.splitext(bam_file)[0] + with open(wig_file, "w") as out_handle: + chr_sizes, wig_valid = write_bam_track(bam_file, regions, config, out_handle) + try: + if wig_valid: + convert_to_bigwig(wig_file, chr_sizes, config, outfile) + finally: + os.remove(wig_file) + +@contextmanager +def indexed_bam(bam_file, config): + if not os.path.exists(bam_file + ".bai"): + pysam.index(bam_file) + sam_reader = pysam.Samfile(bam_file, "rb") + yield sam_reader + sam_reader.close() + +def write_bam_track(bam_file, regions, config, out_handle): + out_handle.write("track %s\n" % " ".join(["type=wiggle_0", + "name=%s" % os.path.splitext(os.path.split(bam_file)[-1])[0], + "visibility=full", + ])) + sizes = [] + is_valid = False + with indexed_bam(bam_file, config) as work_bam: + for ref_info in work_bam.header.get("SQ", []): + sizes.append((ref_info["SN"], ref_info["LN"])) + if len(regions) == 1 and regions[0][0] == "all": + regions = [] + for ref_info in work_bam.header.get("SQ", []): + regions.append((ref_info["SN"], 0, None)) + for chrom, start, end in regions: + if end is None: + for ref_info in work_bam.header.get("SQ", []): + if ref_info["SN"] == chrom: + end = int(ref_info["LN"]) + break + assert end is not None, "Could not find %s in header" % chrom + out_handle.write("variableStep chrom=%s\n" % chrom) + for col in work_bam.pileup(chrom, start, end): + out_handle.write("%s %s\n" % (col.pos+1, col.n)) + is_valid = True + return sizes, is_valid + +def convert_to_bigwig(wig_file, chr_sizes, config, bw_file=None): + if not bw_file: + bw_file = "%s.bigwig" % (os.path.splitext(wig_file)[0]) + size_file = "%s-sizes.txt" % (os.path.splitext(wig_file)[0]) + with open(size_file, "w") as out_handle: + for chrom, size in chr_sizes: + out_handle.write("%s\t%s\n" % (chrom, size)) + try: + cl = [config["program"]["ucsc_bigwig"], wig_file, size_file, bw_file] + subprocess.check_call(cl) + finally: + os.remove(size_file) + return bw_file + +if __name__ == "__main__": + parser = OptionParser() + parser.add_option("-o", "--outfile", dest="outfile") + parser.add_option("-c", "--chrom", dest="chrom") + parser.add_option("-s", "--start", dest="start") + parser.add_option("-e", "--end", dest="end") + (options, args) = parser.parse_args() + if len(args) not in [1, 2]: + print "Incorrect arguments" + print __doc__ + sys.exit() + kwargs = dict( + outfile=options.outfile, + chrom=options.chrom or 'all', + start=options.start or 0, + end=options.end) + main(*args, **kwargs)