Mercurial > repos > devteam > correlation

comparison cor.py @ 0:24e01abf9e34 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Imported from capsule None
author devteam
date Mon, 28 Jul 2014 11:55:23 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:24e01abf9e34
1 #!/usr/bin/env python
2 #Greg Von Kuster
3 """
4 Calculate correlations between numeric columns in a tab delim file.
5 usage: %prog infile output.txt columns method
6 """
7
8 import sys
9 from rpy import *
10
11 def stop_err(msg):
12 sys.stderr.write(msg)
13 sys.exit()
14
15 def main():
16 method = sys.argv[4]
17 assert method in ( "pearson", "kendall", "spearman" )
18
19 try:
20 columns = map( int, sys.argv[3].split( ',' ) )
21 except:
22 stop_err( "Problem determining columns, perhaps your query does not contain a column of numerical data." )
23
24 matrix = []
25 skipped_lines = 0
26 first_invalid_line = 0
27 invalid_value = ''
28 invalid_column = 0
29
30 for i, line in enumerate( file( sys.argv[1] ) ):
31 valid = True
32 line = line.rstrip('\n\r')
33
34 if line and not line.startswith( '#' ):
35 # Extract values and convert to floats
36 row = []
37 for column in columns:
38 column -= 1
39 fields = line.split( "\t" )
40 if len( fields ) <= column:
41 valid = False
42 else:
43 val = fields[column]
44 if val.lower() == "na":
45 row.append( float( "nan" ) )
46 else:
47 try:
48 row.append( float( fields[column] ) )
49 except:
50 valid = False
51 skipped_lines += 1
52 if not first_invalid_line:
53 first_invalid_line = i+1
54 invalid_value = fields[column]
55 invalid_column = column+1
56 else:
57 valid = False
58 skipped_lines += 1
59 if not first_invalid_line:
60 first_invalid_line = i+1
61
62 if valid:
63 matrix.append( row )
64
65 if skipped_lines < i:
66 try:
67 out = open( sys.argv[2], "w" )
68 except:
69 stop_err( "Unable to open output file" )
70
71 # Run correlation
72 try:
73 value = r.cor( array( matrix ), use="pairwise.complete.obs", method=method )
74 except Exception, exc:
75 out.close()
76 stop_err("%s" %str( exc ))
77 for row in value:
78 print >> out, "\t".join( map( str, row ) )
79 out.close()
80
81 if skipped_lines > 0:
82 msg = "..Skipped %d lines starting with line #%d. " %( skipped_lines, first_invalid_line )
83 if invalid_value and invalid_column > 0:
84 msg += "Value '%s' in column %d is not numeric." % ( invalid_value, invalid_column )
85 print msg
86
87 if __name__ == "__main__":
88 main()