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