#!/usr/bin/env python3
"""
07_blast_primers.py  ->  primer_specificity.json

Check the 10 junction colony-PCR primers for cross-reactivity with the CHO-K1
PICRH-1.0 genome (GCF_003668045.3), so they stay construct-specific during later
stable-line genotyping. Two complementary tests:

  (1) in-silico PCR (isPcr) of each pair vs the host genome, run at a SENSITIVE
      setting (minPerfect=8) -> would any pair amplify a host locus?  This is the
      decisive test (requires a real 3' anchor + convergent pair within 4 kb).
  (2) blastn-short with btop -> the longest PERFECT 3'-anchored match of each
      single primer in the host (priming-competence proxy). 'covers 3' end' alone
      is NOT a risk; a long perfect 3' run is.
"""
import json, os, subprocess, tempfile, re
HERE = os.path.dirname(os.path.abspath(__file__))
DB_FA = os.path.expanduser("~/genomes/cho_PICRH1.0_GCF_003668045.3.fna")
DB_BLAST = os.path.expanduser("~/genomes/cho_PICRH1.0")
ISPCR = os.path.expanduser("~/genomes/isPcr")
primers = json.load(open(os.path.join(HERE, "primers.json")))

pairs = {j: (d["forward"]["seq"], d["reverse"]["seq"]) for j, d in primers.items()}
P = {f"{j}__F": d["forward"]["seq"] for j, d in primers.items()}
P.update({f"{j}__R": d["reverse"]["seq"] for j, d in primers.items()})

# ---- (1) isPcr: pairs vs CHO genome (sensitive: minPerfect=8) ----------------
q = tempfile.NamedTemporaryFile("w", suffix=".txt", delete=False)
for j, (f, r) in pairs.items():
    tag = j.replace('|', '_').replace(' ', '_')
    q.write(f"{tag}\t{f}\t{r}\n")            # the designed pair
    q.write(f"{tag}_FF\t{f}\t{f}\n")          # forward self-pair (single-primer convergence)
    q.write(f"{tag}_RR\t{r}\t{r}\n")          # reverse self-pair
q.close()
isp = subprocess.run([ISPCR, DB_FA, q.name, "stdout", "-maxSize=4000", "-minPerfect=8", "-minGood=12"],
                     capture_output=True, text=True, timeout=1200)
ispcr_products = [l for l in isp.stdout.splitlines() if l.startswith(">")]
print(f"isPcr (minPerfect=8, maxSize=4000; 5 designed pairs + 10 self-pairs) host products: {len(ispcr_products)}")
for h in ispcr_products: print("   ", h)

# ---- (2) blastn-short with btop: longest perfect 3' run per primer -----------
fa = tempfile.NamedTemporaryFile("w", suffix=".fa", delete=False)
for n, s in P.items(): fa.write(f">{n}\n{s}\n")
fa.close()
cols = "qseqid pident length qstart qend qlen sstrand btop"
out = subprocess.run(["blastn", "-task", "blastn-short", "-db", DB_BLAST, "-query", fa.name,
                      "-word_size", "7", "-evalue", "1000", "-dust", "no", "-max_target_seqs", "200",
                      "-outfmt", f"6 {cols}"], capture_output=True, text=True, timeout=600)

def perfect_3prime_run(btop, qstart, qend, qlen, strand):
    """length of the contiguous exact-match run ending at the primer's 3' base."""
    if qend < qlen:        # alignment does not include the 3'-terminal base
        return 0
    # btop: numbers = run of matches; letter pairs = mismatch/indel. Walk to find
    # the match run that includes the 3' end. For + strand the 3' end is the last
    # aligned query base, which corresponds to the END of the btop walk.
    toks = re.findall(r"\d+|[A-Z-]{2}", btop)
    # trailing matches in btop = matches at the high-query-coordinate end
    run = 0
    for t in reversed(toks):
        if t.isdigit():
            run += int(t);
            # a mismatch after (i.e., before in 5' direction) breaks it
        else:
            break
    return run

best3p = {n: 0 for n in P}
for line in out.stdout.strip().split("\n"):
    if not line: continue
    qid, pid, ln, qs, qe, ql, strand, btop = line.split("\t")
    qs, qe, ql = int(qs), int(qe), int(ql)
    r = perfect_3prime_run(btop, qs, qe, ql, strand)
    if r > best3p[qid]: best3p[qid] = r

report = {"db": "CHO-K1 PICRH-1.0 (GCF_003668045.3)",
          "ispcr": {"settings": "minPerfect=8 minGood=12 maxSize=4000; 5 designed pairs + 10 self-pairs",
                    "n_host_products": len(ispcr_products), "products": ispcr_products},
          "per_primer_longest_perfect_3prime_run_nt": best3p,
          "note_3prime_runs": ("12-17 nt perfect runs are within chance expectation for a 2.4 Gb genome "
                               "(a 15-mer occurs ~4x by chance) and yield no product without a convergent "
                               "partner; isPcr confirms no such pairing exists.")}
# decisive criterion: does ANY pair (designed or single-primer self-pair) amplify the host?
report["summary"] = {
    "any_host_amplicon_designed_or_self": len(ispcr_products) > 0,
    "max_perfect_3prime_run_nt": max(best3p.values()),
    "verdict": "PASS - construct-specific (no CHO-K1 host amplicon)" if len(ispcr_products) == 0 else "REVIEW",
}
json.dump(report, open(os.path.join(HERE, "primer_specificity.json"), "w"), indent=2)

print("\nLongest perfect 3'-anchored single-primer match in CHO genome (informational):")
for n in P:
    print(f"  {n:22} {best3p[n]:>2}/{len(P[n])} nt")
print("\nSUMMARY:", json.dumps(report["summary"]))
os.unlink(q.name); os.unlink(fa.name)
