From be91b8eed06ecc24e42f4766f6a4a55018906c3b Mon Sep 17 00:00:00 2001 From: Greg Von Kuster Date: Wed, 2 Jan 2008 14:47:30 +0000 Subject: [PATCH] Better error handling for gencode_partition tool, also added functional test for the tool. --- tools/encode/gencode_partition.xml | 13 +- tools/encode/split_by_partitions.py | 206 +++++++++++++++------------- 2 files changed, 117 insertions(+), 102 deletions(-) diff --git a/tools/encode/gencode_partition.xml b/tools/encode/gencode_partition.xml index 040f6fc6b65..1cb80528d7a 100755 --- a/tools/encode/gencode_partition.xml +++ b/tools/encode/gencode_partition.xml @@ -1,5 +1,5 @@ -partition an interval file + an interval file split_by_partitions.py /home/universe/encode_feature_partitions/partition_list.txt $input1 $out_file1 $input1_chromCol $input1_startCol $input1_endCol $input1_strandCol @@ -7,14 +7,20 @@ - + + + + + + + For detailed information about partitioning, click here_. .. _here: http://genome.imim.es/gencode/wiki/index.php/Collecting_Feature_Sets_from_All_Analysis_Groups Datasets are partitioned according to the protocol below: -A partition scheme has been defined that is similar to what has previously been done with TARs/TRANSFRAGs such that any feature can be cla ssified as falling into one of the following 6 categories: +A partition scheme has been defined that is similar to what has previously been done with TARs/TRANSFRAGs such that any feature can be classified as falling into one of the following 6 categories: 1. **Coding** -- coding exons defined from the GENCODE experimentally verified coding set (coding in any transcript) 2. **5UTR** -- 5' UTR exons defined from the GENCODE experimentally verified coding set (5' UTR in some transcript but never coding in any other) 3. **3UTR** -- 3' UTR exons defined from the GENCODE experimentally verified coding set (3' UTR in some transcript but never coding in any other) @@ -29,6 +35,5 @@ A partition scheme has been defined that is similar to what has previously been **Note:** Features overlapping more than one partition will take the identity of the lower-numbered partition. - \ No newline at end of file diff --git a/tools/encode/split_by_partitions.py b/tools/encode/split_by_partitions.py index e11bb6b8330..6c0ecb2ab14 100755 --- a/tools/encode/split_by_partitions.py +++ b/tools/encode/split_by_partitions.py @@ -10,105 +10,115 @@ import pkg_resources; pkg_resources.require( "bx-python" ) from bx.bitset import * from bx.bitset_builders import * +def stop_err( msg ): + sys.stderr.write( msg ) + sys.exit() -partition_index = sys.argv[1] -partition_offset = "/home/universe/encode_feature_partitions/" #should parse perhaps from index filepath - -warnings = [] - -# Load up the partitions -partitions = list() -try: - for line in open( partition_index ): - name, score, filename = line.split() - partitions.append( ( name, score, binned_bitsets_from_file( open( partition_offset+filename ) ) ) ) -except: - print >> sys.stderr, "Error loading partitioning dataset." - sys.exit(0) - -try: - in_file = open(sys.argv[2]) -except: - print >> sys.stderr, "Bad input data." - sys.exit(0) +def main(): + partition_index = sys.argv[1] + partition_offset = "/home/universe/encode_feature_partitions/" #should parse perhaps from index filepath -try: - out_file = open(sys.argv[3], "w") -except: - print >> sys.stderr, "Bad output file." - sys.exit(0) - -try: - chrCol = int (sys.argv[4])-1 -except: - print >> sys.stderr, "Bad chr column." - sys.exit(0) -try: - startCol = int (sys.argv[5])-1 -except: - print >> sys.stderr, "Bad start column." - sys.exit(0) -try: - endCol = int (sys.argv[6])-1 -except: - print >> sys.stderr, "Bad end column." - sys.exit(0) -try: - strandCol = int (sys.argv[7])-1 -except: - print >> sys.stderr, "Bad strand column." - sys.exit(0) - -line_count = 0 -try: - for line in in_file: - line_count+=1 - #ignore comment lines - if line[0:1] == "#": - continue - fields = line.rstrip().split( "\t" ) - try: - chr, start, end = fields[chrCol], int( fields[startCol] ), int( fields[endCol] ) - except: - print >> sys.stderr, "Not enough columns on line ", line_count - continue - label = "input_line_"+str(line_count) #if input file type was known to be bed, then could guess at label column + warnings = [] + + # Load up the partitions + partitions = list() + try: + for line in open( partition_index ): + name, score, filename = line.split() + partitions.append( ( name, score, binned_bitsets_from_file( open( partition_offset+filename ) ) ) ) + except: + stop_err( "Error loading partitioning dataset." ) + + try: + in_file = open( sys.argv[2] ) + except: + stop_err( "Bad input data." ) - if strandCol < 0 : - strand = "+" - else: - try: - strand = fields[strandCol] - except: - strand = "+" - - # Find which partition it overlaps - overlap = 0 - for name, score, bb in partitions: - # Is there at least 1bp overlap? - if chr in bb: - overlap = bb[chr].count_range( start, end-start ) - if overlap > 0: - break - else: - # No overlap with any partition? For now throw this since the - # partitions tile the encode regions completely, indicate an interval - # that does not even overlap an encode region - warning = "warning: Interval (%s, %d, %d) does not overlap any partition" % ( chr, start, end )+", line["+str(line_count)+"]" - warnings.append(warning) - name = "no_overlap" - score = 0 - # Annotate with the name of the partition - frac_overlap = overlap / (end-start) - # BED6 plus? - print >>out_file, "%s\t%d\t%d\t%s\t%s\t%s\t%s\t%0.4f" % ( chr, start, end, label, score, strand, name, frac_overlap ) -except: - print >> sys.stderr, "Unknown error while processing line",line_count -out_file.close() -in_file.close() + try: + out_file = open( sys.argv[3], "w" ) + except: + stop_err( "Bad output file." ) + + try: + chrCol = int( sys.argv[4] ) - 1 + except: + stop_err( "Bad chr column: %s" % ( str( sys.argv[4] ) ) ) + try: + startCol = int( sys.argv[5] ) - 1 + except: + stop_err( "Bad start column: %s" % ( str( sys.argv[5] ) ) ) + try: + endCol = int( sys.argv[6] ) - 1 + except: + stop_err( "Bad end column: %s" % ( str( sys.argv[6] ) ) ) + try: + strandCol = int( sys.argv[7] )-1 + except: + stop_err( "Bad strand column: %s" % ( str( sys.argv[7] ) ) ) + + line_count = 0 + skipped_lines = 0 + first_invalid_line = None + invalid_line = '' + try: + for line in in_file: + line_count += 1 + #ignore comment lines + if line and not line.startswith( '#' ): + fields = line.rstrip( '\r\n' ).split( '\t' ) + try: + chr, start, end = fields[chrCol], int( fields[startCol] ), int( fields[endCol] ) + except: + skipped_lines += 1 + if first_invalid_line is None: + first_invalid_line = line_count + invalid_line = line + continue + label = "input_line_" + str( line_count ) #if input file type was known to be bed, then could guess at label column + + if strandCol < 0: + strand = "+" + else: + try: + strand = fields[strandCol] + except: + strand = "+" + + # Find which partition it overlaps + overlap = 0 + for name, score, bb in partitions: + # Is there at least 1bp overlap? + if chr in bb: + overlap = bb[chr].count_range( start, end-start ) + if overlap > 0: + break + else: + # No overlap with any partition? For now throw this since the + # partitions tile the encode regions completely, indicate an interval + # that does not even overlap an encode region + warning = "warning: Interval (%s, %d, %d) does not overlap any partition" % ( chr, start, end ) + ", line[" + str( line_count ) + "]" + warnings.append( warning ) + name = "no_overlap" + score = 0 + # Annotate with the name of the partition + frac_overlap = overlap / ( end-start ) + # BED6 plus? + print >>out_file, "%s\t%d\t%d\t%s\t%s\t%s\t%s\t%0.4f" % ( chr, start, end, label, score, strand, name, frac_overlap ) + except: + out_file.close() + in_file.close() + stop_err( "Unknown error while processing line # %d: %s" % ( line_count, line ) ) + out_file.close() + in_file.close() -if len(warnings) > 5: - print >> sys.stderr, "There were more than 5 warnings, this tool is useful on ENCODE regions only." -else: - for warning in warnings: - print >> sys.stderr, warning \ No newline at end of file + if warnings: + err_msg = "" + for warning in warnings: + err_msg += warning + "\n" + if len( warnings ) > 5: + err_msg += "There were more than 5 warnings, this tool is useful on ENCODE regions only." + stop_err( err_msg ) + if skipped_lines: + print "Skipped %d invalid lines starting at line # %d: %s" % ( skipped_lines, first_invalid_line, invalid_line ) + +if __name__ == "__main__": main()