diff --git a/tools/annotation_profiler/annotation_profiler.xml b/tools/annotation_profiler/annotation_profiler.xml index 7b918c7f61e..2e8422c837c 100644 --- a/tools/annotation_profiler/annotation_profiler.xml +++ b/tools/annotation_profiler/annotation_profiler.xml @@ -1,6 +1,6 @@ for a set of genomic intervals - 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 + 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 @@ -12,6 +12,10 @@ + + + + @@ -31,15 +35,24 @@ + + - + + + + + + + + **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 + + diff --git a/tools/annotation_profiler/annotation_profiler_for_interval.py b/tools/annotation_profiler/annotation_profiler_for_interval.py index 378ea6b592c..d999464f3c9 100644 --- a/tools/annotation_profiler/annotation_profiler_for_interval.py +++ b/tools/annotation_profiler/annotation_profiler_for_interval.py @@ -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__()