|
6
|
1 # -*- coding: utf-8 -*-
|
|
|
2 """
|
|
|
3 Created on Thu Jul 19 13:45:05 2018
|
|
|
4
|
|
|
5 @author: Marta
|
|
|
6 """
|
|
|
7
|
|
|
8 from Bio import SeqIO
|
|
|
9 import numpy as np
|
|
|
10 import pandas as pd
|
|
|
11 from auxiliar import free_energy,freq_base
|
|
|
12 from Bio.Seq import Seq
|
|
|
13 from Bio.SeqRecord import SeqRecord
|
|
|
14 from Bio.Alphabet import IUPAC
|
|
|
15 from auxiliar import get_bacteria, get_families, get_max_pssm, get_scores, get_lifecycle
|
|
|
16
|
|
|
17 #division of the test genome in sequences of 65 bp
|
|
|
18 def get_testseqs65(form,fic,both=False):
|
|
|
19 ALL = []
|
|
|
20 indexes = []
|
|
|
21 a = 0
|
|
|
22 rec = SeqIO.read(fic,form)
|
|
|
23 genome = rec.seq
|
|
|
24 i = 0
|
|
|
25 j = 65
|
|
|
26 while j < len(genome):
|
|
|
27 s = genome[i:j]
|
|
|
28 ALL.append([1,i,j,s])
|
|
|
29 i += 20
|
|
|
30 j += 20
|
|
|
31 a += 1
|
|
|
32 indexes.append(rec.name+":"+str(a))
|
|
|
33 if both:
|
|
|
34 comp = genome.reverse_complement()
|
|
|
35 size = len(rec.seq)
|
|
|
36 i = 0
|
|
|
37 j = 65
|
|
|
38 while j < len(comp):
|
|
|
39 s = comp[i:j]
|
|
|
40 ALL.append([-1,size-j,size-i,s])
|
|
|
41 i += 20
|
|
|
42 j += 20
|
|
|
43 a += 1
|
|
|
44 indexes.append(rec.name+":"+str(a))
|
|
|
45 df = pd.DataFrame(ALL, index=indexes, columns=['strand','iniprom','endprom','seq'])
|
|
|
46 return df
|
|
|
47
|
|
|
48 #calculate the scores of all sequences (similar to get_posScores and get_negScores)
|
|
|
49 def get_testScores(loc,test):
|
|
|
50 scores = []
|
|
|
51 posis = []
|
|
|
52 sizes = []
|
|
|
53 dic = {}
|
|
|
54 for ind,row in test.iterrows():
|
|
|
55 _,window = ind.split(':')
|
|
|
56 strand = row['strand']
|
|
|
57 ini = row['iniprom']
|
|
|
58 end = row['endprom']
|
|
|
59 seq = row['seq']
|
|
|
60 pos = [ini,end,strand]
|
|
|
61 dic[window] = pos
|
|
|
62 s = seq
|
|
|
63 score10_6,pos10_6 = get_scores(os.path.join(loc,'pssm10_6.txt'), s)
|
|
|
64 maxi10_6 = get_max_pssm(os.path.join(loc,'pssm10_6.txt'))
|
|
|
65 score10_8,pos10_8 = get_scores(os.path.join(loc,'pssm10_8.txt'), s)
|
|
|
66 maxi10_8 = get_max_pssm(os.path.join(loc,'pssm10_8.txt'))
|
|
|
67 scores23,pos23 = get_scores(os.path.join(loc,'pssm_23.txt'), s)
|
|
|
68 maxi23 = get_max_pssm(os.path.join(loc,'pssm_23.txt'))
|
|
|
69 scores21,pos21 = get_scores(os.path.join(loc,'pssm_21.txt'), s)
|
|
|
70 maxi21 = get_max_pssm(os.path.join(loc,'pssm_21.txt'))
|
|
|
71 scores27,pos27 = get_scores(os.path.join(loc,'pssm_27.txt'), s)
|
|
|
72 maxi27 = get_max_pssm(os.path.join(loc,'pssm_27.txt'))
|
|
|
73 scores32,pos32 = get_scores(os.path.join(loc,'pssm_32.txt'), s)
|
|
|
74 maxi32 = get_max_pssm(os.path.join(loc,'pssm_32.txt'))
|
|
|
75 score23 = max(scores23)
|
|
|
76 score21 = max(scores21)
|
|
|
77 score27 = max(scores27)
|
|
|
78 score32 = max(scores32)
|
|
|
79 maxiphage = max(score23,score21,score27,score32)
|
|
|
80 if maxiphage == score23: phage_max = score23*maxi23
|
|
|
81 elif maxiphage == score21: phage_max = score21*maxi21
|
|
|
82 elif maxiphage == score27: phage_max = score27*maxi27
|
|
|
83 elif maxiphage == score32: phage_max = score32*maxi32
|
|
|
84 score35_6,pos35_6 = get_scores(os.path.join(loc,'pssm35_6.txt'), s)
|
|
|
85 maxi35_6 = get_max_pssm(os.path.join(loc,'pssm35_6.txt'))
|
|
|
86 score35_9,pos35_9 = get_scores(os.path.join(loc,'pssm35_9.txt'), s)
|
|
|
87 maxi35_9 = get_max_pssm(os.path.join(loc,'pssm35_9.txt'))
|
|
|
88 score35_t4,pos35_t4 = get_scores(os.path.join(loc,'pssm35_t4.txt'), s)
|
|
|
89 maxi35_t4 = get_max_pssm(os.path.join(loc,'pssm35_t4.txt'))
|
|
|
90 score35_cbb,pos35_cbb = get_scores(os.path.join(loc,'pssm35_cbb.txt'), s)
|
|
|
91 maxi35_cbb = get_max_pssm(os.path.join(loc,'pssm35_cbb.txt'))
|
|
|
92 score35_lb,pos35_lb = get_scores(os.path.join(loc,'pssm35_lb.txt'),s)
|
|
|
93 maxi35_lb = get_max_pssm(os.path.join(loc,'pssm35_lb.txt'))
|
|
|
94 score35_mu, pos35_mu = get_scores(os.path.join(loc,'pssm35_mu.txt'),s)
|
|
|
95 maxi35_mu = get_max_pssm(os.path.join(loc,'pssm35_mu.txt'))
|
|
|
96
|
|
|
97 dists6 = []
|
|
|
98 score6 = []
|
|
|
99 for p in pos10_6:
|
|
|
100 for a in range(14,22):
|
|
|
101 d = p-a-6
|
|
|
102 if d >= 0:
|
|
|
103 s10 = score10_6[p]
|
|
|
104 s35_6 = score35_6[d]
|
|
|
105 score6.append([s35_6,s10])
|
|
|
106 dists6.append([d,p])
|
|
|
107 dists9 = []
|
|
|
108 score9 = []
|
|
|
109 for p in pos10_6:
|
|
|
110 for a in range(11,14):
|
|
|
111 d = p-a-9
|
|
|
112 if d >= 0:
|
|
|
113 s10 = score10_6[p]
|
|
|
114 s35_9 = score35_9[d]
|
|
|
115 score9.append([s35_9,s10])
|
|
|
116 dists9.append([d,p])
|
|
|
117 distst4 = []
|
|
|
118 scoret4 = []
|
|
|
119 distscbb = []
|
|
|
120 scorecbb = []
|
|
|
121 for p in pos10_6:
|
|
|
122 for a in range(16,18):
|
|
|
123 d = p-a-7
|
|
|
124 if d >= 0:
|
|
|
125 s10 = score10_6[p]
|
|
|
126 s35_t4 = score35_t4[d]
|
|
|
127 s35_cbb = score35_cbb[d]
|
|
|
128 scoret4.append([s35_t4,s10])
|
|
|
129 distst4.append([d,p])
|
|
|
130 scorecbb.append([s35_cbb,s10])
|
|
|
131 distscbb.append([d,p])
|
|
|
132
|
|
|
133
|
|
|
134 distsmu = []
|
|
|
135 scoremu = []
|
|
|
136 for p in pos10_6:
|
|
|
137 d = p-16-14
|
|
|
138 if d >= 0:
|
|
|
139 s10 = score10_6[p]
|
|
|
140 s35_mu = score35_mu[d]
|
|
|
141 scoremu.append([s35_mu,s10])
|
|
|
142 distsmu.append([d,p])
|
|
|
143
|
|
|
144 distslb = []
|
|
|
145 scorelb = []
|
|
|
146 for p in pos10_6:
|
|
|
147 d = p-13-14
|
|
|
148 if d >= 0:
|
|
|
149 s10 = score10_6[p]
|
|
|
150 s35_lb = score35_lb[d]
|
|
|
151 scorelb.append([s35_lb,s10])
|
|
|
152 distslb.append([d,p])
|
|
|
153
|
|
|
154 soma6 = [sum(x) for x in score6]
|
|
|
155 soma9 = [sum(x) for x in score9]
|
|
|
156 somat4 = [sum(x) for x in scoret4]
|
|
|
157 somacbb = [sum(x) for x in scorecbb]
|
|
|
158 somamu = [sum(x) for x in scoremu]
|
|
|
159 somalb = [sum(x) for x in scorelb]
|
|
|
160
|
|
|
161 maxi6 = max(soma6)
|
|
|
162 maxi9 = max(soma9)
|
|
|
163 maxit4 = max(somat4)
|
|
|
164 maxicbb = max(somacbb)
|
|
|
165 maximu = max(somamu)
|
|
|
166 maxilb = max(somalb)
|
|
|
167 maxi_elems = max(maxi6,maxi9,maxit4,maxicbb,maxilb,maximu)
|
|
|
168
|
|
|
169 if maxi_elems == maxilb:
|
|
|
170 indmax = somalb.index(maxilb)
|
|
|
171 sc35 = scorelb[indmax][0]
|
|
|
172 sc10 = scorelb[indmax][1]
|
|
|
173 score_elems = [sc35,sc10]
|
|
|
174 posel = distslb[indmax]
|
|
|
175 size35 = 14
|
|
|
176 elems_maxi = sc35*maxi35_lb+sc10*maxi10_6
|
|
|
177 elif maxi_elems == maximu:
|
|
|
178 indmax = somamu.index(maximu)
|
|
|
179 sc35 = scoremu[indmax][0]
|
|
|
180 sc10 = scoremu[indmax][1]
|
|
|
181 score_elems = [sc35,sc10]
|
|
|
182 posel = distsmu[indmax]
|
|
|
183 size35 = 14
|
|
|
184 elems_maxi = sc35*maxi35_mu+sc10*maxi10_6
|
|
|
185 elif maxi_elems == maxi9:
|
|
|
186 indmax = soma9.index(maxi9)
|
|
|
187 sc35 = score9[indmax][0]
|
|
|
188 sc10 = score9[indmax][1]
|
|
|
189 score_elems = [sc35,sc10]
|
|
|
190 posel = dists9[indmax]
|
|
|
191 size35 = 9
|
|
|
192 elems_maxi = sc35*maxi35_9+sc10*maxi10_6
|
|
|
193 elif maxi_elems == maxit4:
|
|
|
194 indmax = somat4.index(maxit4)
|
|
|
195 sc35 = scoret4[indmax][0]
|
|
|
196 sc10 = scoret4[indmax][1]
|
|
|
197 score_elems = [sc35,sc10]
|
|
|
198 posel = distst4[indmax]
|
|
|
199 size35 = 7
|
|
|
200 elems_maxi = sc35*maxi35_t4+sc10*maxi10_6
|
|
|
201 elif maxi_elems == maxicbb:
|
|
|
202 indmax = somacbb.index(maxicbb)
|
|
|
203 sc35 = scorecbb[indmax][0]
|
|
|
204 sc10 = scorecbb[indmax][1]
|
|
|
205 score_elems = [sc35,sc10]
|
|
|
206 posel = distscbb[indmax]
|
|
|
207 size35 = 7
|
|
|
208 elems_maxi = sc35*maxi35_cbb+sc10*maxi10_6
|
|
|
209 else:
|
|
|
210 indmax = soma6.index(maxi6)
|
|
|
211 sc35 = score6[indmax][0]
|
|
|
212 sc10 = score6[indmax][1]
|
|
|
213 score_elems = [sc35,sc10]
|
|
|
214 posel = dists6[indmax]
|
|
|
215 size35 = 6
|
|
|
216 elems_maxi = sc35*maxi35_6+sc10*maxi10_6
|
|
|
217
|
|
|
218 if score23 == maxiphage:
|
|
|
219 phage_score = score23
|
|
|
220 posiphage = scores23.index(score23)
|
|
|
221 sizephage = 23
|
|
|
222 elif score21 == maxiphage:
|
|
|
223 phage_score = score21
|
|
|
224 posiphage = scores21.index(score21)
|
|
|
225 sizephage = 21
|
|
|
226 elif score27 == maxiphage:
|
|
|
227 phage_score = score27
|
|
|
228 posiphage = scores27.index(score27)
|
|
|
229 sizephage = 27
|
|
|
230 else:
|
|
|
231 phage_score = score32
|
|
|
232 posiphage = scores32.index(score32)
|
|
|
233 sizephage = 32
|
|
|
234
|
|
|
235 if elems_maxi >= max(score10_8)*maxi10_8:
|
|
|
236 i = posel[1]
|
|
|
237 ext = s[i-3:i-1]
|
|
|
238 if ext == 'TG': tg = 1
|
|
|
239 else: tg = 0
|
|
|
240 if elems_maxi > phage_max: host = 1
|
|
|
241 else:
|
|
|
242 host = 0
|
|
|
243 tg = 0
|
|
|
244 sc = max(score10_8)
|
|
|
245 end35 = posel[0]+size35
|
|
|
246 dist = posel[1]-end35
|
|
|
247 scores.append([host, score_elems[1],sc,score_elems[0],phage_score,tg,dist,str(seq)])
|
|
|
248 posis.append([posel[1],posel[0],posiphage])
|
|
|
249 sizes.append([6,size35,sizephage])
|
|
|
250 else:
|
|
|
251 host = 1
|
|
|
252 sc = max(score10_8)
|
|
|
253 i = score10_8.index(sc)
|
|
|
254 ext = s[i-3:i-1]
|
|
|
255 if ext == 'TG': tg = 1
|
|
|
256 else: tg = 0
|
|
|
257 if max(score10_8)*maxi10_8 > phage_max: host = 1
|
|
|
258 else:
|
|
|
259 host = 0
|
|
|
260 tg = 0
|
|
|
261 end35 = posel[0]+size35
|
|
|
262 dist = posel[1]-end35
|
|
|
263 scores.append([host,score_elems[1],sc,score_elems[0],phage_score,tg,dist,str(seq)])
|
|
|
264 posis.append([i,posel[0],posiphage])
|
|
|
265 sizes.append([8,size35,sizephage])
|
|
|
266 score = pd.DataFrame(scores, index=test.index, columns=['host','score10','score10_8','score35','score_phage','tg','dist','seq'])
|
|
|
267 posis = pd.DataFrame(posis, index=test.index, columns=['pos10','pos35','posphage'])
|
|
|
268 sizes = pd.DataFrame(sizes, index=test.index, columns=['size10','size35','size_phage'])
|
|
|
269 df_all = pd.concat([score,posis,sizes],axis=1)
|
|
|
270 return df_all,dic
|
|
|
271
|
|
|
272 def create_dftest(scores_test,dic_window,family,bacteria,lifecycle):
|
|
|
273 tudo = []
|
|
|
274 tudo2 = []
|
|
|
275 for ind,row in scores_test.iterrows():
|
|
|
276 _,window = ind.split(':')
|
|
|
277 posis = dic_window[window]
|
|
|
278 strand=posis[2]
|
|
|
279 if strand == 1: ini=posis[0]
|
|
|
280 else: ini=posis[1]
|
|
|
281 seqprom = row['seq']
|
|
|
282 score10 = row['score10']
|
|
|
283 score10_8 = row['score10_8']
|
|
|
284 score35 = row['score35']
|
|
|
285 scorephage = row['score_phage']
|
|
|
286 size10 = row['size10']
|
|
|
287 size35 = row['size35']
|
|
|
288 sizephage = row['size_phage']
|
|
|
289 ini10 = row['pos10']
|
|
|
290 tg = row['tg']
|
|
|
291 host = row['host']
|
|
|
292 ini35 = row['pos35']
|
|
|
293 dist = row['dist']
|
|
|
294 end10=ini10+size10
|
|
|
295 iniphage = row['posphage']
|
|
|
296 endphage = iniphage+sizephage
|
|
|
297 if strand == 1:
|
|
|
298 if host == 0:
|
|
|
299 new_seq = seqprom[iniphage:endphage]
|
|
|
300 new_ini = ini+iniphage+1
|
|
|
301 new_end = ini+endphage
|
|
|
302 else:
|
|
|
303 if size10 == 6:
|
|
|
304 new_seq = seqprom[ini35:end10]
|
|
|
305 new_ini = ini+ini35+1
|
|
|
306 new_end = ini+end10
|
|
|
307 else:
|
|
|
308 new_seq = seqprom[ini10:end10]
|
|
|
309 new_ini = ini+ini10+1
|
|
|
310 new_end = ini+end10
|
|
|
311 new_pos = '('+str(new_ini)+'..'+str(new_end)+')'
|
|
|
312 else:
|
|
|
313 if host == 0:
|
|
|
314 new_seq = seqprom[iniphage:endphage]
|
|
|
315 new_ini = ini-endphage+1
|
|
|
316 new_end = ini-iniphage
|
|
|
317 else:
|
|
|
318 if size10 == 6:
|
|
|
319 new_seq = seqprom[ini35:end10]
|
|
|
320 new_ini = ini-end10+1
|
|
|
321 new_end = ini-ini35
|
|
|
322 else:
|
|
|
323 new_seq = seqprom[ini10:end10]
|
|
|
324 new_ini = ini-end10+1
|
|
|
325 new_end = ini-ini10
|
|
|
326 new_pos = 'complement('+str(new_ini)+'..'+str(new_end)+')'
|
|
|
327 if size10 == 6:
|
|
|
328 size10_6 = 1
|
|
|
329 else:
|
|
|
330 size10_6 = 0
|
|
|
331 if size35 == 6:
|
|
|
332 size35_6 = 1
|
|
|
333 size35_7 = 0
|
|
|
334 size35_9 = 0
|
|
|
335 size35_14 = 0
|
|
|
336 elif size35 == 7:
|
|
|
337 size35_6 = 0
|
|
|
338 size35_7 = 1
|
|
|
339 size35_9 = 0
|
|
|
340 size35_14 = 0
|
|
|
341 elif size35 == 9:
|
|
|
342 size35_6 = 0
|
|
|
343 size35_7 = 0
|
|
|
344 size35_9 = 1
|
|
|
345 size35_14 = 0
|
|
|
346 else:
|
|
|
347 size35_6 = 0
|
|
|
348 size35_7 = 0
|
|
|
349 size35_9 = 0
|
|
|
350 size35_14 = 1
|
|
|
351 if sizephage == 23:
|
|
|
352 sizephage_23 = 1
|
|
|
353 sizephage_21 = 0
|
|
|
354 sizephage_32 = 0
|
|
|
355 elif sizephage == 21:
|
|
|
356 sizephage_23 = 0
|
|
|
357 sizephage_21 = 1
|
|
|
358 sizephage_32 = 0
|
|
|
359 elif sizephage == 32:
|
|
|
360 sizephage_23 = 0
|
|
|
361 sizephage_21 = 0
|
|
|
362 sizephage_32 = 1
|
|
|
363 else:
|
|
|
364 sizephage_23 = 0
|
|
|
365 sizephage_21 = 0
|
|
|
366 sizephage_32 = 0
|
|
|
367 if family == 'Podoviridae':
|
|
|
368 Podo = 1
|
|
|
369 Sipho = 0
|
|
|
370 Myo = 0
|
|
|
371 elif family == 'Siphoviridae':
|
|
|
372 Podo = 0
|
|
|
373 Sipho = 1
|
|
|
374 Myo = 0
|
|
|
375 elif family == 'Myoviridae':
|
|
|
376 Podo = 0
|
|
|
377 Sipho = 0
|
|
|
378 Myo = 1
|
|
|
379 else:
|
|
|
380 Podo = 0
|
|
|
381 Sipho = 0
|
|
|
382 Myo = 0
|
|
|
383 if bacteria == 'Bacillus':
|
|
|
384 bac = [1,0,0,0,0]
|
|
|
385 elif bacteria == 'Escherichia coli':
|
|
|
386 bac = [0,1,0,0,0]
|
|
|
387 elif bacteria == 'Klebsiella':
|
|
|
388 bac = [0,0,1,0,0]
|
|
|
389 elif bacteria == 'Pectobacterium':
|
|
|
390 bac = [0,0,0,1,0]
|
|
|
391 elif bacteria == 'Cronobacter':
|
|
|
392 bac = [0,0,0,0,1]
|
|
|
393 else:
|
|
|
394 bac = [0,0,0,0,0]
|
|
|
395 if lifecycle == 'virulent': tp = 0
|
|
|
396 else: tp = 1
|
|
|
397 fe = free_energy(str(seqprom))
|
|
|
398 AT = freq_base(str(seqprom))
|
|
|
399 linha = [score10, score10_8, score35, dist, scorephage, fe, AT, host,size10_6,
|
|
|
400 size35_6, size35_7, size35_9, size35_14,
|
|
|
401 sizephage_23, sizephage_21,
|
|
|
402 sizephage_32, tg, Podo, Sipho, Myo,tp]
|
|
|
403 linha.extend(bac)
|
|
|
404 tudo.append(linha)
|
|
|
405 linha2 = [new_pos,str(new_seq), host, size10_6,
|
|
|
406 score10, score10_8, size35_6, size35_7, size35_9,size35_14,
|
|
|
407 score35, dist, sizephage_23, sizephage_21,
|
|
|
408 sizephage_32, scorephage, tg, Podo, Sipho, Myo,tp, fe, AT]
|
|
|
409 linha2.extend(bac)
|
|
|
410 tudo2.append(linha2)
|
|
|
411 df_test = pd.DataFrame(tudo, index=scores_test.index,
|
|
|
412 columns = ['score10', 'score10_8','score35', 'dist35_10',
|
|
|
413 'scorephage','fe', 'freqAT',
|
|
|
414 'host','size10',
|
|
|
415 'size35_6', 'size35_7','size35_9','size35_14',
|
|
|
416 'sizephage_23',
|
|
|
417 'sizephage_21', 'sizephage_32',
|
|
|
418 'TG', 'Podo', 'Sipho', 'Myo', 'tp',
|
|
|
419 'Bacillus', 'EColi',
|
|
|
420 'Pectobacterium','Klebsiella', 'Cronobacter'])
|
|
|
421 df_INFO = pd.DataFrame(tudo2, index=scores_test.index,
|
|
|
422 columns = ['Positions','Promoter Sequence','host','size10',
|
|
|
423 'score10', 'score10_8',
|
|
|
424 'size35_6', 'size35_7','size35_9','size35_14',
|
|
|
425 'score35', 'dist35_10','sizephage_23',
|
|
|
426 'sizephage_21', 'sizephage_32',
|
|
|
427 'scorephage', 'TG', 'Podo', 'Sipho', 'Myo', 'tp','fe', 'freqAT',
|
|
|
428 'EColi', 'Salmonella',
|
|
|
429 'Pectobacterium','Cronobacter', 'Streptococcus'])
|
|
|
430 return df_test,df_INFO
|
|
|
431
|
|
|
432
|
|
|
433 def get_predictions(scaler_file,model_file,test,df_testinfo,threshold):
|
|
|
434 from sklearn.externals import joblib
|
|
|
435 scaler = joblib.load(scaler_file)
|
|
|
436 model = joblib.load(model_file)
|
|
|
437 feat_scaled = pd.DataFrame(scaler.transform(test.iloc[:,:7]),index =test.index, columns=test.columns[:7])
|
|
|
438 TEST_scaled = pd.concat([feat_scaled,test.iloc[:,7:]],axis=1)
|
|
|
439 scores = model.predict_proba(TEST_scaled)
|
|
|
440 pos_scores = np.empty((TEST_scaled.shape[0],0), float)
|
|
|
441 for x in scores: pos_scores = np.append(pos_scores,x[1])
|
|
|
442 try: positive_indexes = np.nonzero(pos_scores>float(threshold))[0] #escolher os positivos, podia ser escolher com score > x
|
|
|
443 except ValueError: return 'The threshold value is not a float'
|
|
|
444 else:
|
|
|
445 if len(positive_indexes) == 0: return None
|
|
|
446 else:
|
|
|
447 positive_windows = TEST_scaled.index[positive_indexes]
|
|
|
448 INFO = df_testinfo.loc[positive_windows,['Positions','Promoter Sequence']]
|
|
|
449 promoter_type = []
|
|
|
450 for x in df_testinfo.loc[positive_windows,'host'].tolist():
|
|
|
451 if x == 0: promoter_type.append('phage')
|
|
|
452 else: promoter_type.append('host')
|
|
|
453 INFO['Type'] = promoter_type
|
|
|
454 INFO['Scores'] = np.around(pos_scores[positive_indexes],decimals=3)
|
|
|
455 INFO.index = positive_windows
|
|
|
456 return INFO
|
|
|
457
|
|
|
458
|
|
|
459 def get_finaldf(test,rec):
|
|
|
460 new_df = test.groupby(['Positions'])
|
|
|
461 groups = list(new_df.groups.keys())
|
|
|
462 for i in range(len(groups)-1):
|
|
|
463 for j in range(i, len(groups)):
|
|
|
464 if 'complement' in groups[i]: inii = int(groups[i][11:].split('..')[0])
|
|
|
465 else: inii = int(groups[i][1:].split('..')[0])
|
|
|
466 if 'complement' in groups[j]: inij = int(groups[j][11:].split('..')[0])
|
|
|
467 else: inij = int(groups[j][1:].split('..')[0])
|
|
|
468 if inij < inii:
|
|
|
469 temp = groups[i]
|
|
|
470 groups[i] = groups[j]
|
|
|
471 groups[j] = temp
|
|
|
472 new_inds = []
|
|
|
473 for g in groups:
|
|
|
474 inds = new_df.groups[g]
|
|
|
475 if len(inds) == 1: new_inds.append(inds[0])
|
|
|
476 else:
|
|
|
477 maxi = max(new_df.get_group(g)['Scores'])
|
|
|
478 i = new_df.groups[g][new_df.get_group(g)['Scores']==maxi][0]
|
|
|
479 new_inds.append(i)
|
|
|
480
|
|
|
481 output = test.loc[new_inds,:]
|
|
|
482 strands = []
|
|
|
483 new_pos = []
|
|
|
484 old_pos = output['Positions'].tolist()
|
|
|
485
|
|
|
486 from Bio.SeqFeature import SeqFeature, FeatureLocation
|
|
|
487 feats = rec.features
|
|
|
488 for ind, row in output.iterrows():
|
|
|
489 pos = row['Positions']
|
|
|
490 if 'complement' in pos:
|
|
|
491 strands.append('-')
|
|
|
492 new_pos.append(pos[10:])
|
|
|
493 ini,end= pos[11:-1].split('..')
|
|
|
494 new_loc = FeatureLocation(int(ini),int(end),strand=-1)
|
|
|
495 else:
|
|
|
496 strands.append('+')
|
|
|
497 new_pos.append(pos)
|
|
|
498 ini,end= pos[1:-1].split('..')
|
|
|
499 new_loc = FeatureLocation(int(ini),int(end),strand=1)
|
|
|
500 feat = SeqFeature(new_loc, type='regulatory',qualifiers={'regulatory_class':['promoter'], 'note=':['predicted by PhagePromoter']})
|
|
|
501 feats.append(feat)
|
|
|
502
|
|
|
503 output.insert(loc=0, column='Strand', value=strands)
|
|
|
504 output['Positions'] = new_pos
|
|
|
505 output.to_html('output.html',index=False, justify='center')
|
|
|
506 recs = []
|
|
|
507 i = 0
|
|
|
508 for ind,row in output.iterrows():
|
|
|
509 s = Seq(row['Promoter Sequence'])
|
|
|
510 posis = old_pos[i]
|
|
|
511 typ = row['Type']
|
|
|
512 score = row['Scores']
|
|
|
513 sq = SeqRecord(seq=s, id=ind, description=typ+' '+posis+' score='+str(score))
|
|
|
514 recs.append(sq)
|
|
|
515 i += 1
|
|
|
516 SeqIO.write(recs, 'output.fasta','fasta')
|
|
|
517 new_rec = rec
|
|
|
518 new_rec.seq.alphabet = IUPAC.IUPACAmbiguousDNA()
|
|
|
519 new_feats = sorted(feats, key=lambda x: x.location.start)
|
|
|
520 new_rec.features = new_feats
|
|
|
521 SeqIO.write(new_rec,'output.gb','genbank')
|
|
|
522
|
|
|
523 if __name__== "__main__":
|
|
|
524
|
|
|
525 import sys
|
|
|
526 import os
|
|
|
527 __location__ = os.path.realpath(os.path.join(os.getcwd(), os.path.dirname(__file__)))
|
|
|
528
|
|
|
529 gen_format = sys.argv[1]
|
|
|
530 genome_file = sys.argv[2]
|
|
|
531 both = sys.argv[3]
|
|
|
532 threshold = sys.argv[4]
|
|
|
533 family = sys.argv[5]
|
|
|
534 host = sys.argv[6]
|
|
|
535 phage_type = sys.argv[7]
|
|
|
536 model = sys.argv[8]
|
|
|
537 '''
|
|
|
538 gen_format = 'genbank'
|
|
|
539 genome_file = 'test-data/NC_015264.gb'
|
|
|
540 both = False
|
|
|
541 threshold = '0.50'
|
|
|
542 family = 'Podoviridae'
|
|
|
543 host = 'Pseudomonas'
|
|
|
544 phage_type = 'virulent'
|
|
|
545 model = 'SVM2400'
|
|
|
546 #model = 'ANN1600'
|
|
|
547 '''
|
|
|
548
|
|
|
549 rec = SeqIO.read(genome_file, gen_format)
|
|
|
550 test_windows = get_testseqs65(gen_format, genome_file,both)
|
|
|
551 try: score_test,dic_window = get_testScores(__location__,test_windows)
|
|
|
552 except IndexError: print('Error. Input sequence can only have A,C,G or T')
|
|
|
553 else:
|
|
|
554 df_test,df_testinfo = create_dftest(score_test,dic_window,family,host,phage_type)
|
|
|
555 if model == 'ANN1600':
|
|
|
556 scaler_file = os.path.join(__location__, 'scaler1600.sav')
|
|
|
557 model_file = os.path.join(__location__, 'model1600.sav')
|
|
|
558 preds = get_predictions(scaler_file, model_file, df_test,df_testinfo,threshold)
|
|
|
559 if preds is None: print('There is no sequence with a score value higher or equal to the threshold '+str(threshold))
|
|
|
560 elif type(preds) == str: print(preds)
|
|
|
561 else: output = get_finaldf(preds,rec)
|
|
|
562 else:
|
|
|
563 scaler_file = os.path.join(__location__, 'scaler2400.sav')
|
|
|
564 model_file = os.path.join(__location__, 'model2400.sav')
|
|
|
565 new_df_test = df_test.iloc[:,[0,1,2,3,4,5,6,7,8,9,13,14,16,17,19,20,22,24,25]]
|
|
|
566 preds = get_predictions(scaler_file, model_file, new_df_test,df_testinfo,threshold)
|
|
|
567 if preds is None: print('There is no sequence with a score value higher or equal to the threshold '+str(threshold))
|
|
|
568 elif type(preds) == str: print(preds)
|
|
|
569 else: output = get_finaldf(preds,rec)
|
|
|
570
|
