# HG changeset patch # User g2cmnty@test-web1.g2.bx.psu.edu # Date 1308687934 14400 # Node ID 6ad594db014300c75defdfb4934e5513a37178f2 # Parent 81a79b4a3f6fc7b02ba183a96ae712378913e7f0 Uploaded diff -r 81a79b4a3f6f -r 6ad594db0143 bam_to_bigwig/README.txt --- a/bam_to_bigwig/README.txt Tue Jun 21 09:44:52 2011 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,18 +0,0 @@ -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 - -Place the script in the same directory as the XML configuration file, or -provide a soft link to it. - -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 81a79b4a3f6f -r 6ad594db0143 bam_to_bigwig/bam_to_bigwig.xml --- a/bam_to_bigwig/bam_to_bigwig.xml Tue Jun 21 09:44:52 2011 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,25 +0,0 @@ - - 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 81a79b4a3f6f -r 6ad594db0143 bam_to_bigwig/bam_to_wiggle.py --- a/bam_to_bigwig/bam_to_wiggle.py Tue Jun 21 09:44:52 2011 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,123 +0,0 @@ -#!/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/) -If a configuration file is used, then PyYAML is also required (http://pyyaml.org/) -""" -import os -import sys -import subprocess -import tempfile -from optparse import OptionParser -from contextlib import contextmanager, closing - -import pysam - -def main(bam_file, config_file=None, chrom='all', start=0, end=None, - outfile=None): - if config_file: - import yaml - 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 os.path.abspath(bam_file) == os.path.abspath(outfile): - sys.stderr.write("Bad arguments, input and output files are the same.\n") - sys.exit(1) - if not (os.path.exists(outfile) and os.path.getsize(outfile) > 0): - #Use a temp file to avoid any possiblity of not having write permission - out_handle = tempfile.NamedTemporaryFile(delete=False) - wig_file = out_handle.name - with closing(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", - ])) - is_valid = False - with indexed_bam(bam_file, config) as work_bam: - sizes = zip(work_bam.references, work_bam.lengths) - if len(regions) == 1 and regions[0][0] == "all": - regions = [(name, 0, length) for name, length in sizes] - for chrom, start, end in regions: - if end is None and chrom in work_bam.references: - end = work_bam.lengths[work_bam.references.index(chrom)] - 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) diff -r 81a79b4a3f6f -r 6ad594db0143 filtering.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/filtering.py Tue Jun 21 16:25:34 2011 -0400 @@ -0,0 +1,128 @@ +#!/usr/bin/env python +# This tool takes a tab-delimited text file as input and creates filters on columns based on certain properties. +# The tool will skip over invalid lines within the file, informing the user about the number of lines skipped. + +from __future__ import division +import sys, re, os.path +from galaxy import eggs + +# Older py compatibility +try: + set() +except: + from sets import Set as set + +assert sys.version_info[:2] >= ( 2, 4 ) + +def get_operands( filter_condition ): + # Note that the order of all_operators is important + items_to_strip = ['+', '-', '**', '*', '//', '/', '%', '<<', '>>', '&', '|', '^', '~', '<=', '<', '>=', '>', '==', '!=', '<>', ' and ', ' or ', ' not ', ' is ', ' is not ', ' in ', ' not in '] + for item in items_to_strip: + if filter_condition.find( item ) >= 0: + filter_condition = filter_condition.replace( item, ' ' ) + operands = set( filter_condition.split( ' ' ) ) + return operands + +def stop_err( msg ): + sys.stderr.write( msg ) + sys.exit() + +in_fname = sys.argv[1] +out_fname = sys.argv[2] +cond_text = sys.argv[3] +try: + in_columns = int( sys.argv[4] ) + assert sys.argv[5] #check to see that the column types varaible isn't null + in_column_types = sys.argv[5].split( ',' ) +except: + stop_err( "Data does not appear to be tabular. This tool can only be used with tab-delimited data." ) + +# Unescape if input has been escaped +mapped_str = { + '__lt__': '<', + '__le__': '<=', + '__eq__': '==', + '__ne__': '!=', + '__gt__': '>', + '__ge__': '>=', + '__sq__': '\'', + '__dq__': '"', +} +for key, value in mapped_str.items(): + cond_text = cond_text.replace( key, value ) + +# Attempt to determine if the condition includes executable stuff and, if so, exit +secured = dir() +operands = get_operands(cond_text) +for operand in operands: + try: + check = int( operand ) + except: + if operand in secured: + stop_err( "Illegal value '%s' in condition '%s'" % ( operand, cond_text ) ) + +# Prepare the column variable names and wrappers for column data types +cols, type_casts = [], [] +for col in range( 1, in_columns + 1 ): + col_name = "c%d" % col + cols.append( col_name ) + col_type = in_column_types[ col - 1 ] + type_cast = "%s(%s)" % ( col_type, col_name ) + type_casts.append( type_cast ) + +col_str = ', '.join( cols ) # 'c1, c2, c3, c4' +type_cast_str = ', '.join( type_casts ) # 'str(c1), int(c2), int(c3), str(c4)' +assign = "%s = line.split( '\\t' )" % col_str +wrap = "%s = %s" % ( col_str, type_cast_str ) +skipped_lines = 0 +first_invalid_line = 0 +invalid_line = None +lines_kept = 0 +total_lines = 0 +out = open( out_fname, 'wt' ) + +# Read and filter input file, skipping invalid lines +code = ''' +for i, line in enumerate( file( in_fname ) ): + total_lines += 1 + line = line.rstrip( '\\r\\n' ) + if not line or line.startswith( '#' ): + skipped_lines += 1 + if not invalid_line: + first_invalid_line = i + 1 + invalid_line = line + continue + try: + %s + %s + if %s: + lines_kept += 1 + print >> out, line + except: + skipped_lines += 1 + if not invalid_line: + first_invalid_line = i + 1 + invalid_line = line +''' % ( assign, wrap, cond_text ) + +valid_filter = True +try: + exec code +except Exception, e: + out.close() + if str( e ).startswith( 'invalid syntax' ): + valid_filter = False + stop_err( 'Filter condition "%s" likely invalid. See tool tips, syntax and examples.' % cond_text ) + else: + stop_err( str( e ) ) + +if valid_filter: + out.close() + valid_lines = total_lines - skipped_lines + print 'Filtering with %s, ' % cond_text + if valid_lines > 0: + print 'kept %4.2f%% of %d lines.' % ( 100.0*lines_kept/valid_lines, total_lines ) + else: + print 'Possible invalid filter condition "%s" or non-existent column referenced. See tool tips, syntax and examples.' % cond_text + if skipped_lines > 0: + print 'Skipped %d invalid lines starting at line #%d: "%s"' % ( skipped_lines, first_invalid_line, invalid_line ) diff -r 81a79b4a3f6f -r 6ad594db0143 filtering.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/filtering.xml Tue Jun 21 16:25:34 2011 -0400 @@ -0,0 +1,65 @@ + + data on any column using simple expressions + + filtering.py $input $out_file1 "$cond" ${input.metadata.columns} "${input.metadata.column_types}" + + + + + + + + + + + + + + + + + + + + + + + + +.. class:: warningmark + +Double equal signs, ==, must be used as *"equal to"* (e.g., **c1 == 'chr22'**) + +.. class:: infomark + +**TIP:** Attempting to apply a filtering condition may throw exceptions if the data type (e.g., string, integer) in every line of the columns being filtered is not appropriate for the condition (e.g., attempting certain numerical calculations on strings). If an exception is thrown when applying the condition to a line, that line is skipped as invalid for the filter condition. The number of invalid skipped lines is documented in the resulting history item as a "Condition/data issue". + +.. class:: infomark + +**TIP:** If your data is not TAB delimited, use *Text Manipulation->Convert* + +----- + +**Syntax** + +The filter tool allows you to restrict the dataset using simple conditional statements. + +- Columns are referenced with **c** and a **number**. For example, **c1** refers to the first column of a tab-delimited file +- Make sure that multi-character operators contain no white space ( e.g., **<=** is valid while **< =** is not valid ) +- When using 'equal-to' operator **double equal sign '==' must be used** ( e.g., **c1=='chr1'** ) +- Non-numerical values must be included in single or double quotes ( e.g., **c6=='+'** ) +- Filtering condition can include logical operators, but **make sure operators are all lower case** ( e.g., **(c1!='chrX' and c1!='chrY') or not c6=='+'** ) + +----- + +**Example** + +- **c1=='chr1'** selects lines in which the first column is chr1 +- **c3-c2<100*c4** selects lines where subtracting column 3 from column 2 is less than the value of column 4 times 100 +- **len(c2.split(',')) < 4** will select lines where the second column has less than four comma separated elements +- **c2>=1** selects lines in which the value of column 2 is greater than or equal to 1 +- Numbers should not contain commas - **c2<=44,554,350** will not work, but **c2<=44554350** will +- Some words in the data can be used, but must be single or double quoted ( e.g., **c3=='exon'** ) + + +