changeset 0:88ef41de020d draft default tip

Imported from capsule None
author devteam
date Tue, 01 Apr 2014 10:52:23 -0400
parents
children
files partialR_square.py partialR_square.xml test-data/partialR_result.tabular test-data/regr_inp.tabular tool_dependencies.xml
diffstat 5 files changed, 244 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/partialR_square.py	Tue Apr 01 10:52:23 2014 -0400
@@ -0,0 +1,145 @@
+#!/usr/bin/env python
+
+import sys
+from rpy import *
+import numpy
+
+#export PYTHONPATH=~/galaxy/lib/
+#running command python partialR_square.py reg_inp.tab 4 1,2,3 partialR_result.tabular
+
+def stop_err(msg):
+    sys.stderr.write(msg)
+    sys.exit()
+
+
+def sscombs(s):
+    if len(s) == 1:
+        return [s]
+    else:
+        ssc = sscombs(s[1:])
+        return [s[0]] + [s[0]+comb for comb in ssc] + ssc
+
+
+infile = sys.argv[1]
+y_col = int(sys.argv[2])-1
+x_cols = sys.argv[3].split(',')
+outfile = sys.argv[4]
+
+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([])
+    """
+    try:
+        float( elems[x_cols[k]] )
+    except:
+        try:
+            msg = "This operation cannot be performed on non-numeric column %d containing value '%s'." % ( col, elems[x_cols[k]] )
+        except:
+            msg = "This operation cannot be performed on non-numeric data."
+        stop_err( msg )
+    """
+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')
+                #print >> sys.stderr, "ey = %s" %ey
+            y_vals.append(yval)
+            for k, col in enumerate(x_cols):
+                try:
+                    xval = float(fields[col])
+                except Exception, ex:
+                    xval = r('NA')
+                    #print >> sys.stderr, "ex = %s" %ex
+                x_vals[k].append(xval)
+        except:
+            pass
+
+x_vals1 = numpy.asarray(x_vals).transpose()
+dat = r.list(x=array(x_vals1), y=y_vals)
+
+set_default_mode(NO_CONVERSION)
+try:
+    full = r.lm(r("y ~ x"), data= r.na_exclude(dat))    #full model includes all the predictor variables specified by the user
+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(full)
+fullr2 = summary.get('r.squared','NA')
+
+if fullr2 == 'NA':
+    stop_err("Error in linear regression")
+
+if len(x_vals) < 10:
+    s = ""
+    for ch in range(len(x_vals)):
+        s += str(ch)
+else:
+    stop_err("This tool only works with less than 10 predictors.")
+
+print >> fout, "#Model\tR-sq\tpartial_R_Terms\tpartial_R_Value"
+all_combos = sorted(sscombs(s), key=len)
+all_combos.reverse()
+for j, cols in enumerate(all_combos):
+    #if len(cols) == len(s):    #Same as the full model above
+    #    continue
+    if len(cols) == 1:
+        x_vals1 = x_vals[int(cols)]
+    else:
+        x_v = []
+        for col in cols:
+            x_v.append(x_vals[int(col)])
+        x_vals1 = numpy.asarray(x_v).transpose()
+    dat = r.list(x=array(x_vals1), y=y_vals)
+    set_default_mode(NO_CONVERSION)
+    red = r.lm(r("y ~ x"), data= dat)    #Reduced model
+    set_default_mode(BASIC_CONVERSION)
+    summary = r.summary(red)
+    redr2 = summary.get('r.squared','NA')
+    try:
+        partial_R = (float(fullr2)-float(redr2))/(1-float(redr2))
+    except:
+        partial_R = 'NA'
+    col_str = ""
+    for col in cols:
+        col_str = col_str + str(int(x_cols[int(col)]) + 1) + " "
+    col_str.strip()
+    partial_R_col_str = ""
+    for col in s:
+        if col not in cols:
+            partial_R_col_str = partial_R_col_str + str(int(x_cols[int(col)]) + 1) + " "
+    partial_R_col_str.strip()
+    if len(cols) == len(s):    #full model
+        partial_R_col_str = "-"
+        partial_R = "-"
+    try:
+        redr2 = "%.4f" % (float(redr2))
+    except:
+        pass
+    try:
+        partial_R = "%.4f" % (float(partial_R))
+    except:
+        pass
+    print >> fout, "%s\t%s\t%s\t%s" % ( col_str, redr2, partial_R_col_str, partial_R )
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/partialR_square.xml	Tue Apr 01 10:52:23 2014 -0400
@@ -0,0 +1,73 @@
+<tool id="partialRsq" name="Compute partial R square" version="1.0.0">
+  <description> </description>
+  <requirements>
+    <requirement type="package" version="2.11.0">R</requirement>
+    <requirement type="package" version="1.7.1">numpy</requirement>
+    <requirement type="package" version="1.0.3">rpy</requirement>
+  </requirements>
+  <command interpreter="python">
+    partialR_square.py 
+      $input1
+      $response_col
+      $predictor_cols
+      $out_file1
+      1>/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" />
+  </outputs>
+  <requirements>
+    <requirement type="python-module">rpy</requirement>
+  </requirements>
+  <tests>
+    <!-- Test data with vlid values -->
+  	<test>
+      <param name="input1" value="regr_inp.tabular"/>
+      <param name="response_col" value="3"/>
+      <param name="predictor_cols" value="1,2"/>
+      <output name="out_file1" file="partialR_result.tabular"/>
+    </test>
+    
+  </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 computes the Partial R squared for all possible variable subsets using the following formula:
+
+**Partial R squared = [SSE(without i: 1,2,...,p-1) - SSE (full: 1,2,..,i..,p-1) / SSE(without i: 1,2,...,p-1)]**, which denotes the case where the 'i'th predictor is dropped. 
+
+
+
+In general, **Partial R squared = [SSE(without i: 1,2,...,p-1) - SSE (full: 1,2,..,i..,p-1) / SSE(without i: 1,2,...,p-1)]**, where,
+
+- SSE (full: 1,2,..,i..,p-1) = Sum of Squares left out by the full set of predictors SSE(X1, X2 … Xp)
+- SSE (full: 1,2,..,i..,p-1) = Sum of Squares left out by the set of predictors excluding; for example, if we omit the first predictor, it will be SSE(X2 … Xp).
+
+
+The 4 columns in the output are described below:
+
+- Column 1 (Model): denotes the variables present in the model
+- Column 2 (R-sq): denotes the R-squared value corresponding to the model in Column 1
+- Column 3 (Partial R squared_Terms): denotes the variable/s for which Partial R squared is computed. These are the variables that are absent in the reduced model in Column 1. A '-' in this column indicates that the model in Column 1 is the Full model.
+- Column 4 (Partial R squared): denotes the Partial R squared value corresponding to the variable/s in Column 3. A '-' in this column indicates that the model in Column 1 is the Full model.
+
+*R Development Core Team (2010). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.*  
+  
+  </help>
+</tool>
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/partialR_result.tabular	Tue Apr 01 10:52:23 2014 -0400
@@ -0,0 +1,4 @@
+#Model	R-sq	partial_R_Terms	partial_R_Value
+1 2 	0.9388	-	-
+2 	0.7280	1 	0.7750
+1 	0.9104	2 	0.3167
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/regr_inp.tabular	Tue Apr 01 10:52:23 2014 -0400
@@ -0,0 +1,10 @@
+7	33	42
+4	41	33
+16	7	75
+3	49	28
+21	5	91
+8	31	55
+7	35	52
+5	30	16
+15	10	69
+20	10	94
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tool_dependencies.xml	Tue Apr 01 10:52:23 2014 -0400
@@ -0,0 +1,12 @@
+<?xml version="1.0"?>
+<tool_dependency>
+  <package name="numpy" version="1.7.1">
+      <repository changeset_revision="0c288abd2a1e" name="package_numpy_1_7" owner="devteam" prior_installation_required="False" toolshed="http://toolshed.g2.bx.psu.edu" />
+    </package>
+    <package name="rpy" version="1.0.3">
+      <repository changeset_revision="82170c94ca7c" name="package_rpy_1_0_3" owner="devteam" prior_installation_required="False" toolshed="http://toolshed.g2.bx.psu.edu" />
+    </package>
+    <package name="R" version="2.11.0">
+      <repository changeset_revision="6e1b17857732" name="package_r_2_11_0" owner="devteam" prior_installation_required="False" toolshed="http://toolshed.g2.bx.psu.edu" />
+    </package>
+</tool_dependency>