root/galaxy-central/tools/multivariate_stats/kcca.py

リビジョン 2, 4.9 KB (コミッタ: hatakeyama, 14 年 前)

import galaxy-central

行番号 
1#!/usr/bin/env python
2
3"""
4Run kernel CCA using kcca() from R 'kernlab' package
5
6usage: %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
19usage: %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
22from galaxy import eggs
23import sys, string
24from rpy import *
25import numpy
26import pkg_resources; pkg_resources.require( "bx-python" )
27from bx.cookbook import doc_optparse
28
29
30def stop_err(msg):
31    sys.stderr.write(msg)
32    sys.exit()
33
34#Parse Command Line
35options, 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
38infile = options.input
39x_cols = options.x_cols.split(',')
40y_cols = options.y_cols.split(',')
41kernel = options.kernel
42outfile = options.output1
43ncomps = int(options.features)
44fout = open(outfile,'w')
45
46if 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()
49elems = []
50for 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
58if len( elems )<1:
59    stop_err( "The data in your input dataset is either missing or not formatted properly." )
60
61x_vals = []
62for k,col in enumerate(x_cols):
63    x_cols[k] = int(col)-1
64    x_vals.append([])
65y_vals = []
66for k,col in enumerate(y_cols):
67    y_cols[k] = int(col)-1
68    y_vals.append([])
69NA = 'NA'
70skipped = 0
71for 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
99x_vals1 = numpy.asarray(x_vals).transpose()
100y_vals1 = numpy.asarray(y_vals).transpose()
101
102x_dat= r.list(array(x_vals1))
103y_dat= r.list(array(y_vals1))
104
105try:
106    r.suppressWarnings(r.library('kernlab'))
107except:
108    stop_err('Missing R library kernlab')
109           
110set_default_mode(NO_CONVERSION)
111if kernel=="rbfdot" or kernel=="anovadot":
112    pars = r.list(sigma=float(options.sigma))
113elif kernel=="polydot":
114    pars = r.list(degree=float(options.degree),scale=float(options.scale),offset=float(options.offset))
115elif kernel=="tanhdot":
116    pars = r.list(scale=float(options.scale),offset=float(options.offset))
117elif kernel=="besseldot":
118    pars = r.list(degree=float(options.degree),sigma=float(options.sigma),order=float(options.order))
119elif kernel=="anovadot":
120    pars = r.list(degree=float(options.degree),sigma=float(options.sigma))
121else:
122    pars = rlist()
123   
124try:
125    kcc = r.kcca(x=x_dat, y=y_dat, kernel=kernel, kpar=pars, ncomps=ncomps)
126except RException, rex:
127    stop_err("Encountered error while performing kCCA on the input data: %s" %(rex))
128
129set_default_mode(BASIC_CONVERSION)   
130kcor = r.kcor(kcc)
131if ncomps == 1:
132    kcor = [kcor]
133xcoef = r.xcoef(kcc)
134ycoef = r.ycoef(kcc)
135
136print >>fout, "#Component\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
137
138print >>fout, "#Correlation\t%s" %("\t".join(["%.4g" % el for el in kcor]))
139   
140print >>fout, "#Estimated X-coefficients\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
141for obs,val in enumerate(xcoef):
142    print >>fout, "%s\t%s" %(obs+1, "\t".join(["%.4g" % el for el in val]))
143
144print >>fout, "#Estimated Y-coefficients\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
145for obs,val in enumerate(ycoef):
146    print >>fout, "%s\t%s" %(obs+1, "\t".join(["%.4g" % el for el in val]))
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。