|
0
|
1 #!/usr/bin/env python
|
|
|
2 # -*- coding: utf-8 -*-
|
|
|
3 # vim: set fileencoding=utf-8
|
|
|
4 import os
|
|
|
5 import sys
|
|
|
6 import json
|
|
|
7 import math
|
|
|
8 import numpy
|
|
|
9 import argparse
|
|
|
10 import itertools
|
|
|
11 import logging
|
|
|
12 from gff3 import (
|
|
|
13 feature_lambda,
|
|
|
14 coding_genes,
|
|
|
15 genes,
|
|
|
16 get_gff3_id,
|
|
|
17 feature_test_location,
|
|
|
18 get_rbs_from,
|
|
|
19 nice_name,
|
|
|
20 )
|
|
|
21 from shinefind import NaiveSDCaller
|
|
|
22 from CPT_GFFParser import gffParse, gffWrite, gffSeqFeature
|
|
|
23 from Bio import SeqIO
|
|
|
24 from Bio.SeqRecord import SeqRecord
|
|
|
25 from Bio.SeqFeature import SeqFeature, FeatureLocation
|
|
|
26 from jinja2 import Environment, FileSystemLoader
|
|
|
27 from cpt import MGAFinder
|
|
|
28
|
|
|
29 logging.basicConfig(level=logging.DEBUG)
|
|
|
30 log = logging.getLogger(name="pav")
|
|
|
31
|
|
|
32 # Path to script, required because of Galaxy.
|
|
|
33 SCRIPT_PATH = os.path.dirname(os.path.realpath(__file__))
|
|
|
34 # Path to the HTML template for the report
|
|
|
35
|
|
|
36 ENCOURAGEMENT = (
|
|
|
37 (100, "Perfection itself!"),
|
|
|
38 (90, "Amazing!"),
|
|
|
39 (80, "Not too bad, a few minor things to fix..."),
|
|
|
40 (70, "Some issues to address"),
|
|
|
41 (
|
|
|
42 50,
|
|
|
43 """Issues detected! </p><p class="text-muted">Have you heard of the
|
|
|
44 <a href="https://cpt.tamu.edu">CPT</a>\'s Automated Phage Annotation
|
|
|
45 Pipeline?""",
|
|
|
46 ),
|
|
|
47 (
|
|
|
48 0,
|
|
|
49 """<b>MAJOR</b> issues detected! Please consider using the
|
|
|
50 <a href="https://cpt.tamu.edu">CPT</a>\'s Automated Phage Annotation Pipeline""",
|
|
|
51 ),
|
|
|
52 )
|
|
|
53
|
|
|
54
|
|
|
55 def gen_qc_feature(start, end, message, strand=0, id_src=None, type_src="gene"):
|
|
|
56 kwargs = {"qualifiers": {"note": [message]}}
|
|
|
57 kwargs["type"] = type_src
|
|
|
58 kwargs["strand"] = strand
|
|
|
59 kwargs["phase"]=0
|
|
|
60 kwargs["score"]=0.0
|
|
|
61 kwargs["source"]="feature"
|
|
|
62 if id_src is not None:
|
|
|
63 kwargs["id"] = id_src.id
|
|
|
64 kwargs["qualifiers"]["ID"] = [id_src.id]
|
|
|
65 kwargs["qualifiers"]["Name"] = id_src.qualifiers.get("Name", [])
|
|
|
66
|
|
|
67
|
|
|
68 if end >= start:
|
|
|
69 return gffSeqFeature(FeatureLocation(start, end, strand=strand), **kwargs)
|
|
|
70 else:
|
|
|
71 return gffSeqFeature(FeatureLocation(end, start, strand=strand), **kwargs)
|
|
|
72
|
|
|
73
|
|
|
74 def __ensure_location_in_bounds(start=0, end=0, parent_length=0):
|
|
|
75 # This prevents frameshift errors
|
|
|
76 while start < 0:
|
|
|
77 start += 3
|
|
|
78 while end < 0:
|
|
|
79 end += 3
|
|
|
80 while start > parent_length:
|
|
|
81 start -= 3
|
|
|
82 while end > parent_length:
|
|
|
83 end -= 3
|
|
|
84 return (start, end)
|
|
|
85
|
|
|
86
|
|
|
87 def missing_rbs(record, lookahead_min=5, lookahead_max=15):
|
|
|
88 """
|
|
|
89 Identify gene features with missing RBSs
|
|
|
90
|
|
|
91 This "looks ahead" 5-15 bases ahead of each gene feature, and checks if
|
|
|
92 there's an RBS feature in those bounds.
|
|
|
93
|
|
|
94 The returned data is a set of genes with the RBS sequence in the __upstream
|
|
|
95 attribute, and a message in the __message attribute.
|
|
|
96 """
|
|
|
97 results = []
|
|
|
98 good = 0
|
|
|
99 bad = 0
|
|
|
100 qc_features = []
|
|
|
101 sd_finder = NaiveSDCaller()
|
|
|
102
|
|
|
103 any_rbss = False
|
|
|
104
|
|
|
105 for gene in coding_genes(record.features):
|
|
|
106 # Check if there are RBSs, TODO: make this recursive. Each feature in
|
|
|
107 # gene.sub_features can also have sub_features.
|
|
|
108 rbss = get_rbs_from(gene)
|
|
|
109 # No RBS found
|
|
|
110 if len(rbss) == 0:
|
|
|
111 # Get the sequence lookahead_min to lookahead_max upstream
|
|
|
112 if gene.strand > 0:
|
|
|
113 start = gene.location.start - lookahead_max
|
|
|
114 end = gene.location.start - lookahead_min
|
|
|
115 else:
|
|
|
116 start = gene.location.end + lookahead_min
|
|
|
117 end = gene.location.end + lookahead_max
|
|
|
118 # We have to ensure the feature is ON the genome, otherwise we may
|
|
|
119 # be trying to access a location outside of the length of the
|
|
|
120 # genome, which would be bad.
|
|
|
121 (start, end) = __ensure_location_in_bounds(
|
|
|
122 start=start, end=end, parent_length=len(record)
|
|
|
123 )
|
|
|
124 # Temporary feature to extract sequence
|
|
|
125 tmp = gffSeqFeature(
|
|
|
126 FeatureLocation(start, end, strand=gene.strand), type="domain"
|
|
|
127 )
|
|
|
128 # Get the sequence
|
|
|
129 seq = str(tmp.extract(record.seq))
|
|
|
130 # Set the default properties
|
|
|
131 gene.__upstream = seq.lower()
|
|
|
132 gene.__message = "No RBS annotated, None found"
|
|
|
133
|
|
|
134 # Try and do an automated shinefind call
|
|
|
135 sds = sd_finder.list_sds(seq)
|
|
|
136 if len(sds) > 0:
|
|
|
137 sd = sds[0]
|
|
|
138 gene.__upstream = sd_finder.highlight_sd(
|
|
|
139 seq.lower(), sd["start"], sd["end"]
|
|
|
140 )
|
|
|
141 gene.__message = "Unannotated but valid RBS"
|
|
|
142
|
|
|
143 qc_features.append(
|
|
|
144 gen_qc_feature(
|
|
|
145 start, end, "Missing RBS", strand=gene.strand, id_src=gene, type_src="gene"
|
|
|
146 )
|
|
|
147 )
|
|
|
148
|
|
|
149 bad += 1
|
|
|
150 results.append(gene)
|
|
|
151 results[-1].location = FeatureLocation(results[-1].location.start + 1, results[-1].location.end, results[-1].location.strand)
|
|
|
152 else:
|
|
|
153 if len(rbss) > 1:
|
|
|
154 log.warn("%s RBSs found for gene %s", rbss[0].id, get_gff3_id(gene))
|
|
|
155 any_rbss = True
|
|
|
156 # get first RBS/CDS
|
|
|
157 cds = list(genes(gene.sub_features, feature_type="CDS"))[0]
|
|
|
158 rbs = rbss[0]
|
|
|
159
|
|
|
160 # Get the distance between the two
|
|
|
161 if gene.strand > 0:
|
|
|
162 distance = cds.location.start - rbs.location.end
|
|
|
163 else:
|
|
|
164 distance = rbs.location.start - cds.location.end
|
|
|
165
|
|
|
166 # If the RBS is too far away, annotate that
|
|
|
167 if distance > lookahead_max:
|
|
|
168 gene.__message = "RBS too far away (%s nt)" % distance
|
|
|
169
|
|
|
170 qc_features.append(
|
|
|
171 gen_qc_feature(
|
|
|
172 rbs.location.start,
|
|
|
173 rbs.location.end,
|
|
|
174 gene.__message,
|
|
|
175 strand=gene.strand,
|
|
|
176 id_src=gene,
|
|
|
177 type_src="gene"
|
|
|
178 )
|
|
|
179 )
|
|
|
180
|
|
|
181 bad += 1
|
|
|
182 results.append(gene)
|
|
|
183 results[-1].location = FeatureLocation(results[-1].location.start + 1, results[-1].location.end, results[-1].location.strand)
|
|
|
184 else:
|
|
|
185 good += 1
|
|
|
186
|
|
|
187 return good, bad, results, qc_features, any_rbss
|
|
|
188
|
|
|
189
|
|
|
190 # modified from get_orfs_or_cdss.py
|
|
|
191 # -----------------------------------------------------------
|
|
|
192
|
|
|
193
|
|
|
194 def require_sd(data, record, chrom_start, sd_min, sd_max):
|
|
|
195 sd_finder = NaiveSDCaller()
|
|
|
196 for putative_gene in data:
|
|
|
197 if putative_gene[2] > 0: # strand
|
|
|
198 start = chrom_start + putative_gene[0] - sd_max
|
|
|
199 end = chrom_start + putative_gene[0] - sd_min
|
|
|
200 else:
|
|
|
201 start = chrom_start + putative_gene[1] + sd_min
|
|
|
202 end = chrom_start + putative_gene[1] + sd_max
|
|
|
203
|
|
|
204 (start, end) = __ensure_location_in_bounds(
|
|
|
205 start=start, end=end, parent_length=len(record)
|
|
|
206 )
|
|
|
207 tmp = gffSeqFeature(
|
|
|
208 FeatureLocation(start, end, strand=putative_gene[2]), type="domain"
|
|
|
209 )
|
|
|
210 # Get the sequence
|
|
|
211 seq = str(tmp.extract(record.seq))
|
|
|
212 sds = sd_finder.list_sds(seq)
|
|
|
213 if len(sds) > 0:
|
|
|
214 yield putative_gene + (start, end)
|
|
|
215
|
|
|
216
|
|
|
217 def excessive_gap(
|
|
|
218 record,
|
|
|
219 excess=50,
|
|
|
220 excess_divergent=200,
|
|
|
221 min_gene=30,
|
|
|
222 slop=30,
|
|
|
223 lookahead_min=5,
|
|
|
224 lookahead_max=15,
|
|
|
225 ):
|
|
|
226 """
|
|
|
227 Identify excessive gaps between gene features.
|
|
|
228
|
|
|
229 Default "excessive" gap size is 10, but that should likely be larger.
|
|
|
230 """
|
|
|
231 results = []
|
|
|
232 good = 0
|
|
|
233 bad = 0
|
|
|
234
|
|
|
235 contiguous_regions = []
|
|
|
236
|
|
|
237 sorted_genes = sorted(
|
|
|
238 genes(record.features), key=lambda feature: feature.location.start
|
|
|
239 )
|
|
|
240 if len(sorted_genes) == 0:
|
|
|
241 log.warn("NO GENES FOUND")
|
|
|
242 return good, bad, results, []
|
|
|
243
|
|
|
244 current_gene = None
|
|
|
245 for gene in sorted_genes:
|
|
|
246 # If the gene's start is contiguous to the "current_gene", then we
|
|
|
247 # extend current_gene
|
|
|
248 for cds in genes(gene.sub_features, feature_type="CDS"):
|
|
|
249 if current_gene is None:
|
|
|
250 current_gene = [int(cds.location.start), int(cds.location.end)]
|
|
|
251
|
|
|
252 if cds.location.start <= current_gene[1] + excess:
|
|
|
253 # Don't want to decrease size
|
|
|
254 if int(cds.location.end) >= current_gene[1]:
|
|
|
255 current_gene[1] = int(cds.location.end)
|
|
|
256 else:
|
|
|
257 # If it's discontiguous, we append the region and clear.
|
|
|
258 contiguous_regions.append(current_gene)
|
|
|
259 current_gene = [int(cds.location.start), int(cds.location.end)]
|
|
|
260
|
|
|
261 # This generally expected that annotations would NOT continue unto the end
|
|
|
262 # of the genome, however that's a bug, and we can make it here with an
|
|
|
263 # empty contiguous_regions list
|
|
|
264 contiguous_regions.append(current_gene)
|
|
|
265
|
|
|
266 for i in range(len(contiguous_regions) + 1):
|
|
|
267 if i == 0:
|
|
|
268 a = (1, 1)
|
|
|
269 b = contiguous_regions[i]
|
|
|
270 elif i >= len(contiguous_regions):
|
|
|
271 a = contiguous_regions[i - 1]
|
|
|
272 b = (len(record.seq), None)
|
|
|
273 else:
|
|
|
274 a = contiguous_regions[i - 1]
|
|
|
275 b = contiguous_regions[i]
|
|
|
276
|
|
|
277 gap_size = abs(b[0] - a[1])
|
|
|
278
|
|
|
279 if gap_size > min(excess, excess_divergent):
|
|
|
280 a_feat_l = itertools.islice(
|
|
|
281 feature_lambda(
|
|
|
282 sorted_genes,
|
|
|
283 feature_test_location,
|
|
|
284 {"loc": a[1]},
|
|
|
285 subfeatures=False,
|
|
|
286 ),
|
|
|
287 1,
|
|
|
288 )
|
|
|
289 b_feat_l = itertools.islice(
|
|
|
290 feature_lambda(
|
|
|
291 sorted_genes,
|
|
|
292 feature_test_location,
|
|
|
293 {"loc": b[0]},
|
|
|
294 subfeatures=False,
|
|
|
295 ),
|
|
|
296 1,
|
|
|
297 )
|
|
|
298
|
|
|
299 try:
|
|
|
300 a_feat = next(a_feat_l)
|
|
|
301 except StopIteration:
|
|
|
302 # Triggers on end of genome
|
|
|
303 a_feat = None
|
|
|
304 try:
|
|
|
305 b_feat = next(b_feat_l)
|
|
|
306 except StopIteration:
|
|
|
307 # Triggers on end of genome
|
|
|
308 b_feat = None
|
|
|
309
|
|
|
310 result_obj = [
|
|
|
311 a[1],
|
|
|
312 b[0],
|
|
|
313 None if not a_feat else a_feat.location.strand,
|
|
|
314 None if not b_feat else b_feat.location.strand,
|
|
|
315 ]
|
|
|
316
|
|
|
317 if a_feat is None or b_feat is None:
|
|
|
318 if gap_size > excess_divergent:
|
|
|
319 results.append(result_obj)
|
|
|
320 else:
|
|
|
321 if (
|
|
|
322 a_feat.location.strand == b_feat.location.strand
|
|
|
323 and gap_size > excess
|
|
|
324 ):
|
|
|
325 results.append(result_obj)
|
|
|
326 elif (
|
|
|
327 a_feat.location.strand != b_feat.location.strand
|
|
|
328 and gap_size > excess_divergent
|
|
|
329 ):
|
|
|
330 results.append(result_obj)
|
|
|
331
|
|
|
332 better_results = []
|
|
|
333 qc_features = []
|
|
|
334 of = MGAFinder(11, "CDS", "closed", min_gene)
|
|
|
335 # of = OrfFinder(11, 'CDS', 'closed', min_gene)
|
|
|
336
|
|
|
337 for result_obj in results:
|
|
|
338 start = result_obj[0]
|
|
|
339 end = result_obj[1]
|
|
|
340 f = gen_qc_feature(start, end, "Excessive gap, %s bases" % abs(end - start), type_src="gene")
|
|
|
341 qc_features.append(f)
|
|
|
342 putative_genes = of.putative_genes_in_sequence(
|
|
|
343 str(record[start - slop : end + slop].seq)
|
|
|
344 )
|
|
|
345 putative_genes = list(
|
|
|
346 require_sd(putative_genes, record, start, lookahead_min, lookahead_max)
|
|
|
347 )
|
|
|
348 for putative_gene in putative_genes:
|
|
|
349 # (0, 33, 1, 'ATTATTTTATCAAAACGCTTTACAATCTTTTAG', 'MILSKRFTIF', 123123, 124324)
|
|
|
350 possible_gene_start = start + putative_gene[0]
|
|
|
351 possible_gene_end = start + putative_gene[1]
|
|
|
352
|
|
|
353 if possible_gene_start <= possible_gene_end:
|
|
|
354 possible_cds = gffSeqFeature(
|
|
|
355 FeatureLocation(
|
|
|
356 possible_gene_start, possible_gene_end, strand=putative_gene[2]
|
|
|
357 ),
|
|
|
358 type="CDS",
|
|
|
359 )
|
|
|
360 else:
|
|
|
361 possible_cds = gffSeqFeature(
|
|
|
362 FeatureLocation(
|
|
|
363 possible_gene_end, possible_gene_start, strand=putative_gene[2],
|
|
|
364 ),
|
|
|
365 type="CDS",
|
|
|
366 )
|
|
|
367
|
|
|
368 # Now we adjust our boundaries for the RBS that's required
|
|
|
369 # There are only two cases, the rbs is upstream of it, or downstream
|
|
|
370 if putative_gene[5] < possible_gene_start:
|
|
|
371 possible_gene_start = putative_gene[5]
|
|
|
372 else:
|
|
|
373 possible_gene_end = putative_gene[6]
|
|
|
374
|
|
|
375 if putative_gene[5] <= putative_gene[6]:
|
|
|
376 possible_rbs = gffSeqFeature(
|
|
|
377 FeatureLocation(
|
|
|
378 putative_gene[5], putative_gene[6], strand=putative_gene[2]
|
|
|
379 ),
|
|
|
380 type="Shine_Dalgarno_sequence",
|
|
|
381 )
|
|
|
382 else:
|
|
|
383 possible_rbs = gffSeqFeature(
|
|
|
384 FeatureLocation(
|
|
|
385 putative_gene[6], putative_gene[5], strand=putative_gene[2],
|
|
|
386 ),
|
|
|
387 type="Shine_Dalgarno_sequence",
|
|
|
388 )
|
|
|
389
|
|
|
390 if possible_gene_start <= possible_gene_end:
|
|
|
391 possible_gene = gffSeqFeature(
|
|
|
392 FeatureLocation(
|
|
|
393 possible_gene_start, possible_gene_end, strand=putative_gene[2]
|
|
|
394 ),
|
|
|
395 type="gene",
|
|
|
396 qualifiers={"note": ["Possible gene"]},
|
|
|
397 )
|
|
|
398 else:
|
|
|
399 possible_gene = gffSeqFeature(
|
|
|
400 FeatureLocation(
|
|
|
401 possible_gene_end, possible_gene_start, strand=putative_gene[2],
|
|
|
402 ),
|
|
|
403 type="gene",
|
|
|
404 qualifiers={"note": ["Possible gene"]},
|
|
|
405 )
|
|
|
406 possible_gene.sub_features = [possible_rbs, possible_cds]
|
|
|
407 qc_features.append(possible_gene)
|
|
|
408
|
|
|
409 better_results.append(result_obj + [len(putative_genes)])
|
|
|
410
|
|
|
411 # Bad gaps are those with more than zero possible genes found
|
|
|
412 bad = len([x for x in better_results if x[2] > 0])
|
|
|
413 # Generally taking "good" here as every possible gap in the genome
|
|
|
414 # Thus, good is TOTAL - gaps
|
|
|
415 good = len(sorted_genes) + 1 - bad
|
|
|
416 # and bad is just gaps
|
|
|
417 return good, bad, better_results, qc_features
|
|
|
418
|
|
|
419
|
|
|
420 def phi(x):
|
|
|
421 """Standard phi function used in calculation of normal distribution"""
|
|
|
422 return math.exp(-1 * math.pi * x * x)
|
|
|
423
|
|
|
424
|
|
|
425 def norm(x, mean=0, sd=1):
|
|
|
426 """
|
|
|
427 Normal distribution. Given an x position, a mean, and a standard
|
|
|
428 deviation, calculate the "y" value. Useful for score scaling
|
|
|
429
|
|
|
430 Modified to multiply by SD. This means even at sd=5, norm(x, mean) where x = mean => 1, rather than 1/5.
|
|
|
431 """
|
|
|
432 return (1 / float(sd)) * phi(float(x - mean) / float(sd)) * sd
|
|
|
433
|
|
|
434
|
|
|
435 def coding_density(record, mean=92.5, sd=20):
|
|
|
436 """
|
|
|
437 Find coding density in the genome
|
|
|
438 """
|
|
|
439 feature_lengths = 0
|
|
|
440
|
|
|
441 for gene_a in coding_genes(record.features):
|
|
|
442 feature_lengths += sum(
|
|
|
443 [len(x) for x in genes(gene_a.sub_features, feature_type="CDS")]
|
|
|
444 )
|
|
|
445
|
|
|
446 avgFeatLen = float(feature_lengths) / float(len(record.seq))
|
|
|
447 return int(norm(100 * avgFeatLen, mean=mean, sd=sd) * 100), int(100 * avgFeatLen)
|
|
|
448
|
|
|
449
|
|
|
450 def exact_coding_density(record, mean=92.5, sd=20):
|
|
|
451 """
|
|
|
452 Find exact coding density in the genome
|
|
|
453 """
|
|
|
454 data = numpy.zeros(len(record.seq))
|
|
|
455
|
|
|
456 for gene_a in coding_genes(record.features):
|
|
|
457 for cds in genes(gene_a.sub_features, feature_type="CDS"):
|
|
|
458 for i in range(cds.location.start, cds.location.end + 1):
|
|
|
459 data[i - 1] = 1
|
|
|
460
|
|
|
461 return float(sum(data)) / len(data)
|
|
|
462
|
|
|
463
|
|
|
464 def excessive_overlap(record, excess=15, excess_divergent=30):
|
|
|
465 """
|
|
|
466 Find excessive overlaps in the genome, where excessive is defined as 15
|
|
|
467 bases for same strand, and 30 for divergent translation.
|
|
|
468
|
|
|
469 Does a product of all the top-level features in the genome, and calculates
|
|
|
470 gaps.
|
|
|
471 """
|
|
|
472 results = []
|
|
|
473 bad = 0
|
|
|
474 qc_features = []
|
|
|
475
|
|
|
476 for (gene_a, gene_b) in itertools.combinations(coding_genes(record.features), 2):
|
|
|
477 # Get the CDS from the subfeature list.
|
|
|
478 # TODO: not recursive.
|
|
|
479 cds_a = [x for x in genes(gene_a.sub_features, feature_type="CDS")]
|
|
|
480 cds_b = [x for x in genes(gene_b.sub_features, feature_type="CDS")]
|
|
|
481
|
|
|
482 if len(cds_a) == 0:
|
|
|
483 log.warn("Gene missing subfeatures; %s", get_gff3_id(gene_a))
|
|
|
484 continue
|
|
|
485
|
|
|
486 if len(cds_b) == 0:
|
|
|
487 log.warn("Gene missing subfeatures; %s", get_gff3_id(gene_b))
|
|
|
488 continue
|
|
|
489
|
|
|
490 cds_a = cds_a[0]
|
|
|
491 cds_b = cds_b[0]
|
|
|
492
|
|
|
493 # Set of locations that are included in the CDS of A and the
|
|
|
494 # CDS of B
|
|
|
495 cas = set(range(cds_a.location.start, cds_a.location.end))
|
|
|
496 cbs = set(range(cds_b.location.start, cds_b.location.end))
|
|
|
497
|
|
|
498 # Here we calculate the intersection between the two sets, and
|
|
|
499 # if it's larger than our excessive size, we know that they're
|
|
|
500 # overlapped
|
|
|
501 ix = cas.intersection(cbs)
|
|
|
502
|
|
|
503 if (cds_a.location.strand == cds_b.location.strand and len(ix) >= excess) or (
|
|
|
504 cds_a.location.strand != cds_b.location.strand
|
|
|
505 and len(ix) >= excess_divergent
|
|
|
506 ):
|
|
|
507 bad += float(len(ix)) / float(min(excess, excess_divergent))
|
|
|
508 qc_features.append(
|
|
|
509 gen_qc_feature(min(ix), max(ix), "Excessive Overlap", id_src=gene_a, type_src="gene")
|
|
|
510 )
|
|
|
511 results.append((gene_a, gene_b, min(ix), max(ix)))
|
|
|
512
|
|
|
513 # Good isn't accurate here. It's a triangle number and just ugly, but we
|
|
|
514 # don't care enough to fix it.
|
|
|
515 good = len(list(coding_genes(record.features)))
|
|
|
516 good = int(good - bad)
|
|
|
517 if good < 0:
|
|
|
518 good = 0
|
|
|
519 return good, int(bad), results, qc_features
|
|
|
520
|
|
|
521
|
|
|
522 def get_encouragement(score):
|
|
|
523 """Some text telling the user how they did
|
|
|
524 """
|
|
|
525 for encouragement in ENCOURAGEMENT:
|
|
|
526 if score > encouragement[0]:
|
|
|
527 return encouragement[1]
|
|
|
528 return ENCOURAGEMENT[-1][1]
|
|
|
529
|
|
|
530
|
|
|
531 def genome_overview(record):
|
|
|
532 """Genome overview
|
|
|
533 """
|
|
|
534 data = {
|
|
|
535 "genes": {
|
|
|
536 "count": 0,
|
|
|
537 "bases": len(record.seq),
|
|
|
538 "density": 0, # genes / kb
|
|
|
539 "avg_len": [],
|
|
|
540 "comp": {"A": 0, "C": 0, "G": 0, "T": 0},
|
|
|
541 },
|
|
|
542 "overall": {
|
|
|
543 "comp": {
|
|
|
544 "A": record.seq.count("A") + record.seq.count("a"),
|
|
|
545 "C": record.seq.count("C") + record.seq.count("c"),
|
|
|
546 "G": record.seq.count("G") + record.seq.count("g"),
|
|
|
547 "T": record.seq.count("T") + record.seq.count("t"),
|
|
|
548 },
|
|
|
549 "gc": 0,
|
|
|
550 },
|
|
|
551 }
|
|
|
552 gene_features = list(coding_genes(record.features))
|
|
|
553 data["genes"]["count"] = len(gene_features)
|
|
|
554
|
|
|
555 for feat in gene_features:
|
|
|
556 data["genes"]["comp"]["A"] += feat.extract(record).seq.count("A") + feat.extract(record).seq.count("a")
|
|
|
557 data["genes"]["comp"]["C"] += feat.extract(record).seq.count("C") + feat.extract(record).seq.count("c")
|
|
|
558 data["genes"]["comp"]["T"] += feat.extract(record).seq.count("T") + feat.extract(record).seq.count("t")
|
|
|
559 data["genes"]["comp"]["G"] += feat.extract(record).seq.count("G") + feat.extract(record).seq.count("g")
|
|
|
560 #data["genes"]["bases"] += len(feat)
|
|
|
561 data["genes"]["avg_len"].append(len(feat))
|
|
|
562
|
|
|
563 data["genes"]["avg_len"] = float(sum(data["genes"]["avg_len"])) / len(gene_features)
|
|
|
564 data["overall"]["gc"] = float(
|
|
|
565 data["overall"]["comp"]["G"] + data["overall"]["comp"]["C"]
|
|
|
566 ) / len(record.seq)
|
|
|
567 return data
|
|
|
568
|
|
|
569
|
|
|
570 def find_morons(record):
|
|
|
571 """Locate morons in the genome
|
|
|
572
|
|
|
573 Don't even know why...
|
|
|
574
|
|
|
575 TODO: remove? Idk.
|
|
|
576 """
|
|
|
577 results = []
|
|
|
578 good = 0
|
|
|
579 bad = 0
|
|
|
580
|
|
|
581 gene_features = list(coding_genes(record.features))
|
|
|
582 for i, gene in enumerate(gene_features):
|
|
|
583 two_left = gene_features[i - 2 : i]
|
|
|
584 two_right = gene_features[i + 1 : i + 1 + 2]
|
|
|
585 strands = [x.strand for x in two_left] + [x.strand for x in two_right]
|
|
|
586 anticon = [x for x in strands if x != gene.strand]
|
|
|
587
|
|
|
588 if len(anticon) == 4:
|
|
|
589 has_rbs = [x.type == "Shine_Dalgarno_sequence" for x in gene.sub_features]
|
|
|
590 if any(has_rbs):
|
|
|
591 rbs = [
|
|
|
592 x for x in gene.sub_features if x.type == "Shine_Dalgarno_sequence"
|
|
|
593 ][0]
|
|
|
594 rbs_msg = str(rbs.extract(record.seq))
|
|
|
595 else:
|
|
|
596 rbs_msg = "No RBS Available"
|
|
|
597 results.append((gene, two_left, two_right, rbs_msg))
|
|
|
598 bad += 1
|
|
|
599 else:
|
|
|
600 good += 1
|
|
|
601 return good, bad, results, []
|
|
|
602
|
|
|
603
|
|
|
604 def bad_gene_model(record):
|
|
|
605 """Find features without product
|
|
|
606 """
|
|
|
607 results = []
|
|
|
608 good = 0
|
|
|
609 bad = 0
|
|
|
610 qc_features = []
|
|
|
611
|
|
|
612 for gene in coding_genes(record.features):
|
|
|
613 exons = [
|
|
|
614 x for x in genes(gene.sub_features, feature_type="exon") if len(x) > 10
|
|
|
615 ]
|
|
|
616 CDSs = [x for x in genes(gene.sub_features, feature_type="CDS")]
|
|
|
617 if len(exons) >= 1 and len(CDSs) >= 1:
|
|
|
618 if len(exons) != len(CDSs):
|
|
|
619 results.append(
|
|
|
620 (
|
|
|
621 get_gff3_id(gene),
|
|
|
622 None,
|
|
|
623 None,
|
|
|
624 "Mismatched number of exons and CDSs in gff3 representation",
|
|
|
625 )
|
|
|
626 )
|
|
|
627 qc_features.append(
|
|
|
628 gen_qc_feature(
|
|
|
629 gene.location.start,
|
|
|
630 gene.location.end,
|
|
|
631 "Mismatched number of exons and CDSs in gff3 representation",
|
|
|
632 strand=gene.strand,
|
|
|
633 id_src=gene,
|
|
|
634 type_src="gene"
|
|
|
635 )
|
|
|
636 )
|
|
|
637 bad += 1
|
|
|
638 else:
|
|
|
639 for (exon, cds) in zip(
|
|
|
640 sorted(exons, key=lambda x: x.location.start),
|
|
|
641 sorted(CDSs, key=lambda x: x.location.start),
|
|
|
642 ):
|
|
|
643 if len(exon) != len(cds):
|
|
|
644 results.append(
|
|
|
645 (
|
|
|
646 get_gff3_id(gene),
|
|
|
647 exon,
|
|
|
648 cds,
|
|
|
649 "CDS does not extend to full length of gene",
|
|
|
650 )
|
|
|
651 )
|
|
|
652 qc_features.append(
|
|
|
653 gen_qc_feature(
|
|
|
654 exon.location.start,
|
|
|
655 exon.location.end,
|
|
|
656 "CDS does not extend to full length of gene",
|
|
|
657 strand=exon.strand,
|
|
|
658 id_src=gene,
|
|
|
659 type_src="CDS"
|
|
|
660 )
|
|
|
661 )
|
|
|
662 bad += 1
|
|
|
663 else:
|
|
|
664 good += 1
|
|
|
665 else:
|
|
|
666 log.warn("Could not handle %s, %s", exons, CDSs)
|
|
|
667 results.append(
|
|
|
668 (
|
|
|
669 get_gff3_id(gene),
|
|
|
670 None,
|
|
|
671 None,
|
|
|
672 "{0} exons, {1} CDSs".format(len(exons), len(CDSs)),
|
|
|
673 )
|
|
|
674 )
|
|
|
675
|
|
|
676 return good, len(results) + bad, results, qc_features
|
|
|
677
|
|
|
678
|
|
|
679 def weird_starts(record):
|
|
|
680 """Find features without product
|
|
|
681 """
|
|
|
682 good = 0
|
|
|
683 bad = 0
|
|
|
684 qc_features = []
|
|
|
685 results = []
|
|
|
686
|
|
|
687 overall = {}
|
|
|
688 for gene in coding_genes(record.features):
|
|
|
689 seq = [x for x in genes(gene.sub_features, feature_type="CDS")]
|
|
|
690 if len(seq) == 0:
|
|
|
691 log.warn("No CDS for gene %s", get_gff3_id(gene))
|
|
|
692 continue
|
|
|
693 else:
|
|
|
694 seq = seq[0]
|
|
|
695
|
|
|
696 seq_str = str(seq.extract(record.seq))
|
|
|
697 start_codon = seq_str[0:3]
|
|
|
698 if len(seq_str) < 3:
|
|
|
699 sys.stderr.write("Fatal Error: CDS of length less than 3 at " + str(seq.location) + '\n')
|
|
|
700 exit(2)
|
|
|
701 # if len(seq_str) % 3 != 0:
|
|
|
702 # if len(seq_str) < 3:
|
|
|
703 # stop_codon = seq_str[-(len(seq_str))]
|
|
|
704 # else:
|
|
|
705 # stop_codon = seq_str[-3]
|
|
|
706 #
|
|
|
707 # log.warn("CDS at %s length is not a multiple of three (Length = %d)", get_gff3_id(gene), len(seq_str))
|
|
|
708 # seq.__error = "Bad CDS Length"
|
|
|
709 # results.append(seq)
|
|
|
710 # qc_features.append(
|
|
|
711 # gen_qc_feature(
|
|
|
712 # s, e, "Bad Length", strand=seq.strand, id_src=gene
|
|
|
713 # )
|
|
|
714 # )
|
|
|
715 # bad += 1
|
|
|
716 # seq.__start = start_codon
|
|
|
717 # seq.__stop = stop_codon
|
|
|
718 # continue
|
|
|
719
|
|
|
720 stop_codon = seq_str[-3]
|
|
|
721 seq.__start = start_codon
|
|
|
722 seq.__stop = stop_codon
|
|
|
723 if start_codon not in overall:
|
|
|
724 overall[start_codon] = 1
|
|
|
725 else:
|
|
|
726 overall[start_codon] += 1
|
|
|
727
|
|
|
728 if start_codon not in ("ATG", "TTG", "GTG"):
|
|
|
729 log.warn("Weird start codon (%s) on %s", start_codon, get_gff3_id(gene))
|
|
|
730 seq.__error = "Unusual start codon %s" % start_codon
|
|
|
731
|
|
|
732 s = 0
|
|
|
733 e = 0
|
|
|
734 if seq.strand > 0:
|
|
|
735 s = seq.location.start
|
|
|
736 e = seq.location.start + 3
|
|
|
737 else:
|
|
|
738 s = seq.location.end
|
|
|
739 e = seq.location.end - 3
|
|
|
740
|
|
|
741 results.append(seq)
|
|
|
742 results[-1].location = FeatureLocation(results[-1].location.start + 1, results[-1].location.end, results[-1].location.strand)
|
|
|
743 qc_features.append(
|
|
|
744 gen_qc_feature(
|
|
|
745 s, e, "Weird start codon", strand=seq.strand, id_src=gene, type_src="gene"
|
|
|
746 )
|
|
|
747 )
|
|
|
748 bad += 1
|
|
|
749 else:
|
|
|
750 good += 1
|
|
|
751
|
|
|
752 return good, bad, results, qc_features, overall
|
|
|
753
|
|
|
754
|
|
|
755 def missing_genes(record):
|
|
|
756 """Find features without product
|
|
|
757 """
|
|
|
758 results = []
|
|
|
759 good = 0
|
|
|
760 bad = 0
|
|
|
761 qc_features = []
|
|
|
762
|
|
|
763 for gene in coding_genes(record.features):
|
|
|
764 if gene.qualifiers.get("cpt_source", [None])[0] == "CPT_GENE_MODEL_CORRECTION":
|
|
|
765 results.append(gene)
|
|
|
766 bad += 1
|
|
|
767 else:
|
|
|
768 good += 1
|
|
|
769
|
|
|
770 return good, bad, results, qc_features
|
|
|
771
|
|
|
772
|
|
|
773 def gene_model_correction_issues(record):
|
|
|
774 """Find features that have issues from the gene model correction step.
|
|
|
775 These have qualifiers beginning with CPT_GMS
|
|
|
776 """
|
|
|
777 results = []
|
|
|
778 good = 0
|
|
|
779 bad = 0
|
|
|
780 qc_features = []
|
|
|
781
|
|
|
782 # For each gene
|
|
|
783 for gene in coding_genes(record.features):
|
|
|
784 # Get the list of child CDSs
|
|
|
785 cdss = [x for x in genes(gene.sub_features, feature_type="CDS")]
|
|
|
786 # And our matching qualifiers
|
|
|
787 gene_data = [(k, v) for (k, v) in gene.qualifiers.items() if k == "cpt_gmc"]
|
|
|
788 # If there are problems with ONLY the parent, let's complain
|
|
|
789 local_results = []
|
|
|
790 local_qc_features = []
|
|
|
791 for x in gene_data:
|
|
|
792 if "Missing Locus Tag" in x[1]:
|
|
|
793 # Missing locus tag is an either or thing, if it hits here
|
|
|
794 # there shouldn't be anything else wrong with it.
|
|
|
795
|
|
|
796 # Obviously missing so we remove it
|
|
|
797 gene.qualifiers["locus_tag"] = [""]
|
|
|
798 # Translation from bp_genbank2gff3.py
|
|
|
799 cdss[0].qualifiers["locus_tag"] = cdss[0].qualifiers["Name"]
|
|
|
800 # Append our results
|
|
|
801 local_results.append((gene, cdss[0], "Gene is missing a locus_tag"))
|
|
|
802 local_qc_features.append(
|
|
|
803 gen_qc_feature(
|
|
|
804 gene.location.start,
|
|
|
805 gene.location.end,
|
|
|
806 "Gene is missing a locus_tag",
|
|
|
807 strand=gene.strand,
|
|
|
808 type_src="gene"
|
|
|
809 )
|
|
|
810 )
|
|
|
811
|
|
|
812 # We need to alert on any child issues as well.
|
|
|
813 for cds in cdss:
|
|
|
814 cds_data = [
|
|
|
815 (k, v[0]) for (k, v) in cds.qualifiers.items() if k == "cpt_gmc"
|
|
|
816 ]
|
|
|
817 if len(gene_data) == 0 and len(cds_data) == 0:
|
|
|
818 # Alles gut
|
|
|
819 pass
|
|
|
820 else:
|
|
|
821 for _, problem in cds_data:
|
|
|
822 if problem == "BOTH Missing Locus Tag":
|
|
|
823 gene.qualifiers["locus_tag"] = [""]
|
|
|
824 cds.qualifiers["locus_tag"] = [""]
|
|
|
825 local_results.append(
|
|
|
826 (gene, cds, "Both gene and CDS are missing locus tags")
|
|
|
827 )
|
|
|
828 local_qc_features.append(
|
|
|
829 gen_qc_feature(
|
|
|
830 cds.location.start,
|
|
|
831 cds.location.end,
|
|
|
832 "CDS is missing a locus_tag",
|
|
|
833 strand=cds.strand,
|
|
|
834 type_src="CDS"
|
|
|
835 )
|
|
|
836 )
|
|
|
837 local_qc_features.append(
|
|
|
838 gen_qc_feature(
|
|
|
839 gene.location.start,
|
|
|
840 gene.location.end,
|
|
|
841 "Gene is missing a locus_tag",
|
|
|
842 strand=gene.strand,
|
|
|
843 type_src="gene"
|
|
|
844 )
|
|
|
845 )
|
|
|
846 elif problem == "Different locus tag from associated gene.":
|
|
|
847 gene.qualifiers["locus_tag"] = gene.qualifiers["Name"]
|
|
|
848 cds.qualifiers["locus_tag"] = cds.qualifiers["cpt_gmc_locus"]
|
|
|
849 local_results.append(
|
|
|
850 (gene, cds, "Gene and CDS have differing locus tags")
|
|
|
851 )
|
|
|
852 local_qc_features.append(
|
|
|
853 gen_qc_feature(
|
|
|
854 gene.location.start,
|
|
|
855 gene.location.end,
|
|
|
856 "Gene and CDS have differing locus tags",
|
|
|
857 strand=gene.strand,
|
|
|
858 type_src="gene"
|
|
|
859 )
|
|
|
860 )
|
|
|
861 elif problem == "Missing Locus Tag":
|
|
|
862 # Copy this over
|
|
|
863 gene.qualifiers["locus_tag"] = gene.qualifiers["Name"]
|
|
|
864 # This one is missing
|
|
|
865 cds.qualifiers["locus_tag"] = [""]
|
|
|
866 local_results.append((gene, cds, "CDS is missing a locus_tag"))
|
|
|
867 local_qc_features.append(
|
|
|
868 gen_qc_feature(
|
|
|
869 cds.location.start,
|
|
|
870 cds.location.end,
|
|
|
871 "CDS is missing a locus_tag",
|
|
|
872 strand=cds.strand,
|
|
|
873 type_src="CDS"
|
|
|
874 )
|
|
|
875 )
|
|
|
876 else:
|
|
|
877 log.warn("Cannot handle %s", problem)
|
|
|
878
|
|
|
879 if len(local_results) > 0:
|
|
|
880 bad += 1
|
|
|
881 else:
|
|
|
882 good += 1
|
|
|
883
|
|
|
884 qc_features.extend(local_qc_features)
|
|
|
885 results.extend(local_results)
|
|
|
886 return good, bad, results, qc_features
|
|
|
887
|
|
|
888
|
|
|
889 def missing_tags(record):
|
|
|
890 """Find features without product
|
|
|
891 """
|
|
|
892 results = []
|
|
|
893 good = 0
|
|
|
894 bad = 0
|
|
|
895 qc_features = []
|
|
|
896
|
|
|
897 for gene in coding_genes(record.features):
|
|
|
898 cds = [x for x in genes(gene.sub_features, feature_type="CDS")]
|
|
|
899 if len(cds) == 0:
|
|
|
900 log.warn("Gene missing CDS subfeature %s", get_gff3_id(gene))
|
|
|
901 continue
|
|
|
902
|
|
|
903 cds = cds[0]
|
|
|
904
|
|
|
905 if "product" not in cds.qualifiers:
|
|
|
906 log.info("Missing product tag on %s", get_gff3_id(gene))
|
|
|
907 qc_features.append(
|
|
|
908 gen_qc_feature(
|
|
|
909 cds.location.start,
|
|
|
910 cds.location.end,
|
|
|
911 "Missing product tag",
|
|
|
912 strand=cds.strand,
|
|
|
913 type_src="CDS"
|
|
|
914 )
|
|
|
915 )
|
|
|
916 results.append(cds)
|
|
|
917 bad += 1
|
|
|
918 else:
|
|
|
919 good += 1
|
|
|
920
|
|
|
921 return good, bad, results, qc_features
|
|
|
922
|
|
|
923
|
|
|
924 def evaluate_and_report(
|
|
|
925 annotations,
|
|
|
926 genome,
|
|
|
927 gff3=None,
|
|
|
928 tbl=None,
|
|
|
929 sd_min=5,
|
|
|
930 sd_max=15,
|
|
|
931 min_gene_length=30,
|
|
|
932 excessive_gap_dist=50,
|
|
|
933 excessive_gap_divergent_dist=200,
|
|
|
934 excessive_overlap_dist=25,
|
|
|
935 excessive_overlap_divergent_dist=50,
|
|
|
936 reportTemplateName="phage_annotation_validator.html",
|
|
|
937 ):
|
|
|
938 """
|
|
|
939 Generate our HTML evaluation of the genome
|
|
|
940 """
|
|
|
941 # Get features from GFF file
|
|
|
942 seq_dict = SeqIO.to_dict(SeqIO.parse(genome, "fasta"))
|
|
|
943 # Get the first GFF3 record
|
|
|
944 # TODO: support multiple GFF3 files.
|
|
|
945 mostFeat = 0
|
|
|
946 for rec in list(gffParse(annotations, base_dict=seq_dict)):
|
|
|
947 if len(rec.features) > mostFeat:
|
|
|
948 mostFeat = len(rec.features)
|
|
|
949 record = rec
|
|
|
950
|
|
|
951 gff3_qc_record = SeqRecord(record.id, id=record.id)
|
|
|
952 gff3_qc_record.features = []
|
|
|
953 gff3_qc_features = []
|
|
|
954
|
|
|
955 log.info("Locating missing RBSs")
|
|
|
956 # mb_any = "did they annotate ANY rbss? if so, take off from score."
|
|
|
957 mb_good, mb_bad, mb_results, mb_annotations, mb_any = missing_rbs(
|
|
|
958 record, lookahead_min=sd_min, lookahead_max=sd_max
|
|
|
959 )
|
|
|
960 gff3_qc_features += mb_annotations
|
|
|
961
|
|
|
962 log.info("Locating excessive gaps")
|
|
|
963 eg_good, eg_bad, eg_results, eg_annotations = excessive_gap(
|
|
|
964 record,
|
|
|
965 excess=excessive_gap_dist,
|
|
|
966 excess_divergent=excessive_gap_divergent_dist,
|
|
|
967 min_gene=min_gene_length,
|
|
|
968 slop=excessive_overlap_dist,
|
|
|
969 lookahead_min=sd_min,
|
|
|
970 lookahead_max=sd_max,
|
|
|
971 )
|
|
|
972 gff3_qc_features += eg_annotations
|
|
|
973
|
|
|
974 log.info("Locating excessive overlaps")
|
|
|
975 eo_good, eo_bad, eo_results, eo_annotations = excessive_overlap(
|
|
|
976 record,
|
|
|
977 excess=excessive_overlap_dist,
|
|
|
978 excess_divergent=excessive_overlap_divergent_dist,
|
|
|
979 )
|
|
|
980 gff3_qc_features += eo_annotations
|
|
|
981
|
|
|
982 log.info("Locating morons")
|
|
|
983 mo_good, mo_bad, mo_results, mo_annotations = find_morons(record)
|
|
|
984 gff3_qc_features += mo_annotations
|
|
|
985
|
|
|
986 log.info("Locating missing tags")
|
|
|
987 mt_good, mt_bad, mt_results, mt_annotations = missing_tags(record)
|
|
|
988 gff3_qc_features += mt_annotations
|
|
|
989
|
|
|
990 log.info("Locating missing gene features")
|
|
|
991 mg_good, mg_bad, mg_results, mg_annotations = missing_genes(record)
|
|
|
992 gff3_qc_features += mg_annotations
|
|
|
993
|
|
|
994 log.info("Determining coding density")
|
|
|
995 cd, cd_real = coding_density(record)
|
|
|
996
|
|
|
997 log.info("Locating weird starts")
|
|
|
998 ws_good, ws_bad, ws_results, ws_annotations, ws_overall = weird_starts(record)
|
|
|
999 gff3_qc_features += ws_annotations
|
|
|
1000
|
|
|
1001 log.info("Locating bad gene models")
|
|
|
1002 gm_good, gm_bad, gm_results, gm_annotations = bad_gene_model(record)
|
|
|
1003 if gm_good + gm_bad == 0:
|
|
|
1004 gm_bad = 1
|
|
|
1005
|
|
|
1006 log.info("Locating more bad gene models")
|
|
|
1007 gmc_good, gmc_bad, gmc_results, gmc_annotations = gene_model_correction_issues(
|
|
|
1008 record
|
|
|
1009 )
|
|
|
1010 if gmc_good + gmc_bad == 0:
|
|
|
1011 gmc_bad = 1
|
|
|
1012
|
|
|
1013 good_scores = [eg_good, eo_good, mt_good, ws_good, gm_good, gmc_good]
|
|
|
1014 bad_scores = [eg_bad, eo_bad, mt_bad, ws_bad, gm_bad, gmc_bad]
|
|
|
1015
|
|
|
1016 # Only if they tried to annotate RBSs do we consider them.
|
|
|
1017 if mb_any:
|
|
|
1018 good_scores.append(mb_good)
|
|
|
1019 bad_scores.append(mb_bad)
|
|
|
1020 subscores = []
|
|
|
1021
|
|
|
1022 for (g, b) in zip(good_scores, bad_scores):
|
|
|
1023 if g + b == 0:
|
|
|
1024 s = 0
|
|
|
1025 else:
|
|
|
1026 s = int(100 * float(g) / (float(b) + float(g)))
|
|
|
1027 subscores.append(s)
|
|
|
1028 subscores.append(cd)
|
|
|
1029
|
|
|
1030 score = int(float(sum(subscores)) / float(len(subscores)))
|
|
|
1031
|
|
|
1032 # This is data that will go into our HTML template
|
|
|
1033 kwargs = {
|
|
|
1034 "upstream_min": sd_min,
|
|
|
1035 "upstream_max": sd_max,
|
|
|
1036 "record_name": record.id,
|
|
|
1037 "record_nice_name": nice_name(record),
|
|
|
1038 "params": {
|
|
|
1039 "sd_min": sd_min,
|
|
|
1040 "sd_max": sd_max,
|
|
|
1041 "min_gene_length": min_gene_length,
|
|
|
1042 "excessive_gap_dist": excessive_gap_dist,
|
|
|
1043 "excessive_gap_divergent_dist": excessive_gap_divergent_dist,
|
|
|
1044 "excessive_overlap_dist": excessive_overlap_dist,
|
|
|
1045 "excessive_overlap_divergent_dist": excessive_overlap_divergent_dist,
|
|
|
1046 },
|
|
|
1047 "score": score,
|
|
|
1048 "encouragement": get_encouragement(score),
|
|
|
1049 "genome_overview": genome_overview(record),
|
|
|
1050 "rbss_annotated": mb_any,
|
|
|
1051 "missing_rbs": mb_results,
|
|
|
1052 "missing_rbs_good": mb_good,
|
|
|
1053 "missing_rbs_bad": mb_bad,
|
|
|
1054 "missing_rbs_score": 0
|
|
|
1055 if mb_good + mb_bad == 0
|
|
|
1056 else (100 * mb_good / (mb_good + mb_bad)),
|
|
|
1057 "excessive_gap": eg_results,
|
|
|
1058 "excessive_gap_good": eg_good,
|
|
|
1059 "excessive_gap_bad": eg_bad,
|
|
|
1060 "excessive_gap_score": 0
|
|
|
1061 if eo_good + eo_bad == 0
|
|
|
1062 else (100 * eo_good / (eo_good + eo_bad)),
|
|
|
1063 "excessive_overlap": eo_results,
|
|
|
1064 "excessive_overlap_good": eo_good,
|
|
|
1065 "excessive_overlap_bad": eo_bad,
|
|
|
1066 "excessive_overlap_score": 0
|
|
|
1067 if eo_good + eo_bad == 0
|
|
|
1068 else (100 * eo_good / (eo_good + eo_bad)),
|
|
|
1069 "morons": mo_results,
|
|
|
1070 "morons_good": mo_good,
|
|
|
1071 "morons_bad": mo_bad,
|
|
|
1072 "morons_score": 0
|
|
|
1073 if mo_good + mo_bad == 0
|
|
|
1074 else (100 * mo_good / (mo_good + mo_bad)),
|
|
|
1075 "missing_tags": mt_results,
|
|
|
1076 "missing_tags_good": mt_good,
|
|
|
1077 "missing_tags_bad": mt_bad,
|
|
|
1078 "missing_tags_score": 0
|
|
|
1079 if mt_good + mt_bad == 0
|
|
|
1080 else (100 * mt_good / (mt_good + mt_bad)),
|
|
|
1081 "missing_genes": mg_results,
|
|
|
1082 "missing_genes_good": mg_good,
|
|
|
1083 "missing_genes_bad": mg_bad,
|
|
|
1084 "missing_genes_score": 0
|
|
|
1085 if mg_good + mg_bad == 0
|
|
|
1086 else (100 * mg_good / (mg_good + mg_bad)),
|
|
|
1087 "weird_starts": ws_results,
|
|
|
1088 "weird_starts_good": ws_good,
|
|
|
1089 "weird_starts_bad": ws_bad,
|
|
|
1090 "weird_starts_overall": ws_overall,
|
|
|
1091 "weird_starts_overall_sorted_keys": sorted(
|
|
|
1092 ws_overall, reverse=True, key=lambda x: ws_overall[x]
|
|
|
1093 ),
|
|
|
1094 "weird_starts_score": 0
|
|
|
1095 if ws_good + ws_bad == 0
|
|
|
1096 else (100 * ws_good / (ws_good + ws_bad)),
|
|
|
1097 "gene_model": gm_results,
|
|
|
1098 "gene_model_good": gm_good,
|
|
|
1099 "gene_model_bad": gm_bad,
|
|
|
1100 "gene_model_score": 0
|
|
|
1101 if gm_good + gm_bad == 0
|
|
|
1102 else (100 * gm_good / (gm_good + gm_bad)),
|
|
|
1103 "gene_model_correction": gmc_results,
|
|
|
1104 "gene_model_correction_good": gmc_good,
|
|
|
1105 "gene_model_correction_bad": gmc_bad,
|
|
|
1106 "gene_model_correction_score": 0
|
|
|
1107 if gmc_good + gmc_bad == 0
|
|
|
1108 else (100 * gmc_good / (gmc_good + gmc_bad)),
|
|
|
1109 "coding_density": cd,
|
|
|
1110 "coding_density_exact": exact_coding_density(record),
|
|
|
1111 "coding_density_real": cd_real,
|
|
|
1112 "coding_density_score": cd,
|
|
|
1113 }
|
|
|
1114
|
|
|
1115 with open(tbl, "w") as handle:
|
|
|
1116 kw_subset = {}
|
|
|
1117 for key in kwargs:
|
|
|
1118 if (
|
|
|
1119 key in ("score", "record_name")
|
|
|
1120 or "_good" in key
|
|
|
1121 or "_bad" in key
|
|
|
1122 or "_overall" in key
|
|
|
1123 ):
|
|
|
1124 kw_subset[key] = kwargs[key]
|
|
|
1125 json.dump(kw_subset, handle)
|
|
|
1126
|
|
|
1127 with open(gff3, "w") as handle:
|
|
|
1128 gff3_qc_record.features = gff3_qc_features
|
|
|
1129 gff3_qc_record.annotations = {}
|
|
|
1130 gffWrite([gff3_qc_record], handle)
|
|
|
1131
|
|
|
1132 def nice_strand(direction):
|
|
|
1133 # It is somehow possible for whole gffSeqFeature objects to end up in here, apparently at the gene level
|
|
|
1134 if "SeqFeature" in str(type(direction)):
|
|
|
1135 direction = direction.location.strand
|
|
|
1136 if direction > 0:
|
|
|
1137 return "→"#.decode("utf-8")
|
|
|
1138 else:
|
|
|
1139 return "←"#.decode("utf-8")
|
|
|
1140
|
|
|
1141 def nice_strand_tex(direction):
|
|
|
1142 if "SeqFeature" in str(type(direction)):
|
|
|
1143 direction = direction.location.strand
|
|
|
1144 if direction > 0:
|
|
|
1145 return "$\\rightarrow$"
|
|
|
1146 else:
|
|
|
1147 return "$\\leftarrow$"
|
|
|
1148
|
|
|
1149 def texify(data):
|
|
|
1150 return data.replace("_", "\\_").replace("$", "\\$")
|
|
|
1151
|
|
|
1152 def length(data):
|
|
|
1153 return len(data)
|
|
|
1154
|
|
|
1155 def my_encode(data):
|
|
|
1156 return str(data)#.encode("utf-8")
|
|
|
1157
|
|
|
1158 def my_decode(data):
|
|
|
1159 # For production
|
|
|
1160 return str(data)#.decode("utf-8")
|
|
|
1161 # For local testing. No, I do not understand.
|
|
|
1162 return str(data)#.encode("utf-8")).decode("utf-8")
|
|
|
1163
|
|
|
1164 env = Environment(
|
|
|
1165 loader=FileSystemLoader(SCRIPT_PATH), trim_blocks=True, lstrip_blocks=True
|
|
|
1166 )
|
|
|
1167 env.filters.update(
|
|
|
1168 {
|
|
|
1169 "nice_id": get_gff3_id,
|
|
|
1170 "nice_strand": nice_strand,
|
|
|
1171 "nice_strand_tex": nice_strand_tex,
|
|
|
1172 "texify": texify,
|
|
|
1173 "length": length,
|
|
|
1174 "encode": my_encode,
|
|
|
1175 "decode": my_decode,
|
|
|
1176 }
|
|
|
1177 )
|
|
|
1178 tpl = env.get_template(reportTemplateName)
|
|
|
1179 return tpl.render(**kwargs)#.encode("utf-8")
|
|
|
1180
|
|
|
1181
|
|
|
1182 if __name__ == "__main__":
|
|
|
1183 parser = argparse.ArgumentParser(
|
|
|
1184 description="rebase gff3 features against parent locations", epilog=""
|
|
|
1185 )
|
|
|
1186 parser.add_argument(
|
|
|
1187 "annotations", type=argparse.FileType("r"), help="Parent GFF3 annotations"
|
|
|
1188 )
|
|
|
1189 parser.add_argument("genome", type=argparse.FileType("r"), help="Genome Sequence")
|
|
|
1190 parser.add_argument(
|
|
|
1191 "--gff3", type=str, help="GFF3 Annotations", default="qc_annotations.gff3"
|
|
|
1192 )
|
|
|
1193 parser.add_argument(
|
|
|
1194 "--tbl",
|
|
|
1195 type=str,
|
|
|
1196 help="Table for noninteractive parsing",
|
|
|
1197 default="qc_results.json",
|
|
|
1198 )
|
|
|
1199
|
|
|
1200 parser.add_argument(
|
|
|
1201 "--sd_min",
|
|
|
1202 type=int,
|
|
|
1203 help="Minimum distance from gene start for an SD to be",
|
|
|
1204 default=5,
|
|
|
1205 )
|
|
|
1206 parser.add_argument(
|
|
|
1207 "--sd_max",
|
|
|
1208 type=int,
|
|
|
1209 help="Maximum distance from gene start for an SD to be",
|
|
|
1210 default=15,
|
|
|
1211 )
|
|
|
1212
|
|
|
1213 parser.add_argument(
|
|
|
1214 "--min_gene_length",
|
|
|
1215 type=int,
|
|
|
1216 help="Minimum length for a putative gene call (AAs)",
|
|
|
1217 default=30,
|
|
|
1218 )
|
|
|
1219
|
|
|
1220 parser.add_argument(
|
|
|
1221 "--excessive_overlap_dist",
|
|
|
1222 type=int,
|
|
|
1223 help="Excessive overlap for genes in same direction",
|
|
|
1224 default=25,
|
|
|
1225 )
|
|
|
1226 parser.add_argument(
|
|
|
1227 "--excessive_overlap_divergent_dist",
|
|
|
1228 type=int,
|
|
|
1229 help="Excessive overlap for genes in diff directions",
|
|
|
1230 default=50,
|
|
|
1231 )
|
|
|
1232
|
|
|
1233 parser.add_argument(
|
|
|
1234 "--excessive_gap_dist",
|
|
|
1235 type=int,
|
|
|
1236 help="Maximum distance between two genes",
|
|
|
1237 default=40,
|
|
|
1238 )
|
|
|
1239 parser.add_argument(
|
|
|
1240 "--excessive_gap_divergent_dist",
|
|
|
1241 type=int,
|
|
|
1242 help="Maximum distance between two divergent genes",
|
|
|
1243 default=200,
|
|
|
1244 )
|
|
|
1245
|
|
|
1246 parser.add_argument(
|
|
|
1247 "--reportTemplateName",
|
|
|
1248 help="Report template file name",
|
|
|
1249 default="phageqc_report_full.html",
|
|
|
1250 )
|
|
|
1251
|
|
|
1252 args = parser.parse_args()
|
|
|
1253
|
|
|
1254 sys.stdout.write(evaluate_and_report(**vars(args)))
|
