|
0
|
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 )
|
