Mercurial > repos > tyty > structurefold

comparison reactivity_cal/react_cal.py @ 63:c1f1b552c1b8 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author tyty
date Tue, 18 Nov 2014 15:54:11 -0500
parents
children
comparison
equal deleted inserted replaced
62:58df1060fea1 63:c1f1b552c1b8
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 import random
10 import string
11
12
13 dist_file1 = sys.argv[1] #plus library
14 dist_file2 = sys.argv[2] #minus library
15 seq_file = sys.argv[3] #Reference library(genome/cDNA)
16 nt_spec = sys.argv[4] #only show reactivity for AC or ATCG
17 flag_in = sys.argv[5] # perform 2-8% normalization (1) or not (0)
18 threshold = sys.argv[6] #Threshold to cap the reactivities
19 output_file = sys.argv[7]
20
21
22 distri_p = parse_dist(dist_file1)
23 distri_m = parse_dist(dist_file2)
24 threshold = float(threshold)
25
26
27 syspathrs = os.getcwd()
28
29 h = file(syspathrs+"react.txt",'w')
30 flag_in = int(flag_in)
31
32 seqs = SeqIO.parse(open(seq_file),'fasta');
33 nt_s = set()
34 for i in range(len(nt_spec)):
35 nt_s.add(nt_spec[i])
36
37 flag = 0
38 trans = []
39 distri_p = distri_p[1]
40 distri_m = distri_m[1]
41
42 #thres = int(threshold)
43
44
45 transcripts = {}
46 for seq in seqs:
47 n = seq.id
48 trans.append(n)
49 transcripts[n] = seq.seq.tostring()
50
51
52 #print(distri_p)
53
54
55 for i in range(0, len(trans)):
56 h.write(trans[i])
57 h.write('\n')
58 for j in range(len(distri_p[trans[i]])):
59 distri_p[trans[i]][j] = math.log((int(distri_p[trans[i]][j])+1),math.e)
60 for j in range(len(distri_m[trans[i]])):
61 distri_m[trans[i]][j] = math.log((int(distri_m[trans[i]][j])+1),math.e)
62 s_p = sum(distri_p[trans[i]])
63 s_m = sum(distri_m[trans[i]])
64 length = len(distri_p[trans[i]])
65 if s_p!= 0 and s_m!= 0:
66 r = []
67 for j in range(0, len(distri_p[trans[i]])):
68 f_p = (float(distri_p[trans[i]][j]))/float(s_p)*length
69 f_m = (float(distri_m[trans[i]][j]))/float(s_m)*length
70 raw_react = f_p-f_m
71 r.append(max(0, raw_react))
72
73 if s_p!= 0 and s_m!= 0:
74 for k in range(1,(len(r)-1)):
75 if transcripts[trans[i]][k-1] in nt_s:
76 h.write(str(r[k]))
77 h.write('\t')
78 else:
79 h.write('NA')
80 h.write('\t')
81 k = k+1
82 if transcripts[trans[i]][k-1] in nt_s:
83 h.write(str(r[k]))
84 h.write('\n')
85 else:
86 h.write('NA')
87 h.write('\n')
88
89
90 h.close()
91
92 if flag_in:
93 react_norm((syspathrs+"react.txt"),output_file, threshold)
94 else:
95 h_o = file(output_file, 'w')
96 f_i = open(syspathrs+"react.txt")
97 for aline in f_i.readlines():
98 h_o.write(aline.strip())
99 h_o.write('\n')
100 os.system("rm -f "+syspathrs+"react.txt")
101
102 #os.system("rm -r "+syspathrs)
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135