Mercurial > repos > tyty > structurefold

comparison reactivity_cal/react_cal.py @ 5:7a8ddf1819b1 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author tyty
date Mon, 15 Sep 2014 14:52:52 -0400
parents
children
comparison
equal deleted inserted replaced
4:a292aaf51735 5:7a8ddf1819b1
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