|
0
|
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 numpy
|
|
|
26 import pkg_resources; pkg_resources.require( "bx-python" )
|
|
|
27 from bx.cookbook import doc_optparse
|
|
|
28
|
|
|
29
|
|
|
30 def stop_err(msg):
|
|
|
31 sys.stderr.write(msg)
|
|
|
32 sys.exit()
|
|
|
33
|
|
|
34 #Parse Command Line
|
|
|
35 options, args = doc_optparse.parse( __doc__ )
|
|
|
36 #{'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'}
|
|
|
37
|
|
|
38 infile = options.input
|
|
|
39 x_cols = options.x_cols.split(',')
|
|
|
40 y_cols = options.y_cols.split(',')
|
|
|
41 kernel = options.kernel
|
|
|
42 outfile = options.output1
|
|
|
43 ncomps = int(options.features)
|
|
|
44 fout = open(outfile,'w')
|
|
|
45
|
|
|
46 if ncomps < 1:
|
|
|
47 print "You chose to return '0' canonical components. Please try rerunning the tool with number of components = 1 or more."
|
|
|
48 sys.exit()
|
|
|
49 elems = []
|
|
|
50 for i, line in enumerate( file ( infile )):
|
|
|
51 line = line.rstrip('\r\n')
|
|
|
52 if len( line )>0 and not line.startswith( '#' ):
|
|
|
53 elems = line.split( '\t' )
|
|
|
54 break
|
|
|
55 if i == 30:
|
|
|
56 break # Hopefully we'll never get here...
|
|
|
57
|
|
|
58 if len( elems )<1:
|
|
|
59 stop_err( "The data in your input dataset is either missing or not formatted properly." )
|
|
|
60
|
|
|
61 x_vals = []
|
|
|
62 for k,col in enumerate(x_cols):
|
|
|
63 x_cols[k] = int(col)-1
|
|
|
64 x_vals.append([])
|
|
|
65 y_vals = []
|
|
|
66 for k,col in enumerate(y_cols):
|
|
|
67 y_cols[k] = int(col)-1
|
|
|
68 y_vals.append([])
|
|
|
69 NA = 'NA'
|
|
|
70 skipped = 0
|
|
|
71 for ind,line in enumerate( file( infile )):
|
|
|
72 if line and not line.startswith( '#' ):
|
|
|
73 try:
|
|
|
74 fields = line.strip().split("\t")
|
|
|
75 valid_line = True
|
|
|
76 for col in x_cols+y_cols:
|
|
|
77 try:
|
|
|
78 assert float(fields[col])
|
|
|
79 except:
|
|
|
80 skipped += 1
|
|
|
81 valid_line = False
|
|
|
82 break
|
|
|
83 if valid_line:
|
|
|
84 for k,col in enumerate(x_cols):
|
|
|
85 try:
|
|
|
86 xval = float(fields[col])
|
|
|
87 except:
|
|
|
88 xval = NaN#
|
|
|
89 x_vals[k].append(xval)
|
|
|
90 for k,col in enumerate(y_cols):
|
|
|
91 try:
|
|
|
92 yval = float(fields[col])
|
|
|
93 except:
|
|
|
94 yval = NaN#
|
|
|
95 y_vals[k].append(yval)
|
|
|
96 except:
|
|
|
97 skipped += 1
|
|
|
98
|
|
|
99 x_vals1 = numpy.asarray(x_vals).transpose()
|
|
|
100 y_vals1 = numpy.asarray(y_vals).transpose()
|
|
|
101
|
|
|
102 x_dat= r.list(array(x_vals1))
|
|
|
103 y_dat= r.list(array(y_vals1))
|
|
|
104
|
|
|
105 try:
|
|
|
106 r.suppressWarnings(r.library('kernlab'))
|
|
|
107 except:
|
|
|
108 stop_err('Missing R library kernlab')
|
|
|
109
|
|
|
110 set_default_mode(NO_CONVERSION)
|
|
|
111 if kernel=="rbfdot" or kernel=="anovadot":
|
|
|
112 pars = r.list(sigma=float(options.sigma))
|
|
|
113 elif kernel=="polydot":
|
|
|
114 pars = r.list(degree=float(options.degree),scale=float(options.scale),offset=float(options.offset))
|
|
|
115 elif kernel=="tanhdot":
|
|
|
116 pars = r.list(scale=float(options.scale),offset=float(options.offset))
|
|
|
117 elif kernel=="besseldot":
|
|
|
118 pars = r.list(degree=float(options.degree),sigma=float(options.sigma),order=float(options.order))
|
|
|
119 elif kernel=="anovadot":
|
|
|
120 pars = r.list(degree=float(options.degree),sigma=float(options.sigma))
|
|
|
121 else:
|
|
|
122 pars = rlist()
|
|
|
123
|
|
|
124 try:
|
|
|
125 kcc = r.kcca(x=x_dat, y=y_dat, kernel=kernel, kpar=pars, ncomps=ncomps)
|
|
|
126 except RException, rex:
|
|
|
127 stop_err("Encountered error while performing kCCA on the input data: %s" %(rex))
|
|
|
128
|
|
|
129 set_default_mode(BASIC_CONVERSION)
|
|
|
130 kcor = r.kcor(kcc)
|
|
|
131 if ncomps == 1:
|
|
|
132 kcor = [kcor]
|
|
|
133 xcoef = r.xcoef(kcc)
|
|
|
134 ycoef = r.ycoef(kcc)
|
|
|
135
|
|
|
136 print >>fout, "#Component\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
|
|
|
137
|
|
|
138 print >>fout, "#Correlation\t%s" %("\t".join(["%.4g" % el for el in kcor]))
|
|
|
139
|
|
|
140 print >>fout, "#Estimated X-coefficients\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
|
|
|
141 for obs,val in enumerate(xcoef):
|
|
|
142 print >>fout, "%s\t%s" %(obs+1, "\t".join(["%.4g" % el for el in val]))
|
|
|
143
|
|
|
144 print >>fout, "#Estimated Y-coefficients\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
|
|
|
145 for obs,val in enumerate(ycoef):
|
|
|
146 print >>fout, "%s\t%s" %(obs+1, "\t".join(["%.4g" % el for el in val]))
|
