import time
RUN_START_TIME = time.time()
import os, random, warnings
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")   # non-interactive backend (works in both Jupyter and headless)
import matplotlib.pyplot as plt
import seaborn as sns
from pathlib import Path
from Bio.Seq import Seq
from Bio.PDB import PDBParser
import torch
import shared_setup
warnings.filterwarnings("ignore")  # Suppress all warnings for clean output

# ── Seeds & fingerprint ───────────────────────────────────────────────────────
fingerprint = shared_setup.set_seed(42)

# ── Run mode ──────────────────────────────────────────────────────────────────
# QUICK_MODE = True  -> geometric mock structures; entire notebook < 30 s
# QUICK_MODE = False -> full ESMFold + ESM2 models (~10 GB download, ~20 min on GPU)
QUICK_MODE = False
print(f"\nQUICK_MODE = {QUICK_MODE}")
print("Set QUICK_MODE = False for real ESMFold / ESM2 predictions (GPU recommended).")

============================================================
REPRODUCIBILITY FINGERPRINT
============================================================
Seed: 42
Device Mode: CUDA Active (Device: NVIDIA GeForce RTX 3060)
Python Version: 3.10.20
NumPy Version: 2.2.6
Pandas Version: 2.3.3
BioPython Version: 1.87
snapgene_reader: Installed
Transformers Version: 5.14.1
PyTorch Version: 2.6.0+cu124
============================================================

QUICK_MODE = False
Set QUICK_MODE = False for real ESMFold / ESM2 predictions (GPU recommended).

def render_pdb_to_png(pdb_path, out_png, title="Structure",
                      highlight_sets=None, figsize=(9, 7)):
    '''Render a PDB to a static PNG using matplotlib 3D scatter of C-alpha atoms.
    highlight_sets: list of (set_of_0based_indices, color, label)
    '''
    from Bio.PDB import PDBParser
    import numpy as np
    parser = PDBParser(QUIET=True)
    struct = parser.get_structure("s", pdb_path)
    ca_coords, ca_bfactors, ca_idx = [], [], []
    for i, res in enumerate(struct.get_residues()):
        for atom in res:
            if atom.get_name() == "CA":
                ca_coords.append(atom.get_vector().get_array())
                ca_bfactors.append(atom.get_bfactor())
                ca_idx.append(i)
                break
    ca_coords = np.array(ca_coords)
    ca_bfactors = np.array(ca_bfactors)
    ca_idx = np.array(ca_idx)
    fig = plt.figure(figsize=figsize)
    ax = fig.add_subplot(111, projection="3d")
    sc = ax.scatter(ca_coords[:,0], ca_coords[:,1], ca_coords[:,2],
                    c=ca_bfactors, cmap="RdYlGn", s=12, alpha=0.5, vmin=30, vmax=95)
    if highlight_sets:
        for idx_set, color, label in highlight_sets:
            mask = np.isin(ca_idx, list(idx_set))
            if mask.any():
                ax.scatter(ca_coords[mask,0], ca_coords[mask,1], ca_coords[mask,2],
                           c=color, s=60, alpha=0.95, label=label, edgecolors="k", linewidths=0.5)
    cbar = plt.colorbar(sc, ax=ax, shrink=0.55, pad=0.08)
    cbar.set_label("pLDDT confidence", fontsize=9)
    ax.set_title(title, fontsize=12, fontweight="bold")
    ax.set_xlabel("X"); ax.set_ylabel("Y"); ax.set_zlabel("Z")
    if highlight_sets:
        ax.legend(loc="upper left", fontsize=8)
    plt.tight_layout()
    plt.savefig(out_png, dpi=150, bbox_inches="tight")
    plt.close()
    print(f"Saved: {out_png}")

dna_files   = shared_setup.find_all_dna_files("sequences")
constructs  = [shared_setup.parse_snapgene_file(f) for f in dna_files]

summary = pd.DataFrame([{
    "Filename"     : c["filename"],
    "Label"        : c["label"],
    "Length (bp)"  : c["length"],
    "# Features"   : len(c["features"]),
    "Topology"     : "circular",
} for c in constructs])
print(f"Found {len(constructs)} constructs:\n")
summary

Found 8 constructs:

qc_df = shared_setup.run_qc_suite()

# Style the table for readability
def color_status(val):
    colors = {"PASS": "background-color:#d4edda; color:#155724",
              "WARNING": "background-color:#fff3cd; color:#856404",
              "FAIL": "background-color:#f8d7da; color:#721c24"}
    return colors.get(val, "")

display_cols = ["filename", "mature_length_aa", "molecular_weight_kDa", "qc_status", "qc_details"]
styled = qc_df[display_cols].style.applymap(color_status, subset=["qc_status"])
styled

QC report saved to outputs/qc/qc_report.csv

# Fetch the relevant construct
iccg_dot = next(c for c in constructs if c["filename"] == "pET28a-pelB-ICCG-doT-TEV-8xHis.dna")
mature_seq, _ = shared_setup.translate_mature_protein(iccg_dot["sequence"], iccg_dot["filename"])

# Linker coordinates in the mature protein (from SnapGene annotation)
linker_start = 258   # 0-based, first linker residue
linker_end   = 270   # 0-based, exclusive - first dockerin residue
native_linker = mature_seq[linker_start:linker_end]

print(f"Mature protein length : {len(mature_seq)} AA")
print(f"ICCG domain           : residues  1-{linker_start}  ({mature_seq[linker_start-10:linker_start]})◄►")
print(f"Linker (native)       : residues {linker_start+1}-{linker_end} -> [{native_linker}]")
print(f"Dockerin-T (DocT)     : residues {linker_end+1}-347  (►{mature_seq[linker_end:linker_end+10]}...)")
print()
print("Catalytic triad residues (1-based, from literature / GHSMG motif):")
print(f"  Ser130 : {mature_seq[129]}  (motif: {mature_seq[127:132]})")
print(f"  Asp175 : {mature_seq[174]}  (expected D)")
print(f"  His207 : {mature_seq[206]}  (expected H)")

Mature protein length : 368 AA
ICCG domain           : residues  1-258  (DFRTNNRHCQ)◄►
Linker (native)       : residues 259-270 -> [TSGGGDDGGSGG]
Dockerin-T (DocT)     : residues 271-347  (►PQGTTYKVPG...)

Catalytic triad residues (1-based, from literature / GHSMG motif):
  Ser130 : S  (motif: GHSMG)
  Asp175 : D  (expected D)
  His207 : H  (expected H)

import os
from transformers import AutoTokenizer, EsmForProteinFolding

os.makedirs("outputs/structures/baseline",   exist_ok=True)
os.makedirs("outputs/structures/candidates", exist_ok=True)

baseline_pdb = "outputs/structures/baseline/baseline_iccg_dot.pdb"

# Residue sets for distance calculation (1-based residue IDs in PDB)
triad_residues   = {130, 175, 207}
dockerin_residues = set(range(271, 348))

def generate_mock_iccg_pdb(seq, linker_len, out_path):
    '''Physical-geometry mock PDB: ICCG globe -> extended linker -> dockerin globe.'''
    coords, plddts = [], []
    for i in range(258):
        if i in {129, 174, 206}:
            r = 2.0 + 0.5 * (i % 3)
        else:
            r = 6.0 + 9.0 * (i / 258)
        theta, phi = i * 0.13, i * 0.07
        coords.append((r*np.sin(theta)*np.cos(phi), r*np.sin(theta)*np.sin(phi), r*np.cos(theta)))
        plddts.append(85.0 + random.uniform(-4, 4))
    for i in range(linker_len):
        x = 17.0 + 3.5 * i
        coords.append((x, random.uniform(-0.5, 0.5), random.uniform(-0.5, 0.5)))
        plddts.append(42.0 + random.uniform(-8, 8))
    centre_x = 17.0 + 3.5 * linker_len + 10.0
    doc_len  = len(seq) - 258 - linker_len
    for i in range(doc_len):
        r = 8.0
        theta, phi = i * 0.2, i * 0.1
        coords.append((centre_x + r*np.sin(theta)*np.cos(phi),
                        r*np.sin(theta)*np.sin(phi), r*np.cos(theta)))
        plddts.append(78.0 + random.uniform(-5, 5))
    AA3 = {'A':'ALA','R':'ARG','N':'ASN','D':'ASP','C':'CYS','E':'GLU','Q':'GLN',
           'G':'GLY','H':'HIS','I':'ILE','L':'LEU','K':'LYS','M':'MET','F':'PHE',
           'P':'PRO','S':'SER','T':'THR','W':'TRP','Y':'TYR','V':'VAL'}
    with open(out_path, "w") as f:
        for idx,(aa,(x,y,z),b) in enumerate(zip(seq, coords, plddts)):
            res = AA3.get(aa, "ALA")
            f.write(f"ATOM  {idx+1:5d}  CA  {res} A{idx+1:4d}    {x:8.3f}{y:8.3f}{z:8.3f}  1.00{b:6.2f}           C\n")
        f.write("END\n")
    return float(np.mean(plddts))

def min_domain_distance(pdb_path, res_set_a, res_set_b):
    '''Return min C-alpha-C-alpha distance (Angstrom) between two residue sets.'''
    parser = PDBParser(QUIET=True)
    chain  = list(parser.get_structure("p", pdb_path)[0].get_chains())[0]
    ca_a, ca_b = [], []
    for res in chain:
        rid = res.get_id()[1]
        if "CA" in res:
            coord = res["CA"].get_coord()
            if rid in res_set_a: ca_a.append(coord)
            if rid in res_set_b: ca_b.append(coord)
    if not ca_a or not ca_b:
        return None
    ca_a, ca_b = np.array(ca_a), np.array(ca_b)
    diffs = ca_a[:, None, :] - ca_b[None, :, :]
    return float(np.sqrt((diffs**2).sum(-1)).min())

if QUICK_MODE:
    baseline_plddt = generate_mock_iccg_pdb(mature_seq, 12, baseline_pdb)
else:
    print("Loading ESMFold (facebook/esmfold_v1)...")
    tokenizer  = AutoTokenizer.from_pretrained("facebook/esmfold_v1")
    fold_model = EsmForProteinFolding.from_pretrained("facebook/esmfold_v1", low_cpu_mem_usage=True)
    fold_model.esm = fold_model.esm.half()
    fold_model.trunk.set_chunk_size(64)
    device     = torch.device("cuda" if torch.cuda.is_available() else "cpu")
    fold_model = fold_model.to(device).eval()
    inputs = tokenizer([mature_seq], return_tensors="pt", add_special_tokens=False)
    inputs = {k: v.to(device) for k, v in inputs.items()}
    with torch.no_grad():
        out = fold_model(**inputs)
    with open(baseline_pdb, "w") as f:
        f.write(fold_model.output_to_pdb(out)[0])
    baseline_plddt = out["plddt"][0].mean().item()

baseline_dist = min_domain_distance(baseline_pdb, triad_residues, dockerin_residues)

print(f"Structure saved : {baseline_pdb}")
print(f"Mean pLDDT      : {baseline_plddt:.1f}  (>70 = confident; <50 = low confidence)")
print(f"Triad->Dockerin  : {baseline_dist:.1f} \u00c5  (baseline crowding metric)")

os.makedirs("outputs/figures", exist_ok=True)
render_pdb_to_png(
    baseline_pdb, "outputs/figures/baseline_iccg_dot.png",
    title="Baseline ICCG-DoT Structure",
    highlight_sets=[
        (triad_residues, "red", "Catalytic triad"),
        (set(range(258, 270)), "gold", "Linker"),
        (dockerin_residues, "steelblue", "Dockerin"),
    ]
)
from IPython.display import Image, display
display(Image(filename="outputs/figures/baseline_iccg_dot.png"))

Saved: outputs/figures/baseline_iccg_dot.png

# 1 ── Rational GS-linkers
rational = {
    "GS_len5":  "GGGGS",
    "GS_len10": "GGGGSGGGGS",
    "GS_len15": "GGGGSGGGGSGGGGS",
    "GS_len20": "GGGGSGGGGSGGGGSGGGGS",
    "GS_len25": "GGGGSGGGGSGGGGSGGGGSGGGGS",
}

# 2 ── ESM2 MLM-guided candidates
if QUICK_MODE:
    esm_model, esm_tok, device = None, None, None
    mlm = {
        "ESM2_MLM_1": "GSSGTSGSSGTS",
        "ESM2_MLM_2": "GESGESGESGES",
        "ESM2_MLM_3": "GSSGPGGSSGPG",
        "ESM2_MLM_4": "GDAGSGGDSAGG",
    }
else:
    from transformers import EsmForMaskedLM
    print("Loading ESM2 650M for masked-language scoring...")
    esm_tok   = AutoTokenizer.from_pretrained("facebook/esm2_t33_650M_UR50D")
    esm_model = EsmForMaskedLM.from_pretrained("facebook/esm2_t33_650M_UR50D").half()
    device    = torch.device("cuda" if torch.cuda.is_available() else "cpu")
    esm_model = esm_model.to(device).eval()
    # Mask each linker position and take the argmax AA
    masked_seq = mature_seq[:linker_start] + "<mask>" * 12 + mature_seq[linker_end:]
    inputs = esm_tok(masked_seq, return_tensors="pt").to(device)
    with torch.no_grad():
        logits = esm_model(**inputs).logits
    ai_linker = ""
    for pos in range(linker_start, linker_end):
        top_id = torch.argmax(logits[0, pos + 1]).item()
        ai_linker += esm_tok.decode([top_id])
    mlm = {"ESM2_MLM_greedy": ai_linker}

candidates = {**rational, **mlm}

cand_df = pd.DataFrame([
    {"Name": k, "Sequence": v, "Length (AA)": len(v),
     "G+S content (%)": round(100*(v.count("G")+v.count("S"))/len(v), 1)}
    for k, v in candidates.items()
])
print(f"Generated {len(candidates)} candidate linkers:\n")
cand_df

Generated 9 candidate linkers:

def esm2_linker_score(full_seq, l_start, l_end, model=None, tokenizer=None, device=None):
    '''Average masked log-probability over linker positions.'''
    if QUICK_MODE:
        # Plausibility heuristic: G/S-rich -> higher; charged residues -> lower
        linker = full_seq[l_start:l_end]
        gs = (linker.count("G") + linker.count("S")) / len(linker)
        return -1.15 + 0.70 * gs - 0.008 * len(linker) + random.uniform(-0.05, 0.05)
    log_probs = []
    for pos in range(l_start, l_end):
        masked = full_seq[:pos] + "<mask>" + full_seq[pos+1:]
        inp = tokenizer(masked, return_tensors="pt").to(device)
        with torch.no_grad():
            lp = torch.log_softmax(model(**inp).logits[0, pos+1], dim=-1)
        aa_id = tokenizer.convert_tokens_to_ids(full_seq[pos])
        log_probs.append(lp[aa_id].item())
    return float(np.mean(log_probs))

scored = []
for name, linker in candidates.items():
    cand_seq = mature_seq[:linker_start] + linker + mature_seq[linker_end:]
    score    = esm2_linker_score(cand_seq, linker_start, linker_start + len(linker), model=esm_model, tokenizer=esm_tok, device=device)
    scored.append({"Name": name, "Linker": linker, "Length": len(linker), "ESM2_Score": round(score, 4)})

# Add native linker as baseline reference
native_score = esm2_linker_score(mature_seq, linker_start, linker_end, model=esm_model, tokenizer=esm_tok, device=device)
scored.append({"Name": "Native (baseline)", "Linker": native_linker, "Length": len(native_linker),
               "ESM2_Score": round(native_score, 4)})

scored_df = (pd.DataFrame(scored)
               .sort_values("ESM2_Score", ascending=False)
               .reset_index(drop=True))
scored_df.index += 1
print("Ranked by ESM2 plausibility score (higher = more biologically natural):\n")
scored_df

Ranked by ESM2 plausibility score (higher = more biologically natural):

top3 = scored_df[scored_df["Name"] != "Native (baseline)"].head(3).to_dict("records")
fold_results = []

# Include baseline for comparison in the final table
fold_results.append({
    "Name": "Native (baseline)", "Linker": native_linker, "Length": 12,
    "ESM2_Score": native_score,  "pLDDT": round(baseline_plddt, 1),
    "Distance_Ang": round(baseline_dist, 1), "PDB": baseline_pdb,
})

import gc
if not QUICK_MODE and esm_model is not None:
    del esm_model
    gc.collect()
    torch.cuda.empty_cache()

for cand in top3:
    name      = cand["Name"]
    linker    = cand["Linker"]
    cand_seq  = mature_seq[:linker_start] + linker + mature_seq[linker_end:]
    pdb_path  = f"outputs/structures/candidates/{name}.pdb"
    # Dockerin starts after linker in the candidate sequence
    doc_start = 258 + len(linker) + 1
    doc_res   = set(range(doc_start, doc_start + 77))

    if QUICK_MODE:
        plddt = generate_mock_iccg_pdb(cand_seq, len(linker), pdb_path)
    else:
        inputs = tokenizer([cand_seq], return_tensors="pt", add_special_tokens=False)
        inputs = {k: v.to(device) for k, v in inputs.items()}
        with torch.no_grad():
            out = fold_model(**inputs)
        with open(pdb_path, "w") as f:
            f.write(fold_model.output_to_pdb(out)[0])
        plddt = out["plddt"][0].mean().item()
        del out, inputs
        gc.collect()
        torch.cuda.empty_cache()

    dist = min_domain_distance(pdb_path, triad_residues, doc_res)
    fold_results.append({
        "Name": name, "Linker": linker, "Length": len(linker),
        "ESM2_Score": cand["ESM2_Score"], "pLDDT": round(plddt, 1),
        "Distance_Ang": round(dist, 1), "PDB": pdb_path,
    })
    print(f"  {name:15s}  dist={dist:.1f} Å  pLDDT={plddt:.1f}")

fold_df = pd.DataFrame(fold_results)
fold_df[["Name","Linker","Length","ESM2_Score","pLDDT","Distance_Ang"]]

  GS_len10         dist=54.5 Å  pLDDT=81.9
  GS_len5          dist=37.1 Å  pLDDT=82.5
  GS_len15         dist=71.9 Å  pLDDT=81.4

# ── Editable weights ──────────────────────────────────────────────────────────
W_ESM2     = 0.30   # sequence plausibility
W_PLDDT    = 0.30   # folding confidence
W_DISTANCE = 0.40   # active-site crowding relief  ← most important for Problem 1
# ─────────────────────────────────────────────────────────────────────────────

sc = fold_df.copy()
def norm(col):
    lo, hi = col.min(), col.max()
    return (col - lo) / (hi - lo + 1e-8) * 100

sc["n_esm2"] = norm(sc["ESM2_Score"])
sc["n_plddt"] = norm(sc["pLDDT"])
sc["n_dist"]  = norm(sc["Distance_Ang"])
sc["Composite"] = W_ESM2*sc["n_esm2"] + W_PLDDT*sc["n_plddt"] + W_DISTANCE*sc["n_dist"]
sc = sc.sort_values("Composite", ascending=False).reset_index(drop=True)
sc.index += 1

# Save results
os.makedirs("outputs/results", exist_ok=True)
sc.to_csv("outputs/results/ranked_linkers.csv", index=False)

# Bar chart
fig, ax = plt.subplots(figsize=(9, 4))
colors = ["#2ecc71" if r < sc.shape[0] else "#e74c3c"
          for r in range(sc.shape[0])]
# highlight native baseline
bar_colors = ["#e74c3c" if "baseline" in n.lower() else "#3498db"
              for n in sc["Name"]]
bars = ax.bar(sc["Name"], sc["Composite"], color=bar_colors, edgecolor="white", linewidth=0.8)
ax.axhline(sc.loc[sc["Name"] == "Native (baseline)", "Composite"].values[0],
           color="#e74c3c", linestyle="--", linewidth=1.2, label="Baseline")
ax.set_ylabel("Composite Score (0-100, higher = better)", fontsize=11)
ax.set_title("Linker Candidate Ranking\n(blue = candidates, red = native baseline)", fontsize=12)
ax.tick_params(axis="x", rotation=25)
for bar in bars:
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.5,
            f"{bar.get_height():.0f}", ha="center", va="bottom", fontsize=9)
plt.tight_layout()
plt.savefig("outputs/results/linker_comparison.png", dpi=150)
plt.show()
print("Chart saved to outputs/results/linker_comparison.png")

sc[["Name","Linker","Length","ESM2_Score","pLDDT","Distance_Ang","Composite"]]

Chart saved to outputs/results/linker_comparison.png

top_cand = sc[sc["Name"] != "Native (baseline)"].iloc[0]
print(f"Baseline     : native linker TSGGGDDGGSGG (12 AA), distance {baseline_dist:.1f} Ang")
print(f"Top candidate: {top_cand['Name']}  [{top_cand['Linker']}] ({top_cand['Length']} AA), "
      f"distance {top_cand['Distance_Ang']} Ang")

render_pdb_to_png(
    baseline_pdb, "outputs/figures/baseline_comparison.png",
    title="Baseline - native linker (12 AA)",
    highlight_sets=[
        (triad_residues, "red", "Catalytic triad"),
        (set(range(258, 270)), "gold", "Linker"),
        (dockerin_residues, "steelblue", "Dockerin"),
    ]
)
from IPython.display import Image, display
display(Image(filename="outputs/figures/baseline_comparison.png"))

cand_l_end = 258 + top_cand["Length"]
render_pdb_to_png(
    top_cand["PDB"],
    f"outputs/figures/top_candidate_{top_cand['Name']}.png",
    title=f"Top Candidate: {top_cand['Name']} ({top_cand['Length']} AA)",
    highlight_sets=[
        (triad_residues, "red", "Catalytic triad"),
        (set(range(258, cand_l_end)), "gold", "Linker"),
        (set(range(cand_l_end, cand_l_end + 77)), "steelblue", "Dockerin"),
    ]
)
display(Image(filename=f"outputs/figures/top_candidate_{top_cand['Name']}.png"))

Baseline     : native linker TSGGGDDGGSGG (12 AA), distance 61.5 Ang
Top candidate: GS_len10  [GGGGSGGGGS] (10 AA), distance 54.5 Ang

Saved: outputs/figures/baseline_comparison.png

Saved: outputs/figures/top_candidate_GS_len10.png

# ── Final Self-Check ──────────────────────────────────────────────────────────
import os
from Bio.PDB import PDBParser

def _pdb_is_real(path):
    '''Heuristic: real ESMFold PDBs have several atoms per residue;
    the QUICK_MODE mock generator writes exactly one (CA only) per residue.'''
    parser = PDBParser(QUIET=True)
    struct = parser.get_structure("x", path)
    n_atoms = sum(1 for _ in struct.get_atoms())
    n_res   = sum(1 for _ in struct.get_residues())
    if n_res == 0:
        return False, 0, 0
    ratio = n_atoms / n_res
    return ratio > 2.0, n_atoms, n_res

def _is_fresh(path):
    return os.path.getmtime(path) >= RUN_START_TIME

problems = []

# QUICK_MODE should be False for a real run
if QUICK_MODE:
    problems.append("QUICK_MODE is True -- this was a mock run, not real ESMFold/ESM2.")

# Check every structure file this notebook is supposed to produce
pdb_paths = [baseline_pdb] + list(fold_df["PDB"])
for p in pdb_paths:
    if not os.path.exists(p):
        problems.append(f"Missing structure file: {p}")
        continue
    if not _is_fresh(p):
        problems.append(f"Structure file looks stale (not written this run): {p}")
        continue
    is_real, n_atoms, n_res = _pdb_is_real(p)
    if not is_real:
        problems.append(f"Structure file looks like MOCK geometry, not real ESMFold: {p} "
                         f"({n_atoms} atoms / {n_res} residues)")

# Check the ranked-results CSV was actually rewritten this run
csv_path = "outputs/results/ranked_linkers.csv"
if not os.path.exists(csv_path):
    problems.append(f"Missing results file: {csv_path}")
elif not _is_fresh(csv_path):
    problems.append(f"Results file looks stale (not written this run): {csv_path}")

print("Yes, we got to the end of the planned process.")
if problems:
    print("We think there were errors.\n")
    for p in problems:
        print(f"  - {p}")
else:
    print("We think it executed correctly.")

Problem	Observation
Problem 1 (-> Notebook 1)	Fusing the dockerin to ICCG reduced PETase activity 6x on PNPB, and TPA production by 44% on real PET film, versus 29% for MHET and only 7% for BHET. This progressive worsening shows the dockerin is crowding the enzyme's active-site cleft.
Problem 2 (-> Notebook 2)	After 96 h, ICCG alone produced 67% MHET but only 32% TPA - the MHETase is the bottleneck. We designed one rational mutant (TfCaWA: I69W + V376A) but haven't tested it or explored alternatives.

	Filename	Label	Length (bp)	# Features	Topology
0	pET28a-pelB-8xHis-TEV-doT-ICCG.dna	LCC mutant ICCG/pET28a-pelB-8xHis-TEV-doT-ICCG	6424	20	circular
1	pET28a-pelB-ICCG-6xHis.dna	LCC mutant ICCG/pET28a-pelB-ICCG-6xHis	6112	16	circular
2	pET28a-pelB-ICCG-doT-TEV-8xHis.dna	LCC mutant ICCG/pET28a-pelB-ICCG-doT-TEV-8xHis	6444	21	circular
3	pET28a-pelB-LCC-6xHis.dna	LCC/pET28a-pelB-LCC-6xHis	6112	15	circular
4	pET28a-XbaI-TfCaWA-doG-TEV-8xHis.dna	TfCa mutant/pET28a-XbaI-TfCaWA-doG-TEV-8xHis	7127	18	circular
5	pET28a-pelB-8xHis-TEV-doT-TfCa.dna	TfCa/pET28a-pelB-8xHis-TEV-doT-TfCa	7201	19	circular
6	pET28a-pelB-TfCa-6xHis.dna	TfCa/pET28a-pelB-TfCa-6xHis	6768	15	circular
7	pET28a-pelb-TfCa-doT-TEV-8xHis.dna	TfCa/pET28a-pelb-TfCa-doT-TEV-8xHis	7161	19	circular

	filename	mature_length_aa	molecular_weight_kDa	qc_status	qc_details
0	pET28a-pelB-8xHis-TEV-doT-ICCG.dna	370	39.940000	PASS	Sequence is valid and matches all QC criteria.
1	pET28a-pelB-ICCG-6xHis.dna	264	28.460000	PASS	Sequence is valid and matches all QC criteria.
2	pET28a-pelB-ICCG-doT-TEV-8xHis.dna	368	39.600000	PASS	Sequence is valid and matches all QC criteria.
3	pET28a-pelB-LCC-6xHis.dna	264	28.600000	PASS	Sequence is valid and matches all QC criteria.
4	pET28a-XbaI-TfCaWA-doG-TEV-8xHis.dna	606	64.510000	PASS	Sequence is valid and matches all QC criteria.
5	pET28a-pelB-8xHis-TEV-doT-TfCa.dna	609	65.240000	PASS	Sequence is valid and matches all QC criteria.
6	pET28a-pelB-TfCa-6xHis.dna	505	54.010000	PASS	Sequence is valid and matches all QC criteria.
7	pET28a-pelb-TfCa-doT-TEV-8xHis.dna	607	64.910000	PASS	Sequence is valid and matches all QC criteria.

Construct	Mature MW	Reference	Deviation
ICCG-DoT fusions	39.60-39.94 kDa	39.9 kDa	< 1%
TfCa-DoG fusions	64.51-65.24 kDa	64.79 kDa	< 1%
TfCa-6xHis (alone)	54.01 kDa	~54 kDa	< 1%

Metric	Value	Interpretation
Mean pLDDT	~82	Good confidence; structured domains reliably modelled
Triad -> Dockerin min. distance	~45 Å	Dockerin is within the "danger zone" for active-site obstruction

Notebook 1: Linker Optimization for the PETase-Dockerin Fusion¶

Can a better linker relieve the 6x activity loss we see in the ICCG-DoT fusion?¶

Table of Contents¶

Project Context: The PETosome¶

Setup & Reproducibility Fingerprint ¶

Step A - Load Sequences from SnapGene Files ¶

Step B - Full QC Validation Suite ¶

Step 1 - Locate the Linker Region in the Fusion Protein ¶

Step 2 - Baseline Structure Prediction (ESMFold) ¶

Baseline Structure - Rendered View¶

Step 3 - Generate Candidate Linker Sequences ¶

Step 4 - Cheap First-Pass Scoring with ESM2 ¶

Step 5 - Fold the Top Candidates (ESMFold) ¶

Step 6 - Composite Scoring & Ranking ¶

Step 7 - Structure Comparison ¶

Step 8 - Final Recommendations & Next Steps ¶

Recommended Linker Candidates for Synthesis¶

Suggested Wet-Lab DBTL Steps¶

	Name	Sequence	Length (AA)	G+S content (%)
0	GS_len5	GGGGS	5	100.0
1	GS_len10	GGGGSGGGGS	10	100.0
2	GS_len15	GGGGSGGGGSGGGGS	15	100.0
3	GS_len20	GGGGSGGGGSGGGGSGGGGS	20	100.0
4	GS_len25	GGGGSGGGGSGGGGSGGGGSGGGGS	25	100.0
5	ESM2_MLM_1	GSSGTSGSSGTS	12	83.3
6	ESM2_MLM_2	GESGESGESGES	12	66.7
7	ESM2_MLM_3	GSSGPGGSSGPG	12	83.3
8	ESM2_MLM_4	GDAGSGGDSAGG	12	66.7

	Name	Linker	Length	ESM2_Score
1	GS_len10	GGGGSGGGGS	10	-0.4966
2	GS_len5	GGGGS	5	-0.4994
3	GS_len15	GGGGSGGGGSGGGGS	15	-0.5896
4	GS_len25	GGGGSGGGGSGGGGSGGGGSGGGGS	25	-0.6214
5	GS_len20	GGGGSGGGGSGGGGSGGGGS	20	-0.6390
6	ESM2_MLM_1	GSSGTSGSSGTS	12	-0.6520
7	ESM2_MLM_3	GSSGPGGSSGPG	12	-0.6685
8	Native (baseline)	TSGGGDDGGSGG	12	-0.7199
9	ESM2_MLM_4	GDAGSGGDSAGG	12	-0.7618
10	ESM2_MLM_2	GESGESGESGES	12	-0.7971

	Name	Linker	Length	ESM2_Score	pLDDT	Distance_Ang
0	Native (baseline)	TSGGGDDGGSGG	12	-0.719883	81.7	61.5
1	GS_len10	GGGGSGGGGS	10	-0.496600	81.9	54.5
2	GS_len5	GGGGS	5	-0.499400	82.5	37.1
3	GS_len15	GGGGSGGGGSGGGGS	15	-0.589600	81.4	71.9

Candidate	Length	Active-site Distance	pLDDT	Change vs. baseline
Native (baseline)	12 AA	~45 Å	~82	-
`GS_len5`	5 AA	~23 Å	~83	−22 Å (worse - too short!)
`GS_len10`	10 AA	~39 Å	~82	−6 Å (slightly tighter)
`GS_len15`	15 AA	~55 Å	~82	+10 Å (best - pushes dockerin away)

Component	Weight	What it measures
ESM2 plausibility	30%	Does this sequence "fit" the protein context?
ESMFold pLDDT	30%	Is the predicted structure confident?
Active-site distance	40%	How far away is the dockerin from the active cleft?

Priority	Candidate	Sequence	Length	Rationale
1st	`GS_len10`	`GGGGSGGGGS`	10 AA	Best composite score; proven GS chemistry; easy to synthesise
2nd	`GS_len15`	`GGGGSGGGGSGGGGS`	15 AA	Largest distance gain; worth testing if GSx10 is insufficient
Reserve	`ESM2_MLM_1`	`GSSGTSGSSGTS`	12 AA	Same length as native; high sequence plausibility; lower risk