Mercurial > repos > bgruening > upload_testing

comparison kcca.py @ 80:c4a3a8999945 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author bernhardlutz
date Mon, 20 Jan 2014 14:39:43 -0500
parents
children
comparison
equal deleted inserted replaced
79:dc82017052ac 80:c4a3a8999945
1 #!/usr/bin/env python
2
3 """
4 Run kernel CCA using kcca() from R 'kernlab' package
5
6 usage: %prog [options]
7 -i, --input=i: Input file
8 -o, --output1=o: Summary output
9 -x, --x_cols=x: X-Variable columns
10 -y, --y_cols=y: Y-Variable columns
11 -k, --kernel=k: Kernel function
12 -f, --features=f: Number of canonical components to return
13 -s, --sigma=s: sigma
14 -d, --degree=d: degree
15 -l, --scale=l: scale
16 -t, --offset=t: offset
17 -r, --order=r: order
18
19 usage: %prog input output1 x_cols y_cols kernel features sigma(or_None) degree(or_None) scale(or_None) offset(or_None) order(or_None)
20 """
21
22 from galaxy import eggs
23 import sys, string
24 #from rpy import *
25 import rpy2.robjects as robjects
26 import rpy2.rlike.container as rlc
27 from rpy2.robjects.packages import importr
28 r = robjects.r
29 import numpy
30 import pkg_resources; pkg_resources.require( "bx-python" )
31 from bx.cookbook import doc_optparse
32
33
34 def stop_err(msg):
35 sys.stderr.write(msg)
36 sys.exit()
37
38 #Parse Command Line
39 options, args = doc_optparse.parse( __doc__ )
40 #{'options= kernel': 'rbfdot', 'var_cols': '1,2,3,4', 'degree': 'None', 'output2': '/afs/bx.psu.edu/home/gua110/workspace/galaxy_bitbucket/database/files/000/dataset_260.dat', 'output1': '/afs/bx.psu.edu/home/gua110/workspace/galaxy_bitbucket/database/files/000/dataset_259.dat', 'scale': 'None', 'offset': 'None', 'input': '/afs/bx.psu.edu/home/gua110/workspace/galaxy_bitbucket/database/files/000/dataset_256.dat', 'sigma': '1.0', 'order': 'None'}
41
42 infile = options.input
43 x_cols = options.x_cols.split(',')
44 y_cols = options.y_cols.split(',')
45 kernel = options.kernel
46 outfile = options.output1
47 ncomps = int(options.features)
48 fout = open(outfile,'w')
49
50 if ncomps < 1:
51 print "You chose to return '0' canonical components. Please try rerunning the tool with number of components = 1 or more."
52 sys.exit()
53 elems = []
54 for i, line in enumerate( file ( infile )):
55 line = line.rstrip('\r\n')
56 if len( line )>0 and not line.startswith( '#' ):
57 elems = line.split( '\t' )
58 break
59 if i == 30:
60 break # Hopefully we'll never get here...
61
62 if len( elems )<1:
63 stop_err( "The data in your input dataset is either missing or not formatted properly." )
64
65 x_vals = []
66 for k,col in enumerate(x_cols):
67 x_cols[k] = int(col)-1
68 #x_vals.append([])
69 y_vals = []
70 for k,col in enumerate(y_cols):
71 y_cols[k] = int(col)-1
72 #y_vals.append([])
73 NA = 'NA'
74 skipped = 0
75 for ind,line in enumerate( file( infile )):
76 if line and not line.startswith( '#' ):
77 try:
78 fields = line.strip().split("\t")
79 valid_line = True
80 for col in x_cols+y_cols:
81 try:
82 assert float(fields[col])
83 except:
84 skipped += 1
85 valid_line = False
86 break
87 if valid_line:
88 for k,col in enumerate(x_cols):
89 try:
90 xval = float(fields[col])
91 except:
92 xval = NaN#
93 #x_vals[k].append(xval)
94 x_vals.append(xval)
95 for k,col in enumerate(y_cols):
96 try:
97 yval = float(fields[col])
98 except:
99 yval = NaN#
100 #y_vals[k].append(yval)
101 y_vals.append(yval)
102 except:
103 skipped += 1
104
105 #x_vals1 = numpy.asarray(x_vals).transpose()
106 #y_vals1 = numpy.asarray(y_vals).transpose()
107
108 #x_dat= r.list(array(x_vals1))
109 #y_dat= r.list(array(y_vals1))
110
111 x_dat = r['matrix'](robjects.FloatVector(x_vals),ncol=len(x_cols),byrow=True)
112 y_dat = r['matrix'](robjects.FloatVector(y_vals),ncol=len(y_cols),byrow=True)
113
114 try:
115 r.suppressWarnings(r.library('kernlab'))
116 except:
117 stop_err('Missing R library kernlab')
118
119 #set_default_mode(NO_CONVERSION)
120 if kernel=="rbfdot" or kernel=="anovadot":
121 pars = r.list(sigma=float(options.sigma))
122 elif kernel=="polydot":
123 pars = r.list(degree=float(options.degree),scale=float(options.scale),offset=float(options.offset))
124 elif kernel=="tanhdot":
125 pars = r.list(scale=float(options.scale),offset=float(options.offset))
126 elif kernel=="besseldot":
127 pars = r.list(degree=float(options.degree),sigma=float(options.sigma),order=float(options.order))
128 elif kernel=="anovadot":
129 pars = r.list(degree=float(options.degree),sigma=float(options.sigma))
130 else:
131 pars = rlist()
132
133 try:
134 kcc = r.kcca(x=x_dat, y=y_dat, kernel=kernel, kpar=pars, ncomps=ncomps)
135 except RException, rex:
136 stop_err("Encountered error while performing kCCA on the input data: %s" %(rex))
137
138 #set_default_mode(BASIC_CONVERSION)
139 kcor = r.kcor(kcc)
140 if ncomps == 1:
141 kcor = [kcor]
142 xcoef = r.xcoef(kcc)
143 ycoef = r.ycoef(kcc)
144
145 print >>fout, "#Component\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
146
147 print >>fout, "#Correlation\t%s" %("\t".join(["%.4g" % el for el in kcor]))
148
149 print >>fout, "#Estimated X-coefficients\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
150 #for obs,val in enumerate(xcoef):
151 # print >>fout, "%s\t%s" %(obs+1, "\t".join(["%.4g" % el for el in val]))
152 for i in range(1,xcoef.nrow+1):
153 vals = []
154 for j in range(1,xcoef.ncol+1):
155 vals.append("%.4g" % xcoef.rx2(i,j)[0])
156 print >>fout, "%s\t%s" %(i, "\t".join(vals))
157
158
159 print >>fout, "#Estimated Y-coefficients\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
160 #for obs,val in enumerate(ycoef):
161 # print >>fout, "%s\t%s" %(obs+1, "\t".join(["%.4g" % el for el in val]))
162 for i in range(1,ycoef.nrow+1):
163 vals = []
164 for j in range(1,ycoef.ncol+1):
165 vals.append("%.4g" % ycoef.rx2(i,j)[0])
166 print >>fout, "%s\t%s" %(i, "\t".join(vals))