From 999a5d4b1928e1ac7ed998f74fa3c964c019e0cf Mon Sep 17 00:00:00 2001 From: Daniel Blankenberg Date: Fri, 23 May 2008 20:16:38 +0000 Subject: [PATCH] Add additional output to annotation profiler. --- .../annotation_profiler.xml | 10 ++-- .../annotation_profiler_for_interval.py | 53 ++++++++++++++----- 2 files changed, 45 insertions(+), 18 deletions(-) diff --git a/tools/annotation_profiler/annotation_profiler.xml b/tools/annotation_profiler/annotation_profiler.xml index a920e01c135..95fb2294ca0 100644 --- a/tools/annotation_profiler/annotation_profiler.xml +++ b/tools/annotation_profiler/annotation_profiler.xml @@ -107,19 +107,21 @@ results in:: Where:: tableName is the name of the table - tableChromosomeCoverage is the number of positions existing in the table (only the chromosomes that were referenced by the interval file are included) - tableChromosomeCount is the number of regions existing in the table (only the chromosomes that were referenced by the interval file are included) + tableChromosomeCoverage is the number of positions existing in the table for only the chromosomes that were referenced by the interval file + tableChromosomeCount is the number of regions existing in the table for only the chromosomes that were referenced by the interval file + tableRegionCoverage is the number of positions existing in the table between the minimal and maximal bounding regions that were referenced by the interval file + tableRegionCount is the number of regions existing in the table between the minimal and maximal bounding regions that were referenced by the interval file allIntervalCount is the number of provided intervals allIntervalSize is the sum of the lengths of the provided interval file allCoverage is the sum of the coverage for each provided interval - allTableRegionsOverlaped is the sum of the number of regions of the table that were overlaped for each interval + allTableRegionsOverlaped is the sum of the number of regions of the table (non-unique) that were overlaped for each interval allIntervalsOverlapingTable is the number of provided intervals which overlap the table nrIntervalCount is the number of non-redundant intervals nrIntervalSize is the sum of the lengths of non-redundant intervals nrCoverage is the sum of the coverage of non-redundant intervals - nrTableRegionsOverlaped is the sum of the number of regions of the table that were overlaped for each non-redundant interval + nrTableRegionsOverlaped is the number of regions of the table (unique) that were overlaped by the non-redundant intervals nrIntervalsOverlapingTable is the number of non-redundant intervals which overlap the table diff --git a/tools/annotation_profiler/annotation_profiler_for_interval.py b/tools/annotation_profiler/annotation_profiler_for_interval.py index a67d191d412..f11edb55c6f 100644 --- a/tools/annotation_profiler/annotation_profiler_for_interval.py +++ b/tools/annotation_profiler/annotation_profiler_for_interval.py @@ -69,20 +69,23 @@ class RegionCoverage: def get_coverage( self, start, end ): return self.get_coverage_regions_overlap( start, end )[0] def get_coverage_regions_overlap( self, start, end ): + return self.get_coverage_regions_index_overlap( start, end )[0:2] + def get_coverage_regions_index_overlap( self, start, end ): if len( self._coverage ) < 1 or start > self._coverage[-1][1] or end < self._coverage[0][0]: - return 0, 0 + return 0, 0, 0 if self._total_coverage and start <= self._coverage[0][0] and end >= self._coverage[-1][1]: - return self._total_coverage, len( self._coverage ) + return self._total_coverage, len( self._coverage ), 0 coverage = 0 region_count = 0 - for i in xrange( self.get_start_index( start ), len( self._coverage ) ): + start_index = self.get_start_index( start ) + for i in xrange( start_index, len( self._coverage ) ): c_start, c_end = self._coverage[i] if c_start > end: break if c_start <= end and c_end >= start: coverage += min( end, c_end ) - max( start, c_start ) region_count += 1 - return coverage, region_count + return coverage, region_count, start_index class CachedCoverageReader: def __init__( self, base_file_path, buffer = 10, table_names = None ): @@ -95,14 +98,17 @@ class CachedCoverageReader: for tablename, coverage, regions in self.iter_table_coverage_regions_by_region( chrom, start, end ): yield tablename, coverage def iter_table_coverage_regions_by_region( self, chrom, start, end ): + for tablename, coverage, regions, index in self.iter_table_coverage_regions_index_by_region( chrom, start, end ): + yield tablename, coverage, regions + def iter_table_coverage_regions_index_by_region( self, chrom, start, end ): for tablename, chromosomes in self._coverage.iteritems(): if chrom not in chromosomes: if len( chromosomes ) >= self._buffer: #randomly remove one chromosome from this table del chromosomes[ chromosomes.keys().pop( random.randint( 0, self._buffer - 1 ) ) ] chromosomes[chrom] = RegionCoverage( os.path.join ( self._base_file_path, tablename, chrom ) ) - coverage, regions = chromosomes[chrom].get_coverage_regions_overlap( start, end ) - yield tablename, coverage, regions + coverage, regions, index = chromosomes[chrom].get_coverage_regions_index_overlap( start, end ) + yield tablename, coverage, regions, index class TableCoverageSummary: def __init__( self, coverage_reader ): @@ -144,14 +150,21 @@ class TableCoverageSummary: interval_table_overlap_count = {} table_regions_overlap_count = {} interval_count = 0 + region_start_end = {} for chrom, chromosome_bitset in self.chromosome_coverage.iteritems(): end = 0 + last_end_index = {} while True: start = chromosome_bitset.next_set( end ) if start >= chromosome_bitset.size: break end = chromosome_bitset.next_clear( start ) interval_count += 1 - for table_name, coverage, region_count in self.coverage_reader.iter_table_coverage_regions_by_region( chrom, start, end ): + if chrom not in region_start_end: + region_start_end[chrom] = [start, end] + else: + if start < region_start_end[chrom][0]: region_start_end[chrom][0] = start + if end > region_start_end[chrom][1]: region_start_end[chrom][1] = end + for table_name, coverage, region_count, start_index in self.coverage_reader.iter_table_coverage_regions_index_by_region( chrom, start, end ): if table_name not in table_coverage: table_coverage[table_name] = 0 interval_table_overlap_count[table_name] = 0 @@ -159,8 +172,20 @@ class TableCoverageSummary: table_coverage[table_name] += coverage if coverage: interval_table_overlap_count[table_name] += 1 - table_regions_overlap_count[table_name] += region_count - return interval_count, table_coverage, table_regions_overlap_count, interval_table_overlap_count + table_regions_overlap_count[table_name] += region_count + if table_name in last_end_index and last_end_index[table_name] == start_index: + table_regions_overlap_count[table_name] -= 1 + last_end_index[table_name] = start_index + region_count - 1 + table_region_coverage = {} + table_region_count = {} + for chrom, start_end in region_start_end.items(): + for table_name, coverage, region_count in self.coverage_reader.iter_table_coverage_regions_by_region( chrom, start_end[0], start_end[1] ): + if table_name not in table_region_coverage: + table_region_coverage[table_name] = 0 + table_region_count[table_name] = 0 + table_region_coverage[table_name] += coverage + table_region_count[table_name] += region_count + return table_region_coverage, table_region_count, interval_count, table_coverage, table_regions_overlap_count, interval_table_overlap_count def get_nr_region_size( self ): if self._nr_region_size is None: self._nr_region_size = 0 @@ -168,10 +193,10 @@ class TableCoverageSummary: self._nr_region_size += chromosome_bitset.count_range() return self._nr_region_size def iter_table_coverage( self ): - nr_interval_count, nr_table_coverage, nr_table_regions_overlap_count, nr_interval_table_overlap_count = self.get_nr_coverage() + table_region_coverage, table_region_count, nr_interval_count, nr_table_coverage, nr_table_regions_overlap_count, nr_interval_table_overlap_count = 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, sum( self.table_chromosome_size.get( table_name, [] ).values() ), sum( self.table_chromosome_count.get( table_name, [] ).values() ), self.total_interval_count, self.total_interval_size, self.table_coverage[table_name], self.table_regions_overlaped_count.get( table_name, 0), self.interval_table_overlap_count.get( table_name, 0 ), nr_interval_count, self.get_nr_region_size(), nr_table_coverage[table_name], nr_table_regions_overlap_count.get( table_name, 0 ), nr_interval_table_overlap_count.get( table_name, 0 ) + yield table_name, sum( self.table_chromosome_size.get( table_name, [] ).values() ), sum( self.table_chromosome_count.get( table_name, [] ).values() ), table_region_coverage.get( table_name, 0 ), table_region_count.get( table_name, 0 ), self.total_interval_count, self.total_interval_size, self.table_coverage[table_name], self.table_regions_overlaped_count.get( table_name, 0), self.interval_table_overlap_count.get( table_name, 0 ), nr_interval_count, self.get_nr_region_size(), nr_table_coverage[table_name], nr_table_regions_overlap_count.get( table_name, 0 ), nr_interval_table_overlap_count.get( table_name, 0 ) def profile_per_interval( interval_filename, chrom_col, start_col, end_col, out_filename, keep_empty, coverage_reader ): out = open( out_filename, 'wb' ) @@ -184,15 +209,15 @@ def profile_per_interval( interval_filename, chrom_col, start_col, end_col, out_ 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\ttableChromosomeCoverage\ttableChromosomeCount\tallIntervalCount\tallIntervalSize\tallCoverage\tallTableRegionsOverlaped\tallIntervalsOverlapingTable\tnrIntervalCount\tnrIntervalSize\tnrCoverage\tnrTableRegionsOverlaped\tnrIntervalsOverlapingTable\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_chromosome_size, table_chromosome_count, total_interval_count, total_interval_size, total_coverage, table_regions_overlaped_count, interval_region_overlap_count, nr_interval_count, nr_region_size, nr_coverage, nr_table_regions_overlaped_count, nr_interval_table_overlap_count in table_coverage_summary.iter_table_coverage(): + out.write( "#tableName\ttableChromosomeCoverage\ttableChromosomeCount\ttableRegionCoverage\ttableRegionCount\tallIntervalCount\tallIntervalSize\tallCoverage\tallTableRegionsOverlaped\tallIntervalsOverlapingTable\tnrIntervalCount\tnrIntervalSize\tnrCoverage\tnrTableRegionsOverlaped\tnrIntervalsOverlapingTable\n" )#\tstatistic\n" ) + for table_name, table_chromosome_size, table_chromosome_count, table_region_coverage, table_region_count, total_interval_count, total_interval_size, total_coverage, table_regions_overlaped_count, interval_region_overlap_count, nr_interval_count, nr_region_size, nr_coverage, nr_table_regions_overlaped_count, nr_interval_table_overlap_count 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\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n" % ( table_name, table_chromosome_size, table_chromosome_count, total_interval_count, total_interval_size, total_coverage, table_regions_overlaped_count, interval_region_overlap_count, nr_interval_count, nr_region_size, nr_coverage, nr_table_regions_overlaped_count, nr_interval_table_overlap_count ) ) + out.write( "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n" % ( table_name, table_chromosome_size, table_chromosome_count, table_region_coverage, table_region_count, total_interval_count, total_interval_size, total_coverage, table_regions_overlaped_count, interval_region_overlap_count, nr_interval_count, nr_region_size, nr_coverage, nr_table_regions_overlaped_count, nr_interval_table_overlap_count ) ) out.close() def __main__():