|
0
|
1 #!/usr/bin/env python
|
|
|
2
|
|
|
3 import sys, string
|
|
|
4 from rpy import *
|
|
|
5 import numpy
|
|
|
6
|
|
|
7 def stop_err(msg):
|
|
|
8 sys.stderr.write(msg)
|
|
|
9 sys.exit()
|
|
|
10
|
|
|
11 infile = sys.argv[1]
|
|
|
12 x_cols = sys.argv[2].split(',')
|
|
|
13 y_cols = sys.argv[3].split(',')
|
|
|
14
|
|
|
15 x_scale = x_center = "FALSE"
|
|
|
16 if sys.argv[4] == 'both':
|
|
|
17 x_scale = x_center = "TRUE"
|
|
|
18 elif sys.argv[4] == 'center':
|
|
|
19 x_center = "TRUE"
|
|
|
20 elif sys.argv[4] == 'scale':
|
|
|
21 x_scale = "TRUE"
|
|
|
22
|
|
|
23 y_scale = y_center = "FALSE"
|
|
|
24 if sys.argv[5] == 'both':
|
|
|
25 y_scale = y_center = "TRUE"
|
|
|
26 elif sys.argv[5] == 'center':
|
|
|
27 y_center = "TRUE"
|
|
|
28 elif sys.argv[5] == 'scale':
|
|
|
29 y_scale = "TRUE"
|
|
|
30
|
|
|
31 std_scores = "FALSE"
|
|
|
32 if sys.argv[6] == "yes":
|
|
|
33 std_scores = "TRUE"
|
|
|
34
|
|
|
35 outfile = sys.argv[7]
|
|
|
36 outfile2 = sys.argv[8]
|
|
|
37
|
|
|
38 fout = open(outfile,'w')
|
|
|
39 elems = []
|
|
|
40 for i, line in enumerate( file ( infile )):
|
|
|
41 line = line.rstrip('\r\n')
|
|
|
42 if len( line )>0 and not line.startswith( '#' ):
|
|
|
43 elems = line.split( '\t' )
|
|
|
44 break
|
|
|
45 if i == 30:
|
|
|
46 break # Hopefully we'll never get here...
|
|
|
47
|
|
|
48 if len( elems )<1:
|
|
|
49 stop_err( "The data in your input dataset is either missing or not formatted properly." )
|
|
|
50
|
|
|
51 x_vals = []
|
|
|
52
|
|
|
53 for k,col in enumerate(x_cols):
|
|
|
54 x_cols[k] = int(col)-1
|
|
|
55 x_vals.append([])
|
|
|
56
|
|
|
57 y_vals = []
|
|
|
58
|
|
|
59 for k,col in enumerate(y_cols):
|
|
|
60 y_cols[k] = int(col)-1
|
|
|
61 y_vals.append([])
|
|
|
62
|
|
|
63 skipped = 0
|
|
|
64 for ind,line in enumerate( file( infile )):
|
|
|
65 if line and not line.startswith( '#' ):
|
|
|
66 try:
|
|
|
67 fields = line.strip().split("\t")
|
|
|
68 valid_line = True
|
|
|
69 for col in x_cols+y_cols:
|
|
|
70 try:
|
|
|
71 assert float(fields[col])
|
|
|
72 except:
|
|
|
73 skipped += 1
|
|
|
74 valid_line = False
|
|
|
75 break
|
|
|
76 if valid_line:
|
|
|
77 for k,col in enumerate(x_cols):
|
|
|
78 try:
|
|
|
79 xval = float(fields[col])
|
|
|
80 except:
|
|
|
81 xval = NaN#
|
|
|
82 x_vals[k].append(xval)
|
|
|
83 for k,col in enumerate(y_cols):
|
|
|
84 try:
|
|
|
85 yval = float(fields[col])
|
|
|
86 except:
|
|
|
87 yval = NaN#
|
|
|
88 y_vals[k].append(yval)
|
|
|
89 except:
|
|
|
90 skipped += 1
|
|
|
91
|
|
|
92 x_vals1 = numpy.asarray(x_vals).transpose()
|
|
|
93 y_vals1 = numpy.asarray(y_vals).transpose()
|
|
|
94
|
|
|
95 x_dat= r.list(array(x_vals1))
|
|
|
96 y_dat= r.list(array(y_vals1))
|
|
|
97
|
|
|
98 try:
|
|
|
99 r.suppressWarnings(r.library("yacca"))
|
|
|
100 except:
|
|
|
101 stop_err("Missing R library yacca.")
|
|
|
102
|
|
|
103 set_default_mode(NO_CONVERSION)
|
|
|
104 try:
|
|
|
105 xcolnames = ["c%d" %(el+1) for el in x_cols]
|
|
|
106 ycolnames = ["c%d" %(el+1) for el in y_cols]
|
|
|
107 cc = r.cca(x=x_dat, y=y_dat, xlab=xcolnames, ylab=ycolnames, xcenter=r(x_center), ycenter=r(y_center), xscale=r(x_scale), yscale=r(y_scale), standardize_scores=r(std_scores))
|
|
|
108 ftest = r.F_test_cca(cc)
|
|
|
109 except RException, rex:
|
|
|
110 stop_err("Encountered error while performing CCA on the input data: %s" %(rex))
|
|
|
111
|
|
|
112 set_default_mode(BASIC_CONVERSION)
|
|
|
113 summary = r.summary(cc)
|
|
|
114
|
|
|
115 ncomps = len(summary['corr'])
|
|
|
116 comps = summary['corr'].keys()
|
|
|
117 corr = summary['corr'].values()
|
|
|
118 xlab = summary['xlab']
|
|
|
119 ylab = summary['ylab']
|
|
|
120
|
|
|
121 for i in range(ncomps):
|
|
|
122 corr[comps.index('CV %s' %(i+1))] = summary['corr'].values()[i]
|
|
|
123
|
|
|
124 ftest=ftest.as_py()
|
|
|
125 print >>fout, "#Component\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
|
|
|
126 print >>fout, "#Correlation\t%s" %("\t".join(["%.4g" % el for el in corr]))
|
|
|
127 print >>fout, "#F-statistic\t%s" %("\t".join(["%.4g" % el for el in ftest['statistic']]))
|
|
|
128 print >>fout, "#p-value\t%s" %("\t".join(["%.4g" % el for el in ftest['p.value']]))
|
|
|
129
|
|
|
130 print >>fout, "#X-Coefficients\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
|
|
|
131 for i,val in enumerate(summary['xcoef']):
|
|
|
132 print >>fout, "%s\t%s" %(xlab[i], "\t".join(["%.4g" % el for el in val]))
|
|
|
133
|
|
|
134 print >>fout, "#Y-Coefficients\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
|
|
|
135 for i,val in enumerate(summary['ycoef']):
|
|
|
136 print >>fout, "%s\t%s" %(ylab[i], "\t".join(["%.4g" % el for el in val]))
|
|
|
137
|
|
|
138 print >>fout, "#X-Loadings\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
|
|
|
139 for i,val in enumerate(summary['xstructcorr']):
|
|
|
140 print >>fout, "%s\t%s" %(xlab[i], "\t".join(["%.4g" % el for el in val]))
|
|
|
141
|
|
|
142 print >>fout, "#Y-Loadings\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
|
|
|
143 for i,val in enumerate(summary['ystructcorr']):
|
|
|
144 print >>fout, "%s\t%s" %(ylab[i], "\t".join(["%.4g" % el for el in val]))
|
|
|
145
|
|
|
146 print >>fout, "#X-CrossLoadings\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
|
|
|
147 for i,val in enumerate(summary['xcrosscorr']):
|
|
|
148 print >>fout, "%s\t%s" %(xlab[i], "\t".join(["%.4g" % el for el in val]))
|
|
|
149
|
|
|
150 print >>fout, "#Y-CrossLoadings\t%s" %("\t".join(["%s" % el for el in range(1,ncomps+1)]))
|
|
|
151 for i,val in enumerate(summary['ycrosscorr']):
|
|
|
152 print >>fout, "%s\t%s" %(ylab[i], "\t".join(["%.4g" % el for el in val]))
|
|
|
153
|
|
|
154 r.pdf( outfile2, 8, 8 )
|
|
|
155 #r.plot(cc)
|
|
|
156 for i in range(ncomps):
|
|
|
157 r.helio_plot(cc, cv = i+1, main = r.paste("Explained Variance for CV",i+1), type = "variance")
|
|
|
158 r.dev_off() |
