0
+ − 1 #!/usr/bin/env python
+ − 2
+ − 3 import sys
+ − 4 from rpy import *
+ − 5 import numpy
+ − 6
+ − 7 #export PYTHONPATH=~/galaxy/lib/
+ − 8 #running command python partialR_square.py reg_inp.tab 4 1,2,3 partialR_result.tabular
+ − 9
+ − 10 def stop_err(msg):
+ − 11 sys.stderr.write(msg)
+ − 12 sys.exit()
+ − 13
+ − 14
+ − 15 def sscombs(s):
+ − 16 if len(s) == 1:
+ − 17 return [s]
+ − 18 else:
+ − 19 ssc = sscombs(s[1:])
+ − 20 return [s[0]] + [s[0]+comb for comb in ssc] + ssc
+ − 21
+ − 22
+ − 23 infile = sys.argv[1]
+ − 24 y_col = int(sys.argv[2])-1
+ − 25 x_cols = sys.argv[3].split(',')
+ − 26 outfile = sys.argv[4]
+ − 27
+ − 28 print "Predictor columns: %s; Response column: %d" % ( x_cols, y_col+1 )
+ − 29 fout = open(outfile,'w')
+ − 30
+ − 31 for i, line in enumerate( file ( infile )):
+ − 32 line = line.rstrip('\r\n')
+ − 33 if len( line )>0 and not line.startswith( '#' ):
+ − 34 elems = line.split( '\t' )
+ − 35 break
+ − 36 if i == 30:
+ − 37 break # Hopefully we'll never get here...
+ − 38
+ − 39 if len( elems )<1:
+ − 40 stop_err( "The data in your input dataset is either missing or not formatted properly." )
+ − 41
+ − 42 y_vals = []
+ − 43 x_vals = []
+ − 44
+ − 45 for k, col in enumerate(x_cols):
+ − 46 x_cols[k] = int(col)-1
+ − 47 x_vals.append([])
+ − 48 """
+ − 49 try:
+ − 50 float( elems[x_cols[k]] )
+ − 51 except:
+ − 52 try:
+ − 53 msg = "This operation cannot be performed on non-numeric column %d containing value '%s'." % ( col, elems[x_cols[k]] )
+ − 54 except:
+ − 55 msg = "This operation cannot be performed on non-numeric data."
+ − 56 stop_err( msg )
+ − 57 """
+ − 58 NA = 'NA'
+ − 59 for ind, line in enumerate( file( infile )):
+ − 60 if line and not line.startswith( '#' ):
+ − 61 try:
+ − 62 fields = line.split("\t")
+ − 63 try:
+ − 64 yval = float(fields[y_col])
+ − 65 except Exception, ey:
+ − 66 yval = r('NA')
+ − 67 #print >> sys.stderr, "ey = %s" %ey
+ − 68 y_vals.append(yval)
+ − 69 for k, col in enumerate(x_cols):
+ − 70 try:
+ − 71 xval = float(fields[col])
+ − 72 except Exception, ex:
+ − 73 xval = r('NA')
+ − 74 #print >> sys.stderr, "ex = %s" %ex
+ − 75 x_vals[k].append(xval)
+ − 76 except:
+ − 77 pass
+ − 78
+ − 79 x_vals1 = numpy.asarray(x_vals).transpose()
+ − 80 dat = r.list(x=array(x_vals1), y=y_vals)
+ − 81
+ − 82 set_default_mode(NO_CONVERSION)
+ − 83 try:
+ − 84 full = r.lm(r("y ~ x"), data= r.na_exclude(dat)) #full model includes all the predictor variables specified by the user
+ − 85 except RException, rex:
+ − 86 stop_err("Error performing linear regression on the input data.\nEither the response column or one of the predictor columns contain no numeric values.")
+ − 87 set_default_mode(BASIC_CONVERSION)
+ − 88
+ − 89 summary = r.summary(full)
+ − 90 fullr2 = summary.get('r.squared','NA')
+ − 91
+ − 92 if fullr2 == 'NA':
+ − 93 stop_err("Error in linear regression")
+ − 94
+ − 95 if len(x_vals) < 10:
+ − 96 s = ""
+ − 97 for ch in range(len(x_vals)):
+ − 98 s += str(ch)
+ − 99 else:
+ − 100 stop_err("This tool only works with less than 10 predictors.")
+ − 101
+ − 102 print >> fout, "#Model\tR-sq\tpartial_R_Terms\tpartial_R_Value"
+ − 103 all_combos = sorted(sscombs(s), key=len)
+ − 104 all_combos.reverse()
+ − 105 for j, cols in enumerate(all_combos):
+ − 106 #if len(cols) == len(s): #Same as the full model above
+ − 107 # continue
+ − 108 if len(cols) == 1:
+ − 109 x_vals1 = x_vals[int(cols)]
+ − 110 else:
+ − 111 x_v = []
+ − 112 for col in cols:
+ − 113 x_v.append(x_vals[int(col)])
+ − 114 x_vals1 = numpy.asarray(x_v).transpose()
+ − 115 dat = r.list(x=array(x_vals1), y=y_vals)
+ − 116 set_default_mode(NO_CONVERSION)
+ − 117 red = r.lm(r("y ~ x"), data= dat) #Reduced model
+ − 118 set_default_mode(BASIC_CONVERSION)
+ − 119 summary = r.summary(red)
+ − 120 redr2 = summary.get('r.squared','NA')
+ − 121 try:
+ − 122 partial_R = (float(fullr2)-float(redr2))/(1-float(redr2))
+ − 123 except:
+ − 124 partial_R = 'NA'
+ − 125 col_str = ""
+ − 126 for col in cols:
+ − 127 col_str = col_str + str(int(x_cols[int(col)]) + 1) + " "
+ − 128 col_str.strip()
+ − 129 partial_R_col_str = ""
+ − 130 for col in s:
+ − 131 if col not in cols:
+ − 132 partial_R_col_str = partial_R_col_str + str(int(x_cols[int(col)]) + 1) + " "
+ − 133 partial_R_col_str.strip()
+ − 134 if len(cols) == len(s): #full model
+ − 135 partial_R_col_str = "-"
+ − 136 partial_R = "-"
+ − 137 try:
+ − 138 redr2 = "%.4f" % (float(redr2))
+ − 139 except:
+ − 140 pass
+ − 141 try:
+ − 142 partial_R = "%.4f" % (float(partial_R))
+ − 143 except:
+ − 144 pass
+ − 145 print >> fout, "%s\t%s\t%s\t%s" % ( col_str, redr2, partial_R_col_str, partial_R )
