Restored version 1.1.0 of aggregate scores tool and moved the fixes to version 1.1.1.

This commit is contained in:
Greg Von Kuster
2008-02-28 00:08:05 +00:00
parent 72bfcf8281
commit 335ebef0af
6 changed files with 315 additions and 136 deletions
+1 -66
View File
@@ -65,11 +65,6 @@
</section>
-->
<section name="Fetch Sequences" id="fetchSeq">
<!--
These are deprecated
<tool file="extract/Extract_genomic_DNA_1/1.0.0/fasta-subseq-wrapper.xml" />
<tool file="extract/Extract_genomic_DNA_1/1.0.0/twoBitToFa_wrapper.xml" />
-->
<tool file="extract/Extract_genomic_DNA_1/2.0.0/extract_genomic_dna.xml" />
</section>
<section name="Fetch Alignments" id="fetchAlign">
@@ -89,7 +84,7 @@
</section>
<section name="Get Genomic Scores" id="scores">
<tool file="stats/wiggle2simple1/1.0.0/wiggle_to_simple.xml" />
<tool file="stats/aggregate_scores_in_intervals2/1.1.0/aggregate_binned_scores_in_intervals.xml" />
<tool file="stats/aggregate_scores_in_intervals2/1.1.1/aggregate_binned_scores_in_intervals.xml" />
<tool file="extract/phastOdds/phastOdds_for_intervals/1.0.0/phastOdds_tool.xml" />
</section>
<section name="Operate on Genomic Intervals" id="bxops">
@@ -131,66 +126,6 @@
<tool file="hyphy/hyphy_nj_tree_wrapper1/1.0.0/hyphy_nj_tree_wrapper.xml" />
<tool file="hyphy/hyphy_dnds_wrapper1/1.0.0/hyphy_dnds_wrapper.xml" />
</section>
<!--
These are deprecated
<section name="EMBOSS" id="EMBOSSLite">
<tool file="emboss/emboss_cai.xml" />
<tool file="emboss/emboss_cai_custom.xml" />
<tool file="emboss/emboss_codcmp.xml" />
<tool file="emboss/emboss_compseq.xml" />
<tool file="emboss/emboss_cpgplot.xml" />
<tool file="emboss/emboss_cpgreport.xml" />
<tool file="emboss/emboss_cusp.xml" />
<tool file="emboss/emboss_cutseq.xml" />
<tool file="emboss/emboss_dan.xml" />
<tool file="emboss/emboss_einverted.xml" />
<tool file="emboss/emboss_equicktandem.xml" />
<tool file="emboss/emboss_est2genome.xml" />
<tool file="emboss/emboss_etandem.xml" />
<tool file="emboss/emboss_freak.xml" />
<tool file="emboss/emboss_fuzznuc.xml" />
<tool file="emboss/emboss_fuzztran.xml" />
<tool file="emboss/emboss_getorf.xml" />
<tool file="emboss/emboss_isochore.xml" />
<tool file="emboss/emboss_msbar.xml" />
<tool file="emboss/emboss_needle.xml" />
<tool file="emboss/emboss_newcpgreport.xml" />
<tool file="emboss/emboss_newcpgseek.xml" />
<tool file="emboss/emboss_newseq.xml" />
<tool file="emboss/emboss_notseq.xml" />
<tool file="emboss/emboss_nthseq.xml" />
<tool file="emboss/emboss_palindrome.xml" />
<tool file="emboss/emboss_pasteseq.xml" />
<tool file="emboss/emboss_plotorf.xml" />
<tool file="emboss/emboss_polydot.xml" />
<tool file="emboss/emboss_prettyseq.xml" />
<tool file="emboss/emboss_primersearch.xml" />
<tool file="emboss/emboss_revseq.xml" />
<tool file="emboss/emboss_seqmatchall.xml" />
<tool file="emboss/emboss_showorf.xml" />
<tool file="emboss/emboss_shuffleseq.xml" />
<tool file="emboss/emboss_sirna.xml" />
<tool file="emboss/emboss_sixpack.xml" />
<tool file="emboss/emboss_splitter.xml" />
<tool file="emboss/emboss_stretcher.xml" />
<tool file="emboss/emboss_supermatcher.xml" />
<tool file="emboss/emboss_syco.xml" />
<tool file="emboss/emboss_tranalign.xml" />
<tool file="emboss/emboss_transeq.xml" />
<tool file="emboss/emboss_trimest.xml" />
<tool file="emboss/emboss_trimseq.xml" />
<tool file="emboss/emboss_union.xml" />
<tool file="emboss/emboss_vectorstrip.xml" />
<tool file="emboss/emboss_water.xml" />
<tool file="emboss/emboss_wobble.xml" />
<tool file="emboss/emboss_wordcount.xml" />
<tool file="emboss/emboss_wordmatch.xml" />
<tool file="emboss/emboss_backtranseq.xml" />
<tool file="emboss/emboss_biosed.xml" />
<tool file="emboss/emboss_charge.xml" />
<tool file="emboss/emboss_checktrans.xml" />
</section>
-->
<section name="EMBOSS" id="EMBOSSLite">
<tool file="emboss_5/emboss_antigenic.xml" />
<tool file="emboss_5/emboss_backtranseq.xml" />
+1 -66
View File
@@ -72,11 +72,6 @@
<tool file="patmat/findcluster_mysql.xml" />
</section>
<section name="Fetch Sequences" id="fetchSeq">
<!--
These are deprecated
<tool file="extract/Extract_genomic_DNA_1/1.0.0/fasta-subseq-wrapper.xml" />
<tool file="extract/Extract_genomic_DNA_1/1.0.0/twoBitToFa_wrapper.xml" />
-->
<tool file="extract/Extract_genomic_DNA_1/2.0.0/extract_genomic_dna.xml" />
</section>
<section name="Fetch Alignments" id="fetchAlign">
@@ -94,7 +89,7 @@
</section>
<section name="Get Genomic Scores" id="scores">
<tool file="stats/wiggle2simple1/1.0.0/wiggle_to_simple.xml" />
<tool file="stats/aggregate_scores_in_intervals2/1.1.0/aggregate_binned_scores_in_intervals.xml" />
<tool file="stats/aggregate_scores_in_intervals2/1.1.1/aggregate_binned_scores_in_intervals.xml" />
<tool file="extract/phastOdds/phastOdds_for_intervals/1.0.0/phastOdds_tool.xml" />
</section>
<section name="Operate on Genomic Intervals" id="bxops">
@@ -135,66 +130,6 @@
<tool file="hyphy/hyphy_nj_tree_wrapper1/1.0.0/hyphy_nj_tree_wrapper.xml" />
<tool file="hyphy/hyphy_dnds_wrapper1/1.0.0/hyphy_dnds_wrapper.xml" />
</section>
<!--
These are deprecated
<section name="EMBOSS" id="EMBOSSLite">
<tool file="emboss/emboss_cai.xml" />
<tool file="emboss/emboss_cai_custom.xml" />
<tool file="emboss/emboss_codcmp.xml" />
<tool file="emboss/emboss_compseq.xml" />
<tool file="emboss/emboss_cpgplot.xml" />
<tool file="emboss/emboss_cpgreport.xml" />
<tool file="emboss/emboss_cusp.xml" />
<tool file="emboss/emboss_cutseq.xml" />
<tool file="emboss/emboss_dan.xml" />
<tool file="emboss/emboss_einverted.xml" />
<tool file="emboss/emboss_equicktandem.xml" />
<tool file="emboss/emboss_est2genome.xml" />
<tool file="emboss/emboss_etandem.xml" />
<tool file="emboss/emboss_freak.xml" />
<tool file="emboss/emboss_fuzznuc.xml" />
<tool file="emboss/emboss_fuzztran.xml" />
<tool file="emboss/emboss_getorf.xml" />
<tool file="emboss/emboss_isochore.xml" />
<tool file="emboss/emboss_msbar.xml" />
<tool file="emboss/emboss_needle.xml" />
<tool file="emboss/emboss_newcpgreport.xml" />
<tool file="emboss/emboss_newcpgseek.xml" />
<tool file="emboss/emboss_newseq.xml" />
<tool file="emboss/emboss_notseq.xml" />
<tool file="emboss/emboss_nthseq.xml" />
<tool file="emboss/emboss_palindrome.xml" />
<tool file="emboss/emboss_pasteseq.xml" />
<tool file="emboss/emboss_plotorf.xml" />
<tool file="emboss/emboss_polydot.xml" />
<tool file="emboss/emboss_prettyseq.xml" />
<tool file="emboss/emboss_primersearch.xml" />
<tool file="emboss/emboss_revseq.xml" />
<tool file="emboss/emboss_seqmatchall.xml" />
<tool file="emboss/emboss_showorf.xml" />
<tool file="emboss/emboss_shuffleseq.xml" />
<tool file="emboss/emboss_sirna.xml" />
<tool file="emboss/emboss_sixpack.xml" />
<tool file="emboss/emboss_splitter.xml" />
<tool file="emboss/emboss_stretcher.xml" />
<tool file="emboss/emboss_supermatcher.xml" />
<tool file="emboss/emboss_syco.xml" />
<tool file="emboss/emboss_tranalign.xml" />
<tool file="emboss/emboss_transeq.xml" />
<tool file="emboss/emboss_trimest.xml" />
<tool file="emboss/emboss_trimseq.xml" />
<tool file="emboss/emboss_union.xml" />
<tool file="emboss/emboss_vectorstrip.xml" />
<tool file="emboss/emboss_water.xml" />
<tool file="emboss/emboss_wobble.xml" />
<tool file="emboss/emboss_wordcount.xml" />
<tool file="emboss/emboss_wordmatch.xml" />
<tool file="emboss/emboss_backtranseq.xml" />
<tool file="emboss/emboss_biosed.xml" />
<tool file="emboss/emboss_charge.xml" />
<tool file="emboss/emboss_checktrans.xml" />
</section>
-->
<!--
TODO: uncomment the following EMBOSS section whenever
moving to test, but comment it in .sample to eliminate
@@ -6,9 +6,7 @@
#end if
</command>
<inputs>
<param format="interval" name="input1" type="data" label="Interval file">
<validator type="unspecified_build" message="Unspecified build, this tool works with data from genome builds hg16, hg17 or hg18. Click the pencil icon in your history item to set the genome build."/>
</param>
<param format="interval" name="input1" type="data" label="Interval file"/>
<conditional name="score_source_type">
<param name="score_source" type="select" label="Score Source">
<option value="cached" selected="true">Locally Cached Scores</option>
@@ -62,7 +62,7 @@ def load_scores_wiggle( fname ):
for chrom, pos, val in bx.wiggle.Reader( UCSCOutWrapper(open( fname ) ) ):
if chrom not in scores_by_chrom:
scores_by_chrom[chrom] = BinnedArray()
scores_by_chrom[chrom][pos] = val
scores_by_chrom[chrom][pos] = val
except UCSCLimitException:
# Wiggle data was truncated, at the very least need to warn the user.
print 'Encountered message from UCSC: "Reached output limit of 100000 data values", so be aware your data was truncated.'
@@ -0,0 +1,115 @@
<tool id="aggregate_scores_in_intervals2" description="such as phastCons, GERP, binCons, and others for a set of genomic intervals" name="Aggregate datapoints" version="1.1.1">
<description>Appends the average, min, max of datapoints per interval</description>
<command interpreter="python2.4">
#if $score_source_type.score_source == "user":#aggregate_scores_in_intervals.py $score_source_type.input2 $input1 $input1_chromCol $input1_startCol $input1_endCol $out_file1
#else:#aggregate_scores_in_intervals.py $score_source_type.datasets $input1 $input1_chromCol $input1_startCol $input1_endCol $out_file1 -b
#end if
</command>
<inputs>
<param format="interval" name="input1" type="data" label="Interval file">
<validator type="unspecified_build" message="Unspecified build, this tool works with data from genome builds hg16, hg17 or hg18. Click the pencil icon in your history item to set the genome build."/>
</param>
<conditional name="score_source_type">
<param name="score_source" type="select" label="Score Source">
<option value="cached" selected="true">Locally Cached Scores</option>
<option value="user">Scores in Your History</option>
</param>
<when value="cached">
<param name="datasets" type="select" label="Available datasets" display="radio">
<options from_file="/depot/data2/galaxy/binned_scores.loc" >
<filter type="data_meta" data_ref="input1" key="build" />
<filter type="column" name="build_col" value="0" />
<filter type="column" name="name_col" value="1" />
<filter type="column" name="value_col" value="2" />
</options>
</param>
</when>
<when value="user">
<param format="wig" name="input2" type="data" label="Score file">
<options>
<filter type="data_meta" data_ref="input1" key="build" />
</options>
</param>
</when>
</conditional>
</inputs>
<outputs>
<data format="interval" name="out_file1" metadata_source="input1"/>
</outputs>
<tests>
<test>
<param name="input1" value="6.bed" dbkey="hg17" ftype="bed"/>
<param name="score_source" value="cached"/>
<param name="datasets" value="/depot/data2/galaxy/binned_scores/hg17/phastcons_encode_sep2005_tba" />
<output name="out_file1" file="aggregate_binned_scores_in_intervals.out" />
</test>
<test>
<param name="input1" value="9_hg18.bed" dbkey="hg18" ftype="bed"/>
<param name="score_source" value="cached"/>
<param name="datasets" value="/depot/data2/galaxy/binned_scores/hg18/phastCons17way/ba" />
<output name="out_file1" file="aggregate_binned_scores_in_intervals2.interval" />
</test>
<test>
<param name="input1" value="6.bed" dbkey="hg17" ftype="bed"/>
<param name="score_source" value="user"/>
<param name="input2" value="aggregate_binned_scores_3.wig" dbkey="hg17" ftype="wig"/>
<output name="out_file1" file="aggregate_binned_scores_in_intervals3.out"/>
</test>
</tests>
<help>
.. class:: warningmark
This tool currently only has cached data for genome builds hg16, hg17 and hg18. However, you may use your own data point (wiggle) data, such as is available from UCSC. If you are trying to use your own data point file and it is not appearing as an option, make sure that the builds for your history items are the same.
.. class:: warningmark
This tool assumes that the input dataset is in interval format and contains at least a chrom column, a start column and an end column. These 3 columns can be dispersed throughout any number of other data columns.
-----
.. class:: infomark
**TIP:** Computing summary information may throw exceptions if the data type (e.g., string, integer) in every line of the columns is not appropriate for the computation (e.g., attempting numerical calculations on strings). If an exception is thrown when computing summary information for a line, that line is skipped as invalid for the computation. The number of invalid skipped lines is documented in the resulting history item as a "Data issue".
-----
**Syntax**
This tool appends columns of summary information for each interval matched against a selected dataset. For each interval, the average, minimum and maximum for the data falling within the interval is computed.
- Several quantitative scores are provided for the ENCODE regions.
- Various Scores
- Regulatory Potential
- Neutral rate (Ancestral Repeats)
- GC fraction
- Conservation Scores
- PhastCons
- binCons
- GERP
-----
**Example**
If your original data has the following format:
+------+-----+-----+---+------+
|other1|chrom|start|end|other2|
+------+-----+-----+---+------+
and you choose to aggregate phastCons scores, your output will look like this:
+------+-----+-----+---+------+---+---+---+
|other1|chrom|start|end|other2|avg|min|max|
+------+-----+-----+---+------+---+---+---+
where:
* **avg** - average phastCons score for each region
* **min** - minimum phastCons score for each region
* **max** - maximum phastCons score for each region
</help>
</tool>
@@ -0,0 +1,196 @@
#!/usr/bin/env python2.4
# Greg Von Kuster
"""
usage: %prog score_file interval_file chrom start stop [out_file] [options]
-b, --binned: 'score_file' is actually a directory of binned array files
-m, --mask=FILE: bed file containing regions not to consider valid
"""
from __future__ import division
import pkg_resources
pkg_resources.require( "bx-python" )
pkg_resources.require( "lrucache" )
try:
pkg_resources.require( "python-lzo" )
except:
pass
import psyco_full
import sys
import os, os.path
from UserDict import DictMixin
import bx.wiggle
from bx.binned_array import BinnedArray, FileBinnedArray
from bx.bitset import *
from bx.bitset_builders import *
from fpconst import isNaN
from bx.cookbook import doc_optparse
from galaxy.tools.exception_handling import *
class FileBinnedArrayDir( DictMixin ):
"""
Adapter that makes a directory of FileBinnedArray files look like
a regular dict of BinnedArray objects.
"""
def __init__( self, dir ):
self.dir = dir
self.cache = dict()
def __getitem__( self, key ):
value = None
if key in self.cache:
value = self.cache[key]
else:
fname = os.path.join( self.dir, "%s.ba" % key )
if os.path.exists( fname ):
value = FileBinnedArray( open( fname ) )
self.cache[key] = value
if value is None:
raise KeyError( "File does not exist: " + fname )
return value
def stop_err(msg):
sys.stderr.write(msg)
sys.exit()
def load_scores_wiggle( fname ):
"""
Read a wiggle file and return a dict of BinnedArray objects keyed
by chromosome.
"""
scores_by_chrom = dict()
try:
for chrom, pos, val in bx.wiggle.Reader( UCSCOutWrapper(open( fname ) ) ):
if chrom not in scores_by_chrom:
scores_by_chrom[chrom] = BinnedArray()
scores_by_chrom[chrom][pos] = val
except UCSCLimitException:
# Wiggle data was truncated, at the very least need to warn the user.
print 'Encountered message from UCSC: "Reached output limit of 100000 data values", so be aware your data was truncated.'
except IndexError:
stop_err('Data error: one or more column data values is missing in "%s"' %fname)
except ValueError:
stop_err('Data error: invalid data type for one or more values in "%s".' %fname)
return scores_by_chrom
def load_scores_ba_dir( dir ):
"""
Return a dict-like object (keyed by chromosome) that returns
FileBinnedArray objects created from "key.ba" files in `dir`
"""
return FileBinnedArrayDir( dir )
def main():
# Parse command line
options, args = doc_optparse.parse( __doc__ )
try:
score_fname = args[0]
interval_fname = args[1]
chrom_col = args[2]
start_col = args[3]
stop_col = args[4]
if len( args ) > 5:
out_file = open( args[5], 'w' )
else:
out_file = sys.stdout
binned = bool( options.binned )
mask_fname = options.mask
except:
doc_optparse.exit()
if score_fname == 'None':
stop_err( 'This tool works with data from genome builds hg16, hg17 or hg18. Click the pencil icon in your history item to set the genome build if appropriate.' )
try:
chrom_col = int(chrom_col) - 1
start_col = int(start_col) - 1
stop_col = int(stop_col) - 1
except:
stop_err( 'Chrom, start & end column not properly set, click the pencil icon in your history item to set these values.' )
if chrom_col < 0 or start_col < 0 or stop_col < 0:
stop_err( 'Chrom, start & end column not properly set, click the pencil icon in your history item to set these values.' )
if binned:
scores_by_chrom = load_scores_ba_dir( score_fname )
else:
scores_by_chrom = load_scores_wiggle( score_fname )
if mask_fname:
masks = binned_bitsets_from_file( open( mask_fname ) )
else:
masks = None
skipped_lines = 0
first_invalid_line = 0
invalid_line = ''
for i, line in enumerate( open( interval_fname )):
valid = True
line = line.rstrip('\r\n')
if line and not line.startswith( '#' ):
fields = line.split()
try:
chrom, start, stop = fields[chrom_col], int( fields[start_col] ), int( fields[stop_col] )
except:
valid = False
skipped_lines += 1
if not invalid_line:
first_invalid_line = i + 1
invalid_line = line
if valid:
total = 0
count = 0
min_score = 100000000
max_score = -100000000
for j in range( start, stop ):
if chrom in scores_by_chrom:
try:
# Skip if base is masked
if masks and chrom in masks:
if masks[chrom][j]:
continue
# Get the score, only count if not 'nan'
score = scores_by_chrom[chrom][j]
if not isNaN( score ):
total += score
count += 1
max_score = max( score, max_score )
min_score = min( score, min_score )
except:
continue
if count > 0:
avg = total/count
else:
avg = "nan"
min_score = "nan"
max_score = "nan"
# Build the resulting line of data
out_line = []
for k in range(0, len(fields)):
out_line.append(fields[k])
out_line.append(avg)
out_line.append(min_score)
out_line.append(max_score)
print >> out_file, "\t".join( map( str, out_line ) )
else:
skipped_lines += 1
if not invalid_line:
first_invalid_line = i + 1
invalid_line = line
elif line.startswith( '#' ):
# We'll save the original comments
print >> out_file, line
out_file.close()
if skipped_lines > 0:
print 'Data issue: skipped %d invalid lines starting at line #%d which is "%s"' % ( skipped_lines, first_invalid_line, invalid_line )
if skipped_lines == i:
print 'Consider changing the metadata for the input dataset by clicking on the pencil icon in the history item.'
if __name__ == "__main__": main()