v[A-Z]_[A-Z][a-z]{2}_.*)")
+
+
+def phage_name_parser(name):
+ host = None
+ phage = None
+ name = name.replace(", complete genome.", "")
+ name = name.replace(", complete genome", "")
+
+ m = BACTERIOPHAGE_IN_MIDDLE.match(name)
+ if m:
+ host = m.group("host")
+ phage = m.group("phage")
+ return (host, phage)
+
+ m = PHAGE_IN_MIDDLE.match(name)
+ if m:
+ host = m.group("host")
+ phage = m.group("phage")
+ return (host, phage)
+
+ m = STARTS_WITH_PHAGE.match(name)
+ if m:
+ phage = m.group("phage")
+ return (host, phage)
+
+ m = NEW_STYLE_NAMES.match(name)
+ if m:
+ phage = m.group("phage")
+ return (host, phage)
+
+ return (host, phage)
+
+
+class OrfFinder(object):
+ def __init__(self, table, ftype, ends, min_len, strand):
+ self.table = table
+ self.table_obj = CodonTable.ambiguous_generic_by_id[table]
+ self.ends = ends
+ self.ftype = ftype
+ self.min_len = min_len
+ self.starts = sorted(self.table_obj.start_codons)
+ self.stops = sorted(self.table_obj.stop_codons)
+ self.re_starts = re.compile("|".join(self.starts))
+ self.re_stops = re.compile("|".join(self.stops))
+ self.strand = strand
+
+ def locate(self, fasta_file, out_nuc, out_prot, out_bed, out_gff3):
+ seq_format = "fasta"
+ log.debug("Genetic code table %i" % self.table)
+ log.debug("Minimum length %i aa" % self.min_len)
+
+ out_count = 0
+
+ out_gff3.write("##gff-version 3\n")
+
+ for idx, record in enumerate(SeqIO.parse(fasta_file, seq_format)):
+ for i, (f_start, f_end, f_strand, n, t) in enumerate(
+ self.get_all_peptides(str(record.seq).upper())
+ ):
+ out_count += 1
+
+ descr = "length %i aa, %i bp, from %s..%s[%s] of %s" % (
+ len(t),
+ len(n),
+ f_start,
+ f_end,
+ f_strand,
+ record.description,
+ )
+ fid = record.id + "|%s%i" % (self.ftype, i + 1)
+
+ r = SeqRecord(Seq(n), id=fid, name="", description=descr)
+ t = SeqRecord(Seq(t), id=fid, name="", description=descr)
+
+ SeqIO.write(r, out_nuc, "fasta")
+ SeqIO.write(t, out_prot, "fasta")
+
+ nice_strand = "+" if f_strand == +1 else "-"
+
+ out_bed.write(
+ "\t".join(
+ map(str, [record.id, f_start, f_end, fid, 0, nice_strand])
+ )
+ + "\n"
+ )
+
+ out_gff3.write(
+ "\t".join(
+ map(
+ str,
+ [
+ record.id,
+ "getOrfsOrCds",
+ "CDS",
+ f_start + 1,
+ f_end,
+ ".",
+ nice_strand,
+ 0,
+ "ID=%s.%s.%s" % (self.ftype, idx, i + 1),
+ ],
+ )
+ )
+ + "\n"
+ )
+ log.info("Found %i %ss", out_count, self.ftype)
+
+ def start_chop_and_trans(self, s, strict=True):
+ """Returns offset, trimmed nuc, protein."""
+ if strict:
+ assert s[-3:] in self.stops, s
+ assert len(s) % 3 == 0
+ for match in self.re_starts.finditer(s, overlapped=True):
+ # Must check the start is in frame
+ start = match.start()
+ if start % 3 == 0:
+ n = s[start:]
+ assert len(n) % 3 == 0, "%s is len %i" % (n, len(n))
+ if strict:
+ t = translate(n, self.table)
+ else:
+ # Use when missing stop codon,
+ t = "M" + translate(n[3:], self.table, to_stop=True)
+ yield start, n, t # Edited by CPT to be a generator
+
+ def break_up_frame(self, s):
+ """Returns offset, nuc, protein."""
+ start = 0
+ for match in self.re_stops.finditer(s, overlapped=True):
+ index = match.start() + 3
+ if index % 3 != 0:
+ continue
+ n = s[start:index]
+ for (offset, n, t) in self.start_chop_and_trans(n):
+ if n and len(t) >= self.min_len:
+ yield start + offset, n, t
+ start = index
+
+ def putative_genes_in_sequence(self, nuc_seq):
+ """Returns start, end, strand, nucleotides, protein.
+ Co-ordinates are Python style zero-based.
+ """
+ nuc_seq = nuc_seq.upper()
+ # TODO - Refactor to use a generator function (in start order)
+ # rather than making a list and sorting?
+ answer = []
+ full_len = len(nuc_seq)
+
+ for frame in range(0, 3):
+ for offset, n, t in self.break_up_frame(nuc_seq[frame:]):
+ start = frame + offset # zero based
+ answer.append((start, start + len(n), +1, n, t))
+
+ rc = reverse_complement(nuc_seq)
+ for frame in range(0, 3):
+ for offset, n, t in self.break_up_frame(rc[frame:]):
+ start = full_len - frame - offset # zero based
+ answer.append((start, start - len(n), -1, n, t))
+ answer.sort()
+ return answer
+
+ def get_all_peptides(self, nuc_seq):
+ """Returns start, end, strand, nucleotides, protein.
+
+ Co-ordinates are Python style zero-based.
+ """
+ # Refactored into generator by CPT
+ full_len = len(nuc_seq)
+ if self.strand != "reverse":
+ for frame in range(0, 3):
+ for offset, n, t in self.break_up_frame(nuc_seq[frame:]):
+ start = frame + offset # zero based
+ yield (start, start + len(n), +1, n, t)
+ if self.strand != "forward":
+ rc = reverse_complement(nuc_seq)
+ for frame in range(0, 3):
+ for offset, n, t in self.break_up_frame(rc[frame:]):
+ start = full_len - frame - offset # zero based
+ yield (start - len(n), start, -1, n, t)
+
+
+class MGAFinder(object):
+ def __init__(self, table, ftype, ends, min_len):
+ self.table = table
+ self.table_obj = CodonTable.ambiguous_generic_by_id[table]
+ self.ends = ends
+ self.ftype = ftype
+ self.min_len = min_len
+ self.starts = sorted(self.table_obj.start_codons)
+ self.stops = sorted(self.table_obj.stop_codons)
+ self.re_starts = re.compile("|".join(self.starts))
+ self.re_stops = re.compile("|".join(self.stops))
+
+ def locate(self, fasta_file, out_nuc, out_prot, out_bed, out_gff3):
+ seq_format = "fasta"
+ log.debug("Genetic code table %i" % self.table)
+ log.debug("Minimum length %i aa" % self.min_len)
+
+ out_count = 0
+
+ out_gff3.write("##gff-version 3\n")
+
+ for idx, record in enumerate(SeqIO.parse(fasta_file, seq_format)):
+ for i, (f_start, f_end, f_strand, n, t) in enumerate(
+ self.get_all_peptides(str(record.seq).upper())
+ ):
+ out_count += 1
+
+ descr = "length %i aa, %i bp, from %s..%s[%s] of %s" % (
+ len(t),
+ len(n),
+ f_start,
+ f_end,
+ f_strand,
+ record.description,
+ )
+ fid = record.id + "|%s%i" % (self.ftype, i + 1)
+
+ r = SeqRecord(Seq(n), id=fid, name="", description=descr)
+ t = SeqRecord(Seq(t), id=fid, name="", description=descr)
+
+ SeqIO.write(r, out_nuc, "fasta")
+ SeqIO.write(t, out_prot, "fasta")
+
+ nice_strand = "+" if f_strand == +1 else "-"
+
+ out_bed.write(
+ "\t".join(
+ map(str, [record.id, f_start, f_end, fid, 0, nice_strand])
+ )
+ + "\n"
+ )
+
+ out_gff3.write(
+ "\t".join(
+ map(
+ str,
+ [
+ record.id,
+ "getOrfsOrCds",
+ "CDS",
+ f_start + 1,
+ f_end,
+ ".",
+ nice_strand,
+ 0,
+ "ID=%s.%s.%s" % (self.ftype, idx, i + 1),
+ ],
+ )
+ )
+ + "\n"
+ )
+ log.info("Found %i %ss", out_count, self.ftype)
+
+ def start_chop_and_trans(self, s, strict=True):
+ """Returns offset, trimmed nuc, protein."""
+ if strict:
+ assert s[-3:] in self.stops, s
+ assert len(s) % 3 == 0
+ for match in self.re_starts.finditer(s, overlapped=True):
+ # Must check the start is in frame
+ start = match.start()
+ if start % 3 == 0:
+ n = s[start:]
+ assert len(n) % 3 == 0, "%s is len %i" % (n, len(n))
+ if strict:
+ t = translate(n, self.table)
+ else:
+ # Use when missing stop codon,
+ t = "M" + translate(n[3:], self.table, to_stop=True)
+ yield start, n, t
+
+ def break_up_frame(self, s):
+ """Returns offset, nuc, protein."""
+ start = 0
+ for match in self.re_stops.finditer(s, overlapped=True):
+ index = match.start() + 3
+ if index % 3 != 0:
+ continue
+ n = s[start:index]
+ for (offset, n, t) in self.start_chop_and_trans(n):
+ if n and len(t) >= self.min_len:
+ yield start + offset, n, t
+ start = index
+
+ def putative_genes_in_sequence(self, nuc_seq):
+ """Returns start, end, strand, nucleotides, protein.
+ Co-ordinates are Python style zero-based.
+ """
+ nuc_seq = nuc_seq.upper()
+ # TODO - Refactor to use a generator function (in start order)
+ # rather than making a list and sorting?
+ answer = []
+ full_len = len(nuc_seq)
+
+ for frame in range(0, 3):
+ for offset, n, t in self.break_up_frame(nuc_seq[frame:]):
+ start = frame + offset # zero based
+ answer.append((start, start + len(n), +1, n, t))
+
+ rc = reverse_complement(nuc_seq)
+ for frame in range(0, 3):
+ for offset, n, t in self.break_up_frame(rc[frame:]):
+ start = full_len - frame - offset # zero based
+ answer.append((start, start - len(n), -1, n, t))
+ answer.sort()
+ return answer
+
+ def get_all_peptides(self, nuc_seq):
+ """Returns start, end, strand, nucleotides, protein.
+
+ Co-ordinates are Python style zero-based.
+ """
+ # Refactored into generator by CPT
+
+ full_len = len(nuc_seq)
+ for frame in range(0, 3):
+ for offset, n, t in self.break_up_frame(nuc_seq[frame:]):
+ start = frame + offset # zero based
+ yield (start, start + len(n), +1, n, t)
+ rc = reverse_complement(nuc_seq)
+ for frame in range(0, 3):
+ for offset, n, t in self.break_up_frame(rc[frame:]):
+ start = full_len - frame - offset # zero based
+ yield (start - len(n), start, -1, n, t)
diff -r 000000000000 -r c3140b08d703 cpt_phageqc_annotation/gff3.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/cpt_phageqc_annotation/gff3.py Fri Jun 17 13:00:50 2022 +0000
@@ -0,0 +1,346 @@
+import copy
+import logging
+
+log = logging.getLogger()
+log.setLevel(logging.WARN)
+
+
+def feature_lambda(
+ feature_list,
+ test,
+ test_kwargs,
+ subfeatures=True,
+ parent=None,
+ invert=False,
+ recurse=True,
+):
+ """Recursively search through features, testing each with a test function, yielding matches.
+
+ GFF3 is a hierachical data structure, so we need to be able to recursively
+ search through features. E.g. if you're looking for a feature with
+ ID='bob.42', you can't just do a simple list comprehension with a test
+ case. You don't know how deeply burried bob.42 will be in the feature tree. This is where feature_lambda steps in.
+
+ :type feature_list: list
+ :param feature_list: an iterable of features
+
+ :type test: function reference
+ :param test: a closure with the method signature (feature, **kwargs) where
+ the kwargs are those passed in the next argument. This
+ function should return True or False, True if the feature is
+ to be yielded as part of the main feature_lambda function, or
+ False if it is to be ignored. This function CAN mutate the
+ features passed to it (think "apply").
+
+ :type test_kwargs: dictionary
+ :param test_kwargs: kwargs to pass to your closure when it is called.
+
+ :type subfeatures: boolean
+ :param subfeatures: when a feature is matched, should just that feature be
+ yielded to the caller, or should the entire sub_feature
+ tree for that feature be included? subfeatures=True is
+ useful in cases such as searching for a gene feature,
+ and wanting to know what RBS/Shine_Dalgarno_sequences
+ are in the sub_feature tree (which can be accomplished
+ with two feature_lambda calls). subfeatures=False is
+ useful in cases when you want to process (and possibly
+ return) the entire feature tree, such as applying a
+ qualifier to every single feature.
+
+ :type invert: boolean
+ :param invert: Negate/invert the result of the filter.
+
+ :rtype: yielded list
+ :return: Yields a list of matching features.
+ """
+ # Either the top level set of [features] or the subfeature attribute
+ for feature in feature_list:
+ feature._parent = parent
+ if not parent:
+ # Set to self so we cannot go above root.
+ feature._parent = feature
+ test_result = test(feature, **test_kwargs)
+ # if (not invert and test_result) or (invert and not test_result):
+ if invert ^ test_result:
+ if not subfeatures:
+ feature_copy = copy.deepcopy(feature)
+ feature_copy.sub_features = list()
+ yield feature_copy
+ else:
+ yield feature
+
+ if recurse and hasattr(feature, "sub_features"):
+ for x in feature_lambda(
+ feature.sub_features,
+ test,
+ test_kwargs,
+ subfeatures=subfeatures,
+ parent=feature,
+ invert=invert,
+ recurse=recurse,
+ ):
+ yield x
+
+
+def fetchParent(feature):
+ if not hasattr(feature, "_parent") or feature._parent is None:
+ return feature
+ else:
+ return fetchParent(feature._parent)
+
+
+def feature_test_true(feature, **kwargs):
+ return True
+
+
+def feature_test_type(feature, **kwargs):
+ if "type" in kwargs:
+ return str(feature.type).upper() == str(kwargs["type"]).upper()
+ elif "types" in kwargs:
+ for x in kwargs["types"]:
+ if str(feature.type).upper() == str(x).upper():
+ return True
+ return False
+ raise Exception("Incorrect feature_test_type call, need type or types")
+
+
+def feature_test_qual_value(feature, **kwargs):
+ """Test qualifier values.
+
+ For every feature, check that at least one value in
+ feature.quailfiers(kwargs['qualifier']) is in kwargs['attribute_list']
+ """
+ if isinstance(kwargs["qualifier"], list):
+ for qualifier in kwargs["qualifier"]:
+ for attribute_value in feature.qualifiers.get(qualifier, []):
+ if attribute_value in kwargs["attribute_list"]:
+ return True
+ else:
+ for attribute_value in feature.qualifiers.get(kwargs["qualifier"], []):
+ if attribute_value in kwargs["attribute_list"]:
+ return True
+ return False
+
+
+def feature_test_location(feature, **kwargs):
+ if "strand" in kwargs:
+ if feature.location.strand != kwargs["strand"]:
+ return False
+
+ return feature.location.start <= kwargs["loc"] <= feature.location.end
+
+
+def feature_test_quals(feature, **kwargs):
+ """
+ Example::
+
+ a = Feature(qualifiers={'Note': ['Some notes', 'Aasdf']})
+
+ # Check if a contains a Note
+ feature_test_quals(a, {'Note': None}) # Returns True
+ feature_test_quals(a, {'Product': None}) # Returns False
+
+ # Check if a contains a note with specific value
+ feature_test_quals(a, {'Note': ['ome']}) # Returns True
+
+ # Check if a contains a note with specific value
+ feature_test_quals(a, {'Note': ['other']}) # Returns False
+ """
+ for key in kwargs:
+ if key not in feature.qualifiers:
+ return False
+
+ # Key is present, no value specified
+ if kwargs[key] is None:
+ return True
+
+ # Otherwise there is a key value we're looking for.
+ # so we make a list of matches
+ matches = []
+ # And check all of the feature qualifier valuse
+ for value in feature.qualifiers[key]:
+ # For that kwargs[key] value
+ for x in kwargs[key]:
+ matches.append(x in value)
+
+ # If none matched, then we return false.
+ if not any(matches):
+ return False
+
+ return True
+
+
+def feature_test_contains(feature, **kwargs):
+ if "index" in kwargs:
+ return feature.location.start < kwargs["index"] < feature.location.end
+ elif "range" in kwargs:
+ return (
+ feature.location.start < kwargs["range"]["start"] < feature.location.end
+ and feature.location.start < kwargs["range"]["end"] < feature.location.end
+ )
+ else:
+ raise RuntimeError("Must use index or range keyword")
+
+
+def get_id(feature=None, parent_prefix=None):
+ result = ""
+ if parent_prefix is not None:
+ result += parent_prefix + "|"
+ if "locus_tag" in feature.qualifiers:
+ result += feature.qualifiers["locus_tag"][0]
+ elif "gene" in feature.qualifiers:
+ result += feature.qualifiers["gene"][0]
+ elif "Gene" in feature.qualifiers:
+ result += feature.qualifiers["Gene"][0]
+ elif "product" in feature.qualifiers:
+ result += feature.qualifiers["product"][0]
+ elif "Product" in feature.qualifiers:
+ result += feature.qualifiers["Product"][0]
+ elif "Name" in feature.qualifiers:
+ result += feature.qualifiers["Name"][0]
+ else:
+ return feature.id
+ # Leaving in case bad things happen.
+ # result += '%s_%s_%s_%s' % (
+ # feature.id,
+ # feature.location.start,
+ # feature.location.end,
+ # feature.location.strand
+ # )
+ return result
+
+
+def get_gff3_id(gene):
+ return gene.qualifiers.get("Name", [gene.id])[0]
+
+
+def ensure_location_in_bounds(start=0, end=0, parent_length=0):
+ # This prevents frameshift errors
+ while start < 0:
+ start += 3
+ while end < 0:
+ end += 3
+ while start > parent_length:
+ start -= 3
+ while end > parent_length:
+ end -= 3
+ return (start, end)
+
+
+def coding_genes(feature_list):
+ for x in genes(feature_list):
+ if (
+ len(
+ list(
+ feature_lambda(
+ x.sub_features,
+ feature_test_type,
+ {"type": "CDS"},
+ subfeatures=False,
+ )
+ )
+ )
+ > 0
+ ):
+ yield x
+
+
+def genes(feature_list, feature_type="gene", sort=False):
+ """
+ Simple filter to extract gene features from the feature set.
+ """
+
+ if not sort:
+ for x in feature_lambda(
+ feature_list, feature_test_type, {"type": feature_type}, subfeatures=True
+ ):
+ yield x
+ else:
+ data = list(genes(feature_list, feature_type=feature_type, sort=False))
+ data = sorted(data, key=lambda feature: feature.location.start)
+ for x in data:
+ yield x
+
+
+def wa_unified_product_name(feature):
+ """
+ Try and figure out a name. We gave conflicting instructions, so
+ this isn't as trivial as it should be. Sometimes it will be in
+ 'product' or 'Product', othertimes in 'Name'
+ """
+ # Manually applied tags.
+ protein_product = feature.qualifiers.get(
+ "product", feature.qualifiers.get("Product", [None])
+ )[0]
+
+ # If neither of those are available ...
+ if protein_product is None:
+ # And there's a name...
+ if "Name" in feature.qualifiers:
+ if not is_uuid(feature.qualifiers["Name"][0]):
+ protein_product = feature.qualifiers["Name"][0]
+
+ return protein_product
+
+
+def is_uuid(name):
+ return name.count("-") == 4 and len(name) == 36
+
+
+def get_rbs_from(gene):
+ # Normal RBS annotation types
+ rbs_rbs = list(
+ feature_lambda(
+ gene.sub_features, feature_test_type, {"type": "RBS"}, subfeatures=False
+ )
+ )
+ rbs_sds = list(
+ feature_lambda(
+ gene.sub_features,
+ feature_test_type,
+ {"type": "Shine_Dalgarno_sequence"},
+ subfeatures=False,
+ )
+ )
+ # Fraking apollo
+ apollo_exons = list(
+ feature_lambda(
+ gene.sub_features, feature_test_type, {"type": "exon"}, subfeatures=False
+ )
+ )
+ apollo_exons = [x for x in apollo_exons if len(x) < 10]
+ # These are more NCBI's style
+ regulatory_elements = list(
+ feature_lambda(
+ gene.sub_features,
+ feature_test_type,
+ {"type": "regulatory"},
+ subfeatures=False,
+ )
+ )
+ rbs_regulatory = list(
+ feature_lambda(
+ regulatory_elements,
+ feature_test_quals,
+ {"regulatory_class": ["ribosome_binding_site"]},
+ subfeatures=False,
+ )
+ )
+ # Here's hoping you find just one ;)
+ return rbs_rbs + rbs_sds + rbs_regulatory + apollo_exons
+
+
+def nice_name(record):
+ """
+ get the real name rather than NCBI IDs and so on. If fails, will return record.id
+ """
+ name = record.id
+ likely_parental_contig = list(genes(record.features, feature_type="contig"))
+ if len(likely_parental_contig) == 1:
+ name = likely_parental_contig[0].qualifiers.get("organism", [name])[0]
+ return name
+
+
+def fsort(it):
+ for i in sorted(it, key=lambda x: int(x.location.start)):
+ yield i
diff -r 000000000000 -r c3140b08d703 cpt_phageqc_annotation/macros.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/cpt_phageqc_annotation/macros.xml Fri Jun 17 13:00:50 2022 +0000
@@ -0,0 +1,23 @@
+
+
+
+
+ python
+ biopython
+ cpt_gffparser
+
+
+
+
+
+
+
+
+
+
+ ln -s $genome_fasta genomeref.fa;
+
+
+ genomeref.fa
+
+
diff -r 000000000000 -r c3140b08d703 cpt_phageqc_annotation/phage_annotation_validator.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/cpt_phageqc_annotation/phage_annotation_validator.py Fri Jun 17 13:00:50 2022 +0000
@@ -0,0 +1,1254 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# vim: set fileencoding=utf-8
+import os
+import sys
+import json
+import math
+import numpy
+import argparse
+import itertools
+import logging
+from gff3 import (
+ feature_lambda,
+ coding_genes,
+ genes,
+ get_gff3_id,
+ feature_test_location,
+ get_rbs_from,
+ nice_name,
+)
+from shinefind import NaiveSDCaller
+from CPT_GFFParser import gffParse, gffWrite, gffSeqFeature
+from Bio import SeqIO
+from Bio.SeqRecord import SeqRecord
+from Bio.SeqFeature import SeqFeature, FeatureLocation
+from jinja2 import Environment, FileSystemLoader
+from cpt import MGAFinder
+
+logging.basicConfig(level=logging.DEBUG)
+log = logging.getLogger(name="pav")
+
+# Path to script, required because of Galaxy.
+SCRIPT_PATH = os.path.dirname(os.path.realpath(__file__))
+# Path to the HTML template for the report
+
+ENCOURAGEMENT = (
+ (100, "Perfection itself!"),
+ (90, "Amazing!"),
+ (80, "Not too bad, a few minor things to fix..."),
+ (70, "Some issues to address"),
+ (
+ 50,
+ """Issues detected! Have you heard of the
+ CPT\'s Automated Phage Annotation
+ Pipeline?""",
+ ),
+ (
+ 0,
+ """MAJOR issues detected! Please consider using the
+ CPT\'s Automated Phage Annotation Pipeline""",
+ ),
+)
+
+
+def gen_qc_feature(start, end, message, strand=0, id_src=None, type_src="gene"):
+ kwargs = {"qualifiers": {"note": [message]}}
+ kwargs["type"] = type_src
+ kwargs["strand"] = strand
+ kwargs["phase"]=0
+ kwargs["score"]=0.0
+ kwargs["source"]="feature"
+ if id_src is not None:
+ kwargs["id"] = id_src.id
+ kwargs["qualifiers"]["ID"] = [id_src.id]
+ kwargs["qualifiers"]["Name"] = id_src.qualifiers.get("Name", [])
+
+
+ if end >= start:
+ return gffSeqFeature(FeatureLocation(start, end, strand=strand), **kwargs)
+ else:
+ return gffSeqFeature(FeatureLocation(end, start, strand=strand), **kwargs)
+
+
+def __ensure_location_in_bounds(start=0, end=0, parent_length=0):
+ # This prevents frameshift errors
+ while start < 0:
+ start += 3
+ while end < 0:
+ end += 3
+ while start > parent_length:
+ start -= 3
+ while end > parent_length:
+ end -= 3
+ return (start, end)
+
+
+def missing_rbs(record, lookahead_min=5, lookahead_max=15):
+ """
+ Identify gene features with missing RBSs
+
+ This "looks ahead" 5-15 bases ahead of each gene feature, and checks if
+ there's an RBS feature in those bounds.
+
+ The returned data is a set of genes with the RBS sequence in the __upstream
+ attribute, and a message in the __message attribute.
+ """
+ results = []
+ good = 0
+ bad = 0
+ qc_features = []
+ sd_finder = NaiveSDCaller()
+
+ any_rbss = False
+
+ for gene in coding_genes(record.features):
+ # Check if there are RBSs, TODO: make this recursive. Each feature in
+ # gene.sub_features can also have sub_features.
+ rbss = get_rbs_from(gene)
+ # No RBS found
+ if len(rbss) == 0:
+ # Get the sequence lookahead_min to lookahead_max upstream
+ if gene.strand > 0:
+ start = gene.location.start - lookahead_max
+ end = gene.location.start - lookahead_min
+ else:
+ start = gene.location.end + lookahead_min
+ end = gene.location.end + lookahead_max
+ # We have to ensure the feature is ON the genome, otherwise we may
+ # be trying to access a location outside of the length of the
+ # genome, which would be bad.
+ (start, end) = __ensure_location_in_bounds(
+ start=start, end=end, parent_length=len(record)
+ )
+ # Temporary feature to extract sequence
+ tmp = gffSeqFeature(
+ FeatureLocation(start, end, strand=gene.strand), type="domain"
+ )
+ # Get the sequence
+ seq = str(tmp.extract(record.seq))
+ # Set the default properties
+ gene.__upstream = seq.lower()
+ gene.__message = "No RBS annotated, None found"
+
+ # Try and do an automated shinefind call
+ sds = sd_finder.list_sds(seq)
+ if len(sds) > 0:
+ sd = sds[0]
+ gene.__upstream = sd_finder.highlight_sd(
+ seq.lower(), sd["start"], sd["end"]
+ )
+ gene.__message = "Unannotated but valid RBS"
+
+ qc_features.append(
+ gen_qc_feature(
+ start, end, "Missing RBS", strand=gene.strand, id_src=gene, type_src="gene"
+ )
+ )
+
+ bad += 1
+ results.append(gene)
+ results[-1].location = FeatureLocation(results[-1].location.start + 1, results[-1].location.end, results[-1].location.strand)
+ else:
+ if len(rbss) > 1:
+ log.warn("%s RBSs found for gene %s", rbss[0].id, get_gff3_id(gene))
+ any_rbss = True
+ # get first RBS/CDS
+ cds = list(genes(gene.sub_features, feature_type="CDS"))[0]
+ rbs = rbss[0]
+
+ # Get the distance between the two
+ if gene.strand > 0:
+ distance = cds.location.start - rbs.location.end
+ else:
+ distance = rbs.location.start - cds.location.end
+
+ # If the RBS is too far away, annotate that
+ if distance > lookahead_max:
+ gene.__message = "RBS too far away (%s nt)" % distance
+
+ qc_features.append(
+ gen_qc_feature(
+ rbs.location.start,
+ rbs.location.end,
+ gene.__message,
+ strand=gene.strand,
+ id_src=gene,
+ type_src="gene"
+ )
+ )
+
+ bad += 1
+ results.append(gene)
+ results[-1].location = FeatureLocation(results[-1].location.start + 1, results[-1].location.end, results[-1].location.strand)
+ else:
+ good += 1
+
+ return good, bad, results, qc_features, any_rbss
+
+
+# modified from get_orfs_or_cdss.py
+# -----------------------------------------------------------
+
+
+def require_sd(data, record, chrom_start, sd_min, sd_max):
+ sd_finder = NaiveSDCaller()
+ for putative_gene in data:
+ if putative_gene[2] > 0: # strand
+ start = chrom_start + putative_gene[0] - sd_max
+ end = chrom_start + putative_gene[0] - sd_min
+ else:
+ start = chrom_start + putative_gene[1] + sd_min
+ end = chrom_start + putative_gene[1] + sd_max
+
+ (start, end) = __ensure_location_in_bounds(
+ start=start, end=end, parent_length=len(record)
+ )
+ tmp = gffSeqFeature(
+ FeatureLocation(start, end, strand=putative_gene[2]), type="domain"
+ )
+ # Get the sequence
+ seq = str(tmp.extract(record.seq))
+ sds = sd_finder.list_sds(seq)
+ if len(sds) > 0:
+ yield putative_gene + (start, end)
+
+
+def excessive_gap(
+ record,
+ excess=50,
+ excess_divergent=200,
+ min_gene=30,
+ slop=30,
+ lookahead_min=5,
+ lookahead_max=15,
+):
+ """
+ Identify excessive gaps between gene features.
+
+ Default "excessive" gap size is 10, but that should likely be larger.
+ """
+ results = []
+ good = 0
+ bad = 0
+
+ contiguous_regions = []
+
+ sorted_genes = sorted(
+ genes(record.features), key=lambda feature: feature.location.start
+ )
+ if len(sorted_genes) == 0:
+ log.warn("NO GENES FOUND")
+ return good, bad, results, []
+
+ current_gene = None
+ for gene in sorted_genes:
+ # If the gene's start is contiguous to the "current_gene", then we
+ # extend current_gene
+ for cds in genes(gene.sub_features, feature_type="CDS"):
+ if current_gene is None:
+ current_gene = [int(cds.location.start), int(cds.location.end)]
+
+ if cds.location.start <= current_gene[1] + excess:
+ # Don't want to decrease size
+ if int(cds.location.end) >= current_gene[1]:
+ current_gene[1] = int(cds.location.end)
+ else:
+ # If it's discontiguous, we append the region and clear.
+ contiguous_regions.append(current_gene)
+ current_gene = [int(cds.location.start), int(cds.location.end)]
+
+ # This generally expected that annotations would NOT continue unto the end
+ # of the genome, however that's a bug, and we can make it here with an
+ # empty contiguous_regions list
+ contiguous_regions.append(current_gene)
+
+ for i in range(len(contiguous_regions) + 1):
+ if i == 0:
+ a = (1, 1)
+ b = contiguous_regions[i]
+ elif i >= len(contiguous_regions):
+ a = contiguous_regions[i - 1]
+ b = (len(record.seq), None)
+ else:
+ a = contiguous_regions[i - 1]
+ b = contiguous_regions[i]
+
+ gap_size = abs(b[0] - a[1])
+
+ if gap_size > min(excess, excess_divergent):
+ a_feat_l = itertools.islice(
+ feature_lambda(
+ sorted_genes,
+ feature_test_location,
+ {"loc": a[1]},
+ subfeatures=False,
+ ),
+ 1,
+ )
+ b_feat_l = itertools.islice(
+ feature_lambda(
+ sorted_genes,
+ feature_test_location,
+ {"loc": b[0]},
+ subfeatures=False,
+ ),
+ 1,
+ )
+
+ try:
+ a_feat = next(a_feat_l)
+ except StopIteration:
+ # Triggers on end of genome
+ a_feat = None
+ try:
+ b_feat = next(b_feat_l)
+ except StopIteration:
+ # Triggers on end of genome
+ b_feat = None
+
+ result_obj = [
+ a[1],
+ b[0],
+ None if not a_feat else a_feat.location.strand,
+ None if not b_feat else b_feat.location.strand,
+ ]
+
+ if a_feat is None or b_feat is None:
+ if gap_size > excess_divergent:
+ results.append(result_obj)
+ else:
+ if (
+ a_feat.location.strand == b_feat.location.strand
+ and gap_size > excess
+ ):
+ results.append(result_obj)
+ elif (
+ a_feat.location.strand != b_feat.location.strand
+ and gap_size > excess_divergent
+ ):
+ results.append(result_obj)
+
+ better_results = []
+ qc_features = []
+ of = MGAFinder(11, "CDS", "closed", min_gene)
+ # of = OrfFinder(11, 'CDS', 'closed', min_gene)
+
+ for result_obj in results:
+ start = result_obj[0]
+ end = result_obj[1]
+ f = gen_qc_feature(start, end, "Excessive gap, %s bases" % abs(end - start), type_src="gene")
+ qc_features.append(f)
+ putative_genes = of.putative_genes_in_sequence(
+ str(record[start - slop : end + slop].seq)
+ )
+ putative_genes = list(
+ require_sd(putative_genes, record, start, lookahead_min, lookahead_max)
+ )
+ for putative_gene in putative_genes:
+ # (0, 33, 1, 'ATTATTTTATCAAAACGCTTTACAATCTTTTAG', 'MILSKRFTIF', 123123, 124324)
+ possible_gene_start = start + putative_gene[0]
+ possible_gene_end = start + putative_gene[1]
+
+ if possible_gene_start <= possible_gene_end:
+ possible_cds = gffSeqFeature(
+ FeatureLocation(
+ possible_gene_start, possible_gene_end, strand=putative_gene[2]
+ ),
+ type="CDS",
+ )
+ else:
+ possible_cds = gffSeqFeature(
+ FeatureLocation(
+ possible_gene_end, possible_gene_start, strand=putative_gene[2],
+ ),
+ type="CDS",
+ )
+
+ # Now we adjust our boundaries for the RBS that's required
+ # There are only two cases, the rbs is upstream of it, or downstream
+ if putative_gene[5] < possible_gene_start:
+ possible_gene_start = putative_gene[5]
+ else:
+ possible_gene_end = putative_gene[6]
+
+ if putative_gene[5] <= putative_gene[6]:
+ possible_rbs = gffSeqFeature(
+ FeatureLocation(
+ putative_gene[5], putative_gene[6], strand=putative_gene[2]
+ ),
+ type="Shine_Dalgarno_sequence",
+ )
+ else:
+ possible_rbs = gffSeqFeature(
+ FeatureLocation(
+ putative_gene[6], putative_gene[5], strand=putative_gene[2],
+ ),
+ type="Shine_Dalgarno_sequence",
+ )
+
+ if possible_gene_start <= possible_gene_end:
+ possible_gene = gffSeqFeature(
+ FeatureLocation(
+ possible_gene_start, possible_gene_end, strand=putative_gene[2]
+ ),
+ type="gene",
+ qualifiers={"note": ["Possible gene"]},
+ )
+ else:
+ possible_gene = gffSeqFeature(
+ FeatureLocation(
+ possible_gene_end, possible_gene_start, strand=putative_gene[2],
+ ),
+ type="gene",
+ qualifiers={"note": ["Possible gene"]},
+ )
+ possible_gene.sub_features = [possible_rbs, possible_cds]
+ qc_features.append(possible_gene)
+
+ better_results.append(result_obj + [len(putative_genes)])
+
+ # Bad gaps are those with more than zero possible genes found
+ bad = len([x for x in better_results if x[2] > 0])
+ # Generally taking "good" here as every possible gap in the genome
+ # Thus, good is TOTAL - gaps
+ good = len(sorted_genes) + 1 - bad
+ # and bad is just gaps
+ return good, bad, better_results, qc_features
+
+
+def phi(x):
+ """Standard phi function used in calculation of normal distribution"""
+ return math.exp(-1 * math.pi * x * x)
+
+
+def norm(x, mean=0, sd=1):
+ """
+ Normal distribution. Given an x position, a mean, and a standard
+ deviation, calculate the "y" value. Useful for score scaling
+
+ Modified to multiply by SD. This means even at sd=5, norm(x, mean) where x = mean => 1, rather than 1/5.
+ """
+ return (1 / float(sd)) * phi(float(x - mean) / float(sd)) * sd
+
+
+def coding_density(record, mean=92.5, sd=20):
+ """
+ Find coding density in the genome
+ """
+ feature_lengths = 0
+
+ for gene_a in coding_genes(record.features):
+ feature_lengths += sum(
+ [len(x) for x in genes(gene_a.sub_features, feature_type="CDS")]
+ )
+
+ avgFeatLen = float(feature_lengths) / float(len(record.seq))
+ return int(norm(100 * avgFeatLen, mean=mean, sd=sd) * 100), int(100 * avgFeatLen)
+
+
+def exact_coding_density(record, mean=92.5, sd=20):
+ """
+ Find exact coding density in the genome
+ """
+ data = numpy.zeros(len(record.seq))
+
+ for gene_a in coding_genes(record.features):
+ for cds in genes(gene_a.sub_features, feature_type="CDS"):
+ for i in range(cds.location.start, cds.location.end + 1):
+ data[i - 1] = 1
+
+ return float(sum(data)) / len(data)
+
+
+def excessive_overlap(record, excess=15, excess_divergent=30):
+ """
+ Find excessive overlaps in the genome, where excessive is defined as 15
+ bases for same strand, and 30 for divergent translation.
+
+ Does a product of all the top-level features in the genome, and calculates
+ gaps.
+ """
+ results = []
+ bad = 0
+ qc_features = []
+
+ for (gene_a, gene_b) in itertools.combinations(coding_genes(record.features), 2):
+ # Get the CDS from the subfeature list.
+ # TODO: not recursive.
+ cds_a = [x for x in genes(gene_a.sub_features, feature_type="CDS")]
+ cds_b = [x for x in genes(gene_b.sub_features, feature_type="CDS")]
+
+ if len(cds_a) == 0:
+ log.warn("Gene missing subfeatures; %s", get_gff3_id(gene_a))
+ continue
+
+ if len(cds_b) == 0:
+ log.warn("Gene missing subfeatures; %s", get_gff3_id(gene_b))
+ continue
+
+ cds_a = cds_a[0]
+ cds_b = cds_b[0]
+
+ # Set of locations that are included in the CDS of A and the
+ # CDS of B
+ cas = set(range(cds_a.location.start, cds_a.location.end))
+ cbs = set(range(cds_b.location.start, cds_b.location.end))
+
+ # Here we calculate the intersection between the two sets, and
+ # if it's larger than our excessive size, we know that they're
+ # overlapped
+ ix = cas.intersection(cbs)
+
+ if (cds_a.location.strand == cds_b.location.strand and len(ix) >= excess) or (
+ cds_a.location.strand != cds_b.location.strand
+ and len(ix) >= excess_divergent
+ ):
+ bad += float(len(ix)) / float(min(excess, excess_divergent))
+ qc_features.append(
+ gen_qc_feature(min(ix), max(ix), "Excessive Overlap", id_src=gene_a, type_src="gene")
+ )
+ results.append((gene_a, gene_b, min(ix), max(ix)))
+
+ # Good isn't accurate here. It's a triangle number and just ugly, but we
+ # don't care enough to fix it.
+ good = len(list(coding_genes(record.features)))
+ good = int(good - bad)
+ if good < 0:
+ good = 0
+ return good, int(bad), results, qc_features
+
+
+def get_encouragement(score):
+ """Some text telling the user how they did
+ """
+ for encouragement in ENCOURAGEMENT:
+ if score > encouragement[0]:
+ return encouragement[1]
+ return ENCOURAGEMENT[-1][1]
+
+
+def genome_overview(record):
+ """Genome overview
+ """
+ data = {
+ "genes": {
+ "count": 0,
+ "bases": len(record.seq),
+ "density": 0, # genes / kb
+ "avg_len": [],
+ "comp": {"A": 0, "C": 0, "G": 0, "T": 0},
+ },
+ "overall": {
+ "comp": {
+ "A": record.seq.count("A") + record.seq.count("a"),
+ "C": record.seq.count("C") + record.seq.count("c"),
+ "G": record.seq.count("G") + record.seq.count("g"),
+ "T": record.seq.count("T") + record.seq.count("t"),
+ },
+ "gc": 0,
+ },
+ }
+ gene_features = list(coding_genes(record.features))
+ data["genes"]["count"] = len(gene_features)
+
+ for feat in gene_features:
+ data["genes"]["comp"]["A"] += feat.extract(record).seq.count("A") + feat.extract(record).seq.count("a")
+ data["genes"]["comp"]["C"] += feat.extract(record).seq.count("C") + feat.extract(record).seq.count("c")
+ data["genes"]["comp"]["T"] += feat.extract(record).seq.count("T") + feat.extract(record).seq.count("t")
+ data["genes"]["comp"]["G"] += feat.extract(record).seq.count("G") + feat.extract(record).seq.count("g")
+ #data["genes"]["bases"] += len(feat)
+ data["genes"]["avg_len"].append(len(feat))
+
+ data["genes"]["avg_len"] = float(sum(data["genes"]["avg_len"])) / len(gene_features)
+ data["overall"]["gc"] = float(
+ data["overall"]["comp"]["G"] + data["overall"]["comp"]["C"]
+ ) / len(record.seq)
+ return data
+
+
+def find_morons(record):
+ """Locate morons in the genome
+
+ Don't even know why...
+
+ TODO: remove? Idk.
+ """
+ results = []
+ good = 0
+ bad = 0
+
+ gene_features = list(coding_genes(record.features))
+ for i, gene in enumerate(gene_features):
+ two_left = gene_features[i - 2 : i]
+ two_right = gene_features[i + 1 : i + 1 + 2]
+ strands = [x.strand for x in two_left] + [x.strand for x in two_right]
+ anticon = [x for x in strands if x != gene.strand]
+
+ if len(anticon) == 4:
+ has_rbs = [x.type == "Shine_Dalgarno_sequence" for x in gene.sub_features]
+ if any(has_rbs):
+ rbs = [
+ x for x in gene.sub_features if x.type == "Shine_Dalgarno_sequence"
+ ][0]
+ rbs_msg = str(rbs.extract(record.seq))
+ else:
+ rbs_msg = "No RBS Available"
+ results.append((gene, two_left, two_right, rbs_msg))
+ bad += 1
+ else:
+ good += 1
+ return good, bad, results, []
+
+
+def bad_gene_model(record):
+ """Find features without product
+ """
+ results = []
+ good = 0
+ bad = 0
+ qc_features = []
+
+ for gene in coding_genes(record.features):
+ exons = [
+ x for x in genes(gene.sub_features, feature_type="exon") if len(x) > 10
+ ]
+ CDSs = [x for x in genes(gene.sub_features, feature_type="CDS")]
+ if len(exons) >= 1 and len(CDSs) >= 1:
+ if len(exons) != len(CDSs):
+ results.append(
+ (
+ get_gff3_id(gene),
+ None,
+ None,
+ "Mismatched number of exons and CDSs in gff3 representation",
+ )
+ )
+ qc_features.append(
+ gen_qc_feature(
+ gene.location.start,
+ gene.location.end,
+ "Mismatched number of exons and CDSs in gff3 representation",
+ strand=gene.strand,
+ id_src=gene,
+ type_src="gene"
+ )
+ )
+ bad += 1
+ else:
+ for (exon, cds) in zip(
+ sorted(exons, key=lambda x: x.location.start),
+ sorted(CDSs, key=lambda x: x.location.start),
+ ):
+ if len(exon) != len(cds):
+ results.append(
+ (
+ get_gff3_id(gene),
+ exon,
+ cds,
+ "CDS does not extend to full length of gene",
+ )
+ )
+ qc_features.append(
+ gen_qc_feature(
+ exon.location.start,
+ exon.location.end,
+ "CDS does not extend to full length of gene",
+ strand=exon.strand,
+ id_src=gene,
+ type_src="CDS"
+ )
+ )
+ bad += 1
+ else:
+ good += 1
+ else:
+ log.warn("Could not handle %s, %s", exons, CDSs)
+ results.append(
+ (
+ get_gff3_id(gene),
+ None,
+ None,
+ "{0} exons, {1} CDSs".format(len(exons), len(CDSs)),
+ )
+ )
+
+ return good, len(results) + bad, results, qc_features
+
+
+def weird_starts(record):
+ """Find features without product
+ """
+ good = 0
+ bad = 0
+ qc_features = []
+ results = []
+
+ overall = {}
+ for gene in coding_genes(record.features):
+ seq = [x for x in genes(gene.sub_features, feature_type="CDS")]
+ if len(seq) == 0:
+ log.warn("No CDS for gene %s", get_gff3_id(gene))
+ continue
+ else:
+ seq = seq[0]
+
+ seq_str = str(seq.extract(record.seq))
+ start_codon = seq_str[0:3]
+ if len(seq_str) < 3:
+ sys.stderr.write("Fatal Error: CDS of length less than 3 at " + str(seq.location) + '\n')
+ exit(2)
+# if len(seq_str) % 3 != 0:
+# if len(seq_str) < 3:
+# stop_codon = seq_str[-(len(seq_str))]
+# else:
+# stop_codon = seq_str[-3]
+#
+# log.warn("CDS at %s length is not a multiple of three (Length = %d)", get_gff3_id(gene), len(seq_str))
+# seq.__error = "Bad CDS Length"
+# results.append(seq)
+# qc_features.append(
+# gen_qc_feature(
+# s, e, "Bad Length", strand=seq.strand, id_src=gene
+# )
+# )
+# bad += 1
+# seq.__start = start_codon
+# seq.__stop = stop_codon
+# continue
+
+ stop_codon = seq_str[-3]
+ seq.__start = start_codon
+ seq.__stop = stop_codon
+ if start_codon not in overall:
+ overall[start_codon] = 1
+ else:
+ overall[start_codon] += 1
+
+ if start_codon not in ("ATG", "TTG", "GTG"):
+ log.warn("Weird start codon (%s) on %s", start_codon, get_gff3_id(gene))
+ seq.__error = "Unusual start codon %s" % start_codon
+
+ s = 0
+ e = 0
+ if seq.strand > 0:
+ s = seq.location.start
+ e = seq.location.start + 3
+ else:
+ s = seq.location.end
+ e = seq.location.end - 3
+
+ results.append(seq)
+ results[-1].location = FeatureLocation(results[-1].location.start + 1, results[-1].location.end, results[-1].location.strand)
+ qc_features.append(
+ gen_qc_feature(
+ s, e, "Weird start codon", strand=seq.strand, id_src=gene, type_src="gene"
+ )
+ )
+ bad += 1
+ else:
+ good += 1
+
+ return good, bad, results, qc_features, overall
+
+
+def missing_genes(record):
+ """Find features without product
+ """
+ results = []
+ good = 0
+ bad = 0
+ qc_features = []
+
+ for gene in coding_genes(record.features):
+ if gene.qualifiers.get("cpt_source", [None])[0] == "CPT_GENE_MODEL_CORRECTION":
+ results.append(gene)
+ bad += 1
+ else:
+ good += 1
+
+ return good, bad, results, qc_features
+
+
+def gene_model_correction_issues(record):
+ """Find features that have issues from the gene model correction step.
+ These have qualifiers beginning with CPT_GMS
+ """
+ results = []
+ good = 0
+ bad = 0
+ qc_features = []
+
+ # For each gene
+ for gene in coding_genes(record.features):
+ # Get the list of child CDSs
+ cdss = [x for x in genes(gene.sub_features, feature_type="CDS")]
+ # And our matching qualifiers
+ gene_data = [(k, v) for (k, v) in gene.qualifiers.items() if k == "cpt_gmc"]
+ # If there are problems with ONLY the parent, let's complain
+ local_results = []
+ local_qc_features = []
+ for x in gene_data:
+ if "Missing Locus Tag" in x[1]:
+ # Missing locus tag is an either or thing, if it hits here
+ # there shouldn't be anything else wrong with it.
+
+ # Obviously missing so we remove it
+ gene.qualifiers["locus_tag"] = [""]
+ # Translation from bp_genbank2gff3.py
+ cdss[0].qualifiers["locus_tag"] = cdss[0].qualifiers["Name"]
+ # Append our results
+ local_results.append((gene, cdss[0], "Gene is missing a locus_tag"))
+ local_qc_features.append(
+ gen_qc_feature(
+ gene.location.start,
+ gene.location.end,
+ "Gene is missing a locus_tag",
+ strand=gene.strand,
+ type_src="gene"
+ )
+ )
+
+ # We need to alert on any child issues as well.
+ for cds in cdss:
+ cds_data = [
+ (k, v[0]) for (k, v) in cds.qualifiers.items() if k == "cpt_gmc"
+ ]
+ if len(gene_data) == 0 and len(cds_data) == 0:
+ # Alles gut
+ pass
+ else:
+ for _, problem in cds_data:
+ if problem == "BOTH Missing Locus Tag":
+ gene.qualifiers["locus_tag"] = [""]
+ cds.qualifiers["locus_tag"] = [""]
+ local_results.append(
+ (gene, cds, "Both gene and CDS are missing locus tags")
+ )
+ local_qc_features.append(
+ gen_qc_feature(
+ cds.location.start,
+ cds.location.end,
+ "CDS is missing a locus_tag",
+ strand=cds.strand,
+ type_src="CDS"
+ )
+ )
+ local_qc_features.append(
+ gen_qc_feature(
+ gene.location.start,
+ gene.location.end,
+ "Gene is missing a locus_tag",
+ strand=gene.strand,
+ type_src="gene"
+ )
+ )
+ elif problem == "Different locus tag from associated gene.":
+ gene.qualifiers["locus_tag"] = gene.qualifiers["Name"]
+ cds.qualifiers["locus_tag"] = cds.qualifiers["cpt_gmc_locus"]
+ local_results.append(
+ (gene, cds, "Gene and CDS have differing locus tags")
+ )
+ local_qc_features.append(
+ gen_qc_feature(
+ gene.location.start,
+ gene.location.end,
+ "Gene and CDS have differing locus tags",
+ strand=gene.strand,
+ type_src="gene"
+ )
+ )
+ elif problem == "Missing Locus Tag":
+ # Copy this over
+ gene.qualifiers["locus_tag"] = gene.qualifiers["Name"]
+ # This one is missing
+ cds.qualifiers["locus_tag"] = [""]
+ local_results.append((gene, cds, "CDS is missing a locus_tag"))
+ local_qc_features.append(
+ gen_qc_feature(
+ cds.location.start,
+ cds.location.end,
+ "CDS is missing a locus_tag",
+ strand=cds.strand,
+ type_src="CDS"
+ )
+ )
+ else:
+ log.warn("Cannot handle %s", problem)
+
+ if len(local_results) > 0:
+ bad += 1
+ else:
+ good += 1
+
+ qc_features.extend(local_qc_features)
+ results.extend(local_results)
+ return good, bad, results, qc_features
+
+
+def missing_tags(record):
+ """Find features without product
+ """
+ results = []
+ good = 0
+ bad = 0
+ qc_features = []
+
+ for gene in coding_genes(record.features):
+ cds = [x for x in genes(gene.sub_features, feature_type="CDS")]
+ if len(cds) == 0:
+ log.warn("Gene missing CDS subfeature %s", get_gff3_id(gene))
+ continue
+
+ cds = cds[0]
+
+ if "product" not in cds.qualifiers:
+ log.info("Missing product tag on %s", get_gff3_id(gene))
+ qc_features.append(
+ gen_qc_feature(
+ cds.location.start,
+ cds.location.end,
+ "Missing product tag",
+ strand=cds.strand,
+ type_src="CDS"
+ )
+ )
+ results.append(cds)
+ bad += 1
+ else:
+ good += 1
+
+ return good, bad, results, qc_features
+
+
+def evaluate_and_report(
+ annotations,
+ genome,
+ gff3=None,
+ tbl=None,
+ sd_min=5,
+ sd_max=15,
+ min_gene_length=30,
+ excessive_gap_dist=50,
+ excessive_gap_divergent_dist=200,
+ excessive_overlap_dist=25,
+ excessive_overlap_divergent_dist=50,
+ reportTemplateName="phage_annotation_validator.html",
+):
+ """
+ Generate our HTML evaluation of the genome
+ """
+ # Get features from GFF file
+ seq_dict = SeqIO.to_dict(SeqIO.parse(genome, "fasta"))
+ # Get the first GFF3 record
+ # TODO: support multiple GFF3 files.
+ mostFeat = 0
+ for rec in list(gffParse(annotations, base_dict=seq_dict)):
+ if len(rec.features) > mostFeat:
+ mostFeat = len(rec.features)
+ record = rec
+
+ gff3_qc_record = SeqRecord(record.id, id=record.id)
+ gff3_qc_record.features = []
+ gff3_qc_features = []
+
+ log.info("Locating missing RBSs")
+ # mb_any = "did they annotate ANY rbss? if so, take off from score."
+ mb_good, mb_bad, mb_results, mb_annotations, mb_any = missing_rbs(
+ record, lookahead_min=sd_min, lookahead_max=sd_max
+ )
+ gff3_qc_features += mb_annotations
+
+ log.info("Locating excessive gaps")
+ eg_good, eg_bad, eg_results, eg_annotations = excessive_gap(
+ record,
+ excess=excessive_gap_dist,
+ excess_divergent=excessive_gap_divergent_dist,
+ min_gene=min_gene_length,
+ slop=excessive_overlap_dist,
+ lookahead_min=sd_min,
+ lookahead_max=sd_max,
+ )
+ gff3_qc_features += eg_annotations
+
+ log.info("Locating excessive overlaps")
+ eo_good, eo_bad, eo_results, eo_annotations = excessive_overlap(
+ record,
+ excess=excessive_overlap_dist,
+ excess_divergent=excessive_overlap_divergent_dist,
+ )
+ gff3_qc_features += eo_annotations
+
+ log.info("Locating morons")
+ mo_good, mo_bad, mo_results, mo_annotations = find_morons(record)
+ gff3_qc_features += mo_annotations
+
+ log.info("Locating missing tags")
+ mt_good, mt_bad, mt_results, mt_annotations = missing_tags(record)
+ gff3_qc_features += mt_annotations
+
+ log.info("Locating missing gene features")
+ mg_good, mg_bad, mg_results, mg_annotations = missing_genes(record)
+ gff3_qc_features += mg_annotations
+
+ log.info("Determining coding density")
+ cd, cd_real = coding_density(record)
+
+ log.info("Locating weird starts")
+ ws_good, ws_bad, ws_results, ws_annotations, ws_overall = weird_starts(record)
+ gff3_qc_features += ws_annotations
+
+ log.info("Locating bad gene models")
+ gm_good, gm_bad, gm_results, gm_annotations = bad_gene_model(record)
+ if gm_good + gm_bad == 0:
+ gm_bad = 1
+
+ log.info("Locating more bad gene models")
+ gmc_good, gmc_bad, gmc_results, gmc_annotations = gene_model_correction_issues(
+ record
+ )
+ if gmc_good + gmc_bad == 0:
+ gmc_bad = 1
+
+ good_scores = [eg_good, eo_good, mt_good, ws_good, gm_good, gmc_good]
+ bad_scores = [eg_bad, eo_bad, mt_bad, ws_bad, gm_bad, gmc_bad]
+
+ # Only if they tried to annotate RBSs do we consider them.
+ if mb_any:
+ good_scores.append(mb_good)
+ bad_scores.append(mb_bad)
+ subscores = []
+
+ for (g, b) in zip(good_scores, bad_scores):
+ if g + b == 0:
+ s = 0
+ else:
+ s = int(100 * float(g) / (float(b) + float(g)))
+ subscores.append(s)
+ subscores.append(cd)
+
+ score = int(float(sum(subscores)) / float(len(subscores)))
+
+ # This is data that will go into our HTML template
+ kwargs = {
+ "upstream_min": sd_min,
+ "upstream_max": sd_max,
+ "record_name": record.id,
+ "record_nice_name": nice_name(record),
+ "params": {
+ "sd_min": sd_min,
+ "sd_max": sd_max,
+ "min_gene_length": min_gene_length,
+ "excessive_gap_dist": excessive_gap_dist,
+ "excessive_gap_divergent_dist": excessive_gap_divergent_dist,
+ "excessive_overlap_dist": excessive_overlap_dist,
+ "excessive_overlap_divergent_dist": excessive_overlap_divergent_dist,
+ },
+ "score": score,
+ "encouragement": get_encouragement(score),
+ "genome_overview": genome_overview(record),
+ "rbss_annotated": mb_any,
+ "missing_rbs": mb_results,
+ "missing_rbs_good": mb_good,
+ "missing_rbs_bad": mb_bad,
+ "missing_rbs_score": 0
+ if mb_good + mb_bad == 0
+ else (100 * mb_good / (mb_good + mb_bad)),
+ "excessive_gap": eg_results,
+ "excessive_gap_good": eg_good,
+ "excessive_gap_bad": eg_bad,
+ "excessive_gap_score": 0
+ if eo_good + eo_bad == 0
+ else (100 * eo_good / (eo_good + eo_bad)),
+ "excessive_overlap": eo_results,
+ "excessive_overlap_good": eo_good,
+ "excessive_overlap_bad": eo_bad,
+ "excessive_overlap_score": 0
+ if eo_good + eo_bad == 0
+ else (100 * eo_good / (eo_good + eo_bad)),
+ "morons": mo_results,
+ "morons_good": mo_good,
+ "morons_bad": mo_bad,
+ "morons_score": 0
+ if mo_good + mo_bad == 0
+ else (100 * mo_good / (mo_good + mo_bad)),
+ "missing_tags": mt_results,
+ "missing_tags_good": mt_good,
+ "missing_tags_bad": mt_bad,
+ "missing_tags_score": 0
+ if mt_good + mt_bad == 0
+ else (100 * mt_good / (mt_good + mt_bad)),
+ "missing_genes": mg_results,
+ "missing_genes_good": mg_good,
+ "missing_genes_bad": mg_bad,
+ "missing_genes_score": 0
+ if mg_good + mg_bad == 0
+ else (100 * mg_good / (mg_good + mg_bad)),
+ "weird_starts": ws_results,
+ "weird_starts_good": ws_good,
+ "weird_starts_bad": ws_bad,
+ "weird_starts_overall": ws_overall,
+ "weird_starts_overall_sorted_keys": sorted(
+ ws_overall, reverse=True, key=lambda x: ws_overall[x]
+ ),
+ "weird_starts_score": 0
+ if ws_good + ws_bad == 0
+ else (100 * ws_good / (ws_good + ws_bad)),
+ "gene_model": gm_results,
+ "gene_model_good": gm_good,
+ "gene_model_bad": gm_bad,
+ "gene_model_score": 0
+ if gm_good + gm_bad == 0
+ else (100 * gm_good / (gm_good + gm_bad)),
+ "gene_model_correction": gmc_results,
+ "gene_model_correction_good": gmc_good,
+ "gene_model_correction_bad": gmc_bad,
+ "gene_model_correction_score": 0
+ if gmc_good + gmc_bad == 0
+ else (100 * gmc_good / (gmc_good + gmc_bad)),
+ "coding_density": cd,
+ "coding_density_exact": exact_coding_density(record),
+ "coding_density_real": cd_real,
+ "coding_density_score": cd,
+ }
+
+ with open(tbl, "w") as handle:
+ kw_subset = {}
+ for key in kwargs:
+ if (
+ key in ("score", "record_name")
+ or "_good" in key
+ or "_bad" in key
+ or "_overall" in key
+ ):
+ kw_subset[key] = kwargs[key]
+ json.dump(kw_subset, handle)
+
+ with open(gff3, "w") as handle:
+ gff3_qc_record.features = gff3_qc_features
+ gff3_qc_record.annotations = {}
+ gffWrite([gff3_qc_record], handle)
+
+ def nice_strand(direction):
+ # It is somehow possible for whole gffSeqFeature objects to end up in here, apparently at the gene level
+ if "SeqFeature" in str(type(direction)):
+ direction = direction.location.strand
+ if direction > 0:
+ return "→"#.decode("utf-8")
+ else:
+ return "←"#.decode("utf-8")
+
+ def nice_strand_tex(direction):
+ if "SeqFeature" in str(type(direction)):
+ direction = direction.location.strand
+ if direction > 0:
+ return "$\\rightarrow$"
+ else:
+ return "$\\leftarrow$"
+
+ def texify(data):
+ return data.replace("_", "\\_").replace("$", "\\$")
+
+ def length(data):
+ return len(data)
+
+ def my_encode(data):
+ return str(data)#.encode("utf-8")
+
+ def my_decode(data):
+ # For production
+ return str(data)#.decode("utf-8")
+ # For local testing. No, I do not understand.
+ return str(data)#.encode("utf-8")).decode("utf-8")
+
+ env = Environment(
+ loader=FileSystemLoader(SCRIPT_PATH), trim_blocks=True, lstrip_blocks=True
+ )
+ env.filters.update(
+ {
+ "nice_id": get_gff3_id,
+ "nice_strand": nice_strand,
+ "nice_strand_tex": nice_strand_tex,
+ "texify": texify,
+ "length": length,
+ "encode": my_encode,
+ "decode": my_decode,
+ }
+ )
+ tpl = env.get_template(reportTemplateName)
+ return tpl.render(**kwargs)#.encode("utf-8")
+
+
+if __name__ == "__main__":
+ parser = argparse.ArgumentParser(
+ description="rebase gff3 features against parent locations", epilog=""
+ )
+ parser.add_argument(
+ "annotations", type=argparse.FileType("r"), help="Parent GFF3 annotations"
+ )
+ parser.add_argument("genome", type=argparse.FileType("r"), help="Genome Sequence")
+ parser.add_argument(
+ "--gff3", type=str, help="GFF3 Annotations", default="qc_annotations.gff3"
+ )
+ parser.add_argument(
+ "--tbl",
+ type=str,
+ help="Table for noninteractive parsing",
+ default="qc_results.json",
+ )
+
+ parser.add_argument(
+ "--sd_min",
+ type=int,
+ help="Minimum distance from gene start for an SD to be",
+ default=5,
+ )
+ parser.add_argument(
+ "--sd_max",
+ type=int,
+ help="Maximum distance from gene start for an SD to be",
+ default=15,
+ )
+
+ parser.add_argument(
+ "--min_gene_length",
+ type=int,
+ help="Minimum length for a putative gene call (AAs)",
+ default=30,
+ )
+
+ parser.add_argument(
+ "--excessive_overlap_dist",
+ type=int,
+ help="Excessive overlap for genes in same direction",
+ default=25,
+ )
+ parser.add_argument(
+ "--excessive_overlap_divergent_dist",
+ type=int,
+ help="Excessive overlap for genes in diff directions",
+ default=50,
+ )
+
+ parser.add_argument(
+ "--excessive_gap_dist",
+ type=int,
+ help="Maximum distance between two genes",
+ default=40,
+ )
+ parser.add_argument(
+ "--excessive_gap_divergent_dist",
+ type=int,
+ help="Maximum distance between two divergent genes",
+ default=200,
+ )
+
+ parser.add_argument(
+ "--reportTemplateName",
+ help="Report template file name",
+ default="phageqc_report_full.html",
+ )
+
+ args = parser.parse_args()
+
+ sys.stdout.write(evaluate_and_report(**vars(args)))
diff -r 000000000000 -r c3140b08d703 cpt_phageqc_annotation/phage_annotation_validator.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/cpt_phageqc_annotation/phage_annotation_validator.xml Fri Jun 17 13:00:50 2022 +0000
@@ -0,0 +1,102 @@
+
+
+ validate phage annotations
+
+ macros.xml
+ cpt-macros.xml
+
+
+ python
+ biopython
+ cpt_gffparser
+ python-levenshtein
+ regex
+ metagene_annotator
+ jinja2
+ numpy
+
+ $output;
+
+#if ".tex" in str($report_format):
+ mv $output tmp.tex;
+ docker run --rm -i --user="1002:1002" --net=none -v \$PWD:/data blang/latex pdflatex tmp.tex &&
+ docker run --rm -i --user="1002:1002" --net=none -v \$PWD:/data blang/latex pdflatex tmp.tex &&
+ mv tmp.pdf $output;
+#end if
+]]>
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff -r 000000000000 -r c3140b08d703 cpt_phageqc_annotation/phageqc_report_464.html
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/cpt_phageqc_annotation/phageqc_report_464.html Fri Jun 17 13:00:50 2022 +0000
@@ -0,0 +1,419 @@
+
+
+
+
+
+
+
+
+
+ [BICH464] Phage QC on {{record_name}} - {{score}}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Phage {{record_name}}
+
+
+
+
+
+
Genome Overview
+
Genes
+
+ - Count: {{ genome_overview.genes.count }}
+ - Bases: {{ genome_overview.genes.bases }}
+ - Average Length: {{ genome_overview.genes.avg_len | round | int}}
+ - Coding Density: {{ '%0.2f' % (100 * coding_density_exact) }}%
+ - Composition
+
+ - A {{ genome_overview.genes.comp.A }}
+ - C {{ genome_overview.genes.comp.C }}
+ - T {{ genome_overview.genes.comp.T }}
+ - G {{ genome_overview.genes.comp.G }}
+
+
+
+
Overall
+
+ - %GC: {{ '%0.2f' % (100 * genome_overview.overall.gc) }}%
+ - Composition
+
+ - A {{ genome_overview.overall.comp.A }}
+ - C {{ genome_overview.overall.comp.C }}
+ - T {{ genome_overview.overall.comp.T }}
+ - G {{ genome_overview.overall.comp.G }}
+
+
+
+
+
+
+
+
Genes Missing RBS {{missing_rbs_good}} / {{missing_rbs_good + missing_rbs_bad}}
+
The following genes have issues with their RBS.
+ {% if not rbss_annotated %}
+
+ Since you have not annotated any possible RBSs, this does not count off from your overall score.
+
+ {% endif %}
+
+
+
+
+ | ID |
+ Location |
+ Error |
+ Upstream (-{{upstream_max}} .. -{{upstream_min}}) |
+
+
+
+ {% for row in missing_rbs %}
+ {% if "None found" in row.__message %}
+
+ | {{row | nice_id | decode}} |
+ {{row.location.start}}..{{row.location.end}} [{{row.strand}}] |
+ None found |
+ {{row.__upstream }} |
+
+ {% endif %}
+ {% endfor %}
+
+
+
+
+
Start Codon Usage
+
This section covers genes with unusual start codons
+
+
+
+
+ | Start Codon |
+ Count |
+
+
+
+ {% for codon_key in weird_starts_overall_sorted_keys %}
+ | {{ codon_key }} | {{ weird_starts_overall[codon_key] }} |
+ {% endfor %}
+
+
+
+
+
+
+
+
+ | ID |
+ Location |
+ Error |
+
+
+
+ {% for row in weird_starts %}
+
+ | {{row | nice_id| decode}} |
+ {{row.location.start}}..{{row.location.end}} [{{row.strand}}] |
+ {{row.qualifiers.get('note', [])}} |
+
+ {% endfor %}
+
+
+
+
+
Intergenic Gaps
+
Phage genomes are under pressure to maintain high coding density. Large intergenic gaps may be a sign of incorrect gene starts or missing genes.
+
+
+
+
+ | Region |
+ Size |
+ Bounding Gene Transcription Direction |
+ Message |
+
+
+
+ {% for row in excessive_gap %}
+
+ | {{row[0]}} .. {{row[1]}} |
+ {{row[1] - row[0]}} |
+ {{row[2] | nice_strand}} {{row[3] | nice_strand}} |
+
+ {% if row[4] == 0 %}
+ {% else %}
+ {{row[4]}} possible genes found in this region
+ {% endif %}
+ |
+
+ {% endfor %}
+
+
+
+
+
Overlapping Genes
+
Large gene overlaps may indicate an incorrect gene start or miscalled gene.
+
+
+
+
+ | Feature A |
+ Feature B |
+ Shared Region |
+ Overlap Length |
+
+
+
+ {% for row in excessive_overlap %}
+
+ | {{row[0] | nice_id | decode}} ({{row[0].location}}) |
+ {{row[1] | nice_id | decode}} ({{row[1].location}}) |
+ {{row[2]}}..{{row[3]}} |
+ {{row[3] - row[2]}}bp |
+
+ {% endfor %}
+
+
+
+
+
+
+
These issues are mostly derived from how Apollo handles the gene model. 
+
+
+
+
+ | ID |
+ Exon |
+ CDS |
+ Message |
+
+
+
+ {% for row in gene_model %}
+
+ | {{row[0]}} |
+ {{row[1].location}} |
+ {{row[2].location}} |
+ {{row[3]}} |
+
+ {% endfor %}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Data on each organism will be accessible from the tabs above.
+
+ {% for (record, data) in annotation_table_data %}
+
+
+
+
+ {% for col in annotation_table_col_names %}
+ | {{ col[0] }} |
+ {% endfor %}
+
+
+
+ {% for row in data %}
+
+ {% for col in row %}
+ {% if col is not string %}{% for val in col %}- {{ val }}
{% endfor %} {% else %}{{ col }}{% endif %} |
+ {% endfor %}
+
+ {% endfor %}
+
+
+
+ {% endfor %}
+
+
+
+
+
+
+
+
+
+
+
Phage {{record_name}}
+
+
+
+
+
+
+
+
Genes missing RBS {{missing_rbs_good}} / {{missing_rbs_good + missing_rbs_bad}}
+
The following genes have issues with their RBS.
+ {% if not rbss_annotated %}
+
+ Since you have not annotated any possible RBSs, this does not count off from your overall score.
+
+ {% endif %}
+
+
+
+
+ | Feature Type |
+ ID |
+ Location |
+ Error |
+ Upstream (-{{upstream_max}} .. -{{upstream_min}}) |
+
+
+
+ {% for row in missing_rbs %}
+
+ | {{row.type}} |
+ {{row.id}} |
+ {{row.location.start}}..{{row.location.end}} [{{row.strand}}] |
+ {{row.__message}} |
+ {{row.__upstream}} |
+
+ {% endfor %}
+
+
+
+
+
Start Codon Usage
+
This section covers genes with unusual start codons
+
+
+
+
+ | Start Codon |
+ Count |
+
+
+
+ {% for codon_key in weird_starts_overall_sorted_keys %}
+ | {{ codon_key }} | {{ weird_starts_overall[codon_key] }} |
+ {% endfor %}
+
+
+
+
+
+
+
+
+ | Feature Type |
+ ID |
+ Location |
+ Error |
+
+
+
+ {% for row in weird_starts %}
+
+ | {{row.type}} |
+ {{row.id}} |
+ {{row.location.start}}..{{row.location.end}} [{{row.strand}}] |
+ {{row.__error}} |
+
+ {% endfor %}
+
+
+
+
+
Intergenic Gaps
+
Phage genomes are under pressure to maintain high coding density. Large intergenic gaps may be a sign of incorrect gene starts or missing genes.
+
+
+
+
+ | Region |
+ Size |
+ Bounding Gene Transcription Direction |
+ Message |
+
+
+
+ {% for row in excessive_gap %}
+
+ | {{row[0]}} .. {{row[1]}} |
+ {{row[1] - row[0]}} |
+ {{row[2] | nice_strand}} {{row[3] | nice_strand}} |
+
+ {% if row[4] == 0 %}
+ {% else %}
+ {{row[4]}} ORFs found in this region
+ {% endif %}
+ |
+
+ {% endfor %}
+
+
+
+
+
Overlapping Genes
+
Large gene overlaps may indicate an incorrect gene start or miscalled gene.
+
+
+
+
+ | Feature A |
+ Feature B |
+ Shared Region |
+ Overlap Length |
+
+
+
+ {% for row in excessive_overlap %}
+
+ | {{row[0].id}} ({{row[0].location}}) |
+ {{row[1].id}} ({{row[1].location}}) |
+ {{row[2]}}..{{row[3]}} |
+ {{row[3] - row[2]}} bp |
+
+ {% endfor %}
+
+
+
+
+
+
+
Possible Morons {{morons_good}} / {{morons_good + morons_bad}} (Doesn't count towards score)
+
+
+
+
+ | Feature |
+ RBS |
+ Surrounding Features |
+
+
+
+ {% for row in morons %}
+
+ | {{row[0].id}} |
+ {{row[3]}} |
+
+ {% for x in row[1] %}
+ {{ x | nice_strand }}
+ {% endfor %}
+ {{ row[0].strand | nice_strand }}
+ {% for x in row[2] %}
+ {{ x | nice_strand }}
+ {% endfor %}
+
+ |
+
+ {% endfor %}
+
+
+
+
+
+
+
Missing Product Tags {{missing_tags_good}} / {{missing_tags_good + missing_tags_bad}}
+
+
+
+
+ | Feature |
+ Qualifiers |
+
+
+
+ {% for row in missing_tags %}
+
+ | {{row.id}} |
+
+ {% for key in row.qualifiers %}
+ {{ key }}
+
+ {% for value in row.qualifiers[key] %}
+ - {{value}}
+ {% endfor %}
+
+ {% endfor %}
+ |
+
+ {% endfor %}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Phage AY216660.2
+
+
+
+
+
+
+
+
Genes missing RBS 0 / 79
+
The following genes have issues with their RBS.
+
+ Since you have not annotated any possible RBSs, this does not count off from your overall score.
+
+
+
+
+
+ | Feature Type |
+ ID |
+ Location |
+ Error |
+ Upstream (-15 .. -5) |
+
+
+
+
+ | gene |
+ CPT-T1_001.gene |
+ 40..576 [1] |
+ No RBS annotated, None found |
+ taaatgttaa |
+
+
+ | gene |
+ CPT-T1_002.gene |
+ 589..2184 [1] |
+ Unannotated but valid RBS |
+ taacttt AGG |
+
+
+ | gene |
+ CPT-T1_003.gene |
+ 2230..3522 [1] |
+ No RBS annotated, None found |
+ acggcccttt |
+
+
+ | gene |
+ CPT-T1_004.gene |
+ 3496..4273 [1] |
+ Unannotated but valid RBS |
+ cc GGAG ccgg |
+
+
+ | gene |
+ CPT-T1_005.gene |
+ 4264..5388 [1] |
+ No RBS annotated, None found |
+ tccagacttt |
+
+
+ | gene |
+ CPT-T1_006.gene |
+ 5389..5876 [1] |
+ No RBS annotated, None found |
+ agcaaagtaa |
+
+
+ | gene |
+ CPT-T1_007.gene |
+ 5926..6705 [1] |
+ No RBS annotated, None found |
+ tagtcctttt |
+
+
+ | gene |
+ CPT-T1_008.gene |
+ 6784..7755 [1] |
+ No RBS annotated, None found |
+ ttttttatta |
+
+
+ | gene |
+ CPT-T1_009.gene |
+ 7794..8093 [1] |
+ No RBS annotated, None found |
+ gtcccttttt |
+
+
+ | gene |
+ CPT-T1_010.gene |
+ 8127..8548 [1] |
+ No RBS annotated, None found |
+ agcgcctttt |
+
+
+ | gene |
+ CPT-T1_011.gene |
+ 8534..8919 [1] |
+ Unannotated but valid RBS |
+ atgac AGGA c |
+
+
+ | gene |
+ CPT-T1_012.gene |
+ 8897..9355 [1] |
+ No RBS annotated, None found |
+ atcgtccagt |
+
+
+ | gene |
+ CPT-T1_013.gene |
+ 9332..9743 [1] |
+ No RBS annotated, None found |
+ gcaatcaagc |
+
+
+ | gene |
+ CPT-T1_014.gene |
+ 9733..10414 [1] |
+ Unannotated but valid RBS |
+ tac GAG taga |
+
+
+ | gene |
+ CPT-T1_015.gene |
+ 10516..10845 [1] |
+ No RBS annotated, None found |
+ tgctatcaac |
+
+
+ | gene |
+ CPT-T1_016.gene |
+ 10516..11162 [1] |
+ No RBS annotated, None found |
+ tgctatcaac |
+
+
+ | gene |
+ CPT-T1_017.gene |
+ 11192..14076 [1] |
+ No RBS annotated, None found |
+ gttttttcgt |
+
+
+ | gene |
+ CPT-T1_018.gene |
+ 14067..14432 [1] |
+ Unannotated but valid RBS |
+ cattc AGGA a |
+
+
+ | gene |
+ CPT-T1_019.gene |
+ 14500..15294 [1] |
+ Unannotated but valid RBS |
+ gttgc AGGT a |
+
+
+ | gene |
+ CPT-T1_020.gene |
+ 15279..16025 [1] |
+ No RBS annotated, None found |
+ tgtcgcttca |
+
+
+ | gene |
+ CPT-T1_021.gene |
+ 16009..16621 [1] |
+ No RBS annotated, None found |
+ taatcgttcg |
+
+
+ | gene |
+ CPT-T1_022.gene |
+ 16688..20217 [1] |
+ No RBS annotated, None found |
+ ataaatagca |
+
+
+ | gene |
+ CPT-T1_023.gene |
+ 20251..20568 [1] |
+ No RBS annotated, None found |
+ atgccctttt |
+
+
+ | gene |
+ CPT-T1_024.gene |
+ 20556..21257 [1] |
+ No RBS annotated, None found |
+ cttaatagca |
+
+
+ | gene |
+ CPT-T1_025.gene |
+ 21279..21518 [-1] |
+ No RBS annotated, None found |
+ aatcacacta |
+
+
+ | gene |
+ CPT-T1_026.gene |
+ 21755..21906 [1] |
+ No RBS annotated, None found |
+ tgtatatcgt |
+
+
+ | gene |
+ CPT-T1_027.gene |
+ 22020..23098 [1] |
+ No RBS annotated, None found |
+ attcatcgta |
+
+
+ | gene |
+ CPT-T1_028.gene |
+ 23129..23820 [1] |
+ No RBS annotated, None found |
+ aaacaaataa |
+
+
+ | gene |
+ CPT-T1_029.gene |
+ 23854..24289 [1] |
+ No RBS annotated, None found |
+ tcgcccataa |
+
+
+ | gene |
+ CPT-T1_030.gene |
+ 24351..26550 [-1] |
+ No RBS annotated, None found |
+ cgtgctattt |
+
+
+ | gene |
+ CPT-T1_031.gene |
+ 26447..26600 [-1] |
+ No RBS annotated, None found |
+ tgcccattgc |
+
+
+ | gene |
+ CPT-T1_032.gene |
+ 26638..27585 [-1] |
+ No RBS annotated, None found |
+ cgcgtttttt |
+
+
+ | gene |
+ CPT-T1_033.gene |
+ 27635..28098 [-1] |
+ No RBS annotated, None found |
+ gcaaaaagtg |
+
+
+ | gene |
+ CPT-T1_034.gene |
+ 28168..30199 [1] |
+ No RBS annotated, None found |
+ gttacaacga |
+
+
+ | gene |
+ CPT-T1_035.gene |
+ 30181..30612 [1] |
+ No RBS annotated, None found |
+ caaagacttc |
+
+
+ | gene |
+ CPT-T1_036.gene |
+ 30669..31393 [1] |
+ No RBS annotated, None found |
+ atctcaccaa |
+
+
+ | gene |
+ CPT-T1_037.gene |
+ 31377..31641 [1] |
+ No RBS annotated, None found |
+ aaaacaaacg |
+
+
+ | gene |
+ CPT-T1_038.gene |
+ 31695..31917 [1] |
+ No RBS annotated, None found |
+ tcatagaaac |
+
+
+ | gene |
+ CPT-T1_039.gene |
+ 31881..32167 [1] |
+ No RBS annotated, None found |
+ cttcatctct |
+
+
+ | gene |
+ CPT-T1_040.gene |
+ 32099..32399 [1] |
+ Unannotated but valid RBS |
+ GAG ttttacc |
+
+
+ | gene |
+ CPT-T1_041.gene |
+ 32467..33611 [1] |
+ No RBS annotated, None found |
+ tttcattatc |
+
+
+ | gene |
+ CPT-T1_042.gene |
+ 33673..34166 [1] |
+ No RBS annotated, None found |
+ ttttatagaa |
+
+
+ | gene |
+ CPT-T1_043.gene |
+ 34226..34415 [1] |
+ No RBS annotated, None found |
+ accacatcga |
+
+
+ | gene |
+ CPT-T1_044.gene |
+ 34525..34749 [1] |
+ No RBS annotated, None found |
+ gttcaaaaaa |
+
+
+ | gene |
+ CPT-T1_045.gene |
+ 34739..35237 [1] |
+ Unannotated but valid RBS |
+ AGG cgcttaa |
+
+
+ | gene |
+ CPT-T1_046.gene |
+ 35227..35638 [1] |
+ No RBS annotated, None found |
+ acgcatattg |
+
+
+ | gene |
+ CPT-T1_047.gene |
+ 35761..36179 [-1] |
+ No RBS annotated, None found |
+ caatcctcga |
+
+
+ | gene |
+ CPT-T1_048.gene |
+ 36173..37755 [-1] |
+ No RBS annotated, None found |
+ tcacgccgtc |
+
+
+ | gene |
+ CPT-T1_049.gene |
+ 37741..38294 [-1] |
+ No RBS annotated, None found |
+ gacggcacaa |
+
+
+ | gene |
+ CPT-T1_050.gene |
+ 38277..38706 [-1] |
+ No RBS annotated, None found |
+ tcaagataac |
+
+
+ | gene |
+ CPT-T1_051.gene |
+ 38774..38994 [-1] |
+ No RBS annotated, None found |
+ tttacattaa |
+
+
+ | gene |
+ CPT-T1_052.gene |
+ 38987..39222 [-1] |
+ No RBS annotated, None found |
+ taccaaacaa |
+
+
+ | gene |
+ CPT-T1_053.gene |
+ 39289..39444 [-1] |
+ No RBS annotated, None found |
+ ggcatatcaa |
+
+
+ | gene |
+ CPT-T1_054.gene |
+ 39429..39758 [-1] |
+ No RBS annotated, None found |
+ tatcctgact |
+
+
+ | gene |
+ CPT-T1_055.gene |
+ 39766..39978 [-1] |
+ Unannotated but valid RBS |
+ tagc GGA tcg |
+
+
+ | gene |
+ CPT-T1_056.gene |
+ 39959..40342 [-1] |
+ No RBS annotated, None found |
+ acatcaacag |
+
+
+ | gene |
+ CPT-T1_057.gene |
+ 41020..41606 [1] |
+ No RBS annotated, None found |
+ tgtatattga |
+
+
+ | gene |
+ CPT-T1_058.gene |
+ 41606..41945 [1] |
+ Unannotated but valid RBS |
+ aaaatgt GGA |
+
+
+ | gene |
+ CPT-T1_059.gene |
+ 42013..42254 [1] |
+ Unannotated but valid RBS |
+ ac GAG atacc |
+
+
+ | gene |
+ CPT-T1_060.gene |
+ 42247..42487 [1] |
+ No RBS annotated, None found |
+ tcagcacttt |
+
+
+ | gene |
+ CPT-T1_061.gene |
+ 42473..42594 [1] |
+ No RBS annotated, None found |
+ ggcgacaaca |
+
+
+ | gene |
+ CPT-T1_062.gene |
+ 42580..42764 [1] |
+ No RBS annotated, None found |
+ ttgctactgt |
+
+
+ | gene |
+ CPT-T1_063.gene |
+ 42824..43324 [1] |
+ No RBS annotated, None found |
+ ccatcgacaa |
+
+
+ | gene |
+ CPT-T1_064.gene |
+ 43385..43881 [1] |
+ No RBS annotated, None found |
+ acatcaacca |
+
+
+ | gene |
+ CPT-T1_065.gene |
+ 43870..44030 [1] |
+ No RBS annotated, None found |
+ actcttgtcc |
+
+
+ | gene |
+ CPT-T1_066.gene |
+ 44098..44485 [1] |
+ No RBS annotated, None found |
+ gcacgacaac |
+
+
+ | gene |
+ CPT-T1_067.gene |
+ 44527..44649 [1] |
+ No RBS annotated, None found |
+ gttccccttt |
+
+
+ | gene |
+ CPT-T1_068.gene |
+ 44625..44852 [1] |
+ No RBS annotated, None found |
+ gacattatct |
+
+
+ | gene |
+ CPT-T1_069.gene |
+ 44898..45356 [1] |
+ No RBS annotated, None found |
+ cctacaccaa |
+
+
+ | gene |
+ CPT-T1_070.gene |
+ 45428..45970 [1] |
+ No RBS annotated, None found |
+ ttaagcaacc |
+
+
+ | gene |
+ CPT-T1_071.gene |
+ 45941..46451 [1] |
+ No RBS annotated, None found |
+ gcctgtgata |
+
+
+ | gene |
+ CPT-T1_072.gene |
+ 46444..47026 [1] |
+ No RBS annotated, None found |
+ tcacaaaaag |
+
+
+ | gene |
+ CPT-T1_073.gene |
+ 47085..47308 [1] |
+ No RBS annotated, None found |
+ ctcatcgaca |
+
+
+ | gene |
+ CPT-T1_074.gene |
+ 47293..47649 [1] |
+ No RBS annotated, None found |
+ cccgatcgaa |
+
+
+ | gene |
+ CPT-T1_075.gene |
+ 47637..47879 [1] |
+ No RBS annotated, None found |
+ aagtttttca |
+
+
+ | gene |
+ CPT-T1_076.gene |
+ 47864..48082 [1] |
+ Unannotated but valid RBS |
+ cgt AGGT aac |
+
+
+ | gene |
+ CPT-T1_077.gene |
+ 48235..48408 [1] |
+ No RBS annotated, None found |
+ gcggcaacaa |
+
+
+ | gene |
+ CPT-T1_078.gene |
+ 48377..48574 [1] |
+ No RBS annotated, None found |
+ ctggcgttat |
+
+
+ | gene |
+ CPT-T1_079.gene |
+ 48564..48803 [1] |
+ No RBS annotated, None found |
+ atattgaatt |
+
+
+
+
+
+
Start Codon Usage
+
This section covers genes with unusual start codons
+
+
+
+
+ | Start Codon |
+ Count |
+
+
+
+ | ATG | 74 |
+ | GTG | 3 |
+ | TTG | 2 |
+
+
+
+
+
+
+
+
+ | Feature Type |
+ ID |
+ Location |
+ Error |
+
+
+
+
+
+
+
+
Intergenic Gaps
+
Phage genomes are under pressure to maintain high coding density. Large intergenic gaps may be a sign of incorrect gene starts or missing genes.
+
+
+
+
+ | Region |
+ Size |
+ Bounding Gene Transcription Direction |
+ Message |
+
+
+
+
+ | 2184 .. 2238 |
+ 54 |
+ → → |
+
+ |
+
+
+ | 5876 .. 5937 |
+ 61 |
+ → → |
+
+ |
+
+
+ | 6705 .. 6795 |
+ 90 |
+ → → |
+
+ |
+
+
+ | 10414 .. 10527 |
+ 113 |
+ → → |
+
+ |
+
+
+ | 14432 .. 14511 |
+ 79 |
+ → → |
+
+ 2 ORFs found in this region
+ |
+
+
+ | 16621 .. 16698 |
+ 77 |
+ → → |
+
+ |
+
+
+ | 21506 .. 21765 |
+ 259 |
+ → → |
+
+ |
+
+
+ | 21906 .. 22033 |
+ 127 |
+ → → |
+
+ |
+
+
+ | 24289 .. 24350 |
+ 61 |
+ → → |
+
+ 1 ORFs found in this region
+ |
+
+
+ | 27573 .. 27634 |
+ 61 |
+ → → |
+
+ |
+
+
+ | 28087 .. 28180 |
+ 93 |
+ → → |
+
+ |
+
+
+ | 30612 .. 30679 |
+ 67 |
+ → → |
+
+ 1 ORFs found in this region
+ |
+
+
+ | 31641 .. 31707 |
+ 66 |
+ → → |
+
+ |
+
+
+ | 32399 .. 32477 |
+ 78 |
+ → → |
+
+ 1 ORFs found in this region
+ |
+
+
+ | 33611 .. 33683 |
+ 72 |
+ → → |
+
+ |
+
+
+ | 34166 .. 34238 |
+ 72 |
+ → → |
+
+ |
+
+
+ | 34415 .. 34533 |
+ 118 |
+ → → |
+
+ |
+
+
+ | 35638 .. 35760 |
+ 122 |
+ → → |
+
+ 2 ORFs found in this region
+ |
+
+
+ | 38693 .. 38773 |
+ 80 |
+ → → |
+
+ 2 ORFs found in this region
+ |
+
+
+ | 39211 .. 39288 |
+ 77 |
+ → → |
+
+ |
+
+
+ | 40330 .. 41030 |
+ 700 |
+ → → |
+
+ 7 ORFs found in this region
+ |
+
+
+ | 41945 .. 42023 |
+ 78 |
+ → → |
+
+ |
+
+
+ | 42764 .. 42835 |
+ 71 |
+ → → |
+
+ |
+
+
+ | 43324 .. 43395 |
+ 71 |
+ → → |
+
+ |
+
+
+ | 44030 .. 44110 |
+ 80 |
+ → → |
+
+ |
+
+
+ | 44485 .. 44538 |
+ 53 |
+ → → |
+
+ |
+
+
+ | 44852 .. 44909 |
+ 57 |
+ → → |
+
+ |
+
+
+ | 45356 .. 45439 |
+ 83 |
+ → → |
+
+ 6 ORFs found in this region
+ |
+
+
+ | 47026 .. 47098 |
+ 72 |
+ → → |
+
+ |
+
+
+ | 48082 .. 48246 |
+ 164 |
+ → → |
+
+ 2 ORFs found in this region
+ |
+
+
+
+
+
+
Overlapping Genes
+
Large gene overlaps may indicate an incorrect gene start or miscalled gene.
+
+
+
+
+ | Feature A |
+ Feature B |
+ Shared Region |
+ Overlap Length |
+
+
+
+
+ | CPT-T1_015.gene ([10516:10845](+)) |
+ CPT-T1_016.gene ([10516:11162](+)) |
+ 10527..10844 |
+ 317 bp |
+
+
+ | CPT-T1_030.gene ([24351:26550](-)) |
+ CPT-T1_031.gene ([26447:26600](-)) |
+ 26446..26536 |
+ 90 bp |
+
+
+ | CPT-T1_039.gene ([31881:32167](+)) |
+ CPT-T1_040.gene ([32099:32399](+)) |
+ 32111..32166 |
+ 55 bp |
+
+
+
+
+
+
+
+
Possible Morons 78 / 79 (Doesn't count towards score)
+
+
+
+
+ | Feature |
+ RBS |
+ Surrounding Features |
+
+
+
+
+ | CPT-T1_025.gene |
+ No RBS Available |
+
+ →
+ →
+ ←
+ →
+ →
+
+ |
+
+
+
+
+
+
+
+
Missing Product Tags 79 / 79
+
+
+
+
+ | Feature |
+ Qualifiers |
+
+
+
+
+
+
+
+
+
+
+
+
+
+