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