Better error handling for gencode_partition tool, also added functional test for the tool.

2026-09-24 16:30:27 +08:00 · 2008-01-02 14:47:30 +00:00
parent 6c7570a8cb
commit be91b8eed0
2 changed files with 117 additions and 102 deletions
@@ -1,5 +1,5 @@
 <tool id="gencode_partition1" name="Gencode Partition">
-<description>partition an interval file</description>
+  <description>an interval file</description>
  <command interpreter="python2.4">split_by_partitions.py /home/universe/encode_feature_partitions/partition_list.txt $input1 $out_file1 $input1_chromCol $input1_startCol $input1_endCol $input1_strandCol</command>
  <inputs>
    <param name="input1" type="data" format="interval" label="File to Partition"/>
@@ -7,14 +7,20 @@
  <outputs>
    <data name="out_file1" format="bed"/>
  </outputs>
-<help>
+  <tests>
+    <test>
+      <param name="input1" value="encode_1.bed"/>
+      <output name="out_file1" file="gencode_partition_out.bed"/>
+    </test>
+  </tests>
+  <help>
 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. 

-
  </help>
 </tool>
@@ -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
+    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()