Mercurial > repos > devteam > kernel_canonical_correlation_analysis

comparison kcca.py @ 0:7a092113eb8c draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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]))