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 tool for profiling locally cached annotations from the UCSC table browser.

Browse Source
This commit is contained in:
Daniel Blankenberg
2008-04-09 17:34:35 +00:00
parent 780e2728c5
commit 7f89a78f35
4 changed files with 270 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tool-data/annotation_profiler.loc.sample
+11
View File
@@ -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)
tool_conf.xml.sample
+8 -7
View File
@@ -93,7 +93,7 @@
<tool file="extract/phastOdds/phastOdds_tool.xml" />
</section>
<section name="Operate on Genomic Intervals" id="bxops">
<tool file="new_operations/intersect.xml" />
<tool file="new_operations/intersect.xml" />
<tool file="new_operations/subtract.xml" />
<tool file="new_operations/merge.xml" />
<tool file="new_operations/concat.xml" />
@@ -104,6 +104,7 @@
<tool file="new_operations/join.xml" />
<tool file="new_operations/get_flanks.xml" />
<tool file="new_operations/flanking_features.xml" />
<tool file="annotation_profiler/annotation_profiler.xml" />
</section>
<section name="Statistics" id="stats">
<tool file="stats/gsummary.xml" />
@@ -120,13 +121,13 @@
<tool file="visualization/build_ucsc_custom_track.xml" />
</section>
<section name="Regional Variation" id="regVar">
<tool file="regVariation/windowSplitter.xml" />
<tool file="regVariation/featureCounter.xml" />
<tool file="regVariation/quality_filter.xml" />
<tool file="regVariation/windowSplitter.xml" />
<tool file="regVariation/featureCounter.xml" />
<tool file="regVariation/quality_filter.xml" />
<tool file="regVariation/maf_cpg_filter.xml" />
<tool file="regVariation/getIndels_2way.xml" />
<tool file="regVariation/getIndels_3way.xml" />
<tool file="regVariation/getIndelRates_3way.xml" />
<tool file="regVariation/getIndels_2way.xml" />
<tool file="regVariation/getIndels_3way.xml" />
<tool file="regVariation/getIndelRates_3way.xml" />
</section>
<section name="Evolution: HyPhy" id="hyphy">
<tool file="hyphy/hyphy_branch_lengths_wrapper.xml" />
tools/annotation_profiler/annotation_profiler.xml
+82
View File
@@ -0,0 +1,82 @@
<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
#if $select_tables.select_table == "some":#-t $select_tables.table_names
#end if
</command>
<inputs>
<param format="interval" name="input1" type="data" label="Choose Intervals">
<validator type="dataset_metadata_in_file" filename="annotation_profiler.loc" metadata_name="dbkey" metadata_column="0" message="Profiling is not currently available for this species."/>
</param>
<param name="keep_empty" type="select" label="Keep Region/Table Pairs with 0 Coverage">
<option value="-k">Keep</option>
<option value="" selected="true">Discard</option>
</param>
<conditional name="select_tables">
<param name="select_table" type="select" label="Limit Tables">
<option value="all" selected="true">Use all tables</option>
<option value="some">Select desired tables</option>
</param>
<when value="all">
</when>
<when value="some">
<param name="table_names" type="select" display="checkboxes" multiple="true" label="Choose Tables to Use" help="Selecting no tables will result in using all tables.">
<options from_file="annotation_profiler.loc" name_col="2" value_col="1">
<filter type="data_meta" data_ref="input1" meta_key="dbkey" meta_key_col="0"/>
</options>
</param>
</when>
</conditional>
</inputs>
<outputs>
<data format="input" name="out_file1"/>
</outputs>
<tests>
<test>
<param name="input1" value="4.bed" dbkey="hg18"/>
<param name="keep_empty" 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>
</tests>
<help>
**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
</help>
</tool>
tools/annotation_profiler/annotation_profiler_for_interval.py
+169
View File
@@ -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__()