opensource/galaxy
3
0
Fork 0
mirror of https://github.com/galaxyproject/galaxy.git synced 2026-09-24 16:30:27 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add additional output to annotation profiler.

Browse Source
This commit is contained in:
Daniel Blankenberg
2008-05-23 20:16:38 +00:00
parent 36cec505b1
commit 999a5d4b19
2 changed files with 45 additions and 18 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tools/annotation_profiler/annotation_profiler.xml
+6 -4
View File
@@ -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
tools/annotation_profiler/annotation_profiler_for_interval.py
+39 -14
View File
@@ -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__():