First pass at adding Summary option to Annotation Profiler.

2026-09-24 16:30:27 +08:00 · 2008-05-01 19:38:49 +00:00
parent 4d520da9f8
commit 44c9ccf4f4
2 changed files with 142 additions and 10 deletions
@@ -1,6 +1,6 @@
 <tool id="Annotation_Profiler_0" name="Profile Annotations" Version="1.0.0">
  <description>for a set of genomic intervals</description>
-  <command interpreter="python2.4">annotation_profiler_for_interval.py -i $input1 -c $input1_chromCol -s $input1_startCol -e $input1_endCol -o $out_file1 $keep_empty -p /depot/data2/galaxy/annotation_profiler/$dbkey -b 3
+  <command interpreter="python2.4">annotation_profiler_for_interval.py -i $input1 -c $input1_chromCol -s $input1_startCol -e $input1_endCol -o $out_file1 $keep_empty -p /depot/data2/galaxy/annotation_profiler/$dbkey $summary -b 3
 #if $select_tables.select_table == "some":#-t $select_tables.table_names
 #end if
  </command>
@@ -12,6 +12,10 @@
      <option value="-k">Keep</option>
      <option value="" selected="true">Discard</option>
    </param>
+    <param name="summary" type="select" label="Output per Region/Summary">
+      <option value="-S">Summary</option>
+      <option value="" selected="true">Per Region</option>
+    </param>
    <conditional name="select_tables">
      <param name="select_table" type="select" label="Limit Tables">
        <option value="all" selected="true">Use all tables</option>
@@ -31,15 +35,24 @@
   <outputs>
     <data format="input" name="out_file1"/>
   </outputs>
+   <code file="annotation_profiler_code.py" />  
   <tests>
     <test>
       <param name="input1" value="4.bed" dbkey="hg18"/>
       <param name="keep_empty" value=""/>
+       <param name="summary" value=""/>
       <param name="select_table" value="some"/>
       <param name="table_names" value="acembly,affyGnf1h,affyHuEx1,knownAlt,knownGene,mrna,multiz17way,multiz28way,refFlat,refGene,snp126"/>
-       <param name="keep_empty" value=""/>
       <output name="out_file1" file="annotation_profiler_1.out" />
     </test>
+     <test>
+       <param name="input1" value="4.bed" dbkey="hg18"/>
+       <param name="keep_empty" value=""/>
+       <param name="summary" value="-S"/>
+       <param name="select_table" value="some"/>
+       <param name="table_names" value="acembly,affyGnf1h,affyHuEx1,knownAlt,knownGene,mrna,multiz17way,multiz28way,refFlat,refGene,snp126"/>
+       <output name="out_file1" file="annotation_profiler_2.out" />
+     </test>
   </tests>
   <help>
 **What it does**
@@ -50,6 +63,8 @@ By default, this tool will check the coverage of your intervals against all avai

 In the listing of available tables, the number in parentheses **()** is the total number of bases covered by the feature across all chromosomes.

+You may alternatively choose to recieve a summary across all of the intervals that you provide.
+
 -----

 **Example** 
@@ -78,5 +93,42 @@ results in::
 chr1 4558 14764 uc001aab.1 0 - netGalGal3 3686
 chr1 4558 14764 uc001aab.1 0 - phastCons28wayPlacMammal 10172

+Alternatively, requesting a summary, using the intervals below and selecting several tables::
+
+ chr1 4558 14764 uc001aab.1 0 -
+ chr1 4558 19346 uc001aac.1 0 -
+
+results in::
+
+ #tableName tableSize totalRegionSize totalCoverage nrRegionSize nrCoverage
+ snp126Exceptions 133601 24994 388 14788 237
+ genomicSuperDups 12268847 24994 24994 14788 14788
+ chainOryLat1 70337730 24994 7436 14788 3718
+ affyHuEx1 15703901 24994 7846 14788 4293
+ multiz28way 225928588 24994 24994 14788 14788
+ intronEst 135796064 24994 24994 14788 14788
+ xenoMrna 129031327 24994 20406 14788 10203
+ ctgPos 224999719 24994 24994 14788 14788
+ netXenTro2 111440392 24994 6100 14788 3050
+ clonePos 224999719 24994 24994 14788 14788
+ chainStrPur2Link 7948016 24994 2646 14788 1323
+ affyTxnPhase3HeLaNuclear 136797870 24994 22601 14788 13590
+ snp126orthoPanTro2RheMac2 700436 24994 124 14788 63
+ snp126 956976 24994 498 14788 293
+ chainEquCab1 246306414 24994 24994 14788 14788
+ netGalGal3 203351973 24994 7372 14788 3686
+ phastCons28wayPlacMammal 221017670 24994 24926 14788 14754
+
+Where::
+
+ tableSize is the number of positions existing in the table for only the chromosomes that were referenced by the interval file.
+ totalRegionSize is the sum of the lengths of the provided interval file.
+ totalCoverage is the sum of the coverage for each interval
+ nrRegionSize is the sum of the lengths of non-redundant intervals
+ nrCoverage is the sum of the coverage of non-redundant intervals
+   
+   where non-redundant indicates that input intervals have been collapsed to resolve overlaps
+
+
  </help>
 </tool>
@@ -8,6 +8,7 @@ import sys, struct, optparse, os, random
 from galaxy import eggs
 import pkg_resources; pkg_resources.require( "bx-python" )
 import bx.intervals.io
+import bx.bitset
 try:
    import psyco
    psyco.full()
@@ -95,6 +96,81 @@ class CachedCoverageReader:
                chromosomes[chrom] = RegionCoverage( os.path.join ( self._base_file_path, tablename, chrom ) )
            yield tablename, chromosomes[chrom].get_coverage( start, end )

+class TableCoverageSummary:
+    def __init__( self, coverage_reader ):
+        self.coverage_reader = coverage_reader
+        self.chromosome_coverage = {}
+        self.total_region_size = 0
+        self.table_coverage = {}
+        self.table_size = {}
+        self._nr_region_size = None
+    def add_region( self, chrom, start, end ):
+        self.total_region_size += ( end - start )
+        if chrom not in self.chromosome_coverage:
+            #utilize lengths file here, if possible, if not use 250mb
+            #currently, no valid method to provide location of lengths file by framework:
+            #gops_complement has it hard coded as dbfile = fileinput.FileInput( "static/ucsc/chrom/"+db+".len" )
+            self.chromosome_coverage[chrom] = bx.bitset.BitSet( 250000000 )
+        self.chromosome_coverage[chrom].set_range( start, end - start )
+        for table_name, coverage in self.coverage_reader.iter_table_coverage_by_region( chrom, start, end ):
+            if table_name not in self.table_coverage:
+                self.table_coverage[table_name] = 0
+                self.table_size[table_name] = {}
+            if chrom not in self.table_size[table_name]:
+                self.table_size[table_name][chrom] = self.coverage_reader._coverage[table_name][chrom]._total_coverage
+            self.table_coverage[table_name] += coverage
+    def get_table_size( self, table_name ):
+        if table_name not in self.table_size: return 0
+        size = 0
+        for chrom, chrom_size in self.table_size[table_name].iteritems():
+            size += chrom_size
+        return size
+    def get_nr_coverage( self ):
+        table_coverage = {}
+        for chrom, chromosome_bitset in self.chromosome_coverage.iteritems():
+            end = 0
+            while True:
+                start = chromosome_bitset.next_set( end )
+                if start >= chromosome_bitset.size: break
+                end = chromosome_bitset.next_clear( start )
+                for table_name, coverage in self.coverage_reader.iter_table_coverage_by_region( chrom, start, end ):
+                    if table_name not in table_coverage:
+                        table_coverage[table_name] = 0
+                    table_coverage[table_name] += coverage
+        return table_coverage
+    def get_nr_region_size( self ):
+        if self._nr_region_size is None:
+            self._nr_region_size = 0
+            for chrom, chromosome_bitset in self.chromosome_coverage.iteritems():
+                self._nr_region_size += chromosome_bitset.count_range()
+        return self._nr_region_size
+    def iter_table_coverage( self ):
+        nr_table_coverage = self.get_nr_coverage()
+        for table_name in self.table_coverage:
+            #TODO: determine a type of statistic, then calculate and report here
+            yield table_name, self.get_table_size( table_name ), self.total_region_size, self.table_coverage[table_name], self.get_nr_region_size(), nr_table_coverage[table_name]
+
+def profile_per_interval( interval_filename, chrom_col, start_col, end_col, out_filename, keep_empty, coverage_reader ):
+    out = open( out_filename, 'wb' )
+    for region in bx.intervals.io.NiceReaderWrapper( open( interval_filename, 'rb' ), chrom_col = chrom_col, start_col = start_col, end_col = end_col, fix_strand = True, return_header = False, return_comments = False ):
+        for table_name, coverage in coverage_reader.iter_table_coverage_by_region( region.chrom, region.start, region.end ):
+            if keep_empty or coverage:
+                #only output regions that have atleast 1 base covered unless empty are requested
+                out.write( "%s\t%s\t%s\n" % ( "\t".join( region.fields ), table_name, coverage ) )
+    out.close()
+
+def profile_summary( interval_filename, chrom_col, start_col, end_col, out_filename, keep_empty, coverage_reader ):
+    out = open( out_filename, 'wb' )
+    out.write( "#tableName\ttableSize\ttotalRegionSize\ttotalCoverage\tnrRegionSize\tnrCoverage\n" )#\tstatistic\n" )
+    table_coverage_summary = TableCoverageSummary( coverage_reader )
+    for region in bx.intervals.io.NiceReaderWrapper( open( interval_filename, 'rb' ), chrom_col = chrom_col, start_col = start_col, end_col = end_col, fix_strand = True, return_header = False, return_comments = False ):
+        table_coverage_summary.add_region( region.chrom, region.start, region.end )
+    
+    for table_name, table_size, total_region_size, total_coverage, nr_region_size, nr_coverage in table_coverage_summary.iter_table_coverage():
+        if keep_empty or total_coverage:
+            #only output tables that have atleast 1 base covered unless empty are requested
+            out.write( "%s\t%s\t%s\t%s\t%s\t%s\n" % ( table_name, table_size, total_region_size, total_coverage, nr_region_size, nr_coverage ) )
+    out.close()

 def __main__():
    parser = optparse.OptionParser()
@@ -153,20 +229,24 @@ def __main__():
        type='str',
        help='Input Interval File'
    )
+    parser.add_option(
+        '-S','--summary',
+        action="store_true",
+        dest='summary',
+        default=False,
+        help='Display Summary Results'
+    )
    
    options, args = parser.parse_args()
    
    table_names = options.table_names.split( "," )
-    if "None" in table_names: table_names = None
+    if table_names == ['None']: table_names = None
    coverage_reader = CachedCoverageReader( options.path, buffer = options.buffer, table_names = table_names )
    
-    out = open( options.out_filename, 'wb' )
+    if options.summary:
+        profile_summary( options.interval_filename, options.chrom_col - 1, options.start_col - 1, options.end_col -1, options.out_filename, options.keep_empty, coverage_reader )
+    else:
+        profile_per_interval( options.interval_filename, options.chrom_col - 1, options.start_col - 1, options.end_col -1, options.out_filename, options.keep_empty, coverage_reader )
    
-    for region in bx.intervals.io.NiceReaderWrapper( open( options.interval_filename, 'rb' ), chrom_col = options.chrom_col - 1, start_col = options.start_col - 1, end_col = options.end_col -1 , fix_strand = True, return_header = False, return_comments = False ):
-        for tablename, coverage in coverage_reader.iter_table_coverage_by_region( region.chrom, region.start, region.end ):
-            if options.keep_empty or coverage:
-                #only output regions that have atleast 1 base covered unless empty are requested
-                out.write("%s\t%s\t%s\n" % ( "\t".join( region.fields ), tablename, coverage ) )
-    out.close()

 if __name__ == "__main__": __main__()