# HG changeset patch
# User devteam
# Date 1396363817 14400
# Node ID 54c7a01a2cc77c8fc816d9dc5029d81e8084fcad

Imported from capsule None

diff -r 000000000000 -r 54c7a01a2cc7 best_regression_subsets.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/best_regression_subsets.py	Tue Apr 01 10:50:17 2014 -0400
@@ -0,0 +1,91 @@
+#!/usr/bin/env python
+
+from galaxy import eggs
+
+import sys
+from rpy import *
+import numpy
+
+def stop_err(msg):
+    sys.stderr.write(msg)
+    sys.exit()
+
+
+infile = sys.argv[1]
+y_col = int(sys.argv[2])-1
+x_cols = sys.argv[3].split(',')
+outfile = sys.argv[4]
+outfile2 = sys.argv[5]
+print "Predictor columns: %s; Response column: %d" % ( x_cols, y_col+1 )
+fout = open(outfile,'w')
+
+for i, line in enumerate( file ( infile )):
+    line = line.rstrip('\r\n')
+    if len( line )>0 and not line.startswith( '#' ):
+        elems = line.split( '\t' )
+        break
+    if i == 30:
+        break # Hopefully we'll never get here...
+
+if len( elems )<1:
+    stop_err( "The data in your input dataset is either missing or not formatted properly." )
+
+y_vals = []
+x_vals = []
+
+for k, col in enumerate(x_cols):
+    x_cols[k] = int(col)-1
+    x_vals.append([])
+    
+NA = 'NA'
+for ind, line in enumerate( file( infile ) ):
+    if line and not line.startswith( '#' ):
+        try:
+            fields = line.split("\t")
+            try:
+                yval = float(fields[y_col])
+            except Exception, ey:
+                yval = r('NA')
+            y_vals.append(yval)
+            for k, col in enumerate(x_cols):
+                try:
+                    xval = float(fields[col])
+                except Exception, ex:
+                    xval = r('NA')
+                x_vals[k].append(xval)
+        except:
+            pass
+
+response_term = ""
+
+x_vals1 = numpy.asarray(x_vals).transpose()
+
+dat = r.list(x=array(x_vals1), y=y_vals)
+
+r.library("leaps")
+
+set_default_mode(NO_CONVERSION)
+try:
+    leaps = r.regsubsets(r("y ~ x"), data= r.na_exclude(dat))
+except RException, rex:
+    stop_err("Error performing linear regression on the input data.\nEither the response column or one of the predictor columns contain no numeric values.")
+set_default_mode(BASIC_CONVERSION)
+
+summary = r.summary(leaps)
+tot = len(x_vals)
+pattern = "["
+for i in range(tot):
+    pattern = pattern + 'c' + str(int(x_cols[int(i)]) + 1) + ' '
+pattern = pattern.strip() + ']'
+print >> fout, "#Vars\t%s\tR-sq\tAdj. R-sq\tC-p\tbic" % (pattern)
+for ind, item in enumerate(summary['outmat']):
+    print >> fout, "%s\t%s\t%s\t%s\t%s\t%s" % (str(item).count('*'), item, summary['rsq'][ind], summary['adjr2'][ind], summary['cp'][ind], summary['bic'][ind])
+
+
+r.pdf( outfile2, 8, 8 )
+r.plot(leaps, scale="Cp", main="Best subsets using Cp Criterion")
+r.plot(leaps, scale="r2", main="Best subsets using R-sq Criterion")
+r.plot(leaps, scale="adjr2", main="Best subsets using Adjusted R-sq Criterion")
+r.plot(leaps, scale="bic", main="Best subsets using bic Criterion")
+
+r.dev_off()
diff -r 000000000000 -r 54c7a01a2cc7 best_regression_subsets.xml
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/best_regression_subsets.xml	Tue Apr 01 10:50:17 2014 -0400
@@ -0,0 +1,66 @@
+<tool id="BestSubsetsRegression1" name="Perform Best-subsets Regression" version="0.0.1">
+  <description> </description>
+  <command interpreter="python">
+    best_regression_subsets.py 
+      $input1
+      $response_col
+      $predictor_cols
+      $out_file1
+      $out_file2
+      1>/dev/null
+      2>/dev/null
+  </command>
+  <inputs>
+    <param format="tabular" name="input1" type="data" label="Select data" help="Dataset missing? See TIP below."/>
+    <param name="response_col" label="Response column (Y)" type="data_column" data_ref="input1" />
+    <param name="predictor_cols" label="Predictor columns (X)" type="data_column" data_ref="input1" multiple="true" >
+        <validator type="no_options" message="Please select at least one column."/>
+    </param>
+  </inputs>
+  <outputs>
+    <data format="input" name="out_file1" metadata_source="input1" />
+    <data format="pdf" name="out_file2" />
+  </outputs>
+  <requirements>
+    <requirement type="python-module">rpy</requirement>
+  </requirements>
+  <tests>
+    <!-- Testing this tool will not be possible because this tool produces a pdf output file.
+    -->
+  </tests>
+  <help>
+
+.. class:: infomark
+
+**TIP:** If your data is not TAB delimited, use *Edit Datasets-&gt;Convert characters*
+
+-----
+
+.. class:: infomark
+
+**What it does**
+
+This tool uses the 'regsubsets' function from R statistical package for regression subset selection. It outputs two files, one containing a table with the best subsets and the corresponding summary statistics, and the other containing the graphical representation of the results.  
+
+-----
+
+.. class:: warningmark
+
+**Note**
+
+- This tool currently treats all predictor and response variables as continuous variables. 
+
+- Rows containing non-numeric (or missing) data in any of the chosen columns will be skipped from the analysis.
+
+- The 6 columns in the output are described below:
+
+  - Column 1 (Vars): denotes the number of variables in the model
+  - Column 2 ([c2 c3 c4...]): represents a list of the user-selected predictor variables (full model). An asterix denotes the presence of the corresponding predictor variable in the selected model.
+  - Column 3 (R-sq): the fraction of variance explained by the model
+  - Column 4 (Adj. R-sq): the above R-squared statistic adjusted, penalizing for higher number of predictors (p)
+  - Column 5 (Cp): Mallow's Cp statistics  
+  - Column 6 (bic): Bayesian Information Criterion. 
+
+
+  </help>
+</tool>