Mercurial > repos > xuebing > sharplabtool

comparison tools/regVariation/linear_regression.py @ 0:9071e359b9a3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author xuebing
date Fri, 09 Mar 2012 19:37:19 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:9071e359b9a3
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()