From 7f89a78f35931edc48fe441d1d4a23af7926fefe Mon Sep 17 00:00:00 2001 From: Daniel Blankenberg Date: Wed, 9 Apr 2008 17:34:35 +0000 Subject: [PATCH] Add tool for profiling locally cached annotations from the UCSC table browser. --- tool-data/annotation_profiler.loc.sample | 11 ++ tool_conf.xml.sample | 15 +- .../annotation_profiler.xml | 82 +++++++++ .../annotation_profiler_for_interval.py | 169 ++++++++++++++++++ 4 files changed, 270 insertions(+), 7 deletions(-) create mode 100644 tool-data/annotation_profiler.loc.sample create mode 100644 tools/annotation_profiler/annotation_profiler.xml create mode 100644 tools/annotation_profiler/annotation_profiler_for_interval.py diff --git a/tool-data/annotation_profiler.loc.sample b/tool-data/annotation_profiler.loc.sample new file mode 100644 index 00000000000..510c97de56d --- /dev/null +++ b/tool-data/annotation_profiler.loc.sample @@ -0,0 +1,11 @@ +#build tableName tableDescription +#hg18 acembly acembly (1635047140) +#hg18 acescan acescan (755672) +#hg18 affyGnf1h affyGnf1h (59183600) +#hg18 affyHuEx1 affyHuEx1 (177129194) +#hg18 affyHumanExon affyHumanExon (177129194) +#hg18 affyRatio affyRatio (474485615) +#hg18 affyTxnPhase3FragsHDF affyTxnPhase3FragsHDF (51209878) +#hg18 affyTxnPhase3FragsHeLaBottomStrand affyTxnPhase3FragsHeLaBottomStrand (5245191) +#hg18 affyTxnPhase3FragsHeLaCyto affyTxnPhase3FragsHeLaCyto (67255037) +#hg18 affyTxnPhase3FragsHeLaNuclear affyTxnPhase3FragsHeLaNuclear (87323622) diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample index 1bc90ed044e..4ef0b967658 100644 --- a/tool_conf.xml.sample +++ b/tool_conf.xml.sample @@ -93,7 +93,7 @@
- + @@ -104,6 +104,7 @@ +
@@ -120,13 +121,13 @@
- - - + + + - - - + + +
diff --git a/tools/annotation_profiler/annotation_profiler.xml b/tools/annotation_profiler/annotation_profiler.xml new file mode 100644 index 00000000000..7b918c7f61e --- /dev/null +++ b/tools/annotation_profiler/annotation_profiler.xml @@ -0,0 +1,82 @@ + + 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 +#if $select_tables.select_table == "some":#-t $select_tables.table_names +#end if + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +**What it does** + +Takes an input set of intervals and for each interval determines the base coverage of the interval by a set of features (tables) available from UCSC. + +By default, this tool will check the coverage of your intervals against all available features; you may, however, choose to select only those tables that you want to include. + +In the listing of available tables, the number in parentheses **()** is the total number of bases covered by the feature across all chromosomes. + +----- + +**Example** + +Using the interval below and selecting several tables:: + + chr1 4558 14764 uc001aab.1 0 - + +results in:: + + chr1 4558 14764 uc001aab.1 0 - xenoMrna 10203 + chr1 4558 14764 uc001aab.1 0 - snp126 205 + chr1 4558 14764 uc001aab.1 0 - chainStrPur2Link 1323 + chr1 4558 14764 uc001aab.1 0 - netXenTro2 3050 + chr1 4558 14764 uc001aab.1 0 - intronEst 10206 + chr1 4558 14764 uc001aab.1 0 - snp126orthoPanTro2RheMac2 61 + chr1 4558 14764 uc001aab.1 0 - multiz28way 10206 + chr1 4558 14764 uc001aab.1 0 - chainEquCab1 10206 + chr1 4558 14764 uc001aab.1 0 - affyTxnPhase3HeLaNuclear 9011 + chr1 4558 14764 uc001aab.1 0 - affyHuEx1 3553 + chr1 4558 14764 uc001aab.1 0 - clonePos 10206 + chr1 4558 14764 uc001aab.1 0 - ctgPos 10206 + chr1 4558 14764 uc001aab.1 0 - snp126Exceptions 151 + chr1 4558 14764 uc001aab.1 0 - chainOryLat1 3718 + chr1 4558 14764 uc001aab.1 0 - genomicSuperDups 10206 + chr1 4558 14764 uc001aab.1 0 - netGalGal3 3686 + chr1 4558 14764 uc001aab.1 0 - phastCons28wayPlacMammal 10172 + + + diff --git a/tools/annotation_profiler/annotation_profiler_for_interval.py b/tools/annotation_profiler/annotation_profiler_for_interval.py new file mode 100644 index 00000000000..1a63cbd0f9f --- /dev/null +++ b/tools/annotation_profiler/annotation_profiler_for_interval.py @@ -0,0 +1,169 @@ +#!/usr/bin/env python2.4 +#Dan Blankenberg +#For a set of intervals, this tool returns the same set of intervals +#with 2 additional fields: the name of a Table/Feature and the number of +#bases covered. The original intervals are repeated for each Table/Feature. + +import struct, optparse, os, random +import pkg_resources; pkg_resources.require( "bx-python" ) +import bx.intervals.io +try: + import psyco + psyco.full() +except: + pass + +class CachedRangesInFile: + fmt = 'I' + fmt_size = struct.calcsize( fmt ) + def __init__( self, filename ): + self.file_size = os.stat( filename ).st_size + self.file = open( filename, 'rb' ) + self.length = int( self.file_size / self.fmt_size / 2 ) + self._cached_ranges = [ None for i in xrange( self.length ) ] + def __getitem__( self, i ): + if self._cached_ranges[i] is not None: + return self._cached_ranges[i] + if i < 0: i = self.length + i + offset = i * self.fmt_size * 2 + self.file.seek( offset ) + try: + start = struct.unpack( self.fmt, self.file.read( self.fmt_size ) )[0] + end = struct.unpack( self.fmt, self.file.read( self.fmt_size ) )[0] + except Exception, e: + raise IndexError, e + self._cached_ranges[i] = ( start, end ) + return start, end + def __len__( self ): + return self.length + +class RegionCoverage: + def __init__( self, filename_base ): + try: + self._coverage = CachedRangesInFile( "%s.covered" % filename_base ) + except Exception, e: + #print "Error loading coverage file %s: %s" % ( "%s.covered" % filename_base, e ) + self._coverage = [] + try: + self._total_coverage = int( open( "%s.total_coverage" % filename_base ).read() ) + except Exception, e: + #print "Error loading total coverage file %s: %s" % ( "%s.total_coverage" % filename_base, e ) + self._total_coverage = 0 + def get_start_index( self, start ): + #binary search: returns index of range closest to start + if start > self._coverage[-1][1]: + return len( self._coverage ) - 1 + i = 0 + j = len( self._coverage) - 1 + while i < j: + k = ( i + j ) / 2 + if start <= self._coverage[k][1]: + j = k + else: + i = k + 1 + return i + def get_coverage( self, start, end ): + if len( self._coverage ) < 1 or start > self._coverage[-1][1] or end < self._coverage[0][0]: + return 0 + if self._total_coverage and start <= self._coverage[0][0] and end >= self._coverage[-1][1]: + return self._total_coverage + coverage = 0 + for i in xrange( self.get_start_index( start ), 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 ) + return coverage + +class CachedCoverageReader: + def __init__( self, base_file_path, buffer = 10, table_names = None ): + self._base_file_path = base_file_path + self._buffer = buffer #number of chromosomes to keep in memory at a time + self._coverage = {} + if table_names is None: table_names = os.listdir( self._base_file_path ) + for tablename in table_names: self._coverage[tablename] = {} + def iter_table_coverage_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 ) ) + yield tablename, chromosomes[chrom].get_coverage( start, end ) + + +def __main__(): + parser = optparse.OptionParser() + parser.add_option( + '-k','--keep_empty', + action="store_true", + dest='keep_empty', + default=False, + help='Keep tables with 0 coverage' + ) + parser.add_option( + '-b','--buffer', + dest='buffer', + type='int',default=10, + help='Number of Chromosomes to keep buffered' + ) + parser.add_option( + '-c','--chrom_col', + dest='chrom_col', + type='int',default=1, + help='Chromosome column' + ) + parser.add_option( + '-s','--start_col', + dest='start_col', + type='int',default=2, + help='Start Column' + ) + parser.add_option( + '-e','--end_col', + dest='end_col', + type='int',default=3, + help='End Column' + ) + parser.add_option( + '-p','--path', + dest='path', + type='str',default='/depot/data2/galaxy/annotation_profiler/hg18', + help='Path to profiled data for this organism' + ) + parser.add_option( + '-t','--table_names', + dest='table_names', + type='str',default='None', + help='Path to profiled data for this organism' + ) + parser.add_option( + '-i','--input', + dest='interval_filename', + type='str', + help='Input Interval File' + ) + parser.add_option( + '-o','--output', + dest='out_filename', + type='str', + help='Input Interval File' + ) + + options, args = parser.parse_args() + + table_names = options.table_names.split( "," ) + if "None" in table_names: table_names = None + coverage_reader = CachedCoverageReader( options.path, buffer = options.buffer, table_names = table_names ) + + out = open( options.out_filename, 'wb' ) + + 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__()