--- /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! </p><p class="text-muted">Have you heard of the
+     <a href="https://cpt.tamu.edu">CPT</a>\'s Automated Phage Annotation
+     Pipeline?""",
+    ),
+    (
+        0,
+        """<b>MAJOR</b> issues detected! Please consider using the
+     <a href="https://cpt.tamu.edu">CPT</a>\'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)))