|
0
|
1 #! /usr/bin/env python
|
|
|
2 # -*- coding: utf-8 -*-
|
|
|
3 ## @file functions_PhageTerm.py
|
|
|
4 #
|
|
|
5 # This file is a part of PhageTerm software
|
|
|
6 # A tool to determine phage termini and packaging strategy
|
|
|
7 # and other useful informations using raw sequencing reads.
|
|
|
8 # (This programs works with sequencing reads from a randomly
|
|
|
9 # sheared DNA library preparations as Illumina TruSeq paired-end or similar)
|
|
|
10 #
|
|
|
11 # ----------------------------------------------------------------------
|
|
|
12 # Copyright (C) 2017 Julian Garneau
|
|
|
13 #
|
|
|
14 # This program is free software; you can redistribute it and/or modify
|
|
|
15 # it under the terms of the GNU General Public License as published by
|
|
|
16 # the Free Software Foundation; either version 3 of the License, or
|
|
|
17 # (at your option) any later version.
|
|
|
18 # <http://www.gnu.org/licenses/gpl-3.0.html>
|
|
|
19 #
|
|
|
20 # This program is distributed in the hope that it will be useful,
|
|
|
21 # but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
22 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
23 # GNU General Public License for more details.
|
|
|
24 # ----------------------------------------------------------------------
|
|
|
25 #
|
|
|
26 # @author Julian Garneau <julian.garneau@usherbrooke.ca>
|
|
|
27 # @author Marc Monot <marc.monot@pasteur.fr>
|
|
|
28 # @author David Bikard <david.bikard@pasteur.fr>
|
|
|
29
|
|
|
30
|
|
|
31 ### PYTHON Module
|
|
|
32 # Base
|
|
|
33 from __future__ import print_function
|
|
|
34
|
|
|
35 import sys
|
|
|
36
|
|
|
37 import os
|
|
|
38
|
|
|
39
|
|
|
40 import matplotlib
|
|
|
41 matplotlib.use('Agg')
|
|
|
42 import matplotlib.pyplot as plt
|
|
|
43 from matplotlib import patches
|
|
|
44 from matplotlib.path import Path
|
|
|
45
|
|
|
46 import numpy as np
|
|
|
47 import pandas as pd
|
|
|
48
|
|
|
49 # String
|
|
|
50 #import cStringIO
|
|
|
51 import io
|
|
|
52
|
|
|
53 # PDF report building
|
|
|
54 import time
|
|
|
55 from reportlab.lib.enums import TA_JUSTIFY, TA_CENTER, TA_LEFT, TA_RIGHT
|
|
|
56 from reportlab.lib.pagesizes import letter, landscape
|
|
|
57 from reportlab.platypus import SimpleDocTemplate, Paragraph, Spacer, Image, Table, TableStyle, PageBreak
|
|
|
58 from reportlab.lib.styles import getSampleStyleSheet, ParagraphStyle
|
|
|
59 from reportlab.lib.units import inch
|
|
|
60 from reportlab.lib import colors
|
|
|
61 from reportlab.lib.utils import ImageReader
|
|
|
62
|
|
|
63 from _modules.utilities import reverseComplement,hybridCoverage,applyCoverage,correctEdge
|
|
|
64 from _modules.common_readsCoverage_processing import picMax
|
|
|
65 from _modules.readsCoverage_res import RCRes, RCCheckpoint_handler,RCWorkingS
|
|
|
66
|
|
|
67 ### UTILITY function
|
|
|
68 def chunks(l, n):
|
|
|
69 """Yield n successive chunks from l."""
|
|
|
70 newn = int(1.0 * len(l) / n + 0.5)
|
|
|
71 for i in range(0, n-1):
|
|
|
72 yield l[i*newn:i*newn+newn]
|
|
|
73 yield l[n*newn-newn:]
|
|
|
74
|
|
|
75 ##
|
|
|
76 # Initializes working structure for readsCoverage
|
|
|
77 def init_ws(p_res,refseq,hostseq):
|
|
|
78 gen_len = len(refseq)
|
|
|
79 host_len = len(hostseq)
|
|
|
80 k = count_line = 0
|
|
|
81 if p_res==None:
|
|
|
82 termini_coverage = np.array([gen_len*[0], gen_len*[0]])
|
|
|
83 whole_coverage = np.array([gen_len*[0], gen_len*[0]])
|
|
|
84 paired_whole_coverage = np.array([gen_len*[0], gen_len*[0]])
|
|
|
85 phage_hybrid_coverage = np.array([gen_len*[0], gen_len*[0]])
|
|
|
86 host_hybrid_coverage = np.array([host_len*[0], host_len*[0]])
|
|
|
87 host_whole_coverage = np.array([host_len*[0], host_len*[0]])
|
|
|
88 list_hybrid = np.array([0,0])
|
|
|
89 insert = []
|
|
|
90 paired_missmatch = 0
|
|
|
91 read_match = 0
|
|
|
92 else:
|
|
|
93 termini_coverage=p_res.interm_res.termini_coverage
|
|
|
94 whole_coverage=p_res.interm_res.whole_coverage
|
|
|
95 paired_whole_coverage=p_res.interm_res.paired_whole_coverage
|
|
|
96 phage_hybrid_coverage=p_res.interm_res.phage_hybrid_coverage
|
|
|
97 host_hybrid_coverage=p_res.interm_res.host_hybrid_coverage
|
|
|
98 host_whole_coverage=p_res.interm_res.host_whole_coverage
|
|
|
99 list_hybrid=p_res.interm_res.list_hybrid
|
|
|
100 insert=p_res.interm_res.insert
|
|
|
101 paired_missmatch=p_res.interm_res.paired_mismatch
|
|
|
102 k=int(p_res.interm_res.reads_tested)
|
|
|
103 #count_line=p_res.count_line-1 # do that because readsCoverage will start by incrementing it of 1
|
|
|
104 read_match=p_res.read_match
|
|
|
105 return gen_len,host_len,termini_coverage,whole_coverage,paired_whole_coverage,phage_hybrid_coverage,host_hybrid_coverage, \
|
|
|
106 host_whole_coverage,list_hybrid,insert,paired_missmatch,k,count_line,read_match #TODO refactor that.
|
|
|
107
|
|
|
108
|
|
|
109
|
|
|
110 ## COVERAGE Starting and Whole function
|
|
|
111 #
|
|
|
112 # VL: use debug mode to keep track of what reads matched and what reads didn't. For those who matched, want to know if it is at the beginning of the read or at the end or if it is its reverse complement.
|
|
|
113 # My aim is to compare the results with those of the GPU version.
|
|
|
114 def readsCoverage(inRawDArgs,refseq,inDArgs,fParms,return_dict, core_id,line_start,line_end,tParms,\
|
|
|
115 chk_handler,idx_refseq,logger=None):
|
|
|
116 """Calculate whole coverage and first base coverage. """
|
|
|
117
|
|
|
118 p_res=chk_handler.load(core_id,idx_refseq)
|
|
|
119 gen_len,host_len,termini_coverage, whole_coverage, paired_whole_coverage, phage_hybrid_coverage, host_hybrid_coverage,\
|
|
|
120 host_whole_coverage, list_hybrid, insert, paired_missmatch, k, count_line, read_match=init_ws(p_res, refseq, inDArgs.hostseq)
|
|
|
121 if logger!=None:
|
|
|
122 logger.add_rw(p_res)
|
|
|
123 test_read_seq = match = 0
|
|
|
124 # Timer
|
|
|
125 if core_id == (tParms.core-1):
|
|
|
126 sys.stdout.write(" 0.0 %")
|
|
|
127 sys.stdout.flush()
|
|
|
128
|
|
|
129 # Mapping
|
|
|
130 filin = open(inRawDArgs.fastq)
|
|
|
131 line = filin.readline()
|
|
|
132
|
|
|
133 if inRawDArgs.paired != "":
|
|
|
134 filin_paired = open(inRawDArgs.paired)
|
|
|
135 line_paired = filin_paired.readline()
|
|
|
136 count_line_tmp=0
|
|
|
137 while line and count_line!=count_line_tmp:
|
|
|
138 count_line_tmp += 1
|
|
|
139 line = filin.readline()
|
|
|
140 while line:
|
|
|
141 count_line+=1
|
|
|
142 if count_line//4 <= line_start:
|
|
|
143 test_read_seq = 0
|
|
|
144 if count_line//4 > line_end:
|
|
|
145 break
|
|
|
146
|
|
|
147 if test_read_seq:
|
|
|
148 k += 1
|
|
|
149 # Read sequence
|
|
|
150 read = line.split("\n")[0].split("\r")[0]
|
|
|
151 line = filin.readline()
|
|
|
152
|
|
|
153 if inRawDArgs.paired != "":
|
|
|
154 read_paired = line_paired.split("\n")[0].split("\r")[0]
|
|
|
155 line_paired = filin_paired.readline()
|
|
|
156
|
|
|
157 readlen = len(read)
|
|
|
158 if readlen < fParms.seed:
|
|
|
159 if logger!=None:
|
|
|
160 print("CPU rejecting read",k)
|
|
|
161 continue
|
|
|
162 corlen = readlen-fParms.seed
|
|
|
163
|
|
|
164 if logger!=None:
|
|
|
165 print("CPU processing read: ",k,read, reverseComplement(read))
|
|
|
166 logger.newRmInfo(k)
|
|
|
167
|
|
|
168 ### Match sense + (multiple, random pick one)
|
|
|
169 #print("read[:fParms.seed]=",read[:fParms.seed])
|
|
|
170 matchPplus_start, matchPplus_end = applyCoverage(read[:fParms.seed], refseq)
|
|
|
171
|
|
|
172 ## Phage
|
|
|
173 if matchPplus_start != -1:
|
|
|
174 if logger!=None:
|
|
|
175 print("CPU found: ",read[:fParms.seed])
|
|
|
176 logger.rMatch("mstart")
|
|
|
177 match = 1
|
|
|
178 termini_coverage[0][matchPplus_start]+=1
|
|
|
179 position_end = matchPplus_end+corlen
|
|
|
180
|
|
|
181 # whole coverage
|
|
|
182 for i in range(matchPplus_start, min(gen_len,position_end)):
|
|
|
183 whole_coverage[0][i]+=1
|
|
|
184
|
|
|
185 # Paired-read
|
|
|
186 if inRawDArgs.paired != "":
|
|
|
187 matchPplus_start_paired, matchPplus_end_paired = applyCoverage(reverseComplement(read_paired)[-fParms.seed:], refseq)
|
|
|
188 insert_length = matchPplus_end_paired - matchPplus_start
|
|
|
189 if insert_length > 0 and insert_length < fParms.insert_max:
|
|
|
190 position_end = matchPplus_end_paired
|
|
|
191 insert.append(insert_length)
|
|
|
192 else:
|
|
|
193 paired_missmatch += 1
|
|
|
194 # Paired hybrid
|
|
|
195 if inDArgs.hostseq != "" and matchPplus_start_paired == -1:
|
|
|
196 matchHplus_start, matchHplus_end = applyCoverage(read_paired[:fParms.seed], inDArgs.hostseq)
|
|
|
197 if matchHplus_start != -1:
|
|
|
198 list_hybrid[0] += 1
|
|
|
199 phage_hybrid_coverage[0] = hybridCoverage(read, refseq, phage_hybrid_coverage[0], matchPplus_start, min(gen_len,matchPplus_end+corlen) )
|
|
|
200 host_hybrid_coverage[0] = hybridCoverage(read_paired, inDArgs.hostseq, host_hybrid_coverage[0], matchHplus_start, min(host_len,matchHplus_end+corlen) )
|
|
|
201 else:
|
|
|
202 matchHminus_start, matchHminus_end = applyCoverage(reverseComplement(read_paired)[:fParms.seed], inDArgs.hostseq)
|
|
|
203 if matchHminus_start != -1:
|
|
|
204 list_hybrid[0] += 1
|
|
|
205 phage_hybrid_coverage[0] = hybridCoverage(read, refseq, phage_hybrid_coverage[0], matchPplus_start, min(gen_len,matchPplus_end+corlen) )
|
|
|
206 host_hybrid_coverage[1] = hybridCoverage(reverseComplement(read_paired), inDArgs.hostseq, host_hybrid_coverage[1], matchHminus_start, min(host_len,matchHminus_end+corlen) )
|
|
|
207
|
|
|
208 # Single hybrid
|
|
|
209 elif inDArgs.hostseq != "":
|
|
|
210 matchPFplus_start, matchPFplus_end = applyCoverage(read[-fParms.seed:], refseq)
|
|
|
211 if matchPFplus_start == -1:
|
|
|
212 matchHplus_start, matchHplus_end = applyCoverage(read[-fParms.seed:], inDArgs.hostseq)
|
|
|
213 if matchHplus_start != -1:
|
|
|
214 list_hybrid[0] += 1
|
|
|
215 phage_hybrid_coverage[0] = hybridCoverage(read, refseq, phage_hybrid_coverage[0], matchPplus_start, min(gen_len,matchPplus_end+corlen) )
|
|
|
216 host_hybrid_coverage[0] = hybridCoverage(read, inDArgs.hostseq, host_hybrid_coverage[0], matchHplus_start, min(host_len,matchHplus_end+corlen) )
|
|
|
217 else:
|
|
|
218 matchHminus_start, matchHminus_end = applyCoverage(reverseComplement(read)[-fParms.seed:], inDArgs.hostseq)
|
|
|
219 if matchHminus_start != -1:
|
|
|
220 list_hybrid[0] += 1
|
|
|
221 phage_hybrid_coverage[0] = hybridCoverage(read, refseq, phage_hybrid_coverage[0], matchPplus_start, min(gen_len,matchPplus_end+corlen) )
|
|
|
222 host_hybrid_coverage[1] = hybridCoverage(reverseComplement(read), inDArgs.hostseq, host_hybrid_coverage[1], matchHminus_start, min(host_len,matchHminus_end+corlen) )
|
|
|
223
|
|
|
224 # whole coverage
|
|
|
225 for i in range(matchPplus_start, min(gen_len,position_end)):
|
|
|
226 paired_whole_coverage[0][i]+=1
|
|
|
227
|
|
|
228 ### if no match sense +, then test sense -
|
|
|
229 if not match:
|
|
|
230 matchPminus_start, matchPminus_end = applyCoverage(reverseComplement(read)[-fParms.seed:], refseq)
|
|
|
231
|
|
|
232 ## Phage
|
|
|
233 if matchPminus_end != -1:
|
|
|
234 if logger != None:
|
|
|
235 print("CPU found: ",reverseComplement(read)[-fParms.seed:]," from ",reverseComplement(read))
|
|
|
236 logger.rMatch("mrcplstart")
|
|
|
237 match = 1
|
|
|
238 termini_coverage[1][matchPminus_end-1]+=1
|
|
|
239 position_start = matchPminus_start-corlen
|
|
|
240
|
|
|
241 # whole coverage
|
|
|
242 for i in range(max(0,position_start), matchPminus_end):
|
|
|
243 whole_coverage[1][i]+=1
|
|
|
244
|
|
|
245 # Paired-read
|
|
|
246 if inRawDArgs.paired != "":
|
|
|
247 matchPminus_start_paired, matchPminus_end_paired = applyCoverage(read_paired[:fParms.seed], refseq)
|
|
|
248 insert_length = matchPminus_end - matchPminus_start_paired
|
|
|
249 if insert_length > 0 and insert_length < fParms.insert_max:
|
|
|
250 position_start = matchPminus_start_paired
|
|
|
251 insert.append(insert_length)
|
|
|
252 else:
|
|
|
253 paired_missmatch += 1
|
|
|
254 # Paired hybrid
|
|
|
255 if inDArgs.hostseq != "" and matchPminus_start_paired == -1:
|
|
|
256 matchHplus_start, matchHplus_end = applyCoverage(read_paired[:fParms.seed], inDArgs.hostseq)
|
|
|
257 if matchHplus_start != -1:
|
|
|
258 list_hybrid[1] += 1
|
|
|
259 phage_hybrid_coverage[1] = hybridCoverage(reverseComplement(read), refseq, phage_hybrid_coverage[1], matchPminus_start, min(gen_len,matchPminus_end+corlen) )
|
|
|
260 host_hybrid_coverage[0] = hybridCoverage(read_paired, inDArgs.hostseq, host_hybrid_coverage[0], matchHplus_start, min(host_len,matchHplus_end+corlen) )
|
|
|
261
|
|
|
262 else:
|
|
|
263 matchHminus_start, matchHminus_end = applyCoverage(reverseComplement(read_paired)[-fParms.seed:], inDArgs.hostseq)
|
|
|
264 if matchHminus_start != -1:
|
|
|
265 list_hybrid[1] += 1
|
|
|
266 phage_hybrid_coverage[1] = hybridCoverage(reverseComplement(read), refseq, phage_hybrid_coverage[1], matchPminus_start, min(gen_len,matchPminus_end+corlen) )
|
|
|
267 host_hybrid_coverage[1] = hybridCoverage(reverseComplement(read_paired), inDArgs.hostseq, host_hybrid_coverage[1], matchHminus_start, min(host_len,matchHminus_end+corlen) )
|
|
|
268
|
|
|
269 # Single hybrid
|
|
|
270 elif inDArgs.hostseq != "":
|
|
|
271 matchPRplus_start, matchPRplus_end = applyCoverage(reverseComplement(read)[:fParms.seed], refseq)
|
|
|
272 if matchPRplus_start == -1:
|
|
|
273 matchHplus_start, matchHplus_end = applyCoverage(read[:fParms.seed], inDArgs.hostseq)
|
|
|
274 if matchHplus_start != -1:
|
|
|
275 list_hybrid[1] += 1
|
|
|
276 phage_hybrid_coverage[1] = hybridCoverage(reverseComplement(read), refseq, phage_hybrid_coverage[1], matchPminus_start, min(gen_len,matchPminus_end+corlen) )
|
|
|
277 host_hybrid_coverage[0] = hybridCoverage(read, inDArgs.hostseq, host_hybrid_coverage[0], matchHplus_start, min(host_len,matchHplus_end+corlen) )
|
|
|
278 else:
|
|
|
279 matchHminus_start, matchHminus_end = applyCoverage(reverseComplement(read)[:fParms.seed], inDArgs.hostseq)
|
|
|
280 if matchHminus_start != -1:
|
|
|
281 list_hybrid[1] += 1
|
|
|
282 phage_hybrid_coverage[1] = hybridCoverage(reverseComplement(read), refseq, phage_hybrid_coverage[1], matchPminus_start, min(gen_len,matchPminus_end+corlen) )
|
|
|
283 host_hybrid_coverage[1] = hybridCoverage(reverseComplement(read), inDArgs.hostseq, host_hybrid_coverage[1], matchHminus_start, min(host_len,matchHminus_end+corlen) )
|
|
|
284
|
|
|
285 # whole coverage
|
|
|
286 for i in range(max(0,position_start), matchPminus_end):
|
|
|
287 paired_whole_coverage[1][i]+=1
|
|
|
288
|
|
|
289 ### if no match on Phage, test Host
|
|
|
290 if not match:
|
|
|
291 matchHplus_start, matchHplus_end = applyCoverage(read[:fParms.seed], inDArgs.hostseq)
|
|
|
292 if matchHplus_start != -1:
|
|
|
293 for i in range(matchHplus_start, min(host_len,matchHplus_end+corlen)):
|
|
|
294 host_whole_coverage[0][i]+=1
|
|
|
295 else:
|
|
|
296 matchHminus_start, matchHminus_end = applyCoverage(reverseComplement(read)[-fParms.seed:], inDArgs.hostseq)
|
|
|
297 if matchHminus_end != -1:
|
|
|
298 for i in range(max(0,matchHminus_start-corlen), matchHminus_end):
|
|
|
299 host_whole_coverage[1][i]+=1
|
|
|
300
|
|
|
301 # TEST limit_coverage
|
|
|
302 read_match += match*readlen
|
|
|
303
|
|
|
304 match = test_read_seq = 0
|
|
|
305 # Timer
|
|
|
306 if core_id == (tParms.core-1):
|
|
|
307 if k%1000 == 0:
|
|
|
308 sys.stdout.write("\b\b\b\b\b\b\b\b\b%3.1f %%" %( float(read_match/gen_len) / tParms.limit_coverage * 100 ))
|
|
|
309 sys.stdout.flush()
|
|
|
310
|
|
|
311 chk_handler.check(count_line,core_id,idx_refseq,termini_coverage,whole_coverage,paired_whole_coverage,\
|
|
|
312 phage_hybrid_coverage, host_hybrid_coverage, \
|
|
|
313 host_whole_coverage,list_hybrid,insert,paired_missmatch,k,read_match) # maybe time to create checkpoint
|
|
|
314
|
|
|
315 else:
|
|
|
316 if line[0] == "@":
|
|
|
317 test_read_seq = 1
|
|
|
318
|
|
|
319 line = filin.readline()
|
|
|
320 if inRawDArgs.paired != "":
|
|
|
321 line_paired = filin_paired.readline()
|
|
|
322
|
|
|
323 # TEST limit_coverage
|
|
|
324 if (read_match/gen_len) > tParms.limit_coverage:
|
|
|
325 line = 0
|
|
|
326
|
|
|
327
|
|
|
328 if core_id == (tParms.core-1):
|
|
|
329 sys.stdout.write("\b\b\b\b\b\b\b\b\b%3.1f %%" %( 100 ))
|
|
|
330 sys.stdout.flush()
|
|
|
331
|
|
|
332 # Close file
|
|
|
333 filin.close()
|
|
|
334 if inRawDArgs.paired != "":
|
|
|
335 filin_paired.close()
|
|
|
336
|
|
|
337
|
|
|
338 # Correct EDGE coverage
|
|
|
339 if sum(termini_coverage[0]) + sum(termini_coverage[1]) == 0 and not fParms.virome:
|
|
|
340 print("WARNING: No read Match, please check your fastq file")
|
|
|
341
|
|
|
342 termini_coverage = correctEdge(termini_coverage, fParms.edge)
|
|
|
343 #paired_whole_coverage = correctEdge(whole_coverage, fParms.edge) #TODO: discuss with Julian and Max about the PE issue that Max reported.
|
|
|
344 whole_coverage = correctEdge(whole_coverage, fParms.edge)
|
|
|
345 phage_hybrid_coverage = correctEdge(phage_hybrid_coverage, fParms.edge)
|
|
|
346 if inDArgs.hostseq != "":
|
|
|
347 host_whole_coverage = correctEdge(host_whole_coverage, fParms.edge)
|
|
|
348 host_hybrid_coverage = correctEdge(host_hybrid_coverage, fParms.edge)
|
|
|
349
|
|
|
350 if return_dict!=None and tParms.dir_cov_mm==None:
|
|
|
351 return_dict[core_id] = [termini_coverage, whole_coverage, paired_whole_coverage, phage_hybrid_coverage, host_hybrid_coverage, host_whole_coverage, list_hybrid, np.array(insert), paired_missmatch, k]
|
|
|
352 elif return_dict==None and tParms.dir_cov_mm!=None:
|
|
|
353 insert = np.array(insert)
|
|
|
354 fic_name = os.path.join(tParms.dir_cov_mm, "coverage" + str(idx_refseq) + "_" + str(core_id))
|
|
|
355 res=RCRes(termini_coverage,whole_coverage,paired_whole_coverage,\
|
|
|
356 phage_hybrid_coverage, host_hybrid_coverage, \
|
|
|
357 host_whole_coverage,list_hybrid,insert,paired_missmatch,k)
|
|
|
358 res.save(fic_name)
|
|
|
359 else:
|
|
|
360 print("Error: readsCoverage must be used either with directory name or return_dict")
|
|
|
361 chk_handler.end(core_id)
|
|
|
362
|
|
|
363 return
|
|
|
364
|
|
|
365
|
|
|
366
|
|
|
367 ### IMAGE Functions
|
|
|
368 def GraphCov(termini_coverage, picMaxPlus, picMaxMinus, phagename, norm, draw, hybrid = 0):
|
|
|
369 """Produces a plot with termini coverage values."""
|
|
|
370 # Remove old plot
|
|
|
371 plt.clf()
|
|
|
372 plt.cla()
|
|
|
373 plt.close()
|
|
|
374 # Create figure
|
|
|
375 plt.figure(1)
|
|
|
376 term_len = len(termini_coverage[0])
|
|
|
377 term_range = list(range(term_len))
|
|
|
378 # MOP: not totally sure what's going on here with the plot formatting
|
|
|
379 # but I refactored this out as it was getting complicated.
|
|
|
380 # Someone who understands the code better might put in more informative var names.
|
|
|
381 zipper = list(zip(*picMaxPlus))
|
|
|
382 max_first_zipper = max(np.array(zipper[0]))
|
|
|
383 if norm == 1:
|
|
|
384 ylim = 1.2
|
|
|
385 elif hybrid == 1:
|
|
|
386 offset = 0.2*(max_first_zipper) + 1
|
|
|
387 ylim = max_first_zipper + offset
|
|
|
388 else:
|
|
|
389 offset = 0.2*(max(max(np.array(list(zip(*picMaxPlus))[0])), max(np.array(list(zip(*picMaxMinus))[0]))))
|
|
|
390 ylim = max(max(np.array(list(zip(*picMaxPlus))[0])), max(np.array(list(zip(*picMaxMinus))[0]))) + offset
|
|
|
391 # Strand (+)
|
|
|
392 plt.subplot(211)
|
|
|
393 if norm == 1:
|
|
|
394 plt.scatter(term_range,termini_coverage[0])
|
|
|
395 else:
|
|
|
396 plt.plot(termini_coverage[0],linewidth=2)
|
|
|
397 plt.title('strand (+)')
|
|
|
398 plt.ylabel('')
|
|
|
399 # Axes
|
|
|
400 axes = plt.gca()
|
|
|
401 axes.set_ylim([0,ylim])
|
|
|
402 # Maximum
|
|
|
403 x_strandplus = np.array(list(zip(*picMaxPlus))[1])
|
|
|
404 y_strandplus = np.array(list(zip(*picMaxPlus))[0])
|
|
|
405 # Plot
|
|
|
406 plt.plot(x_strandplus, y_strandplus, 'ro')
|
|
|
407 if norm == 1:
|
|
|
408 axes.axhline(y=0.5, xmin=0, xmax=x_strandplus, color='grey', linestyle='dashed', linewidth=1.5)
|
|
|
409 axes.axhline(y=1.0, xmin=0, xmax=x_strandplus, color='grey', linestyle='dashed', linewidth=1.5)
|
|
|
410 # Annotation
|
|
|
411 for i,j in zip(x_strandplus,y_strandplus):
|
|
|
412 plt.text(i+(term_len*0.03), j, str(i+1), fontsize=15, bbox=dict(boxstyle='round', facecolor='white', alpha=1))
|
|
|
413 # Plot Option
|
|
|
414 plt.margins(0.1)
|
|
|
415 plt.locator_params(axis = 'x', nbins = 10)
|
|
|
416 plt.locator_params(axis = 'y', nbins = 3)
|
|
|
417 plt.xticks(rotation=75)
|
|
|
418 # Strand (-)
|
|
|
419 plt.subplot(212)
|
|
|
420 if norm == 1:
|
|
|
421 plt.scatter(term_range,termini_coverage[1])
|
|
|
422 else:
|
|
|
423 plt.plot(termini_coverage[1],linewidth=2)
|
|
|
424 plt.title('strand (-)')
|
|
|
425 plt.ylabel('')
|
|
|
426 # Axes
|
|
|
427 if hybrid == 1:
|
|
|
428 offset = 0.2*(max_first_zipper) + 1
|
|
|
429 ylim = max_first_zipper + offset
|
|
|
430 axes = plt.gca()
|
|
|
431 axes.set_ylim([0,ylim])
|
|
|
432 # Maximum
|
|
|
433 x_strandminus = np.array(list(zip(*picMaxMinus))[1])
|
|
|
434 y_strandminus = np.array(list(zip(*picMaxMinus))[0])
|
|
|
435 # Plot
|
|
|
436 plt.plot(x_strandminus, y_strandminus, 'ro')
|
|
|
437 if norm == 1:
|
|
|
438 axes.axhline(y=0.5, xmin=0, xmax=x_strandplus, color='grey', linestyle='dashed', linewidth=1.5)
|
|
|
439 axes.axhline(y=1.0, xmin=0, xmax=x_strandplus, color='grey', linestyle='dashed', linewidth=1.5)
|
|
|
440 # Annotation
|
|
|
441 for i,j in zip(x_strandminus,y_strandminus):
|
|
|
442 plt.text(i+(term_len*0.03), j, str(i+1), fontsize=15, bbox=dict(boxstyle='round', facecolor='white', alpha=1))
|
|
|
443 # Plot Option
|
|
|
444 plt.margins(0.1)
|
|
|
445 plt.locator_params(axis = 'x', nbins = 10)
|
|
|
446 plt.locator_params(axis = 'y', nbins = 3)
|
|
|
447 plt.xticks(rotation=75)
|
|
|
448 # Plot Adjustments
|
|
|
449 plt.tight_layout()
|
|
|
450 # Draw graph
|
|
|
451 if draw:
|
|
|
452 plt.savefig("%s_TCov.png" % phagename, dpi=200)
|
|
|
453 fig = plt.figure(1)
|
|
|
454 return fig
|
|
|
455
|
|
|
456 def GraphWholeCov(added_whole_coverage, phagename, draw, P_left = "", P_right = "", pos_left = 0, pos_right = 0, graphZoom = 0, title = "WHOLE COVERAGE"):
|
|
|
457 """Produces a plot with whole coverage values."""
|
|
|
458 # Remove old plot
|
|
|
459 plt.clf()
|
|
|
460 plt.cla()
|
|
|
461 plt.close()
|
|
|
462 # Create figure
|
|
|
463 offset = 0.2*(max(added_whole_coverage))
|
|
|
464 ylim = max(added_whole_coverage) + offset
|
|
|
465 # Cumulative both strands
|
|
|
466 plt.figure(figsize=(15,8))
|
|
|
467 plt.plot(added_whole_coverage,linewidth=2)
|
|
|
468 plt.title(title)
|
|
|
469 # Axes
|
|
|
470 axes = plt.gca()
|
|
|
471 axes.set_ylim([0,ylim])
|
|
|
472 # Plot Option
|
|
|
473 plt.margins(0.1)
|
|
|
474 plt.locator_params(axis = 'x', nbins = 10)
|
|
|
475 plt.xticks(rotation=75)
|
|
|
476 # Termini vertical dashed line display
|
|
|
477 if graphZoom and isinstance(P_left, np.integer):
|
|
|
478 plt.axvline(x=pos_left, ymin=0, ymax=ylim, color='red', zorder=2, linestyle='dashed', linewidth=1)
|
|
|
479 if graphZoom and isinstance(P_right, np.integer):
|
|
|
480 plt.axvline(x=pos_right, ymin=0, ymax=ylim, color='green', zorder=2, linestyle='dashed', linewidth=1)
|
|
|
481 # Draw graph
|
|
|
482 if draw:
|
|
|
483 plt.savefig("%s_plot_WCov.png" % phagename, dpi=200)
|
|
|
484 fig = plt.figure(1)
|
|
|
485 return fig
|
|
|
486
|
|
|
487 def GraphLogo(P_class, P_left, P_right, draw):
|
|
|
488 """Produce logo."""
|
|
|
489 # Remove old plot
|
|
|
490 plt.clf()
|
|
|
491 plt.cla()
|
|
|
492 plt.close()
|
|
|
493 # Create figure
|
|
|
494 plt.figure(figsize=(10,10))
|
|
|
495 #axes = plt.add_subplot(111)
|
|
|
496 axes = plt.gca()
|
|
|
497 axes.set_frame_on(False)
|
|
|
498 axes.xaxis.set_visible(False)
|
|
|
499 axes.yaxis.set_visible(False)
|
|
|
500 # Cadre
|
|
|
501 axes.add_artist(patches.Rectangle((0.1, 0.1), 0.8, 0.8, edgecolor = 'black', fill = False, linewidth = 15))
|
|
|
502
|
|
|
503 if P_class == "Headful (pac)":
|
|
|
504 # Texte
|
|
|
505 axes.text(0.17, 0.7, r"Headful (pac)", fontsize=50, fontweight='bold')
|
|
|
506 # PAC (blue line)
|
|
|
507 axes.axhline(y=0.35, xmin=0.2, xmax=0.8, color='blue', linewidth=15)
|
|
|
508 # PAC (red line)
|
|
|
509 axes.axvline(x=0.19, ymin=0.30, ymax=0.40, color='red', linewidth=10)
|
|
|
510 axes.axvline(x=0.42, ymin=0.30, ymax=0.40, color='red', linewidth=10)
|
|
|
511 axes.axvline(x=0.65, ymin=0.30, ymax=0.40, color='red', linewidth=10)
|
|
|
512 # PAC (Arrow)
|
|
|
513 axes.axvline(x=0.19, ymin=0.45, ymax=0.55, color='red', linewidth=15)
|
|
|
514 axes.arrow(0.19, 0.55, 0.07, 0, color='red', linewidth=15, head_width=0.07, head_length=0.1)
|
|
|
515
|
|
|
516 elif P_class == "COS (5')":
|
|
|
517 # Texte
|
|
|
518 axes.text(0.3, 0.7, r"COS (5')", fontsize=50, fontweight='bold')
|
|
|
519 axes.add_artist(patches.Ellipse(xy=(0.5,0.4), width=0.5, height=0.35 , edgecolor = 'blue', fill=False, lw=15))
|
|
|
520 axes.add_artist(patches.Ellipse(xy=(0.5,0.4), width=0.58, height=0.43 , edgecolor = 'blue', fill=False, lw=15))
|
|
|
521 axes.add_artist(patches.Rectangle((0.4, 0.5), 0.20, 0.20, edgecolor = 'white', facecolor = 'white', fill = True))
|
|
|
522 axes.axhline(y=0.56, xmin=0.415, xmax=0.48, color='red', linewidth=16)
|
|
|
523 axes.axhline(y=0.601, xmin=0.52, xmax=0.585, color='red', linewidth=16)
|
|
|
524
|
|
|
525 elif P_class == "COS (3')":
|
|
|
526 # Texte
|
|
|
527 axes.text(0.3, 0.7, r"COS (3')", fontsize=50, fontweight='bold')
|
|
|
528 axes.add_artist(patches.Ellipse(xy=(0.5,0.4), width=0.5, height=0.35 , edgecolor = 'blue', fill=False, lw=15))
|
|
|
529 axes.add_artist(patches.Ellipse(xy=(0.5,0.4), width=0.58, height=0.43 , edgecolor = 'blue', fill=False, lw=15))
|
|
|
530 axes.add_artist(patches.Rectangle((0.4, 0.5), 0.20, 0.20, edgecolor = 'white', facecolor = 'white', fill = True))
|
|
|
531 axes.axhline(y=0.601, xmin=0.415, xmax=0.48, color='red', linewidth=16)
|
|
|
532 axes.axhline(y=0.56, xmin=0.52, xmax=0.585, color='red', linewidth=16)
|
|
|
533
|
|
|
534 elif P_class == "COS":
|
|
|
535 # Texte
|
|
|
536 axes.text(0.4, 0.7, r"COS", fontsize=50, fontweight='bold')
|
|
|
537 axes.add_artist(patches.Ellipse(xy=(0.5,0.4), width=0.5, height=0.35 , edgecolor = 'blue', fill=False, lw=15))
|
|
|
538 axes.add_artist(patches.Ellipse(xy=(0.5,0.4), width=0.58, height=0.43 , edgecolor = 'blue', fill=False, lw=15))
|
|
|
539 axes.add_artist(patches.Rectangle((0.4, 0.5), 0.20, 0.20, edgecolor = 'white', facecolor = 'white', fill = True))
|
|
|
540
|
|
|
541 elif P_class == "DTR (short)":
|
|
|
542 # Texte
|
|
|
543 axes.text(0.22, 0.7, r"DTR (short)", fontsize=50, fontweight='bold')
|
|
|
544
|
|
|
545 verts = [(0.5, 0.5), (0.9, 0.4), (0.9, 0.3), (0.5,0.2)]
|
|
|
546 codes = [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]
|
|
|
547 path = Path(verts, codes)
|
|
|
548 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'blue', lw=15)
|
|
|
549 axes.add_patch(patch)
|
|
|
550
|
|
|
551 verts = [(0.5, 0.2), (0.1, 0.30), (0.1, 0.45), (0.5,0.55)]
|
|
|
552 codes = [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]
|
|
|
553 path = Path(verts, codes)
|
|
|
554 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'blue', lw=15)
|
|
|
555 axes.add_patch(patch)
|
|
|
556
|
|
|
557 verts = [(0.5, 0.55), (0.52, 0.545), (0, 0)]
|
|
|
558 codes = [Path.MOVETO, Path.LINETO, Path.CLOSEPOLY]
|
|
|
559 path = Path(verts, codes)
|
|
|
560 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'red', lw=15)
|
|
|
561 axes.add_patch(patch)
|
|
|
562
|
|
|
563 verts = [(0.56, 0.536), (0.58, 0.530), (0, 0)]
|
|
|
564 codes = [Path.MOVETO, Path.LINETO, Path.CLOSEPOLY]
|
|
|
565 path = Path(verts, codes)
|
|
|
566 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'red', lw=15)
|
|
|
567 axes.add_patch(patch)
|
|
|
568
|
|
|
569 verts = [(0.5, 0.50), (0.56, 0.480), (0, 0)]
|
|
|
570 codes = [Path.MOVETO, Path.LINETO, Path.CLOSEPOLY]
|
|
|
571 path = Path(verts, codes)
|
|
|
572 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'white', lw=20)
|
|
|
573 axes.add_patch(patch)
|
|
|
574
|
|
|
575 verts = [(0.5, 0.50), (0.52, 0.495), (0, 0)]
|
|
|
576 codes = [Path.MOVETO, Path.LINETO, Path.CLOSEPOLY]
|
|
|
577 path = Path(verts, codes)
|
|
|
578 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'red', lw=15)
|
|
|
579 axes.add_patch(patch)
|
|
|
580
|
|
|
581 verts = [(0.56, 0.486), (0.58, 0.480), (0, 0)]
|
|
|
582 codes = [Path.MOVETO, Path.LINETO, Path.CLOSEPOLY]
|
|
|
583 path = Path(verts, codes)
|
|
|
584 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'red', lw=15)
|
|
|
585 axes.add_patch(patch)
|
|
|
586
|
|
|
587 elif P_class == "DTR (long)":
|
|
|
588 # Texte
|
|
|
589 axes.text(0.25, 0.7, r"DTR (long)", fontsize=50, fontweight='bold')
|
|
|
590 verts = [(0.5, 0.5), (0.9, 0.4), (0.9, 0.3), (0.5,0.2)]
|
|
|
591 codes = [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]
|
|
|
592 path = Path(verts, codes)
|
|
|
593 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'blue', lw=15)
|
|
|
594 axes.add_patch(patch)
|
|
|
595
|
|
|
596 verts = [(0.5, 0.2), (0.1, 0.30), (0.1, 0.45), (0.5,0.55)]
|
|
|
597 codes = [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]
|
|
|
598 path = Path(verts, codes)
|
|
|
599 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'blue', lw=15)
|
|
|
600 axes.add_patch(patch)
|
|
|
601
|
|
|
602 verts = [(0.5, 0.55), (0.52, 0.545), (0, 0)]
|
|
|
603 codes = [Path.MOVETO, Path.LINETO, Path.CLOSEPOLY]
|
|
|
604 path = Path(verts, codes)
|
|
|
605 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'red', lw=15)
|
|
|
606 axes.add_patch(patch)
|
|
|
607
|
|
|
608 verts = [(0.56, 0.536), (0.58, 0.530), (0, 0)]
|
|
|
609 codes = [Path.MOVETO, Path.LINETO, Path.CLOSEPOLY]
|
|
|
610 path = Path(verts, codes)
|
|
|
611 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'red', lw=15)
|
|
|
612 axes.add_patch(patch)
|
|
|
613
|
|
|
614 verts = [(0.62, 0.521), (0.64, 0.516), (0, 0)]
|
|
|
615 codes = [Path.MOVETO, Path.LINETO, Path.CLOSEPOLY]
|
|
|
616 path = Path(verts, codes)
|
|
|
617 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'red', lw=15)
|
|
|
618 axes.add_patch(patch)
|
|
|
619
|
|
|
620 verts = [(0.68, 0.507), (0.70, 0.501), (0, 0)]
|
|
|
621 codes = [Path.MOVETO, Path.LINETO, Path.CLOSEPOLY]
|
|
|
622 path = Path(verts, codes)
|
|
|
623 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'red', lw=15)
|
|
|
624 axes.add_patch(patch)
|
|
|
625
|
|
|
626 verts = [(0.5, 0.50), (0.65, 0.460), (0, 0)]
|
|
|
627 codes = [Path.MOVETO, Path.LINETO, Path.CLOSEPOLY]
|
|
|
628 path = Path(verts, codes)
|
|
|
629 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'white', lw=25)
|
|
|
630 axes.add_patch(patch)
|
|
|
631
|
|
|
632 verts = [(0.5, 0.50), (0.52, 0.495), (0, 0)]
|
|
|
633 codes = [Path.MOVETO, Path.LINETO, Path.CLOSEPOLY]
|
|
|
634 path = Path(verts, codes)
|
|
|
635 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'red', lw=15)
|
|
|
636 axes.add_patch(patch)
|
|
|
637
|
|
|
638 verts = [(0.56, 0.486), (0.58, 0.480), (0, 0)]
|
|
|
639 codes = [Path.MOVETO, Path.LINETO, Path.CLOSEPOLY]
|
|
|
640 path = Path(verts, codes)
|
|
|
641 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'red', lw=15)
|
|
|
642 axes.add_patch(patch)
|
|
|
643
|
|
|
644 verts = [(0.62, 0.471), (0.64, 0.465), (0, 0)]
|
|
|
645 codes = [Path.MOVETO, Path.LINETO, Path.CLOSEPOLY]
|
|
|
646 path = Path(verts, codes)
|
|
|
647 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'red', lw=15)
|
|
|
648 axes.add_patch(patch)
|
|
|
649
|
|
|
650 verts = [(0.68, 0.456), (0.70, 0.450), (0, 0)]
|
|
|
651 codes = [Path.MOVETO, Path.LINETO, Path.CLOSEPOLY]
|
|
|
652 path = Path(verts, codes)
|
|
|
653 patch = patches.PathPatch(path, facecolor='none', edgecolor = 'red', lw=15)
|
|
|
654 axes.add_patch(patch)
|
|
|
655
|
|
|
656 elif P_class == "Mu-like":
|
|
|
657 # Texte
|
|
|
658 axes.text(0.33, 0.7, r"Mu-like", fontsize=50, fontweight='bold')
|
|
|
659 axes.axhline(y=0.43, xmin=0.3, xmax=0.7, color='blue', linewidth=15)
|
|
|
660 axes.axhline(y=0.47, xmin=0.3, xmax=0.7, color='blue', linewidth=15)
|
|
|
661 axes.axhline(y=0.43, xmin=0.7, xmax=0.8, color='green', linewidth=15)
|
|
|
662 axes.axhline(y=0.47, xmin=0.7, xmax=0.8, color='green', linewidth=15)
|
|
|
663 axes.axhline(y=0.43, xmin=0.2, xmax=0.3, color='green', linewidth=15)
|
|
|
664 axes.axhline(y=0.47, xmin=0.2, xmax=0.3, color='green', linewidth=15)
|
|
|
665
|
|
|
666 elif P_left == "Random" and P_right == "Random":
|
|
|
667 # Texte
|
|
|
668 axes.text(0.25, 0.7, r"UNKNOWN", fontsize=50, fontweight='bold')
|
|
|
669 axes.text(0.44, 0.3, r"?", fontsize=200, fontweight='bold')
|
|
|
670 else:
|
|
|
671 # Texte
|
|
|
672 axes.text(0.4, 0.7, r"NEW", fontsize=50, fontweight='bold')
|
|
|
673 axes.text(0.44, 0.3, r"!", fontsize=200, fontweight='bold')
|
|
|
674
|
|
|
675 # Draw graph
|
|
|
676 if draw:
|
|
|
677 plt.savefig("%s_logo.png" % phagename, dpi=200)
|
|
|
678 fig = plt.figure(1)
|
|
|
679 return fig
|
|
|
680
|
|
|
681
|
|
|
682 ### OUTPUT Result files
|
|
|
683 def exportDataSplit(sequence, split):
|
|
|
684 """Export sequence with split line length."""
|
|
|
685 seq = ""
|
|
|
686 for i in range((len(sequence)//split)+1):
|
|
|
687 seq += "".join(map(str,sequence[i*split:(i+1)*split])) + '\n'
|
|
|
688 return seq
|
|
|
689
|
|
|
690 def ExportStatistics(phagename, whole_coverage, paired_whole_coverage, termini_coverage, phage_plus_norm, phage_minus_norm, paired, test_run):
|
|
|
691 """Export peaks statistics."""
|
|
|
692 if test_run:
|
|
|
693 return
|
|
|
694 export = pd.DataFrame()
|
|
|
695 # ORGANIZE Column
|
|
|
696 export["Position"] = list(phage_plus_norm.sort_values("Position")["Position"])
|
|
|
697 if paired != "":
|
|
|
698 export["Coverage +"] = paired_whole_coverage[0]
|
|
|
699 else:
|
|
|
700 export["Coverage +"] = whole_coverage[0]
|
|
|
701 export["SPC +"] = termini_coverage[0]
|
|
|
702 export["T +"] = [format(x/100.0,'0.2') for x in list(phage_plus_norm.sort_values("Position")["SPC_std"])]
|
|
|
703 export["T + (close)"] = [format(x/100.0,'0.2') for x in list(phage_plus_norm.sort_values("Position")["SPC"])]
|
|
|
704 export["pvalue +"] = [format(x,'0.2e') for x in list(phage_plus_norm.sort_values("Position")["pval_gamma"])]
|
|
|
705 export["padj +"] = [format(x,'0.2e') for x in list(phage_plus_norm.sort_values("Position")["pval_gamma_adj"])]
|
|
|
706 if paired != "":
|
|
|
707 export["Coverage -"] = whole_coverage[1]
|
|
|
708 else:
|
|
|
709 export["Coverage -"] = paired_whole_coverage[1]
|
|
|
710 export["SPC -"] = termini_coverage[1]
|
|
|
711 export["T -"] = [format(x/100.0,'0.2') for x in list(phage_minus_norm.sort_values("Position")["SPC_std"])]
|
|
|
712 export["T - (close)"] = [format(x/100.0,'0.2') for x in list(phage_minus_norm.sort_values("Position")["SPC"])]
|
|
|
713 export["pvalue -"] = [format(x,'0.2e') for x in list(phage_minus_norm.sort_values("Position")["pval_gamma"])]
|
|
|
714 export["padj -"] = [format(x,'0.2e') for x in list(phage_minus_norm.sort_values("Position")["pval_gamma_adj"])]
|
|
|
715 filout = open(phagename + "_statistics.csv", "w")
|
|
|
716 filout.write(export.to_csv(index=False))
|
|
|
717 filout.close()
|
|
|
718 return
|
|
|
719
|
|
|
720 def ExportCohesiveSeq(phagename, ArtcohesiveSeq, P_seqcoh, test_run, multi = 0):
|
|
|
721 """Export cohesive sequence of COS phages."""
|
|
|
722 if test_run:
|
|
|
723 return ""
|
|
|
724 if len(ArtcohesiveSeq) < 3 and len(P_seqcoh) < 3:
|
|
|
725 return ""
|
|
|
726 if len(ArtcohesiveSeq) < 20 and len(P_seqcoh) < 20:
|
|
|
727 export_text = "cohesive sequence"
|
|
|
728 if not multi:
|
|
|
729 filout = open(phagename + "_cohesive-sequence.fasta", "w")
|
|
|
730 else:
|
|
|
731 export_text = "direct terminal repeats sequence"
|
|
|
732 if not multi:
|
|
|
733 filout = open(phagename + "_direct-term-repeats.fasta", "w")
|
|
|
734 if P_seqcoh != '':
|
|
|
735 if not multi:
|
|
|
736 filout.write(">" + phagename + " " + export_text + " (Analysis: Statistics)\n" + exportDataSplit(P_seqcoh, 60))
|
|
|
737 else:
|
|
|
738 return ">" + phagename + " " + export_text + " (Analysis: Statistics)\n" + exportDataSplit(P_seqcoh, 60)
|
|
|
739 if ArtcohesiveSeq != '':
|
|
|
740 if not multi:
|
|
|
741 filout.write(">" + phagename + " " + export_text + " (Analysis: Li)\n" + exportDataSplit(ArtcohesiveSeq, 60))
|
|
|
742 filout.close()
|
|
|
743 return ""
|
|
|
744
|
|
|
745 def ExportPhageSequence(phagename, P_left, P_right, refseq, P_orient, Redundant, Mu_like, P_class, P_seqcoh, test_run, multi = 0):
|
|
|
746 """Export the phage sequence reorganized and completed if needed."""
|
|
|
747 if test_run:
|
|
|
748 return ""
|
|
|
749 seq_out = ""
|
|
|
750 # Mu-like
|
|
|
751 if Mu_like:
|
|
|
752 if P_orient == "Forward":
|
|
|
753 if P_right != "Random":
|
|
|
754 if P_left > P_right:
|
|
|
755 seq_out = refseq[P_right-1:P_left-1]
|
|
|
756 else:
|
|
|
757 seq_out = refseq[P_right-1:] + refseq[:P_left-1]
|
|
|
758 else:
|
|
|
759 seq_out = refseq[P_left-1:] + refseq[:P_left-1]
|
|
|
760 elif P_orient == "Reverse":
|
|
|
761 if P_left != "Random":
|
|
|
762 if P_left > P_right:
|
|
|
763 seq_out = reverseComplement(refseq[P_right-1:P_left-1])
|
|
|
764 else:
|
|
|
765 seq_out = reverseComplement(refseq[P_right-1:] + reverseComplement(refseq[:P_left-1]))
|
|
|
766 else:
|
|
|
767 seq_out = reverseComplement(refseq[P_right-1:] + reverseComplement(refseq[:P_right-1]) )
|
|
|
768 # COS
|
|
|
769 elif isinstance(P_left, np.integer) and isinstance(P_right, np.integer):
|
|
|
770 # Cos or DTR
|
|
|
771 if P_class == "COS (3')":
|
|
|
772 if abs(P_left-P_right) > len(refseq)/2:
|
|
|
773 seq_out = refseq[min(P_left,P_right)-1:max(P_left,P_right)]
|
|
|
774 else:
|
|
|
775 seq_out = refseq[max(P_left,P_right)-1:] + refseq[:min(P_left,P_right)]
|
|
|
776 seq_out = seq_out + P_seqcoh
|
|
|
777 else:
|
|
|
778 # Genome
|
|
|
779 if abs(P_left-P_right) > len(refseq)/2:
|
|
|
780 seq_out = refseq[min(P_left,P_right)-1:max(P_left,P_right)]
|
|
|
781 else:
|
|
|
782 seq_out = refseq[max(P_left,P_right):] + refseq[:min(P_left,P_right)-1]
|
|
|
783 # COS 5'
|
|
|
784 if P_class == "COS (5')":
|
|
|
785 seq_out = P_seqcoh + seq_out
|
|
|
786 # DTR
|
|
|
787 else:
|
|
|
788 seq_out = P_seqcoh + seq_out + P_seqcoh
|
|
|
789 # PAC
|
|
|
790 elif isinstance(P_left, np.integer) or isinstance(P_right, np.integer):
|
|
|
791 if P_orient == "Reverse":
|
|
|
792 seq_out = reverseComplement(refseq[:P_right]) + reverseComplement(refseq[P_right:])
|
|
|
793 else:
|
|
|
794 seq_out = refseq[P_left-1:] + refseq[:P_left-1]
|
|
|
795 # Write Sequence
|
|
|
796 if multi:
|
|
|
797 return ">" + phagename + " sequence re-organized\n" + exportDataSplit(seq_out, 60)
|
|
|
798 else:
|
|
|
799 filout = open(phagename + "_sequence.fasta", "w")
|
|
|
800 filout.write(">" + phagename + " sequence re-organized\n" + exportDataSplit(seq_out, 60))
|
|
|
801 filout.close()
|
|
|
802 return ""
|
|
|
803
|
|
|
804 def CreateReport(phagename, seed, added_whole_coverage, draw, Redundant, P_left, P_right, Permuted, P_orient, termini_coverage_norm_close, picMaxPlus_norm_close, picMaxMinus_norm_close, gen_len, tot_reads, P_seqcoh, phage_plus_norm, phage_minus_norm, ArtPackmode, termini, forward, reverse, ArtOrient, ArtcohesiveSeq, termini_coverage_close, picMaxPlus_close, picMaxMinus_close, picOUT_norm_forw, picOUT_norm_rev, picOUT_forw, picOUT_rev, lost_perc, ave_whole_cov, R1, R2, R3, host, host_len, host_whole_coverage, picMaxPlus_host, picMaxMinus_host, surrounding, drop_cov, paired, insert, phage_hybrid_coverage, host_hybrid_coverage, added_paired_whole_coverage, Mu_like, test_run, P_class, P_type, P_concat, multi = 0, multiReport = 0, *args, **kwargs):
|
|
|
805 """Produce a PDF report."""
|
|
|
806 if not multi:
|
|
|
807 doc = SimpleDocTemplate("%s_PhageTerm_report.pdf" % phagename, pagesize=letter, rightMargin=10,leftMargin=10, topMargin=5, bottomMargin=10)
|
|
|
808 report=[]
|
|
|
809 else:
|
|
|
810 report = multiReport
|
|
|
811
|
|
|
812 styles=getSampleStyleSheet()
|
|
|
813 styles.add(ParagraphStyle(name='Justify', alignment=TA_JUSTIFY))
|
|
|
814 styles.add(ParagraphStyle(name='Center', alignment=TA_CENTER))
|
|
|
815 styles.add(ParagraphStyle(name='Right', alignment=TA_RIGHT))
|
|
|
816 styles.add(ParagraphStyle(name='Left', alignment=TA_LEFT))
|
|
|
817
|
|
|
818 ### GENERAL INFORMATION
|
|
|
819
|
|
|
820 # TITLE
|
|
|
821 ptext = '<b><font size=16>' + phagename + ' PhageTerm Analysis</font></b>'
|
|
|
822 report.append(Paragraph(ptext, styles["Center"]))
|
|
|
823 report.append(Spacer(1, 15))
|
|
|
824
|
|
|
825 ## ZOOMED TERMINI GRAPH AND LOGO RESULT
|
|
|
826
|
|
|
827 # LOGO SLECTION
|
|
|
828
|
|
|
829 imgdata = io.BytesIO()
|
|
|
830 fig_logo = GraphLogo(P_class, P_left, P_right, draw)
|
|
|
831 fig_logo.savefig(imgdata, format='png')
|
|
|
832 imgdata.seek(0)
|
|
|
833 IMG = ImageReader(imgdata)
|
|
|
834 IMAGE_2 = Image(IMG.fileName, width=150, height=150, kind='proportional')
|
|
|
835 IMAGE_2.hAlign = 'CENTER'
|
|
|
836
|
|
|
837 # Zoom on inter-termini seq
|
|
|
838 if isinstance(P_left, np.integer) and isinstance(P_right, np.integer) and not Mu_like:
|
|
|
839 Zoom_left = min(P_left-1000, P_right-1000)
|
|
|
840 Zoom_right = max(P_left+1000, P_right+1000)
|
|
|
841 imgdata = io.BytesIO()
|
|
|
842 if P_orient == "Reverse":
|
|
|
843 zoom_pos_left = P_right-max(0,Zoom_left)
|
|
|
844 zoom_pos_right = P_left-max(0,Zoom_left)
|
|
|
845 else:
|
|
|
846 zoom_pos_left = P_left-max(0,Zoom_left)
|
|
|
847 zoom_pos_right = P_right-max(0,Zoom_left)
|
|
|
848
|
|
|
849 figZ_whole = GraphWholeCov(added_whole_coverage[max(0,Zoom_left):min(gen_len,Zoom_right)], phagename + "-zoom", draw, P_left, P_right, zoom_pos_left, zoom_pos_right, 1, "Zoom Termini")
|
|
|
850 figZ_whole.savefig(imgdata, format='png')
|
|
|
851 imgdata.seek(0)
|
|
|
852 IMG = ImageReader(imgdata)
|
|
|
853 IMAGE = Image(IMG.fileName, width=275, height=340, kind='proportional')
|
|
|
854 IMAGE.hAlign = 'CENTER'
|
|
|
855
|
|
|
856 data = [[IMAGE, IMAGE_2]]
|
|
|
857 t=Table(data, 1*[4*inch]+1*[3*inch], 1*[2*inch], hAlign='CENTER', style=[('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),10), ('ALIGN',(0,0),(-1,-1),'LEFT'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
858 report.append(t)
|
|
|
859 report.append(Spacer(1, 5))
|
|
|
860
|
|
|
861 elif isinstance(P_left, np.integer) and P_orient == "Forward":
|
|
|
862 imgdata = io.BytesIO()
|
|
|
863
|
|
|
864 if Mu_like:
|
|
|
865 figZL_whole = GraphWholeCov(phage_hybrid_coverage[0][max(0,P_left-1000):min(gen_len,P_left+1000)], phagename + "-zoom-left", draw, P_left, "", P_left-max(0,P_left-1000), 0, 1, "Zoom Termini")
|
|
|
866 else:
|
|
|
867 figZL_whole = GraphWholeCov(added_whole_coverage[max(0,P_left-1000):min(gen_len,P_left+1000)], phagename + "-zoom-left", draw, P_left, P_right, P_left-max(0,P_left-1000), 0, 1, "Zoom Termini")
|
|
|
868 figZL_whole.savefig(imgdata, format='png')
|
|
|
869 imgdata.seek(0)
|
|
|
870 IMG = ImageReader(imgdata)
|
|
|
871 IMAGE = Image(IMG.fileName, width=275, height=340, kind='proportional')
|
|
|
872 IMAGE.hAlign = 'CENTER'
|
|
|
873
|
|
|
874 data = [[IMAGE, IMAGE_2]]
|
|
|
875 t=Table(data, 1*[5*inch]+1*[3*inch], 1*[2*inch], hAlign='CENTER', style=[('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),10), ('ALIGN',(0,0),(-1,-1),'LEFT'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
876 report.append(t)
|
|
|
877
|
|
|
878 elif isinstance(P_right, np.integer) and P_orient == "Reverse":
|
|
|
879 imgdata = io.BytesIO()
|
|
|
880
|
|
|
881 if Mu_like:
|
|
|
882 figZR_whole = GraphWholeCov(phage_hybrid_coverage[1][max(0,P_right-1000):min(gen_len,P_right+1000)], phagename + "-zoom-right", draw, "", P_right, 0, P_right-max(0,P_right-1000), 1, "Zoom Termini")
|
|
|
883 else:
|
|
|
884 figZR_whole = GraphWholeCov(added_whole_coverage[max(0,P_right-1000):min(gen_len,P_right+1000)], phagename + "-zoom-right", draw, P_left, P_right, 0, P_right-max(0,P_right-1000), 1, "Zoom Termini")
|
|
|
885 figZR_whole.savefig(imgdata, format='png')
|
|
|
886 imgdata.seek(0)
|
|
|
887 IMG = ImageReader(imgdata)
|
|
|
888 IMAGE = Image(IMG.fileName, width=275, height=340, kind='proportional')
|
|
|
889 IMAGE.hAlign = 'CENTER'
|
|
|
890
|
|
|
891 data = [[IMAGE, IMAGE_2]]
|
|
|
892 t=Table(data, 1*[5*inch]+1*[3*inch], 1*[2*inch], hAlign='CENTER', style=[('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),10), ('ALIGN',(0,0),(-1,-1),'LEFT'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
893 report.append(t)
|
|
|
894 report.append(Spacer(1, 5))
|
|
|
895 else:
|
|
|
896 data = [[IMAGE_2]]
|
|
|
897 t=Table(data, 1*[1.5*inch], 1*[2*inch], hAlign='CENTER', style=[('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),10), ('ALIGN',(0,0),(-1,-1),'LEFT'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
898 report.append(t)
|
|
|
899
|
|
|
900 # Warning coverage message
|
|
|
901 if ave_whole_cov < 50 and test_run == 0:
|
|
|
902 ptextawc = "- - - - - - - - - WARNING: Coverage (" + str(int(ave_whole_cov)) + ") is under the limit of the software, Please consider results carrefuly. - - - - - - - - -"
|
|
|
903 data = [[ptextawc]]
|
|
|
904 t=Table(data, 1*[5*inch], hAlign='LEFT', style=[('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('TEXTCOLOR',(0,0),(-1,-1),'RED'), ('FONTSIZE',(0,0),(-1,-1),10), ('ALIGN',(0,0),(-1,-1),'LEFT'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
905 report.append(t)
|
|
|
906
|
|
|
907 ## Statistics
|
|
|
908 ptext = '<u><font size=14>PhageTerm Method</font></u>'
|
|
|
909 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
910 report.append(Spacer(1, 10))
|
|
|
911
|
|
|
912 if Redundant:
|
|
|
913 Ends = "Redundant"
|
|
|
914 else:
|
|
|
915 Ends = "Non Red."
|
|
|
916
|
|
|
917 data = [["Ends", "Left (red)", "Right (green)", "Permuted", "Orientation", "Class", "Type"], [Ends, P_left, P_right, Permuted, P_orient, P_class, P_type]]
|
|
|
918 t=Table(data, 7*[1.10*inch], 2*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-2),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
919 report.append(t)
|
|
|
920 report.append(Spacer(1, 5))
|
|
|
921
|
|
|
922 # Seq cohesive or Direct terminal repeats
|
|
|
923 if P_seqcoh != "":
|
|
|
924 if len(P_seqcoh) < 20:
|
|
|
925 ptext = '<i><font size=12>*Sequence cohesive: ' + P_seqcoh + '</font></i>'
|
|
|
926 else:
|
|
|
927 ptext = '<i><font size=12>*Direct Terminal Repeats: ' + str(len(P_seqcoh)) + ' bp</font></i>'
|
|
|
928 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
929
|
|
|
930 # Multiple / Multiple (Nextera)
|
|
|
931 if P_left == "Multiple" and P_right == "Multiple":
|
|
|
932 ptext = '<i><font size=12>*This results could be due to a non-random fragmented sequence (e.g. Nextera)</font></i>'
|
|
|
933 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
934
|
|
|
935
|
|
|
936 # Concatermer
|
|
|
937 elif P_class[:7] == "Headful" and paired != "":
|
|
|
938 ptext = '<i><font size=12>*concatemer estimation: ' + str(P_concat) + '</font></i>'
|
|
|
939 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
940
|
|
|
941 # Mu hybrid
|
|
|
942 elif Mu_like:
|
|
|
943 if P_orient == "Forward":
|
|
|
944 Mu_termini = P_left
|
|
|
945 else:
|
|
|
946 Mu_termini = P_right
|
|
|
947 ptext = '<i><font size=12>*Mu estimated termini position with hybrid fragments: ' + str(Mu_termini) + '</font></i>'
|
|
|
948 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
949
|
|
|
950 report.append(Spacer(1, 10))
|
|
|
951
|
|
|
952 # Results
|
|
|
953 imgdata = io.BytesIO()
|
|
|
954 figP_norm = GraphCov(termini_coverage_norm_close, picMaxPlus_norm_close[:1], picMaxMinus_norm_close[:1], phagename + "-norm", 1, draw)
|
|
|
955 figP_norm.savefig(imgdata, format='png')
|
|
|
956 imgdata.seek(0)
|
|
|
957 IMG = ImageReader(imgdata)
|
|
|
958 IMAGE = Image(IMG.fileName, width=240, height=340, kind='proportional')
|
|
|
959 IMAGE.hAlign = 'CENTER'
|
|
|
960
|
|
|
961 data = [["Strand", "Location", "T", "pvalue", "T (Start. Pos. Cov. / Whole Cov.)"], ["+",phage_plus_norm["Position"][0],format(phage_plus_norm["SPC"][0]/100.0, '0.2f'),format(phage_plus_norm["pval_gamma_adj"][0], '0.2e'),IMAGE], ["",phage_plus_norm["Position"][1],format(phage_plus_norm["SPC"][1]/100.0, '0.2f'),format(phage_plus_norm["pval_gamma_adj"][1], '0.2e'),""], ["",phage_plus_norm["Position"][2],format(phage_plus_norm["SPC"][2]/100.0, '0.2f'),format(phage_plus_norm["pval_gamma_adj"][2], '0.2e'),""], ["",phage_plus_norm["Position"][3],format(phage_plus_norm["SPC"][3]/100.0, '0.2f'),format(phage_plus_norm["pval_gamma_adj"][3], '0.2e'),""], ["",phage_plus_norm["Position"][4],format(phage_plus_norm["SPC"][4]/100.0, '0.2f'),format(phage_plus_norm["pval_gamma_adj"][4], '0.2e'),""], ["-",phage_minus_norm["Position"][0],format(phage_minus_norm["SPC"][0]/100.0, '0.2f'),format(phage_minus_norm["pval_gamma_adj"][0], '0.2e'),""], ["",phage_minus_norm["Position"][1],format(phage_minus_norm["SPC"][1]/100.0, '0.2f'),format(phage_minus_norm["pval_gamma_adj"][1], '0.2e'),""], ["",phage_minus_norm["Position"][2],format(phage_minus_norm["SPC"][2]/100.0, '0.2f'),format(phage_minus_norm["pval_gamma_adj"][2], '0.2e'),""], ["",phage_minus_norm["Position"][3],format(phage_minus_norm["SPC"][3]/100.0, '0.2f'),format(phage_minus_norm["pval_gamma_adj"][3], '0.2e'),""], ["",phage_minus_norm["Position"][4],format(phage_minus_norm["SPC"][4]/100.0, '0.2f'),format(phage_minus_norm["pval_gamma_adj"][4], '0.2e'),""]]
|
|
|
962 t=Table(data, 4*[1*inch]+1*[4*inch], 11*[0.25*inch], hAlign='CENTER', style=[('SPAN',(0,1),(0,5)), ('SPAN',(0,6),(0,10)), ('SPAN',(4,1),(4,10)), ('LINEABOVE',(0,1),(4,1),1.5,colors.black), ('LINEABOVE',(0,6),(4,6),1.5,colors.grey), ('FONT',(0,0),(-1,0),'Helvetica-Bold'), ('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),12), ('FONTSIZE',(0,1),(0,-1),16), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
963
|
|
|
964 report.append(t)
|
|
|
965 report.append(Spacer(1, 5))
|
|
|
966
|
|
|
967 ## Li's Analysis
|
|
|
968 ptext = '<u><font size=14>Li\'s Method</font></u>'
|
|
|
969 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
970 report.append(Spacer(1, 10))
|
|
|
971
|
|
|
972 data = [["Packaging", "Termini", "Forward", "Reverse", "Orientation"], [ArtPackmode, termini, forward, reverse, ArtOrient]]
|
|
|
973 t=Table(data, 2*[1*inch] + 2*[2*inch] + 1*[1*inch], 2*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-2),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
974
|
|
|
975 report.append(t)
|
|
|
976 report.append(Spacer(1, 5))
|
|
|
977
|
|
|
978 # Seq cohesive or Direct terminal repeats
|
|
|
979 if len(ArtcohesiveSeq) > 2:
|
|
|
980 if len(ArtcohesiveSeq) < 20:
|
|
|
981 ptext = '<i><font size=12>*Sequence cohesive: ' + ArtcohesiveSeq + '</font></i>'
|
|
|
982 else:
|
|
|
983 ptext = '<i><font size=12>*Direct Terminal Repeats: ' + str(len(ArtcohesiveSeq)) + ' bp</font></i>'
|
|
|
984 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
985 report.append(Spacer(1, 10))
|
|
|
986
|
|
|
987 # Results
|
|
|
988 imgdata = io.BytesIO()
|
|
|
989 figP = GraphCov(termini_coverage_close, picMaxPlus_close[:1], picMaxMinus_close[:1], phagename, 0, draw)
|
|
|
990 figP.savefig(imgdata, format='png')
|
|
|
991 imgdata.seek(0)
|
|
|
992 IMG = ImageReader(imgdata)
|
|
|
993 IMAGE = Image(IMG.fileName, width=240, height=340, kind='proportional')
|
|
|
994 IMAGE.hAlign = 'CENTER'
|
|
|
995
|
|
|
996 data = [["Strand", "Location", "SPC", "R", "SPC"],["+",picMaxPlus_close[0][1]+1,picMaxPlus_close[0][0],R2,IMAGE],["",picMaxPlus_close[1][1]+1,picMaxPlus_close[1][0],"-",""],["",picMaxPlus_close[2][1]+1,picMaxPlus_close[2][0],"-",""],["",picMaxPlus_close[3][1]+1,picMaxPlus_close[3][0],"-",""],["",picMaxPlus_close[4][1]+1,picMaxPlus_close[4][0],"-",""],["-",picMaxMinus_close[0][1]+1,picMaxMinus_close[0][0],R3,""], ["",picMaxMinus_close[1][1]+1,picMaxMinus_close[1][0],"-",""], ["",picMaxMinus_close[2][1]+1,picMaxMinus_close[2][0],"-",""], ["",picMaxMinus_close[3][1]+1,picMaxMinus_close[3][0],"-",""], ["",picMaxMinus_close[4][1]+1,picMaxMinus_close[4][0],"-",""]]
|
|
|
997 t=Table(data, 4*[1*inch]+1*[4*inch], 11*[0.25*inch], hAlign='CENTER', style=[('SPAN',(0,1),(0,5)), ('SPAN',(0,6),(0,10)), ('SPAN',(4,1),(4,10)), ('LINEABOVE',(0,1),(4,1),1.5,colors.black), ('LINEABOVE',(0,6),(4,6),1.5,colors.grey), ('FONT',(0,0),(-1,0),'Helvetica-Bold'), ('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),12), ('FONTSIZE',(0,1),(0,-1),16), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
998
|
|
|
999 report.append(t)
|
|
|
1000
|
|
|
1001
|
|
|
1002 # NEW PAGE
|
|
|
1003 report.append(PageBreak())
|
|
|
1004
|
|
|
1005 # HOST RESULTS
|
|
|
1006 if host != "":
|
|
|
1007 # Host coverage
|
|
|
1008 ptext = '<u><font size=14>Host Analysis</font></u>'
|
|
|
1009 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
1010 report.append(Spacer(1, 10))
|
|
|
1011
|
|
|
1012 ptext = '<i><font size=10></font>Reads that does not match on the phage genome are tested on the host genome. These reads could come from Phage transduction but also Host DNA contamination.</i>'
|
|
|
1013 report.append(Paragraph(ptext, styles["Justify"]))
|
|
|
1014 report.append(Spacer(1, 5))
|
|
|
1015
|
|
|
1016 data = [["Host Genome", str(host_len) + " bp"]]
|
|
|
1017 t=Table(data, 2*[2.25*inch], hAlign='LEFT', style=[('FONT',(0,0),(0,0),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'LEFT') ,('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1018
|
|
|
1019 report.append(t)
|
|
|
1020 report.append(Spacer(1, 5))
|
|
|
1021
|
|
|
1022 imgdata = io.BytesIO()
|
|
|
1023
|
|
|
1024 figH = GraphCov(host_whole_coverage, picMaxPlus_host[:1], picMaxMinus_host[:1], "", 0, draw)
|
|
|
1025 figH.savefig(imgdata, format='png')
|
|
|
1026 imgdata.seek(0)
|
|
|
1027 IMG = ImageReader(imgdata)
|
|
|
1028 IMAGE = Image(IMG.fileName, width=240, height=340, kind='proportional')
|
|
|
1029 IMAGE.hAlign = 'CENTER'
|
|
|
1030
|
|
|
1031 data = [["Strand", "Location", "Coverage", "-", "Whole Coverage"],["+",picMaxPlus_host[0][1]+1,picMaxPlus_host[0][0],"-",IMAGE],["",picMaxPlus_host[1][1]+1,picMaxPlus_host[1][0],"-",""],["",picMaxPlus_host[2][1]+1,picMaxPlus_host[2][0],"-",""],["",picMaxPlus_host[3][1]+1,picMaxPlus_host[3][0],"-",""],["",picMaxPlus_host[4][1]+1,picMaxPlus_host[4][0],"-",""],["-",picMaxMinus_host[0][1]+1,picMaxMinus_host[0][0],"-",""], ["",picMaxMinus_host[1][1]+1,picMaxMinus_host[1][0],"-",""], ["",picMaxMinus_host[2][1]+1,picMaxMinus_host[2][0],"-",""], ["",picMaxMinus_host[3][1]+1,picMaxMinus_host[3][0],"-",""], ["",picMaxMinus_host[4][1]+1,picMaxMinus_host[4][0],"-",""]]
|
|
|
1032 t=Table(data, 4*[1*inch]+1*[4*inch], 11*[0.25*inch], hAlign='CENTER', style=[('SPAN',(0,1),(0,5)), ('SPAN',(0,6),(0,10)), ('SPAN',(4,1),(4,10)), ('LINEABOVE',(0,1),(4,1),1.5,colors.black), ('LINEABOVE',(0,6),(4,6),1.5,colors.grey), ('FONT',(0,0),(-1,0),'Helvetica-Bold'), ('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),12), ('FONTSIZE',(0,1),(0,-1),16), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1033
|
|
|
1034 report.append(t)
|
|
|
1035 report.append(Spacer(1, 10))
|
|
|
1036
|
|
|
1037 # Hybrid coverage
|
|
|
1038 ptext = '<u><font size=14>Hybrid Analysis</font></u>'
|
|
|
1039 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
1040 report.append(Spacer(1, 10))
|
|
|
1041
|
|
|
1042 ptext = '<i><font size=10></font>Hybrid reads with one edge on the phage genome and the other edge on the host genome are detected. Phage Hybrid Coverages are used to detect Mu-like packaging mode. Host Hybrid Coverages could be used to detect Phage Transduction but also genome integration location of prophages.</i>'
|
|
|
1043 report.append(Paragraph(ptext, styles["Justify"]))
|
|
|
1044 report.append(Spacer(1, 5))
|
|
|
1045
|
|
|
1046 picMaxPlus_phage_hybrid, picMaxMinus_phage_hybrid, TopFreqH_phage_hybrid = picMax(phage_hybrid_coverage, 5)
|
|
|
1047 picMaxPlus_host_hybrid, picMaxMinus_host_hybrid, TopFreqH_host_hybrid = picMax(host_hybrid_coverage, 5)
|
|
|
1048
|
|
|
1049 imgdataPH = io.BytesIO()
|
|
|
1050 figPH = GraphCov(phage_hybrid_coverage, picMaxPlus_phage_hybrid[:1], picMaxMinus_phage_hybrid[:1], "", 0, draw, 1)
|
|
|
1051 figPH.savefig(imgdataPH, format='png')
|
|
|
1052 imgdataPH.seek(0)
|
|
|
1053 IMGPH = ImageReader(imgdataPH)
|
|
|
1054 IMAGEPH = Image(IMGPH.fileName, width=240, height=340, kind='proportional')
|
|
|
1055 IMAGEPH.hAlign = 'CENTER'
|
|
|
1056
|
|
|
1057
|
|
|
1058 imgdataHH = io.BytesIO()
|
|
|
1059 figHH = GraphCov(host_hybrid_coverage, picMaxPlus_host_hybrid[:1], picMaxMinus_host_hybrid[:1], "", 0, draw, 1)
|
|
|
1060 figHH.savefig(imgdataHH, format='png')
|
|
|
1061 imgdataHH.seek(0)
|
|
|
1062 IMGHH = ImageReader(imgdataHH)
|
|
|
1063 IMAGEHH = Image(IMGHH.fileName, width=240, height=340, kind='proportional')
|
|
|
1064 IMAGEHH.hAlign = 'CENTER'
|
|
|
1065
|
|
|
1066 data = [["Phage Hybrid Coverage", "Host Hybrid Coverage"],[IMAGEPH,IMAGEHH]]
|
|
|
1067 t=Table(data, 2*[4*inch], 1*[0.25*inch]+1*[2.5*inch], hAlign='CENTER', style=[('LINEABOVE',(0,1),(1,1),1.5,colors.black),('FONT',(0,0),(-1,-1),'Helvetica-Bold'),('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1068
|
|
|
1069 report.append(t)
|
|
|
1070 report.append(Spacer(1, 10))
|
|
|
1071
|
|
|
1072 # NEW PAGE
|
|
|
1073 report.append(PageBreak())
|
|
|
1074
|
|
|
1075
|
|
|
1076 # DETAILED RESULTS
|
|
|
1077 ptext = '<u><font size=14>Analysis Methodology</font></u>'
|
|
|
1078 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
1079 report.append(Spacer(1, 10))
|
|
|
1080
|
|
|
1081 ptext = '<i><font size=10>PhageTerm software uses raw reads of a phage sequenced with a sequencing technology using random fragmentation and its genomic reference sequence to determine the termini position. The process starts with the alignment of NGS reads to the phage genome in order to calculate the starting position coverage (SPC), where a hit is given only to the position of the first base in a successfully aligned read (the alignment algorithm uses the lenght of the seed (default: 20) for mapping and does not accept gap or missmatch to speed up the process). Then the program apply 2 distinct scoring methods: i) a statistical approach based on the Gamma law; and ii) a method derived from LI and al. 2014 paper.</font></i>'
|
|
|
1082 report.append(Paragraph(ptext, styles["Justify"]))
|
|
|
1083 report.append(Spacer(1, 5))
|
|
|
1084
|
|
|
1085
|
|
|
1086 # INFORMATION
|
|
|
1087 ptext = '<u><font size=12>General set-up and mapping informations</font></u>'
|
|
|
1088 report.append(Paragraph(ptext, styles["Justify"]))
|
|
|
1089 report.append(Spacer(1, 5))
|
|
|
1090
|
|
|
1091
|
|
|
1092 imgdata = io.BytesIO()
|
|
|
1093
|
|
|
1094 if paired != "":
|
|
|
1095 figP_whole = GraphWholeCov(added_paired_whole_coverage, phagename, draw)
|
|
|
1096 else:
|
|
|
1097 figP_whole = GraphWholeCov(added_whole_coverage, phagename, draw)
|
|
|
1098 figP_whole.savefig(imgdata, format='png')
|
|
|
1099 imgdata.seek(0)
|
|
|
1100 IMG = ImageReader(imgdata)
|
|
|
1101 IMAGE = Image(IMG.fileName, width=275, height=340, kind='proportional')
|
|
|
1102 IMAGE.hAlign = 'CENTER'
|
|
|
1103
|
|
|
1104 if host == "":
|
|
|
1105 host_analysis = "No"
|
|
|
1106 else:
|
|
|
1107 host_analysis = "Yes"
|
|
|
1108
|
|
|
1109 if paired == "":
|
|
|
1110 sequencing_reads = "Single-ends Reads"
|
|
|
1111 else:
|
|
|
1112 sequencing_reads = "Paired-ends Reads"
|
|
|
1113
|
|
|
1114 data = [["Phage Genome ", str(gen_len) + " bp",IMAGE], ["Sequencing Reads", int(tot_reads),""], ["Mapping Reads", str(int(100 - lost_perc)) + " %",""], ["OPTIONS","",""], ["Mapping Seed",seed,""], ["Surrounding",surrounding,""], ["Host Analysis ", host_analysis,""], ["","",""]]
|
|
|
1115 t=Table(data, 1*[2.25*inch]+1*[1.80*inch]+1*[4*inch], 8*[0.25*inch], hAlign='LEFT', style=[('SPAN',(2,0),(2,-1)), ('FONT',(0,0),(0,2),'Helvetica-Bold'), ('FONT',(0,3),(0,3),'Helvetica-Oblique'), ('FONT',(0,4),(1,-1),'Helvetica-Oblique'), ('FONT',(2,0),(2,0),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-2),'LEFT'), ('ALIGN',(2,0),(2,-1),'CENTER') ,('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1116
|
|
|
1117 report.append(t)
|
|
|
1118 report.append(Spacer(1, 5))
|
|
|
1119
|
|
|
1120
|
|
|
1121 # Img highest peaks of each side even if no significative
|
|
|
1122 ptext = '<u><font size=12>Highest peak of each side coverage graphics</font></u>'
|
|
|
1123 report.append(Paragraph(ptext, styles["Justify"]))
|
|
|
1124 report.append(Spacer(1, 5))
|
|
|
1125
|
|
|
1126
|
|
|
1127 imgdata = io.BytesIO()
|
|
|
1128
|
|
|
1129 if Mu_like and isinstance(P_left, np.integer):
|
|
|
1130 figHL_whole = GraphWholeCov(phage_hybrid_coverage[0][max(0,P_left-1000):min(gen_len,P_left+1000)], phagename + "-zoom-left", draw, P_left, "", P_left-max(0,P_left-1000), 0, 1, "Zoom Termini")
|
|
|
1131 else:
|
|
|
1132 P_left = phage_plus_norm["Position"][0]
|
|
|
1133 figHL_whole = GraphWholeCov(added_whole_coverage[max(0,P_left-1000):min(gen_len,P_left+1000)], phagename + "-zoom-left", draw, P_left, "", P_left-max(0,P_left-1000), 0, 1, "Zoom Termini")
|
|
|
1134 figHL_whole.savefig(imgdata, format='png')
|
|
|
1135 imgdata.seek(0)
|
|
|
1136 IMG = ImageReader(imgdata)
|
|
|
1137 IMAGE = Image(IMG.fileName, width=275, height=340, kind='proportional')
|
|
|
1138 IMAGE.hAlign = 'CENTER'
|
|
|
1139
|
|
|
1140 imgdata2 = io.BytesIO()
|
|
|
1141
|
|
|
1142 if Mu_like and isinstance(P_right, np.integer):
|
|
|
1143 figHR_whole = GraphWholeCov(phage_hybrid_coverage[1][max(0,P_right-1000):min(gen_len,P_right+1000)], phagename + "-zoom-right", draw, "", P_right, 0, P_right-max(0,P_right-1000), 1, "Zoom Termini")
|
|
|
1144 else:
|
|
|
1145 P_right = phage_minus_norm["Position"][0]
|
|
|
1146 figHR_whole = GraphWholeCov(added_whole_coverage[max(0,P_right-1000):min(gen_len,P_right+1000)], phagename + "-zoom-right", draw, "", P_right, 0, P_right-max(0,P_right-1000), 1, "Zoom Termini")
|
|
|
1147 figHR_whole.savefig(imgdata2, format='png')
|
|
|
1148 imgdata2.seek(0)
|
|
|
1149 IMG2 = ImageReader(imgdata2)
|
|
|
1150 IMAGE2 = Image(IMG2.fileName, width=275, height=340, kind='proportional')
|
|
|
1151 IMAGE2.hAlign = 'CENTER'
|
|
|
1152
|
|
|
1153 if Mu_like:
|
|
|
1154 data = [["Hybrid Coverage Zoom (Left)", "Hybrid Coverage Zoom (Right)"],[IMAGE,IMAGE2]]
|
|
|
1155 else:
|
|
|
1156 data = [["Whole Coverage Zoom (Left)", "Whole Coverage Zoom (Right)"],[IMAGE,IMAGE2]]
|
|
|
1157 t=Table(data, 2*[4*inch], 1*[0.25*inch]+1*[2*inch], hAlign='CENTER', style=[('LINEABOVE',(0,1),(1,1),1.5,colors.black),('FONT',(0,0),(-1,-1),'Helvetica-Bold'),('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1158 report.append(t)
|
|
|
1159
|
|
|
1160 # Controls
|
|
|
1161 ptext = '<u><font size=12>General controls information</font></u>'
|
|
|
1162 report.append(Paragraph(ptext, styles["Justify"]))
|
|
|
1163 report.append(Spacer(1, 5))
|
|
|
1164
|
|
|
1165 if ave_whole_cov < 50:
|
|
|
1166 ptextawc = "WARNING: Under the limit of the software (50)"
|
|
|
1167 elif ave_whole_cov < 200:
|
|
|
1168 ptextawc = "WARNING: Low (<200), Li's method could not be reliable"
|
|
|
1169 else:
|
|
|
1170 ptextawc = "OK"
|
|
|
1171
|
|
|
1172 data = [["Whole genome coverage", int(ave_whole_cov), ptextawc]]
|
|
|
1173 t=Table(data, 1*[3.5*inch]+1*[1*inch]+1*[3.5*inch], hAlign='LEFT', style=[('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),10), ('ALIGN',(0,0),(-1,-1),'LEFT'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1174 report.append(t)
|
|
|
1175
|
|
|
1176 drop_perc = len([i for i in added_whole_coverage if i < (ave_whole_cov/2)]) / float(len(added_whole_coverage))
|
|
|
1177 if drop_perc < 1:
|
|
|
1178 ptextdp = "OK"
|
|
|
1179 else:
|
|
|
1180 ptextdp = "Check your genome reference"
|
|
|
1181
|
|
|
1182 data = [["Weak genome coverage", "%.1f %%" %drop_perc, ptextdp]]
|
|
|
1183 t=Table(data, 1*[3.5*inch]+1*[1*inch]+1*[4*inch], hAlign='LEFT', style=[('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),10), ('ALIGN',(0,0),(-1,-1),'LEFT'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1184 report.append(t)
|
|
|
1185
|
|
|
1186 if paired != "":
|
|
|
1187 if len(insert) != 0:
|
|
|
1188 insert_mean = sum(insert)/len(insert)
|
|
|
1189 else:
|
|
|
1190 insert_mean = "-"
|
|
|
1191 data = [["Insert mean size", insert_mean, "Mean insert estimated from paired-end reads"]]
|
|
|
1192 t=Table(data, 1*[3.5*inch]+1*[1*inch]+1*[4*inch], hAlign='LEFT', style=[('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),10), ('ALIGN',(0,0),(-1,-1),'LEFT'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1193 report.append(t)
|
|
|
1194
|
|
|
1195 if lost_perc > 25:
|
|
|
1196 ptextlp = "Warning: high percentage of reads lost"
|
|
|
1197 else:
|
|
|
1198 ptextlp = "OK"
|
|
|
1199
|
|
|
1200 data = [["Reads lost during alignment", "%.1f %%" %lost_perc, ptextlp]]
|
|
|
1201 t=Table(data, 1*[3.5*inch]+1*[1*inch]+1*[4*inch], hAlign='LEFT', style=[('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),10), ('ALIGN',(0,0),(-1,-1),'LEFT'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1202 report.append(t)
|
|
|
1203 report.append(Spacer(1, 5))
|
|
|
1204
|
|
|
1205 # DETAILED SCORE
|
|
|
1206 ptext = '<b><font size=14>i) PhageTerm method</font></b>'
|
|
|
1207 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
1208 report.append(Spacer(1, 10))
|
|
|
1209
|
|
|
1210 ptext = '<i><font size=10>Reads are mapped on the reference to determine the starting position coverage (SPC) as well as the coverage (COV) in each orientation. These values are then used to compute the variable T = SPC/COV. The average value of T at positions along the genome that are not termini is expected to be 1/F, where F is the average fragment size. For the termini that depends of the packaging mode. Cos Phages: no reads should start before the terminus and therefore T=1. DTR phages: for N phages present in the sample, there should be N fragments that start at the terminus and N fragments that cover the edge of the repeat on the other side of the genome as a results T is expected to be 0.5. Pac phages: for N phages in the sample, there should be N/C fragments starting at the pac site, where C is the number of phage genome copies per concatemer. In the same sample N fragments should cover the pac site position, T is expected to be (N/C)/(N+N/C) = 1/(1+C). To assess whether the number of reads starting at a given position along the genome can be considered a significant outlier, PhageTerm first segments the genome according to coverage using a regression tree. A gamma distribution is fitted to SPC for each segment and an adjusted p-value is computed for each position. If several significant peaks are detected within a small sequence window (default: 20bp), their X values are merged.</font></i>'
|
|
|
1211 report.append(Paragraph(ptext, styles["Justify"]))
|
|
|
1212 report.append(Spacer(1, 5))
|
|
|
1213
|
|
|
1214 # surrounding
|
|
|
1215 if surrounding > 0:
|
|
|
1216 data = [["Nearby Termini (Forward / Reverse)", str(len(picOUT_norm_forw)-1) + " / " + str(len(picOUT_norm_rev)-1), "Peaks localized %s bases around the maximum" %surrounding]]
|
|
|
1217 t=Table(data, 1*[3.5*inch]+1*[1*inch]+1*[4*inch], 1*[0.25*inch], hAlign='LEFT', style=[('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),10), ('ALIGN',(0,0),(-1,-1),'LEFT'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1218 report.append(t)
|
|
|
1219
|
|
|
1220 report.append(Spacer(1, 10))
|
|
|
1221
|
|
|
1222 # Li's Method
|
|
|
1223 if not multi:
|
|
|
1224 ptext = '<b><font size=14>ii) Li\'s method</font></b>'
|
|
|
1225 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
1226 report.append(Spacer(1, 10))
|
|
|
1227
|
|
|
1228 ptext = '<i><font size=10>The second approach is based on the calculation and interpretation of three specific ratios R1, R2 and R3 as suggested in a previous publication from Li et al. 2014. The first ratio, is calculated as follow: the highest starting frequency found on either the forward or reverse strands is divided by the average starting frequency, R1 = (highest frequency/average frequency). Li’s et al. have proposed three possible interpretation of the R1 ratio. First, if R1 < 30, the phage genome does not have any termini, and is either circular or completely permuted and terminally redundant. The second interpretation for R1 is when 30 ≤ R1 ≥ 100, suggesting the presence of preferred termini with terminal redundancy and apparition of partially circular permutations. At last if R1 > 100 that is an indication that at least one fixed termini is present with terminase recognizing a specific site. The two other ratios are R2 and R3 and the calculation is done in a similar manner. R2 is calculated using the highest two frequencies (T1-F and T2-F) found on the forward strand and R3 is calculated using the highest two frequencies (T1-R and T2-R) found on the reverse strand. To calculate these two ratios, we divide the highest frequency by the second highest frequency T2. So R2 = (T1-F / T2-F) and R3 = (T1-R / T2-R). These two ratios are used to analyze termini characteristics on each strand taken individually. Li et al. suggested two possible interpretations for R2 and R3 ratios combine to R1. When R1 < 30 and R2 < 3, we either have no obvious termini on the forward strand, or we have multiple preferred termini on the forward strand, if 30 ≤ R1 ≤ 100. If R2 > 3, it is suggested that there is an obvious unique termini on the forward strand. The same reasoning is applicable for the result of R3. Combining the results for ratios found with this approach, it is possible to make the first prediction for the viral packaging mode of the analyzed phage. A unique obvious termini present at both ends (both R2 and R3 > 3) reveals the presence of a COS mode of packaging. The headful mode of packaging PAC is concluded when we have a single obvious termini only on one strand.</font></i>'
|
|
|
1229 report.append(Paragraph(ptext, styles["Justify"]))
|
|
|
1230 report.append(Spacer(1, 5))
|
|
|
1231
|
|
|
1232 if surrounding > 0:
|
|
|
1233 data = [["Nearby Termini (Forward / Reverse)", str(len(picOUT_forw)-1) + " / " + str(len(picOUT_rev)-1), "Peaks localized %s bases around the maximum" %surrounding]]
|
|
|
1234 t=Table(data, 1*[3.5*inch]+1*[1*inch]+1*[3.5*inch], 1*[0.25*inch], hAlign='LEFT', style=[('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),10), ('ALIGN',(0,0),(-1,-1),'LEFT'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1235 report.append(t)
|
|
|
1236 report.append(Spacer(1, 5))
|
|
|
1237
|
|
|
1238 if R1 > 100:
|
|
|
1239 ptextR1 = "At least one fixed termini is present with terminase recognizing a specific site."
|
|
|
1240 elif R1 > 30:
|
|
|
1241 ptextR1 = "Presence of preferred termini with terminal redundancy and apparition of partially circular permutations."
|
|
|
1242 else:
|
|
|
1243 ptextR1 = "Phage genome does not have any termini, and is either circular or completely permuted and terminally redundant."
|
|
|
1244
|
|
|
1245 data = [["R1 - highest freq./average freq.", int(R1), Paragraph(ptextR1, styles["Justify"])]]
|
|
|
1246 t=Table(data, 1*[3.5*inch]+1*[1*inch]+1*[3.5*inch], 1*[0.25*inch], hAlign='LEFT', style=[('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),10), ('ALIGN',(0,0),(-1,-1),'LEFT'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1247 report.append(t)
|
|
|
1248 report.append(Spacer(1, 5))
|
|
|
1249
|
|
|
1250 if R2 < 3 and R1 < 30:
|
|
|
1251 ptextR2 = "No obvious termini on the forward strand."
|
|
|
1252 elif R2 < 3 :
|
|
|
1253 ptextR2 = "Multiple preferred termini on the forward strand."
|
|
|
1254 elif R2 >= 3:
|
|
|
1255 ptextR2 = "Unique termini on the forward strand."
|
|
|
1256
|
|
|
1257 data = [["R2 Forw - highest freq./second freq.", int(R2), Paragraph(ptextR2, styles["Justify"])]]
|
|
|
1258 t=Table(data, 1*[3.5*inch]+1*[1*inch]+1*[3.5*inch], 1*[0.25*inch], hAlign='LEFT', style=[('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),10), ('ALIGN',(0,0),(-1,-1),'LEFT'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1259 report.append(t)
|
|
|
1260 report.append(Spacer(1, 5))
|
|
|
1261
|
|
|
1262 if R3 < 3 and R1 < 30:
|
|
|
1263 ptextR3 = "No obvious termini on the reverse strand."
|
|
|
1264 elif R3 < 3 :
|
|
|
1265 ptextR3 = "Multiple preferred termini on the reverse strand."
|
|
|
1266 elif R3 >= 3:
|
|
|
1267 ptextR3 = "Unique termini on the reverse strand."
|
|
|
1268
|
|
|
1269 data = [["R3 Rev - highest freq./second freq.", int(R3), Paragraph(ptextR3, styles["Justify"])]]
|
|
|
1270 t=Table(data, 1*[3.5*inch]+1*[1*inch]+1*[3.5*inch], 1*[0.25*inch], hAlign='LEFT', style=[('FONT',(0,0),(0,-1),'Helvetica-Bold'), ('FONTSIZE',(0,0),(-1,-1),10), ('ALIGN',(0,0),(-1,-1),'LEFT'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1271 report.append(t)
|
|
|
1272
|
|
|
1273 # CREDITS and TIME
|
|
|
1274 ptext = '<font size=8>%s</font>' % "Please cite: Sci. Rep. DOI 10.1038/s41598-017-07910-5"
|
|
|
1275 report.append(Paragraph(ptext, styles["Center"]))
|
|
|
1276 ptext = '<font size=8>%s</font>' % "Garneau, Depardieu, Fortier, Bikard and Monot. PhageTerm: Determining Bacteriophage Termini and Packaging using NGS data."
|
|
|
1277 report.append(Paragraph(ptext, styles["Center"]))
|
|
|
1278 ptext = '<font size=8>Report generated : %s</font>' % time.ctime()
|
|
|
1279 report.append(Paragraph(ptext, styles["Center"]))
|
|
|
1280
|
|
|
1281 # CREATE PDF
|
|
|
1282 if not multi:
|
|
|
1283 doc.build(report)
|
|
|
1284 else:
|
|
|
1285 report.append(PageBreak())
|
|
|
1286 return report
|
|
|
1287 return
|
|
|
1288
|
|
|
1289 def SummaryReport(phagename, DR, no_match):
|
|
|
1290 """ Create first page of multi reports."""
|
|
|
1291 report=[]
|
|
|
1292 styles=getSampleStyleSheet()
|
|
|
1293 styles.add(ParagraphStyle(name='Justify', alignment=TA_JUSTIFY))
|
|
|
1294 styles.add(ParagraphStyle(name='Center', alignment=TA_CENTER))
|
|
|
1295 styles.add(ParagraphStyle(name='Right', alignment=TA_RIGHT))
|
|
|
1296 styles.add(ParagraphStyle(name='Left', alignment=TA_LEFT))
|
|
|
1297
|
|
|
1298 ### GENERAL INFORMATION
|
|
|
1299
|
|
|
1300 # TITLE
|
|
|
1301 ptext = '<b><font size=16>' + phagename + ' PhageTerm Analysis</font></b>'
|
|
|
1302 report.append(Paragraph(ptext, styles["Center"]))
|
|
|
1303 report.append(Spacer(1, 15))
|
|
|
1304
|
|
|
1305 # No Match
|
|
|
1306 if len(no_match) > 0:
|
|
|
1307 ptext = '<u><font size=14>No Match ('+ str(len(no_match)) +')</font></u>'
|
|
|
1308 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
1309 report.append(Spacer(1, 10))
|
|
|
1310
|
|
|
1311 data = [["Name", "Class", "Left", "Right", "Type", "Orient", "Coverage", "Comments"]]
|
|
|
1312 t=Table(data, 2*[1.50*inch]+5*[0.80*inch]+1*[1.25*inch], 1*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-1),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1313 report.append(t)
|
|
|
1314
|
|
|
1315 for contig in no_match:
|
|
|
1316 P_comments = "No read match"
|
|
|
1317
|
|
|
1318 data = [[contig, "-", "-", "-", "-", "-", 0, P_comments]]
|
|
|
1319 t=Table(data, 2*[1.50*inch]+5*[0.80*inch]+1*[1.25*inch], 1*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-2),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1320 report.append(t)
|
|
|
1321
|
|
|
1322 # COS Phages
|
|
|
1323 count_COS = len(DR["COS (3')"]) + len(DR["COS (5')"]) + len(DR["COS"])
|
|
|
1324 ptext = '<u><font size=14>COS Phages ('+ str(count_COS) +')</font></u>'
|
|
|
1325 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
1326 report.append(Spacer(1, 10))
|
|
|
1327
|
|
|
1328 if count_COS != 0:
|
|
|
1329
|
|
|
1330 data = [["Name", "Class", "Left", "Right", "Type", "Orient", "Coverage", "Comments"]]
|
|
|
1331 t=Table(data, 2*[1.50*inch]+5*[0.80*inch]+1*[1.25*inch], 1*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-1),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1332 report.append(t)
|
|
|
1333
|
|
|
1334 for DC in DR["COS (3')"]:
|
|
|
1335 P_comments = ""
|
|
|
1336 if int(DR["COS (3')"][DC]["ave_whole_cov"]) < 50:
|
|
|
1337 P_comments = "Low coverage"
|
|
|
1338
|
|
|
1339 data = [[DC, DR["COS (3')"][DC]["P_class"], DR["COS (3')"][DC]["P_left"], DR["COS (3')"][DC]["P_right"], DR["COS (3')"][DC]["P_type"], DR["COS (3')"][DC]["P_orient"], int(DR["COS (3')"][DC]["ave_whole_cov"]), P_comments]]
|
|
|
1340 t=Table(data, 2*[1.50*inch]+5*[0.80*inch]+1*[1.25*inch], 1*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-2),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1341 report.append(t)
|
|
|
1342
|
|
|
1343 for DC in DR["COS (5')"]:
|
|
|
1344 P_comments = ""
|
|
|
1345 if int(DR["COS (5')"][DC]["ave_whole_cov"]) < 50:
|
|
|
1346 P_comments = "Low coverage"
|
|
|
1347
|
|
|
1348 data = [[DC, DR["COS (5')"][DC]["P_class"], DR["COS (5')"][DC]["P_left"], DR["COS (5')"][DC]["P_right"], DR["COS (5')"][DC]["P_type"], DR["COS (5')"][DC]["P_orient"], int(DR["COS (5')"][DC]["ave_whole_cov"]), P_comments]]
|
|
|
1349 t=Table(data, 2*[1.50*inch]+5*[0.80*inch]+1*[1.25*inch], 1*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-2),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1350 report.append(t)
|
|
|
1351
|
|
|
1352 for DC in DR["COS"]:
|
|
|
1353 P_comments = ""
|
|
|
1354 if int(DR["COS"][DC]["ave_whole_cov"]) < 50:
|
|
|
1355 P_comments = "Low coverage"
|
|
|
1356
|
|
|
1357 data = [[DC, DR["COS"][DC]["P_class"], DR["COS"][DC]["P_left"], DR["COS"][DC]["P_right"], DR["COS"][DC]["P_type"], DR["COS"][DC]["P_orient"], int(DR["COS"][DC]["ave_whole_cov"]), P_comments]]
|
|
|
1358 t=Table(data, 2*[1.50*inch]+5*[0.80*inch]+1*[1.25*inch], 1*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-2),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1359 report.append(t)
|
|
|
1360
|
|
|
1361 report.append(Spacer(1, 5))
|
|
|
1362
|
|
|
1363 # DTR Phages
|
|
|
1364 count_DTR = len(DR["DTR (short)"]) + len(DR["DTR (long)"])
|
|
|
1365 ptext = '<u><font size=14>DTR Phages ('+ str(count_DTR) +')</font></u>'
|
|
|
1366 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
1367 report.append(Spacer(1, 10))
|
|
|
1368
|
|
|
1369 if count_DTR != 0:
|
|
|
1370
|
|
|
1371 data = [["Name", "Class", "Left", "Right", "Type", "Orient", "Coverage", "Comments"]]
|
|
|
1372 t=Table(data, 2*[1.50*inch]+5*[0.80*inch]+1*[1.25*inch], 1*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-1),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1373 report.append(t)
|
|
|
1374
|
|
|
1375 for DC in DR["DTR (short)"]:
|
|
|
1376 P_comments = ""
|
|
|
1377 if int(DR["DTR (short)"][DC]["ave_whole_cov"]) < 50:
|
|
|
1378 P_comments = "Low coverage"
|
|
|
1379
|
|
|
1380 data = [[DC, DR["DTR (short)"][DC]["P_class"], DR["DTR (short)"][DC]["P_left"], DR["DTR (short)"][DC]["P_right"], DR["DTR (short)"][DC]["P_type"], DR["DTR (short)"][DC]["P_orient"], int(DR["DTR (short)"][DC]["ave_whole_cov"]), P_comments]]
|
|
|
1381 t=Table(data, 2*[1.50*inch]+5*[0.80*inch]+1*[1.25*inch], 1*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-2),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1382 report.append(t)
|
|
|
1383
|
|
|
1384 for DC in DR["DTR (long)"]:
|
|
|
1385 P_comments = ""
|
|
|
1386 if int(DR["DTR (long)"][DC]["ave_whole_cov"]) < 50:
|
|
|
1387 P_comments = "Low coverage"
|
|
|
1388
|
|
|
1389 data = [[DC, DR["DTR (long)"][DC]["P_class"], DR["DTR (long)"][DC]["P_left"], DR["DTR (long)"][DC]["P_right"], DR["DTR (long)"][DC]["P_type"], DR["DTR (long)"][DC]["P_orient"], int(DR["DTR (long)"][DC]["ave_whole_cov"]), P_comments]]
|
|
|
1390 t=Table(data, 2*[1.50*inch]+5*[0.80*inch]+1*[1.25*inch], 1*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-2),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1391 report.append(t)
|
|
|
1392
|
|
|
1393 report.append(Spacer(1, 5))
|
|
|
1394
|
|
|
1395 # Headful Phages
|
|
|
1396 count_Headful = len(DR["Headful (pac)"])
|
|
|
1397 ptext = '<u><font size=14>Headful Phages ('+ str(count_Headful) +')</font></u>'
|
|
|
1398 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
1399 report.append(Spacer(1, 10))
|
|
|
1400
|
|
|
1401 if count_Headful != 0:
|
|
|
1402
|
|
|
1403 data = [["Name", "Class", "Left", "Right", "Type", "Orient", "Coverage", "Comments"]]
|
|
|
1404 t=Table(data, 2*[1.50*inch]+5*[0.80*inch]+1*[1.25*inch], 1*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-1),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1405 report.append(t)
|
|
|
1406
|
|
|
1407 for DC in DR["Headful (pac)"]:
|
|
|
1408 P_comments = ""
|
|
|
1409 if int(DR["Headful (pac)"][DC]["ave_whole_cov"]) < 50:
|
|
|
1410 P_comments = "Low coverage"
|
|
|
1411
|
|
|
1412 data = [[DC, DR["Headful (pac)"][DC]["P_class"], DR["Headful (pac)"][DC]["P_left"], DR["Headful (pac)"][DC]["P_right"], DR["Headful (pac)"][DC]["P_type"], DR["Headful (pac)"][DC]["P_orient"], int(DR["Headful (pac)"][DC]["ave_whole_cov"]), P_comments]]
|
|
|
1413 t=Table(data, 2*[1.50*inch]+5*[0.80*inch]+1*[1.25*inch], 1*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-2),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1414 report.append(t)
|
|
|
1415
|
|
|
1416 report.append(Spacer(1, 5))
|
|
|
1417
|
|
|
1418 # OTHERS Phages
|
|
|
1419 count_Others = len(DR["Mu-like"]) + len(DR["UNKNOWN"]) + len(DR["NEW"])
|
|
|
1420 ptext = '<u><font size=14>Others Phages ('+ str(count_Others) +')</font></u>'
|
|
|
1421 report.append(Paragraph(ptext, styles["Left"]))
|
|
|
1422 report.append(Spacer(1, 10))
|
|
|
1423
|
|
|
1424 if count_Others != 0:
|
|
|
1425
|
|
|
1426 data = [["Name", "Class", "Left", "Right", "Type", "Orient", "Coverage", "Comments"]]
|
|
|
1427 t=Table(data, 2*[1.50*inch]+5*[0.80*inch]+1*[1.25*inch], 1*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-1),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1428 report.append(t)
|
|
|
1429
|
|
|
1430 for DC in DR["Mu-like"]:
|
|
|
1431 P_comments = ""
|
|
|
1432 if int(DR["Mu-like"][DC]["ave_whole_cov"]) < 50:
|
|
|
1433 P_comments = "Low coverage"
|
|
|
1434
|
|
|
1435 data = [[DC, DR["Mu-like"][DC]["P_class"], DR["Mu-like"][DC]["P_left"], DR["Mu-like"][DC]["P_right"], DR["Mu-like"][DC]["P_type"], DR["Mu-like"][DC]["P_orient"], int(DR["Mu-like"][DC]["ave_whole_cov"]), P_comments]]
|
|
|
1436 t=Table(data, 2*[1.50*inch]+5*[0.80*inch]+1*[1.25*inch], 1*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-2),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1437 report.append(t)
|
|
|
1438
|
|
|
1439 for DC in DR["NEW"]:
|
|
|
1440 P_comments = ""
|
|
|
1441 if int(DR["NEW"][DC]["ave_whole_cov"]) < 50:
|
|
|
1442 P_comments = "Low coverage"
|
|
|
1443
|
|
|
1444 data = [[DC, DR["NEW"][DC]["P_class"], DR["NEW"][DC]["P_left"], DR["NEW"][DC]["P_right"], DR["NEW"][DC]["P_type"], DR["NEW"][DC]["P_orient"], int(DR["NEW"][DC]["ave_whole_cov"]), P_comments]]
|
|
|
1445 t=Table(data, 2*[1.50*inch]+5*[0.80*inch]+1*[1.25*inch], 1*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-2),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1446 report.append(t)
|
|
|
1447
|
|
|
1448 for DC in DR["UNKNOWN"]:
|
|
|
1449 P_comments = ""
|
|
|
1450 if int(DR["UNKNOWN"][DC]["ave_whole_cov"]) < 50:
|
|
|
1451 P_comments = "Low coverage"
|
|
|
1452
|
|
|
1453 data = [[DC, DR["UNKNOWN"][DC]["P_class"], DR["UNKNOWN"][DC]["P_left"], DR["UNKNOWN"][DC]["P_right"], DR["UNKNOWN"][DC]["P_type"], DR["UNKNOWN"][DC]["P_orient"], int(DR["UNKNOWN"][DC]["ave_whole_cov"]), P_comments]]
|
|
|
1454 t=Table(data, 2*[1.50*inch]+5*[0.80*inch]+1*[1.25*inch], 1*[0.25*inch], hAlign='CENTER', style=[('FONT',(0,0),(-1,-2),'Helvetica-Bold'), ('GRID',(0,0),(-1,-1),0.5,colors.grey), ('FONTSIZE',(0,0),(-1,-1),12), ('ALIGN',(0,0),(-1,-1),'CENTER'),('VALIGN',(0,0),(-1,-1),'MIDDLE')])
|
|
|
1455 report.append(t)
|
|
|
1456
|
|
|
1457 report.append(Spacer(1, 5))
|
|
|
1458
|
|
|
1459 report.append(PageBreak())
|
|
|
1460
|
|
|
1461 return report
|
|
|
1462
|
|
|
1463 def WorkflowReport(phagename, P_class, P_left, P_right, P_type, P_orient, ave_whole_cov, multi = 0, phage_plus_norm=None, phage_minus_norm=None,*args, **kwargs):
|
|
|
1464 """ Text report for each phage."""
|
|
|
1465
|
|
|
1466 P_comments = ""
|
|
|
1467 if ave_whole_cov < 50:
|
|
|
1468 P_comments = "WARNING: Low coverage"
|
|
|
1469
|
|
|
1470 if ave_whole_cov == 0:
|
|
|
1471 P_comments = "No read match"
|
|
|
1472
|
|
|
1473 if not multi:
|
|
|
1474 filoutWorkflow = open(phagename + "_workflow.txt", "w")
|
|
|
1475 filoutWorkflow.write("#phagename\tClass\tLeft\tRight\tType\tOrient\tCoverage\tComments\n")
|
|
|
1476 filoutWorkflow.write(phagename + "\t" + P_class + "\t" + str(P_left) + "\t" + str(P_right) + "\t" + P_type + "\t" + P_orient + "\t" + str(ave_whole_cov) + "\t" + P_comments + "\n")
|
|
|
1477 filoutWorkflow.close()
|
|
|
1478 else:
|
|
|
1479 pval_left_peak="-"
|
|
|
1480 pval_adj_left_peak="-"
|
|
|
1481 pval_right_peak="-"
|
|
|
1482 pval_adj_right_peak="-"
|
|
|
1483 if isinstance(P_left,np.int64):
|
|
|
1484 # get pvalue and adjusted pvalue for this + peak
|
|
|
1485 left_peak_infos=phage_plus_norm.loc[phage_plus_norm['Position']==P_left]
|
|
|
1486 pval_left_peak=left_peak_infos["pval_gamma"]
|
|
|
1487 pval_left_peak=pval_left_peak.values[0]
|
|
|
1488 pval_adj_left_peak=left_peak_infos["pval_gamma_adj"]
|
|
|
1489 pval_adj_left_peak =pval_adj_left_peak.values[0]
|
|
|
1490 if isinstance(P_right,np.int64):
|
|
|
1491 # get pvalue and adjusted pvalue for this + peak
|
|
|
1492 right_peak_infos=phage_minus_norm.loc[phage_minus_norm['Position']==P_right]
|
|
|
1493 pval_right_peak=right_peak_infos["pval_gamma"]
|
|
|
1494 pval_right_peak=pval_right_peak.values[0]
|
|
|
1495 pval_adj_right_peak=right_peak_infos["pval_gamma_adj"]
|
|
|
1496 pval_adj_right_peak=pval_adj_right_peak.values[0]
|
|
|
1497 return phagename + "\t" + P_class + "\t" + str(P_left) + "\t" +str(pval_left_peak)+ "\t" +str(pval_adj_left_peak)+\
|
|
|
1498 "\t" + str(P_right) + "\t" + str(pval_right_peak) + "\t" + str(pval_adj_right_peak)+ "\t" + P_type +\
|
|
|
1499 "\t" + P_orient + "\t" + str(ave_whole_cov) + "\t" + P_comments + "\n"
|
|
|
1500 return
|
|
|
1501
|
|
|
1502 def EstimateTime(secondes):
|
|
|
1503 """ Convert secondes into time."""
|
|
|
1504 conv = (86400,3600,60,1)
|
|
|
1505 result = [0,0,0,0]
|
|
|
1506 i=0
|
|
|
1507 while secondes>0:
|
|
|
1508 result[i]= secondes/conv[i]
|
|
|
1509 secondes=secondes-result[i]*conv[i]
|
|
|
1510 i+=1
|
|
|
1511 return str(result[0]) + " Days " + str(result[1]) + " Hrs " + str(result[2]) + " Min " + str(result[3]) + " Sec"
|
|
|
1512
|
|
|
1513
|
|
|
1514
|
|
|
1515
|
|
|
1516
|
|
|
1517
|
