Mercurial > repos > devteam > partialr_square

comparison partialR_square.py @ 0:88ef41de020d draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Imported from capsule None
author devteam
date Tue, 01 Apr 2014 10:52:23 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:88ef41de020d
1 #!/usr/bin/env python
2
3 import sys
4 from rpy import *
5 import numpy
6
7 #export PYTHONPATH=~/galaxy/lib/
8 #running command python partialR_square.py reg_inp.tab 4 1,2,3 partialR_result.tabular
9
10 def stop_err(msg):
11 sys.stderr.write(msg)
12 sys.exit()
13
14
15 def sscombs(s):
16 if len(s) == 1:
17 return [s]
18 else:
19 ssc = sscombs(s[1:])
20 return [s[0]] + [s[0]+comb for comb in ssc] + ssc
21
22
23 infile = sys.argv[1]
24 y_col = int(sys.argv[2])-1
25 x_cols = sys.argv[3].split(',')
26 outfile = sys.argv[4]
27
28 print "Predictor columns: %s; Response column: %d" % ( x_cols, y_col+1 )
29 fout = open(outfile,'w')
30
31 for i, line in enumerate( file ( infile )):
32 line = line.rstrip('\r\n')
33 if len( line )>0 and not line.startswith( '#' ):
34 elems = line.split( '\t' )
35 break
36 if i == 30:
37 break # Hopefully we'll never get here...
38
39 if len( elems )<1:
40 stop_err( "The data in your input dataset is either missing or not formatted properly." )
41
42 y_vals = []
43 x_vals = []
44
45 for k, col in enumerate(x_cols):
46 x_cols[k] = int(col)-1
47 x_vals.append([])
48 """
49 try:
50 float( elems[x_cols[k]] )
51 except:
52 try:
53 msg = "This operation cannot be performed on non-numeric column %d containing value '%s'." % ( col, elems[x_cols[k]] )
54 except:
55 msg = "This operation cannot be performed on non-numeric data."
56 stop_err( msg )
57 """
58 NA = 'NA'
59 for ind, line in enumerate( file( infile )):
60 if line and not line.startswith( '#' ):
61 try:
62 fields = line.split("\t")
63 try:
64 yval = float(fields[y_col])
65 except Exception, ey:
66 yval = r('NA')
67 #print >> sys.stderr, "ey = %s" %ey
68 y_vals.append(yval)
69 for k, col in enumerate(x_cols):
70 try:
71 xval = float(fields[col])
72 except Exception, ex:
73 xval = r('NA')
74 #print >> sys.stderr, "ex = %s" %ex
75 x_vals[k].append(xval)
76 except:
77 pass
78
79 x_vals1 = numpy.asarray(x_vals).transpose()
80 dat = r.list(x=array(x_vals1), y=y_vals)
81
82 set_default_mode(NO_CONVERSION)
83 try:
84 full = r.lm(r("y ~ x"), data= r.na_exclude(dat)) #full model includes all the predictor variables specified by the user
85 except RException, rex:
86 stop_err("Error performing linear regression on the input data.\nEither the response column or one of the predictor columns contain no numeric values.")
87 set_default_mode(BASIC_CONVERSION)
88
89 summary = r.summary(full)
90 fullr2 = summary.get('r.squared','NA')
91
92 if fullr2 == 'NA':
93 stop_err("Error in linear regression")
94
95 if len(x_vals) < 10:
96 s = ""
97 for ch in range(len(x_vals)):
98 s += str(ch)
99 else:
100 stop_err("This tool only works with less than 10 predictors.")
101
102 print >> fout, "#Model\tR-sq\tpartial_R_Terms\tpartial_R_Value"
103 all_combos = sorted(sscombs(s), key=len)
104 all_combos.reverse()
105 for j, cols in enumerate(all_combos):
106 #if len(cols) == len(s): #Same as the full model above
107 # continue
108 if len(cols) == 1:
109 x_vals1 = x_vals[int(cols)]
110 else:
111 x_v = []
112 for col in cols:
113 x_v.append(x_vals[int(col)])
114 x_vals1 = numpy.asarray(x_v).transpose()
115 dat = r.list(x=array(x_vals1), y=y_vals)
116 set_default_mode(NO_CONVERSION)
117 red = r.lm(r("y ~ x"), data= dat) #Reduced model
118 set_default_mode(BASIC_CONVERSION)
119 summary = r.summary(red)
120 redr2 = summary.get('r.squared','NA')
121 try:
122 partial_R = (float(fullr2)-float(redr2))/(1-float(redr2))
123 except:
124 partial_R = 'NA'
125 col_str = ""
126 for col in cols:
127 col_str = col_str + str(int(x_cols[int(col)]) + 1) + " "
128 col_str.strip()
129 partial_R_col_str = ""
130 for col in s:
131 if col not in cols:
132 partial_R_col_str = partial_R_col_str + str(int(x_cols[int(col)]) + 1) + " "
133 partial_R_col_str.strip()
134 if len(cols) == len(s): #full model
135 partial_R_col_str = "-"
136 partial_R = "-"
137 try:
138 redr2 = "%.4f" % (float(redr2))
139 except:
140 pass
141 try:
142 partial_R = "%.4f" % (float(partial_R))
143 except:
144 pass
145 print >> fout, "%s\t%s\t%s\t%s" % ( col_str, redr2, partial_R_col_str, partial_R )