[2] | 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])) |
---|