From 335ebef0afb66c7d4d2168f6935eec81588bb8ff Mon Sep 17 00:00:00 2001
From: Greg Von Kuster <greg@bx.psu.edu>
Date: Thu, 28 Feb 2008 00:08:05 +0000
Subject: [PATCH] Restored version 1.1.0 of aggregate scores tool and moved the
 fixes to version 1.1.1.

---
 tool_conf.xml.main                            |  67 +-----
 tool_conf.xml.sample                          |  67 +-----
 .../aggregate_binned_scores_in_intervals.xml  |   4 +-
 .../1.1.0/aggregate_scores_in_intervals.py    |   2 +-
 .../aggregate_binned_scores_in_intervals.xml  | 115 ++++++++++
 .../1.1.1/aggregate_scores_in_intervals.py    | 196 ++++++++++++++++++
 6 files changed, 315 insertions(+), 136 deletions(-)
 create mode 100644 tools/stats/aggregate_scores_in_intervals2/1.1.1/aggregate_binned_scores_in_intervals.xml
 create mode 100755 tools/stats/aggregate_scores_in_intervals2/1.1.1/aggregate_scores_in_intervals.py

diff --git a/tool_conf.xml.main b/tool_conf.xml.main
index 5fe6ade6020..26c1e9c211e 100644
--- a/tool_conf.xml.main
+++ b/tool_conf.xml.main
@@ -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" />
diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample
index 266041134c5..59d2adc3172 100644
--- a/tool_conf.xml.sample
+++ b/tool_conf.xml.sample
@@ -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 
diff --git a/tools/stats/aggregate_scores_in_intervals2/1.1.0/aggregate_binned_scores_in_intervals.xml b/tools/stats/aggregate_scores_in_intervals2/1.1.0/aggregate_binned_scores_in_intervals.xml
index 9d1ec13e971..221c551d988 100644
--- a/tools/stats/aggregate_scores_in_intervals2/1.1.0/aggregate_binned_scores_in_intervals.xml
+++ b/tools/stats/aggregate_scores_in_intervals2/1.1.0/aggregate_binned_scores_in_intervals.xml
@@ -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>
diff --git a/tools/stats/aggregate_scores_in_intervals2/1.1.0/aggregate_scores_in_intervals.py b/tools/stats/aggregate_scores_in_intervals2/1.1.0/aggregate_scores_in_intervals.py
index faaa7935d25..49514c33f80 100755
--- a/tools/stats/aggregate_scores_in_intervals2/1.1.0/aggregate_scores_in_intervals.py
+++ b/tools/stats/aggregate_scores_in_intervals2/1.1.0/aggregate_scores_in_intervals.py
@@ -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.'
diff --git a/tools/stats/aggregate_scores_in_intervals2/1.1.1/aggregate_binned_scores_in_intervals.xml b/tools/stats/aggregate_scores_in_intervals2/1.1.1/aggregate_binned_scores_in_intervals.xml
new file mode 100644
index 00000000000..3f9c95f464f
--- /dev/null
+++ b/tools/stats/aggregate_scores_in_intervals2/1.1.1/aggregate_binned_scores_in_intervals.xml
@@ -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>
diff --git a/tools/stats/aggregate_scores_in_intervals2/1.1.1/aggregate_scores_in_intervals.py b/tools/stats/aggregate_scores_in_intervals2/1.1.1/aggregate_scores_in_intervals.py
new file mode 100755
index 00000000000..faaa7935d25
--- /dev/null
+++ b/tools/stats/aggregate_scores_in_intervals2/1.1.1/aggregate_scores_in_intervals.py
@@ -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()