#!/usr/bin/env python3
"""Verify supplied RootGuard files without contacting Kernel or changing DNA.
Uses the Python standard library. Run from any directory; input files must be
next to this script. Prints JSON; does not overwrite archived QC files.
"""
from pathlib import Path
import hashlib, json, re
from itertools import groupby
WEB=Path(__file__).resolve().parent

def fasta(path):
    return ''.join(x.strip() for x in path.read_text().splitlines() if not x.startswith('>')).upper()

gb_text = (WEB / 'RootGuard_pHT01_Pveg_TasA_v2_Kernel_export_2026-10-01.gb').read_text()
v2 = re.sub('[^a-zA-Z]', '', gb_text.split('\nORIGIN\n', 1)[1].split('//', 1)[0]).upper()
native = fasta(WEB / 'RootGuard_TasA_CDS_786bp.fasta')
candidate = fasta(WEB / 'RootGuard_TasA_Bsub_DNAChisel_candidate_v1_2026-10-01.fasta')
protein = fasta(WEB / 'RootGuard_TasA_protein_261aa.faa')
backbone = fasta(WEB / 'RootGuard_pHT01_Addgene_reference_7956bp.fasta')
cassette = fasta(WEB / 'RootGuard_Pveg_RBS_TasA_Term_v1.fasta')

# Standard internal codon assignments, also shared by bacterial code 11;
# the experimentally untested sequence explicitly retains the ATG start.
bases = 'TCAG'
aa = ('FFLLSSSSYY**CC*W' 'LLLLPPPPHHQQRRRR'
      'IIIMTTTTNNKKSSRR' 'VVVVAAAADDEEGGGG')
codons = [x+y+z for x in bases for y in bases for z in bases]
code = dict(zip(codons, aa))
translate = lambda seq: ''.join(code[seq[i:i+3]] for i in range(0, len(seq), 3))
parts = [
 ('Pveg_BBa_K143012', 'promoter', 2684, 2780, 'https://registry.igem.org/parts/BBa_K143012'),
 ('RBS_Candidate_BBa_K090505', 'ribosome-binding-site', 2781, 2791, 'https://registry.igem.org/parts/BBa_K090505'),
 ('Spacer_RBS_tasA_5nt_v1', '5-utr', 2792, 2796, 'Local designed spacer; biological spacing unvalidated'),
 ('TasA_BBa_K1701005_local', 'cds-yes-start-yes-stop', 2797, 3582, 'https://www.ncbi.nlm.nih.gov/nuccore/NZ_CP020102.1; WP_004398632.1; https://registry.igem.org/parts/BBa_K1701005'),
 ('Terminator_BBa_K780000_local', 'terminator', 3583, 3636, 'https://registry.igem.org/parts/BBa_K780000'),
]
for seq in [native, candidate, v2, cassette, backbone]:
    assert set(seq) <= set('ACGT')
checks = {
 'v2_length_8734_bp': len(v2) == 8734,
 'v2_DNA_hash_matches_previous_QC': hashlib.sha256(v2.encode()).hexdigest() == '23c372fad614d39adb670cf5291f002524fd956d518e39186b48047c7d81eba1',
 'v2_is_circular_in_GenBank_metadata': 'circular' in gb_text.splitlines()[0],
 'native_length_786_bp': len(native) == 786,
 'candidate_length_786_bp': len(candidate) == 786,
 'native_v2_CDS_exact_match': v2[2796:3582] == native,
 'v2_cassette_matches_v1': v2[2683:3636] == cassette,
 'backbone_length_7956_bp': len(backbone) == 7956,
 'replacement_arithmetic_7956_minus_175_plus_953': len(v2) == len(backbone)-175+len(cassette),
 'unchanged_backbone_left_2683_bp': v2[:2683] == backbone[:2683],
 'unchanged_backbone_right_5098_bp': v2[3636:] == backbone[2858:],
 'native_start_ATG_stop_TAA': native[:3] == 'ATG' and native[-3:] == 'TAA',
 'candidate_start_ATG_stop_TAA': candidate[:3] == 'ATG' and candidate[-3:] == 'TAA',
 'native_no_internal_stops': '*' not in translate(native)[:-1],
 'candidate_no_internal_stops': '*' not in translate(candidate)[:-1],
 'native_translation_matches_261aa_reference': translate(native) == protein+'*' and len(protein)==261,
 'candidate_translation_matches_261aa_reference': translate(candidate) == protein+'*',
 '52_nucleotide_substitutions': sum(a != b for a,b in zip(native,candidate)) == 52,
 '47_synonymous_codons_changed': sum(native[i:i+3] != candidate[i:i+3] for i in range(0,786,3)) == 47,
 'RBS_spacer_ATG_junction': v2[2780:2799] == 'AAAGGAGGTGTAATTAATG',
}
for name, typ, start, end, source in parts:
    checks['typed_forward_annotation_'+typ] = bool(re.search(
        rf'(?m)^     \S+\s+{start}\.\.{end}\s*\n(?:(?!^     \S).)*?/kernel_part_type="{re.escape(typ)}"',
        gb_text, re.S | re.M))
assert all(checks.values()), {k:v for k,v in checks.items() if not v}

def metrics(seq):
    return {'length_bp': len(seq), 'GC_percent': 100*sum(seq.count(b) for b in 'GC')/len(seq),
            'longest_homopolymer_nt': max(len(list(g)) for _,g in groupby(seq)),
            'DNA_sha256': hashlib.sha256(seq.encode()).hexdigest()}

result = {
 'reverification_date_WIB': '2026-10-02',
 'execution': 'Read-only local file comparisons; no new Kernel tool run and no new biological data',
 'baseline': metrics(v2), 'native_CDS': metrics(native), 'external_candidate_CDS': metrics(candidate),
 'protein_length_without_terminal_stop_aa': len(protein),
 'checks': checks, 'checks_passed': sum(checks.values()), 'checks_total': len(checks),
 'nucleotide_identity_percent_native_vs_candidate': 100*734/786,
 'native_CDS_EcoRI_sites_1based_inclusive': [[m.start()+1,m.end()] for m in re.finditer('GAATTC', native)],
 'candidate_Kernel_full_reexport_hash_verified': False,
 'screening_input_full_hash_or_independent_count_verified': False,
 'Kernel_host_specific_Bacillus_optimization_completed': False,
 'instructor_acceptance_external_alternative_confirmed': False,
 'biological_function_verified': False,
 'source_generation_QC_kept_unchanged': True,
 'historical_QC_note': 'October 1 generation/QC pending fields are time-specific snapshots; use this record and report for later import, screening, schematic and search evidence.',
}

print(json.dumps(result, indent=2))