From 240fb5933e2eedbfa87ee7d3501b1a7f94bfebb4 Mon Sep 17 00:00:00 2001 From: Jeremy Goecks Date: Mon, 9 Aug 2010 22:37:43 -0400 Subject: [PATCH] Enhance GFF-to-BED converter to output block data. Added test for new functionality. --- lib/galaxy/tools/util/gff_util.py | 20 +++++- tools/filters/gff2bed.xml | 6 +- tools/filters/gff_to_bed_converter.py | 98 +++++++++++++++++++++++++-- 3 files changed, 118 insertions(+), 6 deletions(-) diff --git a/lib/galaxy/tools/util/gff_util.py b/lib/galaxy/tools/util/gff_util.py index 85843fa4ebf..a60714f1535 100644 --- a/lib/galaxy/tools/util/gff_util.py +++ b/lib/galaxy/tools/util/gff_util.py @@ -39,4 +39,22 @@ def convert_gff_coords_to_bed( interval ): elif type ( interval ) is list: interval[ 0 ] -= 1 return interval - \ No newline at end of file + +def parse_gff_attributes( attr_str ): + """ + Parses a GFF attribute string and returns a dictionary of name-value pairs. + The general format for a GFF attribute string is name1 "value1" ; name2 "value2" + """ + attributes_list = attr_str.split(";") + attributes = {} + for name_value_pair in attributes_list: + pair = name_value_pair.strip().split(" ") + if pair == '': + continue + name = pair[0].strip() + if name == '': + continue + # Need to strip double quote from values + value = pair[1].strip(" \"") + attributes[ name ] = value + return attributes \ No newline at end of file diff --git a/tools/filters/gff2bed.xml b/tools/filters/gff2bed.xml index 40a80cc0541..5536638f260 100644 --- a/tools/filters/gff2bed.xml +++ b/tools/filters/gff2bed.xml @@ -1,4 +1,4 @@ - + converter gff_to_bed_converter.py $input $out_file1 @@ -12,6 +12,10 @@ + + + + diff --git a/tools/filters/gff_to_bed_converter.py b/tools/filters/gff_to_bed_converter.py index 3851d417845..98e9e01a994 100644 --- a/tools/filters/gff_to_bed_converter.py +++ b/tools/filters/gff_to_bed_converter.py @@ -1,8 +1,50 @@ #!/usr/bin/env python import sys +from galaxy import eggs +from galaxy.tools.util.gff_util import parse_gff_attributes assert sys.version_info[:2] >= ( 2, 4 ) +def get_bed_line( chrom, name, strand, blocks ): + """ Returns a BED line for given data. """ + + + if len( blocks ) == 1: + # Use simple BED format if there is only a single block: + # chrom, chromStart, chromEnd, name, score, strand + # + start, end = blocks[0] + return "%s\t%i\t%i\t%s\t0\t%s\n" % ( chrom, start, end, name, strand ) + + # + # Build lists for transcript blocks' starts, sizes. + # + + # Get transcript start, end. + t_start = sys.maxint + t_end = -1 + for block_start, block_end in blocks: + if block_start < t_start: + t_start = block_start + if block_end > t_end: + t_end = block_end + + # Get block starts, sizes. + block_starts = [] + block_sizes = [] + for block_start, block_end in blocks: + block_starts.append( str( block_start - t_start ) ) + block_sizes.append( str( block_end - block_start ) ) + + # + # Create BED entry. + # Bed format: chrom, chromStart, chromEnd, name, score, strand, \ + # thickStart, thickEnd, itemRgb, blockCount, blockSizes, blockStarts + # + return "%s\t%i\t%i\t%s\t0\t%s\t%i\t%i\t0\t%i\t%s\t%s\n" % \ + ( chrom, t_start, t_end, name, strand, t_start, t_start, len( block_starts ), + ",".join( block_sizes ), ",".join( block_starts ) ) + def __main__(): input_name = sys.argv[1] output_name = sys.argv[2] @@ -10,18 +52,61 @@ def __main__(): first_skipped_line = 0 out = open( output_name, 'w' ) i = 0 + cur_transcript_chrom = None + cur_transcript_id = None + cur_transcript_strand = None + cur_transcripts_blocks = [] # (start, end) for each block. for i, line in enumerate( file( input_name ) ): line = line.rstrip( '\r\n' ) if line and not line.startswith( '#' ): try: + # GFF format: chrom source, name, chromStart, chromEnd, score, strand, attributes elems = line.split( '\t' ) - start = str( int( elems[3] ) - 1 ) + start = str( long( elems[3] ) - 1 ) + coords = [ long( start ), long( elems[4] ) ] strand = elems[6] if strand not in ['+', '-']: strand = '+' - # GFF format: chrom source, name, chromStart, chromEnd, score, strand - # Bed format: chrom, chromStart, chromEnd, name, score, strand - out.write( "%s\t%s\t%s\t%s\t0\t%s\n" %( elems[0], start, elems[4], elems[2], strand ) ) + attributes = parse_gff_attributes( elems[8] ) + t_id = attributes.get( "transcript_id", None ) + + if not t_id: + # + # No transcript ID, so write last transcript and write current line as its own line. + # + + # Write previous transcript. + if cur_transcript_id: + # Write BED entry. + out.write( get_bed_line( cur_transcript_chrome, cur_transcript_id, cur_transcript_strand, cur_transcripts_blocks ) ) + + # Replace any spaces in the name with underscores so UCSC will not complain. + name = elems[2].replace(" ", "_") + out.write( get_bed_line( elems[0], name, strand, [ coords ] ) ) + continue + + # There is a transcript ID, so process line at transcript level. + if t_id == cur_transcript_id: + # Line is element of transcript and will be a block in the BED entry. + cur_transcripts_blocks.append( coords ) + continue + + # + # Line is part of new transcript; write previous transcript and start + # new transcript. + # + + # Write previous transcript. + if cur_transcript_id: + # Write BED entry. + out.write( get_bed_line( cur_transcript_chrome, cur_transcript_id, cur_transcript_strand, cur_transcripts_blocks ) ) + + # Start new transcript. + cur_transcript_chrome = elems[0] + cur_transcript_id = t_id + cur_transcript_strand = strand + cur_transcripts_blocks = [] + cur_transcripts_blocks.append( coords ) except: skipped_lines += 1 if not first_skipped_line: @@ -30,6 +115,11 @@ def __main__(): skipped_lines += 1 if not first_skipped_line: first_skipped_line = i + 1 + + # Write last transcript. + if cur_transcript_id: + # Write BED entry. + out.write( get_bed_line( cur_transcript_chrome, cur_transcript_id, cur_transcript_strand, cur_transcripts_blocks ) ) out.close() info_msg = "%i lines converted to BED. " % ( i + 1 - skipped_lines ) if skipped_lines > 0: