Mercurial > repos > tyty > structurefold
comparison reactivity_cal/react_cal.py @ 22:ec9348a4c42f draft
Uploaded
author | tyty |
---|---|
date | Mon, 20 Oct 2014 14:32:13 -0400 |
parents | 7a8ddf1819b1 |
children |
comparison
equal
deleted
inserted
replaced
21:362993bbcf2b | 22:ec9348a4c42f |
---|---|
1 #!/usr/bin/env python | |
2 # -*- coding: utf-8 -*- | |
3 import sys | |
4 from Bio import SeqIO | |
5 import math | |
6 from parse_dis_react import * | |
7 from react_norm_function import * | |
8 import os | |
9 | |
10 | |
11 dist_file1 = sys.argv[1] #plus library | |
12 dist_file2 = sys.argv[2] #minus library | |
13 seq_file = sys.argv[3] | |
14 nt_spec = sys.argv[4] | |
15 flag_in = sys.argv[5] | |
16 threshold = sys.argv[6] | |
17 output_file = sys.argv[7] | |
18 | |
19 | |
20 distri_p = parse_dist(dist_file1) | |
21 distri_m = parse_dist(dist_file2) | |
22 threshold = float(threshold) | |
23 | |
24 print(flag_in) | |
25 | |
26 ospath = os.path.realpath(sys.argv[0]) | |
27 ost = ospath.split('/') | |
28 syspath = "" | |
29 for i in range(len(ost)-1): | |
30 syspath = syspath+ost[i].strip() | |
31 syspath = syspath+'/' | |
32 | |
33 h = file(syspath+"react.txt",'w') | |
34 flag_in = int(flag_in) | |
35 | |
36 seqs = SeqIO.parse(open(seq_file),'fasta'); | |
37 nt_s = set() | |
38 for i in range(len(nt_spec)): | |
39 nt_s.add(nt_spec[i]) | |
40 | |
41 flag = 0 | |
42 trans = [] | |
43 distri_p = distri_p[1] | |
44 distri_m = distri_m[1] | |
45 | |
46 #thres = int(threshold) | |
47 | |
48 | |
49 transcripts = {} | |
50 for seq in seqs: | |
51 n = seq.id | |
52 trans.append(n) | |
53 transcripts[n] = seq.seq.tostring() | |
54 | |
55 | |
56 #print(distri_p) | |
57 | |
58 | |
59 for i in range(0, len(trans)): | |
60 h.write(trans[i]) | |
61 h.write('\n') | |
62 if (trans[i].find('AT1G29930')==-1) and (trans[i].find('At1g29930')==-1): | |
63 for j in range(len(distri_p[trans[i]])): | |
64 distri_p[trans[i]][j] = math.log((int(distri_p[trans[i]][j])+1),math.e) | |
65 for j in range(len(distri_m[trans[i]])): | |
66 distri_m[trans[i]][j] = math.log((int(distri_m[trans[i]][j])+1),math.e) | |
67 s_p = sum(distri_p[trans[i]]) | |
68 s_m = sum(distri_m[trans[i]]) | |
69 length = len(distri_p[trans[i]]) | |
70 if s_p!= 0 and s_m!= 0: | |
71 r = [] | |
72 for j in range(0, len(distri_p[trans[i]])): | |
73 f_p = (float(distri_p[trans[i]][j]))/float(s_p)*length | |
74 f_m = (float(distri_m[trans[i]][j]))/float(s_m)*length | |
75 raw_react = f_p-f_m | |
76 r.append(max(0, raw_react)) | |
77 else: | |
78 for j in range(len(distri_p[trans[i]])): | |
79 distri_p[trans[i]][j] = int(distri_p[trans[i]][j]) | |
80 for j in range(len(distri_m[trans[i]])): | |
81 distri_m[trans[i]][j] = int(distri_m[trans[i]][j]) | |
82 s_p = sum(distri_p[trans[i]]) | |
83 s_m = sum(distri_m[trans[i]]) | |
84 if s_p!= 0 and s_m!= 0: | |
85 r = [] | |
86 for j in range(0, len(distri_p[trans[i]])): | |
87 f_p = float(distri_p[trans[i]][j])/float(s_p) | |
88 f_m = float(distri_m[trans[i]][j])/float(s_m) | |
89 r.append((max(0,(f_p-f_m)))*100) | |
90 | |
91 if s_p!= 0 and s_m!= 0: | |
92 for k in range(1,(len(r)-1)): | |
93 if transcripts[trans[i]][k-1] in nt_s: | |
94 h.write(str(r[k])) | |
95 h.write('\t') | |
96 else: | |
97 h.write('NA') | |
98 h.write('\t') | |
99 k = k+1 | |
100 if transcripts[trans[i]][k-1] in nt_s: | |
101 h.write(str(r[k])) | |
102 h.write('\n') | |
103 else: | |
104 h.write('NA') | |
105 h.write('\n') | |
106 | |
107 | |
108 h.close() | |
109 | |
110 if flag_in: | |
111 react_norm((syspath+"react.txt"),output_file, threshold) | |
112 else: | |
113 h_o = file(output_file, 'w') | |
114 f_i = open(syspath+"react.txt") | |
115 for aline in f_i.readlines(): | |
116 h_o.write(aline.strip()) | |
117 h_o.write('\n') | |
118 os.system("rm -f "+syspath+"react.txt") | |
119 | |
120 | |
121 | |
122 | |
123 | |
124 | |
125 | |
126 | |
127 | |
128 | |
129 | |
130 | |
131 | |
132 | |
133 | |
134 | |
135 | |
136 | |
137 | |
138 | |
139 | |
140 | |
141 | |
142 | |
143 | |
144 | |
145 | |
146 | |
147 | |
148 | |
149 | |
150 | |
151 |