|
0
|
1 #!/usr/bin/env python
|
|
|
2
|
|
|
3 from galaxy import eggs
|
|
|
4
|
|
|
5 import sys
|
|
|
6 from rpy import *
|
|
|
7 import numpy
|
|
|
8
|
|
|
9 def stop_err(msg):
|
|
|
10 sys.stderr.write(msg)
|
|
|
11 sys.exit()
|
|
|
12
|
|
|
13
|
|
|
14 infile = sys.argv[1]
|
|
|
15 y_col = int(sys.argv[2])-1
|
|
|
16 x_cols = sys.argv[3].split(',')
|
|
|
17 outfile = sys.argv[4]
|
|
|
18 outfile2 = sys.argv[5]
|
|
|
19 print "Predictor columns: %s; Response column: %d" % ( x_cols, y_col+1 )
|
|
|
20 fout = open(outfile,'w')
|
|
|
21
|
|
|
22 for i, line in enumerate( file ( infile )):
|
|
|
23 line = line.rstrip('\r\n')
|
|
|
24 if len( line )>0 and not line.startswith( '#' ):
|
|
|
25 elems = line.split( '\t' )
|
|
|
26 break
|
|
|
27 if i == 30:
|
|
|
28 break # Hopefully we'll never get here...
|
|
|
29
|
|
|
30 if len( elems )<1:
|
|
|
31 stop_err( "The data in your input dataset is either missing or not formatted properly." )
|
|
|
32
|
|
|
33 y_vals = []
|
|
|
34 x_vals = []
|
|
|
35
|
|
|
36 for k, col in enumerate(x_cols):
|
|
|
37 x_cols[k] = int(col)-1
|
|
|
38 x_vals.append([])
|
|
|
39
|
|
|
40 NA = 'NA'
|
|
|
41 for ind, line in enumerate( file( infile ) ):
|
|
|
42 if line and not line.startswith( '#' ):
|
|
|
43 try:
|
|
|
44 fields = line.split("\t")
|
|
|
45 try:
|
|
|
46 yval = float(fields[y_col])
|
|
|
47 except Exception, ey:
|
|
|
48 yval = r('NA')
|
|
|
49 y_vals.append(yval)
|
|
|
50 for k, col in enumerate(x_cols):
|
|
|
51 try:
|
|
|
52 xval = float(fields[col])
|
|
|
53 except Exception, ex:
|
|
|
54 xval = r('NA')
|
|
|
55 x_vals[k].append(xval)
|
|
|
56 except:
|
|
|
57 pass
|
|
|
58
|
|
|
59 response_term = ""
|
|
|
60
|
|
|
61 x_vals1 = numpy.asarray(x_vals).transpose()
|
|
|
62
|
|
|
63 dat = r.list(x=array(x_vals1), y=y_vals)
|
|
|
64
|
|
|
65 r.library("leaps")
|
|
|
66
|
|
|
67 set_default_mode(NO_CONVERSION)
|
|
|
68 try:
|
|
|
69 leaps = r.regsubsets(r("y ~ x"), data= r.na_exclude(dat))
|
|
|
70 except RException, rex:
|
|
|
71 stop_err("Error performing linear regression on the input data.\nEither the response column or one of the predictor columns contain no numeric values.")
|
|
|
72 set_default_mode(BASIC_CONVERSION)
|
|
|
73
|
|
|
74 summary = r.summary(leaps)
|
|
|
75 tot = len(x_vals)
|
|
|
76 pattern = "["
|
|
|
77 for i in range(tot):
|
|
|
78 pattern = pattern + 'c' + str(int(x_cols[int(i)]) + 1) + ' '
|
|
|
79 pattern = pattern.strip() + ']'
|
|
|
80 print >> fout, "#Vars\t%s\tR-sq\tAdj. R-sq\tC-p\tbic" % (pattern)
|
|
|
81 for ind, item in enumerate(summary['outmat']):
|
|
|
82 print >> fout, "%s\t%s\t%s\t%s\t%s\t%s" % (str(item).count('*'), item, summary['rsq'][ind], summary['adjr2'][ind], summary['cp'][ind], summary['bic'][ind])
|
|
|
83
|
|
|
84
|
|
|
85 r.pdf( outfile2, 8, 8 )
|
|
|
86 r.plot(leaps, scale="Cp", main="Best subsets using Cp Criterion")
|
|
|
87 r.plot(leaps, scale="r2", main="Best subsets using R-sq Criterion")
|
|
|
88 r.plot(leaps, scale="adjr2", main="Best subsets using Adjusted R-sq Criterion")
|
|
|
89 r.plot(leaps, scale="bic", main="Best subsets using bic Criterion")
|
|
|
90
|
|
|
91 r.dev_off()
|
