PageRenderTime 28ms CodeModel.GetById 13ms app.highlight 12ms RepoModel.GetById 1ms app.codeStats 0ms

/tools/regVariation/best_regression_subsets.py

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