Better memory management, extensive code cleanup for Summary Statistics tool. Added new EmptyTextField validator. Deleted ExtractAxt_wrapper tool from code base since it is no longer used.

This commit is contained in:
Greg Von Kuster
2008-04-24 20:29:47 +00:00
parent ea8a46baaa
commit 5c1cae4e57
5 changed files with 105 additions and 389 deletions
+14
View File
@@ -205,6 +205,19 @@ class NoOptionsValidator( Validator ):
self.message = "No options available for selection"
raise ValueError( self.message )
class EmptyTextfieldValidator( Validator ):
"""Validator that checks for empty text field"""
def __init__( self, message=None ):
self.message = message
@classmethod
def from_element( cls, param, elem ):
return cls( elem.get( 'message', None ) )
def validate( self, value, history=None ):
if value == '':
if self.message is None:
self.message = "Field requires a value"
raise ValueError( self.message )
class MetadataInFileColumnValidator( Validator ):
"""
Validator that checks if the value for a dataset's metadata item exists in a file.
@@ -247,6 +260,7 @@ validator_types = dict( expression=ExpressionValidator,
metadata=MetadataValidator,
unspecified_build=UnspecifiedBuildValidator,
no_options=NoOptionsValidator,
empty_field=EmptyTextfieldValidator,
dataset_metadata_in_file=MetadataInFileColumnValidator,
dataset_ok_validator=DatasetOkValidator )
-66
View File
@@ -1,66 +0,0 @@
#! /usr/bin/perl -w
# Galaxy (universe) wrapper for Rico's extractAxt
# For this to work extractAxt should be intalled in bins/
# directory of universe's tools section
# Takes the following parameters:
# extractorAxt_wrapper.pl -i $inp_file1 -o $out_file1 --species $species -g $dbkey $chroCol $startCol $endCol $strandCol
# Location of alignment files is taken from $GALAXY_DATA_INDEX_DIR/alignseq.loc (to change -> edit line 19)
use strict;
use warnings;
use File::Temp "tempfile";
die "Your query genome, $ARGV[7], is the same as your target genome, $ARGV[5]. Please go back and select different target genome\n" if ($ARGV[5] eq $ARGV[7]);
die "Not enough params -> check\n" unless @ARGV == 12;
my $alignseqLoc = "$GALAXY_DATA_INDEX_DIR/alignseq.loc";
my %alignLocation = ();
my @locFields = ();
my $extractAxtStatus = 0;
my $alignDir = "";
my @bed = ();
$ARGV[5] = "musMus$1" if ($ARGV[5]=~ m/^mm(\d)$/);
$ARGV[7] = "musMus$1" if ($ARGV[7]=~ m/^mm(\d)$/);
$ARGV[5] = "ratNor$1" if ($ARGV[5]=~ m/^rn(\d)$/);
$ARGV[7] = "ratNor$1" if ($ARGV[7]=~ m/^rn(\d)$/);
open (LOC, "<$alignseqLoc") or die "Cannot open $alignseqLoc:$!\n";
while (<LOC>) {
if (m/^align/) {
chop;
@locFields = split " ";
$alignLocation{$locFields[1]."-".$locFields[2]} = "$locFields[3]/";
}
}
close LOC;
$alignDir = $alignLocation{"$ARGV[7]-$ARGV[5]"};
die "No alignments between $ARGV[7] and $ARGV[5] are presently stored on Galaxy site. E-mail to galaxy-bugs\@bx.psu.edu to request them\n" if !defined($alignDir);
if ($ARGV[8] == 1 and $ARGV[9] == 2 and $ARGV[10] == 3 and $ARGV[11] == 6) {
$extractAxtStatus = system("extractAxt -b $ARGV[1] -c -n -o $ARGV[3] $alignDir");
} else {
my ($fh, $filename) = tempfile();
open (DATA, "<$ARGV[1]") or die "Cannot open $ARGV[1]:$!\n";
while (<DATA>) {
chop;
my @line = split /\t/;
$line[$ARGV[11]-1] = "+" if !defined$line[$ARGV[11]-1];
my $nameLine = substr(join("-",@line),0,50);
my $bedLine = "$line[$ARGV[8]-1]\t$line[$ARGV[9]-1]\t$line[$ARGV[10]-1]\t\t$nameLine\t$line[$ARGV[11]-1]\n";
print $fh $bedLine if !m/^#/;
}
close DATA;
$extractAxtStatus = system("extractAxt -b $filename -c -n -o $ARGV[3] $alignDir");
`rm -f $filename`;
}
die "axt extractor exited abnormally: $?. E-mail to galaxy-bugs\@bx.psu.edu to report this problem\n" unless $extractAxtStatus == 0;
-65
View File
@@ -1,65 +0,0 @@
<tool id="Extract blastz alignments1" name="Extract blastz alignments">
<description>between query genome and another genome</description>
<command interpreter="perl">extractAxt_wrapper.pl -i $input -o $out_file1 --species $species -g $dbkey $input_chromCol $input_startCol $input_endCol $input_strandCol</command>
<inputs>
<param format="interval" name="input" type="data" label="Between regions of Query"/>
<param name="species" type="select" label="and one of these genomes">
<validator type="unspecified_build" />
<options from_file="alignseq.loc" name_col="1" value_col="2">
<filter type="data_meta" data_ref="input" meta_key="dbkey" meta_key_col="0" />
</options>
</param>
</inputs>
<outputs>
<data format="axt" name="out_file1" />
</outputs>
<tests>
<test>
<param name="input" value="1.bed" dbkey="hg17" ftype="bed" />
<param name="species" value="musMus6"/>
<output name="out_file1" file="fsa_extract_blastz_alignments.dat" />
</test>
</tests>
<help>
.. class:: warningmark
**IMPORTANT**: AXT formatted alignments will be phased out from Galaxy in the coming weeks. They will be replaced with pairwise MAF alignments, which are already available. To try pairwise MAF alignments use "Extract Pairwise MAF blocks" tool in *Fetch Sequences and Alignments* section.
--------
.. class:: infomark
**TIP:** The last field of the axt header is the same as the column 4 (name) of the corresponding BED file. *This is very useful* for establishing correspondence between alignments and the original BED file.
-----
**Syntax**
This tool uses coordinates specified in the query to extract alignments, which is pre-computed blastZ output. Alignments are in the axt format as shown below.
- **blastz alignments** Blastz alignment program is available from Webb Miller's lab at Penn State University (http://www.bx.psu.edu/miller_lab/).
- **AXT format** The alignments are produced from Blastz. The lav format Blastz output, which does not include the sequence, was converted to AXT format with lavToAxt. Each alignment block in an AXT file contains three lines: a summary line and 2 sequence lines. Blocks are separated from one another by blank lines.
-----
**Example**
- Input file::
chr7 127486022 127486166 NM_000230 0 +
chr7 127486011 127486166 D49487 0 +
- Extract the hg17 and mm5 assemblies blastZ alignments of the above file::
0 chr7 127475282 127475310 chr6 28912408 28912436 + 53096 NM_000230
GTAGGAATCGCAGCGCCAGCGGTTGCAAG
GGAGGGATCCCTGCTCCAGCAGCTGCAAG
1 chr7 127486012 127486166 chr6 28921129 28921283 + 61453 D49487
TGGGAAGGAAAATGCATTGGGGAACCCTGTGCGGATTCTTGTGGCTTTGGCCCTATCTTTTCTATGTCCAAGCTGTGCCCATCCAAAAAGTCCAAGATGACACCAAAACCCTCATCAAGACAATTGTCACCAGGATCAATGACATTTCACACACG
CAGGGAGGAAAATGTGCTGGAGACCCCTGTGTCGGTTCCTGTGGCTTTGGTCCTATCTGTCTTATGTTCAAGCAGTGCCTATCCAGAAAGTCCAGGATGACACCAAAACCCTCATCAAGACCATTGTCACCAGGATCAATGACATTTCACACACG
</help>
</tool>
+85 -254
View File
@@ -1,277 +1,108 @@
#!/usr/bin/python
# Greg Von Kuster
import sys
import sys, sets, re, tempfile
from rpy import *
import sets,re
def stop_err(msg):
sys.stderr.write(msg)
assert sys.version_info[:2] >= ( 2, 4 )
def stop_err( msg ):
sys.stderr.write( msg )
sys.exit()
def mode_func(c):
try:
check = float(c)
return "r.as_numeric"
except:
return "r.as_factor"
def order_for_display(l):
if l[0].startswith("chr"):
l.sort(byChr)
else:
l.sort() # alphanumerically
return l
def byChr(a,b):
fa = a.split()
fb = b.split()
if chrint(fa[0]) < chrint(fb[0]):
return -1
elif chrint(fa[0]) > chrint(fb[0]):
return 1
else:
if len(fa) > 1 and len(fb) > 1:
if int(fa[1]) < int(fb[1]):
return -1
elif int(fa[1]) > int(fb[1]):
return 1
else:
if int(fa[2]) < int(fb[2]):
return -1
elif int(fa[2]) > int(fb[2]):
return 1
return 0
def chrint( x ):
i = x.replace("chr","")
if (i == "X"):
return 23
elif (i == "Y"):
return 24
elif (i == "Un"):
return 25
# for randoms and whatever
elif ( i.find("_") > -1):
parts = i.split("_")
i = parts[0]
if (i == "X"):
return 23 + 25
elif (i == "Y"):
return 24 + 25
elif (i == "Un"):
return 24 + 25
try:
return 24 + int(i)
except ValueError:
return 0
else:
try:
return int(i)
except ValueError:
return 0
def S3_METHODS(all="key"):
# See the R documentation for groupGeneric
# help(.Method)
Group_Math = [ "abs", "sign", "sqrt",
"floor", "ceiling", "trunc",
"round", "signif",
"exp", "log",
"cos", "sin", "tan",
"acos", "asin", "atan",
"cosh", "sinh", "tanh",
"acosh", "asinh", "atanh",
"lgamma", "gamma", "gammaCody",
"digamma", "trigamma",
"cumsum", "cumprod", "cummax", "cummin"]
Group_Ops = [ "+", "-", "*", "/", "^", "%%", "%/%", "&", "|", "!", "==", "!=", "<", "<=", ">=", ">"]
Group_Summary = [ "all", "any", "sum", "prod", "min", "max", "range" ]
def S3_METHODS( all="key" ):
Group_Math = [ "abs", "sign", "sqrt", "floor", "ceiling", "trunc", "round", "signif",
"exp", "log", "cos", "sin", "tan", "acos", "asin", "atan", "cosh", "sinh", "tanh",
"acosh", "asinh", "atanh", "lgamma", "gamma", "gammaCody", "digamma", "trigamma",
"cumsum", "cumprod", "cummax", "cummin", "c" ]
Group_Ops = [ "+", "-", "*", "/", "^", "%%", "%/%", "&", "|", "!", "==", "!=", "<", "<=", ">=", ">", "(", ")", "~", "," ]
if all is "key":
return { 'Math' : Group_Math, 'Ops' : Group_Ops, 'Summary' : Group_Summary }
def read_table(datafile, cols):
table = {}
width = 0
skipped_lines = 0
first_invalid_line = 0
for i, line in enumerate(file(datafile)):
valid = True
line = line.rstrip('\r\n')
if line and not line.startswith( '#' ):
f = line.split( '\t' )
for col in cols:
if col > len(f):
valid = False
skipped_lines += 1
if not first_invalid_line:
first_invalid_line = i+1
break
if valid:
# Make sure the column value is numeric
try:
check = float(f[col-1])
except:
valid = False
skipped_lines += 1
if not first_invalid_line:
first_invalid_line = i+1
break
if valid:
for col_i, val in enumerate(f):
# Create column names c1..cn
colname = "c" + str(col_i + 1)
if not table.has_key( colname ):
table[colname] = []
table[colname].append( val )
if not width:
width = len(f)
if len(table) > 0:
# terms will look like this: c7 = r.as_numeric(table["c7"]), c8 = r.as_numeric(table["c8"]), ...
terms = ["%s = %s(table[\"%s\"])" % (x, mode_func(table[x][0]), x) for x in ["c" + str(col_i + 1) for col_i in range(0, width)]]
code = "d = r.data_frame(%s)" % ",".join(terms)
try:
exec code
except Exception, e:
stop_err(str(e))
return (skipped_lines, first_invalid_line, d)
else:
return (skipped_lines, first_invalid_line, None)
return { 'Math' : Group_Math, 'Ops' : Group_Ops }
def main():
if len(sys.argv) >= 4:
try:
datafile = sys.argv[1]
outfile = sys.argv[2]
outfile_name = sys.argv[2]
expression = sys.argv[3]
else:
print sys.argv
stop_err('Usage: python gsummary.py input_file ouput_file expression')
except:
stop_err( 'Usage: python gsummary.py input_file ouput_file expression' )
if len(sys.argv) == 5:
if sys.argv[4].find('none') < 0:
tmp_rhs = sys.argv[4]
tmp_rhs = tmp_rhs.replace(',','/')
group_terms = re.compile('c[0-9]+').findall(tmp_rhs)
if len(group_terms) > 0:
dep_var = group_terms[0]
tmp_rhs = "|".join([ dep_var, tmp_rhs])
expression = '~'.join([expression,tmp_rhs])
else:
stop_err("%s unrecognized for groups" % tmp_rhs)
elif sys.argv[4] is 'none':
pass
math_allowed = S3_METHODS()[ 'Math' ]
ops_allowed = S3_METHODS()[ 'Ops' ]
# summary function and return labels
f = r("function(x) { c(sum=sum(x,na.rm=T),mean=mean(x,na.rm=T),stdev=sd(x,na.rm=T),quantile(x,na.rm=TRUE))}")
returns = ['sum', 'mean','stdev','0%', '25%', '50%', '75%', '100%']
lhs = ""
rhs = ""
lhs_allowed = S3_METHODS()['Math']
lhs_allowed.append( 'c' )
ops_allowed = S3_METHODS()['Ops']
ops_allowed = ops_allowed + ['(', ')', '~', ',']
of = open(outfile,'w')
if expression.find("~") > 0:
lhs,rhs = expression.split('~')
else:
lhs = expression
for word in re.compile('[a-zA-Z]+').findall(expression):
if word and not word in lhs_allowed:
of.close()
# Check for invalid expressions
for word in re.compile( '[a-zA-Z]+' ).findall( expression ):
if word and not word in math_allowed:
stop_err( "Invalid expression '%s': term '%s' is not recognized or allowed" %( expression, word ) )
"""
Users sometimes want statistics for more than 1 column, so they enter a comma-separated
string of columns in the free text field. This tool only handles a single column or an
expression (computed for 1 or more columns), so we'll use the following hack to provide a
useful response for multiple column entries.
"""
symbols = sets.Set()
for symbol in re.compile('[^a-z0-9\s]+').findall(expression):
for symbol in re.compile( '[^a-z0-9\s]+' ).findall( expression ):
if symbol and not symbol in ops_allowed:
of.close()
stop_err( "Invalid expression '%s': operator '%s' is not recognized or allowed" %(expression, symbol ) )
stop_err( "Invalid expression '%s': operator '%s' is not recognized or allowed" % ( expression, symbol ) )
else:
symbols.add(symbol)
if len(symbols) == 1 and ',' in symbols:
of.close()
stop_err( "Invalid columns '%s': this tool requires a single column or expression" %expression )
symbols.add( symbol )
if len( symbols ) == 1 and ',' in symbols:
# User may have entered a comma-separated list r_data_frame columns
stop_err( "Invalid columns '%s': this tool requires a single column or expression" % expression )
# Find all column references in the expression
cols = []
if lhs:
for col in re.compile('c[0-9]+').findall(lhs):
try:
cols.append(int(col[1:]))
except:
# This is a weakness, but hopefully we won't arrive here
pass
for col in re.compile( 'c[0-9]+' ).findall( expression ):
try:
cols.append( int( col[1:] ) - 1 )
except:
pass
tmp_file = tempfile.NamedTemporaryFile( 'w+b' )
# Write the R header row to the temporary file
hdr_str = "\t".join( "c%s" % str( col+1 ) for col in cols )
tmp_file.write( "%s\n" % hdr_str )
skipped_lines = 0
first_invalid_line = 0
for i, line in enumerate( file( datafile ) ):
line = line.rstrip( '\r\n' )
if line and not line.startswith( '#' ):
valid = True
fields = line.split( '\t' )
# Write the R data row to the temporary file
for col in cols:
try:
float( fields[ col ] )
except:
skipped_lines += 1
if not first_invalid_line:
first_invalid_line = i + 1
valid = False
break
if valid:
data_str = "\t".join( fields[ col ] for col in cols )
tmp_file.write( "%s\n" % data_str )
tmp_file.flush()
if rhs:
for col in re.compile('c[0-9]+').findall(rhs):
try:
cols.append(int(col[1:]))
except:
# This is a weakness, but hopefully we won't arrive here
pass
set_default_mode(NO_CONVERSION)
skipped_lines, first_invalid_line, df = read_table(datafile, cols)
if df is not None:
if not rhs:
for col in re.compile('c[0-9]+').findall(lhs):
r.assign(col,r["$"](df,col))
try:
summary = f(r(lhs))
except RException, s:
# Due to previous checking, this should not occur
of.close()
stop_err("Computation attempted on invalid data in column %s. Exception: %s" %(lhs, s))
summary = summary.as_py(BASIC_CONVERSION)
print >>of,"#%s" % "\t".join(returns)
print >>of,"\t".join([ "%.3f" % (summary[k]) for k in returns])
else:
# Prepare R structures
r.library("nlme",warn_conflicts=r.FALSE)
r.library("lattice",warn_conflicts=r.FALSE)
set_default_mode(NO_CONVERSION)
try:
df_g = r.groupedData(r.expression(expression), df)
df_r = r.groupedData(r.expression(expression), r.data_frame(df_g, response=r.getResponse(df_g)))
except RException, s:
stop-err("Computation attempted on invalid data in column on the left hand side of expression. Exception:\n\t%s" % s)
# Try some plotting stuff
if (0):
outfile = "plots.pdf"
# r.pdf(file=outfile,width=6,height=6)
# r.histogram(df_g)
r.plot_default(df_g)
r.dev_off()
# Apply summary function and returns
summary_obj = r.gsummary(df_r,FUN=f)
summary_response = r["$"](summary_obj,"response").as_py(BASIC_CONVERSION).tolist()
summary_labels = r.rownames(summary_obj).as_py(BASIC_CONVERSION)
print >>of,"%s\t%s" % ("#group","\t".join(returns))
for index,row in enumerate( summary_labels ):
print >>of,"%s: " % row,
print >>of,"\t".join([ "%.3f" % (summary_response[index][k]) for k in returns])
if skipped_lines > 0:
print "..Skipped %d lines in query beginning with line #%d due to data issues. See tool tips for data requirements." % (skipped_lines, first_invalid_line)
if skipped_lines == i + 1:
stop_err( "Invalid column or column data values invalid for computation. See tool tips and syntax for data requirements." )
else:
stop_err( "Entire data column consisting of %d lines invalid for computation. See tool tips for data requirements." %skipped_lines )
# summary function and return labels
summary_func = r( "function( x ) { c( sum=sum( as.numeric( x ), na.rm=T ), mean=mean( as.numeric( x ), na.rm=T ), stdev=sd( as.numeric( x ), na.rm=T ), quantile( as.numeric( x ), na.rm=TRUE ) ) }" )
headings = [ 'sum', 'mean', 'stdev', '0%', '25%', '50%', '75%', '100%' ]
headings_str = "\t".join( headings )
set_default_mode( NO_CONVERSION )
r_data_frame = r.read_table( tmp_file.name, header=True, sep="\t" )
outfile = open( outfile_name, 'w' )
for col in re.compile( 'c[0-9]+' ).findall( expression ):
r.assign( col, r[ "$" ]( r_data_frame, col ) )
try:
summary = summary_func( r( expression ) )
except RException, s:
outfile.close()
stop_err( "Computation resulted in the following error: %s" % str( s ) )
summary = summary.as_py( BASIC_CONVERSION )
outfile.write( "#%s\n" % headings_str )
outfile.write( "%s\n" % "\t".join( [ "%.3f" % ( summary[ k ] ) for k in headings ] ) )
outfile.close()
if skipped_lines:
print "Skipped %d invalid lines beginning with line #%d. See tool tips for data requirements." % ( skipped_lines, first_invalid_line )
if __name__ == "__main__": main()
+6 -4
View File
@@ -1,9 +1,11 @@
<tool id="Summary_Statistics1" name="Summary Statistics">
<description>for any numerical column</description>
<command interpreter="python">gsummary.py $input $out_file1 "$cond" "none"</command>
<command interpreter="python">gsummary.py $input $out_file1 "$cond"</command>
<inputs>
<param format="tabular" name="input" type="data" label="Summary statistics on" help="Query missing? See TIP below"/>
<param name="cond" size="40" type="text" value="c5" label="Column or expression" help="See tool syntax below" />
<param name="cond" size="30" type="text" value="c5" label="Column or expression" help="See syntax below">
<validator type="empty_field" message="Enter a valid column or expression, see syntax below for examples"/>
</param>
</inputs>
<outputs>
<data format="tabular" name="out_file1" />
@@ -11,7 +13,7 @@
<tests>
<test>
<param name="input" value="1.bed"/>
<output name="out_file1" file="sta_summary.dat"/>
<output name="out_file1" file="gsummary_out1.tabular"/>
<param name="cond" value="c2"/>
</test>
</tests>
@@ -23,7 +25,7 @@ This tool expects input datasets to consist of tab-delimited columns (blank or c
.. class:: infomark
**TIP:** If your data is not TAB delimited, use *Edit Queries-&gt;Convert delimiters to TAB*
**TIP:** If your data is not TAB delimited, use *Text Manipulation-&gt;Convert delimiters to TAB*
.. class:: infomark