root/galaxy-central/tools/stats/cor.py

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

import galaxy-central

行番号 
1#!/usr/bin/env python
2#Greg Von Kuster
3"""
4Calculate correlations between numeric columns in a tab delim file.
5usage: %prog infile output.txt columns method
6"""
7
8import sys
9from rpy import *
10
11def stop_err(msg):
12    sys.stderr.write(msg)
13    sys.exit()
14   
15def 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
87if __name__ == "__main__":
88    main()
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。