annotate poisson2test.py @ 1:5ff18ec88181 draft default tip

planemo upload for repository https://github.com/galaxyproject/tools-devteam/tree/master/tool_collections/taxonomy/poisson2test commit a1517c9d22029095120643bbe2c8fa53754dd2b7
author devteam
date Wed, 11 Nov 2015 12:56:24 -0500
parents 8cd5945559b8
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
1 #!/usr/local/bin/python
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
2
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
3 import sys
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
4 from math import *
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
5 from rpy import *
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
6
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
7
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
8 if ((len(sys.argv)-1) != 6):
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
9 print 'too few parameters'
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
10 print 'usage: inputfile, col1, col2, d-value(not 0), p-val correction method(0 or 1)'
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
11 sys.exit()
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
12
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
13 try:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
14 lines_arr = open(sys.argv[1]).readlines()
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
15 except IOError:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
16 print'cannot open',sys.argv[1]
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
17 sys.exit()
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
18
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
19 try:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
20 i = int(sys.argv[2]) #first column to compare
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
21 j = int(sys.argv[3]) #second colum to compare
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
22 d = float(sys.argv[4]) #correction factor
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
23 k = int(sys.argv[5]) #p-val correction method
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
24 outfile = open(sys.argv[6],'w') # output data
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
25
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
26 if (i>j):
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
27 print 'column order not correct col1 < col2'
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
28 print 'usage: inputfile, col1, col2, d-value, p-val correction method'
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
29 sys.exit()
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
30
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
31 try:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
32 a = 1 / d
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
33 assert k in [0,1]
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
34 except ZeroDivisionError:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
35 print 'd cannot be 0'
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
36 print 'usage: inputfile, col1, col2, d-value, p-val correction method'
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
37 sys.exit()
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
38 except:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
39 print ' p-val correction should be 0 or 1 (0 = "bonferroni", 1 = "fdr")'
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
40 print 'usage: inputfile, col1, col2, d-value, p-val correction method'
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
41 sys.exit()
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
42 except ValueError:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
43 print 'parameters are not integers'
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
44 print 'usage: inputfile, col1, col2, d-value, p-val correction method'
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
45 sys.exit()
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
46
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
47
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
48 fsize = len(lines_arr)
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
49
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
50 z1 = []
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
51 z2 = []
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
52 pz1 = []
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
53 pz2 = []
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
54 field = []
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
55
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
56 if d<1: # Z score calculation
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
57 for line in lines_arr:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
58 line.strip()
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
59 field = line.split('\t')
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
60
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
61 x = int(field[j-1]) #input column 2
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
62 y = int(field[i-1]) #input column 1
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
63 if y>x:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
64 z1.append(float((y - ((1/d)*x))/sqrt((1/d)*(x + y))))
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
65 z2.append(float((2*(sqrt(y+(3/8))-sqrt((1/d)*(x+(3/8)))))/sqrt(1+(1/d))))
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
66 else:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
67 tmp_var1 = x
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
68 x = y
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
69 y = tmp_var1
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
70 z1.append(float((y - (d*x))/sqrt(d*(x + y))))
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
71 z2.append(float((2*(sqrt(y+(3/8))-sqrt(d*(x+(3/8)))))/sqrt(1+d)))
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
72
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
73 else: #d>1 Z score calculation
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
74 for line in lines_arr:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
75 line.strip()
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
76 field = line.split('\t')
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
77 x = int(field[i-1]) #input column 1
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
78 y = int(field[j-1]) #input column 2
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
79
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
80 if y>x:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
81 z1.append(float((y - (d*x))/sqrt(d*(x + y))))
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
82 z2.append(float((2*(sqrt(y+(3/8))-sqrt(d*(x+(3/8)))))/sqrt(1+d)))
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
83 else:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
84 tmp_var2 = x
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
85 x = y
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
86 y = tmp_var2
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
87 z1.append(float((y - ((1/d)*x))/sqrt((1/d)*(x + y))))
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
88 z2.append(float((2*(sqrt(y+(3/8))-sqrt((1/d)*(x+(3/8)))))/sqrt(1+(1/d))))
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
89
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
90
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
91
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
92
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
93
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
94 # P-value caluculation for z1 and z2
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
95 for p in z1:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
96
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
97 pz1.append(float(r.pnorm(-abs(float(p)))))
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
98
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
99 for q in z2:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
100
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
101 pz2.append(float(r.pnorm(-abs(float(q)))))
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
102
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
103 # P-value correction for pz1 and pz2
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
104
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
105 if k == 0:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
106 corrz1 = r.p_adjust(pz1,"bonferroni",fsize)
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
107 corrz2 = r.p_adjust(pz2,"bonferroni",fsize)
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
108
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
109
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
110 else:
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
111
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
112 corrz1 = r.p_adjust(pz1,"fdr",fsize)
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
113 corrz2 = r.p_adjust(pz2,"fdr",fsize)
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
114
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
115
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
116 #printing all columns
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
117 for n in range(fsize):
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
118 print >> outfile, "%s\t%4.3f\t%4.3f\t%8.6f\t%8.6f\t%8.6f\t%8.6f" %(lines_arr[n].strip(),z1[n],z2[n],pz1[n],pz2[n],corrz1[n],corrz2[n])
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
119
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
120
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
121
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
122
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
123
8cd5945559b8 Imported from capsule None
devteam
parents:
diff changeset
124