annotate reactivity_cal/react_cal.py @ 58:971ef91f6b68 draft

Uploaded
author tyty
date Tue, 18 Nov 2014 01:01:38 -0500
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
58
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
1 #!/usr/bin/env python
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
2 # -*- coding: utf-8 -*-
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
3 import sys
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
4 from Bio import SeqIO
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
5 import math
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
6 from parse_dis_react import *
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
7 from react_norm_function import *
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
8 import os
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
9 import random
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
10 import string
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
11
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
12
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
13 dist_file1 = sys.argv[1] #plus library
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
14 dist_file2 = sys.argv[2] #minus library
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
15 seq_file = sys.argv[3] #Reference library(genome/cDNA)
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
16 nt_spec = sys.argv[4] #only show reactivity for AC or ATCG
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
17 flag_in = sys.argv[5] # perform 2-8% normalization (1) or not (0)
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
18 threshold = sys.argv[6] #Threshold to cap the reactivities
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
19 output_file = sys.argv[7]
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
20
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
21
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
22 distri_p = parse_dist(dist_file1)
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
23 distri_m = parse_dist(dist_file2)
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
24 threshold = float(threshold)
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
25
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
26
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
27 ospath = os.path.realpath(sys.argv[0])
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
28 ost = ospath.split('/')
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
29 syspath = ""
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
30 for i in range(len(ost)-1):
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
31 syspath = syspath+ost[i].strip()
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
32 syspath = syspath+'/'
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
33
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
34 rs = ''.join(random.sample(string.ascii_letters + string.digits, 8))
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
35 os.system("mkdir "+syspath+"output_"+rs)
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
36
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
37 syspathrs = syspath+"output_"+rs+'/'
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
38
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
39 h = file(syspathrs+"react.txt",'w')
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
40 flag_in = int(flag_in)
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
41
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
42 seqs = SeqIO.parse(open(seq_file),'fasta');
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
43 nt_s = set()
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
44 for i in range(len(nt_spec)):
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
45 nt_s.add(nt_spec[i])
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
46
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
47 flag = 0
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
48 trans = []
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
49 distri_p = distri_p[1]
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
50 distri_m = distri_m[1]
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
51
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
52 #thres = int(threshold)
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
53
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
54
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
55 transcripts = {}
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
56 for seq in seqs:
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
57 n = seq.id
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
58 trans.append(n)
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
59 transcripts[n] = seq.seq.tostring()
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
60
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
61
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
62 #print(distri_p)
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
63
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
64
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
65 for i in range(0, len(trans)):
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
66 h.write(trans[i])
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
67 h.write('\n')
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
68 for j in range(len(distri_p[trans[i]])):
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
69 distri_p[trans[i]][j] = math.log((int(distri_p[trans[i]][j])+1),math.e)
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
70 for j in range(len(distri_m[trans[i]])):
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
71 distri_m[trans[i]][j] = math.log((int(distri_m[trans[i]][j])+1),math.e)
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
72 s_p = sum(distri_p[trans[i]])
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
73 s_m = sum(distri_m[trans[i]])
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
74 length = len(distri_p[trans[i]])
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
75 if s_p!= 0 and s_m!= 0:
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
76 r = []
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
77 for j in range(0, len(distri_p[trans[i]])):
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
78 f_p = (float(distri_p[trans[i]][j]))/float(s_p)*length
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
79 f_m = (float(distri_m[trans[i]][j]))/float(s_m)*length
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
80 raw_react = f_p-f_m
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
81 r.append(max(0, raw_react))
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
82
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
83 if s_p!= 0 and s_m!= 0:
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
84 for k in range(1,(len(r)-1)):
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
85 if transcripts[trans[i]][k-1] in nt_s:
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
86 h.write(str(r[k]))
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
87 h.write('\t')
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
88 else:
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
89 h.write('NA')
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
90 h.write('\t')
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
91 k = k+1
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
92 if transcripts[trans[i]][k-1] in nt_s:
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
93 h.write(str(r[k]))
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
94 h.write('\n')
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
95 else:
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
96 h.write('NA')
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
97 h.write('\n')
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
98
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
99
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
100 h.close()
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
101
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
102 if flag_in:
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
103 react_norm((syspathrs+"react.txt"),output_file, threshold)
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
104 else:
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
105 h_o = file(output_file, 'w')
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
106 f_i = open(syspathrs+"react.txt")
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
107 for aline in f_i.readlines():
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
108 h_o.write(aline.strip())
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
109 h_o.write('\n')
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
110 os.system("rm -f "+syspathrs+"react.txt")
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
111
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
112 os.system("rm -r "+syspathrs)
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
113
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
114
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
115
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
116
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
117
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
118
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
119
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
120
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
121
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
122
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
123
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
124
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
125
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
126
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
127
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
128
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
129
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
130
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
131
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
132
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
133
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
134
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
135
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
136
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
137
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
138
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
139
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
140
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
141
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
142
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
143
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
144
971ef91f6b68 Uploaded
tyty
parents:
diff changeset
145