|
0
|
1 #!/usr/bin/env python
|
|
|
2
|
|
|
3 from galaxy import eggs
|
|
|
4 import sys, string
|
|
|
5 from rpy import *
|
|
|
6 import numpy
|
|
|
7
|
|
|
8 def stop_err(msg):
|
|
|
9 sys.stderr.write(msg)
|
|
|
10 sys.exit()
|
|
|
11
|
|
|
12 infile = sys.argv[1]
|
|
|
13 y_col = int(sys.argv[2])-1
|
|
|
14 x_cols = sys.argv[3].split(',')
|
|
|
15 outfile = sys.argv[4]
|
|
|
16 outfile2 = sys.argv[5]
|
|
|
17
|
|
|
18 print "Predictor columns: %s; Response column: %d" %(x_cols,y_col+1)
|
|
|
19 fout = open(outfile,'w')
|
|
|
20 elems = []
|
|
|
21 for i, line in enumerate( file ( infile )):
|
|
|
22 line = line.rstrip('\r\n')
|
|
|
23 if len( line )>0 and not line.startswith( '#' ):
|
|
|
24 elems = line.split( '\t' )
|
|
|
25 break
|
|
|
26 if i == 30:
|
|
|
27 break # Hopefully we'll never get here...
|
|
|
28
|
|
|
29 if len( elems )<1:
|
|
|
30 stop_err( "The data in your input dataset is either missing or not formatted properly." )
|
|
|
31
|
|
|
32 y_vals = []
|
|
|
33 x_vals = []
|
|
|
34
|
|
|
35 for k,col in enumerate(x_cols):
|
|
|
36 x_cols[k] = int(col)-1
|
|
|
37 x_vals.append([])
|
|
|
38
|
|
|
39 NA = 'NA'
|
|
|
40 for ind,line in enumerate( file( infile )):
|
|
|
41 if line and not line.startswith( '#' ):
|
|
|
42 try:
|
|
|
43 fields = line.split("\t")
|
|
|
44 try:
|
|
|
45 yval = float(fields[y_col])
|
|
|
46 except:
|
|
|
47 yval = r('NA')
|
|
|
48 y_vals.append(yval)
|
|
|
49 for k,col in enumerate(x_cols):
|
|
|
50 try:
|
|
|
51 xval = float(fields[col])
|
|
|
52 except:
|
|
|
53 xval = r('NA')
|
|
|
54 x_vals[k].append(xval)
|
|
|
55 except:
|
|
|
56 pass
|
|
|
57
|
|
|
58 x_vals1 = numpy.asarray(x_vals).transpose()
|
|
|
59
|
|
|
60 dat= r.list(x=array(x_vals1), y=y_vals)
|
|
|
61
|
|
|
62 set_default_mode(NO_CONVERSION)
|
|
|
63 try:
|
|
|
64 linear_model = r.lm(r("y ~ x"), data = r.na_exclude(dat))
|
|
|
65 except RException, rex:
|
|
|
66 stop_err("Error performing linear regression on the input data.\nEither the response column or one of the predictor columns contain only non-numeric or invalid values.")
|
|
|
67 set_default_mode(BASIC_CONVERSION)
|
|
|
68
|
|
|
69 coeffs=linear_model.as_py()['coefficients']
|
|
|
70 yintercept= coeffs['(Intercept)']
|
|
|
71 summary = r.summary(linear_model)
|
|
|
72
|
|
|
73 co = summary.get('coefficients', 'NA')
|
|
|
74 """
|
|
|
75 if len(co) != len(x_vals)+1:
|
|
|
76 stop_err("Stopped performing linear regression on the input data, since one of the predictor columns contains only non-numeric or invalid values.")
|
|
|
77 """
|
|
|
78
|
|
|
79 try:
|
|
|
80 yintercept = r.round(float(yintercept), digits=10)
|
|
|
81 pvaly = r.round(float(co[0][3]), digits=10)
|
|
|
82 except:
|
|
|
83 pass
|
|
|
84
|
|
|
85 print >>fout, "Y-intercept\t%s" %(yintercept)
|
|
|
86 print >>fout, "p-value (Y-intercept)\t%s" %(pvaly)
|
|
|
87
|
|
|
88 if len(x_vals) == 1: #Simple linear regression case with 1 predictor variable
|
|
|
89 try:
|
|
|
90 slope = r.round(float(coeffs['x']), digits=10)
|
|
|
91 except:
|
|
|
92 slope = 'NA'
|
|
|
93 try:
|
|
|
94 pval = r.round(float(co[1][3]), digits=10)
|
|
|
95 except:
|
|
|
96 pval = 'NA'
|
|
|
97 print >>fout, "Slope (c%d)\t%s" %(x_cols[0]+1,slope)
|
|
|
98 print >>fout, "p-value (c%d)\t%s" %(x_cols[0]+1,pval)
|
|
|
99 else: #Multiple regression case with >1 predictors
|
|
|
100 ind=1
|
|
|
101 while ind < len(coeffs.keys()):
|
|
|
102 try:
|
|
|
103 slope = r.round(float(coeffs['x'+str(ind)]), digits=10)
|
|
|
104 except:
|
|
|
105 slope = 'NA'
|
|
|
106 print >>fout, "Slope (c%d)\t%s" %(x_cols[ind-1]+1,slope)
|
|
|
107 try:
|
|
|
108 pval = r.round(float(co[ind][3]), digits=10)
|
|
|
109 except:
|
|
|
110 pval = 'NA'
|
|
|
111 print >>fout, "p-value (c%d)\t%s" %(x_cols[ind-1]+1,pval)
|
|
|
112 ind+=1
|
|
|
113
|
|
|
114 rsq = summary.get('r.squared','NA')
|
|
|
115 adjrsq = summary.get('adj.r.squared','NA')
|
|
|
116 fstat = summary.get('fstatistic','NA')
|
|
|
117 sigma = summary.get('sigma','NA')
|
|
|
118
|
|
|
119 try:
|
|
|
120 rsq = r.round(float(rsq), digits=5)
|
|
|
121 adjrsq = r.round(float(adjrsq), digits=5)
|
|
|
122 fval = r.round(fstat['value'], digits=5)
|
|
|
123 fstat['value'] = str(fval)
|
|
|
124 sigma = r.round(float(sigma), digits=10)
|
|
|
125 except:
|
|
|
126 pass
|
|
|
127
|
|
|
128 print >>fout, "R-squared\t%s" %(rsq)
|
|
|
129 print >>fout, "Adjusted R-squared\t%s" %(adjrsq)
|
|
|
130 print >>fout, "F-statistic\t%s" %(fstat)
|
|
|
131 print >>fout, "Sigma\t%s" %(sigma)
|
|
|
132
|
|
|
133 r.pdf( outfile2, 8, 8 )
|
|
|
134 if len(x_vals) == 1: #Simple linear regression case with 1 predictor variable
|
|
|
135 sub_title = "Slope = %s; Y-int = %s" %(slope,yintercept)
|
|
|
136 try:
|
|
|
137 r.plot(x=x_vals[0], y=y_vals, xlab="X", ylab="Y", sub=sub_title, main="Scatterplot with regression")
|
|
|
138 r.abline(a=yintercept, b=slope, col="red")
|
|
|
139 except:
|
|
|
140 pass
|
|
|
141 else:
|
|
|
142 r.pairs(dat, main="Scatterplot Matrix", col="blue")
|
|
|
143 try:
|
|
|
144 r.plot(linear_model)
|
|
|
145 except:
|
|
|
146 pass
|
|
|
147 r.dev_off()
|
