kbank.biosynthesis¶

class TargetType(StrEnum):
    """Structural target category in a bipartite glycan graph for a biosynthetic gene role."""

    NONE = "NONE"
    REDUCING_END = "REDUCING_END"
    LINKAGE = "LINKAGE"
    MODIFIED_MONO = "MODIFIED_MONO"
    MONOSACCHARIDE = "MONOSACCHARIDE"
    ANY = "ANY"

class AssignmentEvidence(IntEnum):
    """Provenance and evidence level for a gene-to-glycan component assignment."""

    RULE_INFERRED = 0
    GRAPH_ALIGNED = 1
    LITERATURE = 2
    LOCUS_COMPLETED = 3
    PAN_GRAPH_CONSENSUS = 4

@dataclass(slots=True, frozen=True)
class BiosyntheticRule:
    """Declarative rule defining what structural elements a `GeneRole` can synthesize."""

    role_idx: int
    role_name: str
    target_type: TargetType
    compatible_sugars: frozenset[int] = frozenset()
    target_modifiers: frozenset[int] = frozenset()
    matches_linkages: bool = False
    is_polymerase: bool = False

class BiosyntheticRuleRegistry:
    """Builds and caches biochemical compatibility rules from `GeneRole` and `Sugars`."""

    def __init__(self, sugars_db: type[Sugars] | Sugars = Sugars) -> None:
        """Initialize the registry and precompute rules for all `GeneRole` members."""
        self.sugars_db = sugars_db
        self.rules: dict[int, BiosyntheticRule] = {}
        self._build_rules()

    def _build_rules(self) -> None:
        names = self.sugars_db.names
        classes = self.sugars_db.classes
        name_to_idx = self.sugars_db.name_to_index

        for role in GeneRole:
            spec = _ROLE_GLYCAN_SPECS[role]
            mod_indices = frozenset(name_to_idx[m] for m in spec.target_modifiers if m in name_to_idx)

            compat_indices: set[int] = set()
            for s_idx in range(len(names)):
                s_name = names[s_idx]
                s_cls = classes[s_idx]

                if spec.disallowed_classes and s_cls in spec.disallowed_classes:
                    continue

                if spec.allowed_names or spec.allowed_classes:
                    if s_name in spec.allowed_names or s_cls in spec.allowed_classes:
                        compat_indices.add(s_idx)
                elif spec.disallowed_classes:
                    compat_indices.add(s_idx)

            self.rules[int(role)] = BiosyntheticRule(
                role_idx=int(role),
                role_name=role.name,
                target_type=spec.target_type,
                compatible_sugars=frozenset(compat_indices),
                target_modifiers=mod_indices,
                matches_linkages=spec.matches_linkages,
                is_polymerase=(role is GeneRole.POLYMERASE),
            )

    def get_rule(self, role_idx: int | GeneRole) -> BiosyntheticRule:
        """Retrieve the biochemical rule for a given `GeneRole` index."""
        return self.rules.get(int(role_idx), self.rules[int(GeneRole.UNKNOWN)])

class PanGraphConsensusPrior:
    """Archetype- and cross-species pan-graph consensus model for `Init`, `Wzy`, and `GT` -> linkage resolution.

    Pools node-level homology (`GeneNode.cluster_id`) and edge-level syntenic transitions
    (`Loci.as_pan_graph()` and `Glycans.as_pan_graph()`) across paired records to resolve:
    1. `INITIATOR` homology `cluster_id` -> biological Und-PP initiating monosaccharide (`WbaP` -> `D-Galp`,
       `WcaJ` -> `D-Glcp`, `WecA` -> `D-GlcNAc`/`D-GalNAc`, `Itr`/`PglC` -> `D-QuiNAc`/`D-FucNAc`/`D-Gal`).
    2. `GLYCOSYLTRANSFERASE` & `POLYMERASE` homology `cluster_id` -> `(donor, anomer, acceptor_carbon, acceptor)`
       and retaining vs. inverting stereochemical mechanism.
    3. Consecutive syntenic GT transitions `(cluster_a -> cluster_b)` in the pan-locus network paired with
       consecutive glycosidic linkages `(e_1 -> e_2)` in the pan-glycan network.
    """

    def __init__(self) -> None:
        """Initialize an empty `PanGraphConsensusPrior`."""
        self.init_cluster_sugar: dict[int, dict[int, float]] = {}
        self.cluster_record_counts: dict[int, int] = {}
        self.cluster_sig_full: dict[int, dict[tuple[int, ...], float]] = {}
        self.cluster_sig_dl: dict[int, dict[tuple[int, ...], float]] = {}
        self.cluster_sig_da: dict[int, dict[tuple[int, ...], float]] = {}
        self.cluster_mech: dict[int, dict[int, float]] = {}
        self.syntenic_pair_dl: dict[tuple[int, int], dict[tuple[tuple[int, ...], tuple[int, ...]], float]] = {}
        self.init_adj_gt_step1: dict[tuple[int, int], dict[tuple[int, ...], float]] = {}
        self.precursor_gt_donor: dict[tuple[int, int], dict[int, float]] = {}
        self.precursor_cluster_sugar: dict[int, dict[int, float]] = {}
        self.syntenic_pair_record_counts: dict[tuple[int, int], int] = {}
        self.pan_locus_edge_counts: dict[tuple[int, int], int] = {}
        self.pan_glycan_node_counts: dict[int, int] = {}
        self.pan_glycan_bigram_counts: dict[tuple[tuple[int, ...], tuple[int, ...]], int] = {}
        self.modifier_cluster_spec: dict[int, dict[tuple[int, int, int], float]] = {}
        self.fine_to_med: dict[int, int] = {}
        self.rescued_loci: list[Locus] = []
        self.classifier_suite: Any | None = None

    def preferred_init_sugar(
        self,
        cluster_id: int,
        candidate_sugars: set[int],
        aa_seq: bytes | None = None,
    ) -> int | None:
        """Return the highest-scoring initiating sugar index in `candidate_sugars` for `cluster_id`, or `None`."""
        if (
            cluster_id > 0
            and cluster_id in self.init_cluster_sugar
            and self.cluster_record_counts.get(cluster_id, 0) >= 2
        ):
            s_map = self.init_cluster_sugar[cluster_id]
            best_s: int | None = None
            best_val = 0.0
            for s_idx in candidate_sugars:
                val = s_map.get(s_idx, 0.0)
                if val > best_val:
                    best_val = val
                    best_s = s_idx
            if best_s is not None and best_val >= 0.5:
                return best_s
        if self.classifier_suite is not None:
            clf_s, clf_p = self.classifier_suite.preferred_init_sugar(cluster_id, candidate_sugars, aa_seq=aa_seq)
            if clf_s is not None and clf_p >= 0.40:
                return clf_s
        if cluster_id > 0 and cluster_id in self.init_cluster_sugar:
            s_map = self.init_cluster_sugar[cluster_id]
            best_s = None
            best_val = 0.0
            for s_idx in candidate_sugars:
                val = s_map.get(s_idx, 0.0)
                if val > best_val:
                    best_val = val
                    best_s = s_idx
            return best_s if best_val >= 0.5 else None
        return None

    def preferred_precursor_sugar(self, cluster_id: int, candidate_sugars: set[int]) -> int | None:
        """Return the multi-record consensus sugar synthesized by precursor `cluster_id` among `candidate_sugars`."""
        for c_tier in (cluster_id, self.fine_to_med.get(cluster_id, 0)):
            if c_tier <= 0 or self.cluster_record_counts.get(c_tier, 0) < 2:
                continue
            s_map = self.precursor_cluster_sugar.get(c_tier, {})
            if not s_map:
                continue
            cand_w = {s_idx: s_map.get(s_idx, 0.0) for s_idx in candidate_sugars if s_map.get(s_idx, 0.0) > 0.0}
            if not cand_w:
                continue
            tot_cand = sum(cand_w.values())
            best_s, best_w = max(cand_w.items(), key=lambda kv: kv[1])
            if best_w >= 1.5 and best_w >= 0.55 * tot_cand:
                return best_s
        return None

    def has_init_consensus(self, cluster_id: int, sugar_idx: int, aa_seq: bytes | None = None) -> bool:
        """Return True if `cluster_id` has multi-record pan-graph or classifier consensus for initiating on `sugar_idx`."""
        if cluster_id > 0 and self.cluster_record_counts.get(cluster_id, 0) >= 2:
            s_map = self.init_cluster_sugar.get(cluster_id, {})
            tot = sum(s_map.values())
            if tot > 0 and s_map.get(sugar_idx, 0.0) >= 0.45 * tot:
                return True
        if self.classifier_suite is not None:
            clf_s, clf_p = self.classifier_suite.preferred_init_sugar(cluster_id, {sugar_idx}, aa_seq=aa_seq)
            if clf_s == sugar_idx and clf_p >= 0.65:
                return True
        return False

    def gt_edge_bonus(
        self,
        cluster_id: int,
        glycan: Glycan,
        edge_idx: int,
        adj_clusters: Sequence[tuple[int, int]] = (),
        adj_init_clusters: Sequence[int] = (),
        adj_precursor_clusters: Sequence[int] = (),
        adj_precursor_sugars: set[int] | frozenset[int] = frozenset(),
        aa_seq: bytes | None = None,
    ) -> tuple[float, bool]:
        """Compute cross-archetype pan-graph homology & synteny cost bonus and whether pan-consensus holds."""
        sig_full, sig_dl, sig_da, mech = _edge_signatures(glycan, edge_idx)
        bonus = 0.0
        has_consensus = False

        red_idx = glycan.reducing_end_index
        wzy_k, _, _ = glycan.polymer_edge_and_assembly_depths(red_idx)
        all_core = [
            k
            for k, e in enumerate(glycan.edges)
            if not glycan.nodes[e.source].is_modifier and not glycan.nodes[e.target].is_modifier
        ]
        intra_core = [k for k in all_core if k != wzy_k]
        cand_edges = intra_core if (intra_core and edge_idx != wzy_k) else all_core

        for c_tier, tier_scale in ((cluster_id, 1.0), (self.fine_to_med.get(cluster_id, 0), 0.72)):
            if c_tier <= 0 or self.cluster_record_counts.get(c_tier, 0) < 2:
                continue

            full_map = self.cluster_sig_full.get(c_tier, {})
            dl_map = self.cluster_sig_dl.get(c_tier, {})
            da_map = self.cluster_sig_da.get(c_tier, {})

            gly_full_w = {
                _edge_signatures(glycan, k)[0]: full_map.get(_edge_signatures(glycan, k)[0], 0.0) for k in cand_edges
            }
            gly_dl_w = {
                _edge_signatures(glycan, k)[1]: dl_map.get(_edge_signatures(glycan, k)[1], 0.0) for k in cand_edges
            }
            gly_da_w = {
                _edge_signatures(glycan, k)[2]: da_map.get(_edge_signatures(glycan, k)[2], 0.0) for k in cand_edges
            }

            w_full = gly_full_w.get(sig_full, 0.0)
            max_gly_full = max(gly_full_w.values(), default=0.0)
            other_full = max((v for k_s, v in gly_full_w.items() if k_s != sig_full), default=0.0)

            w_dl = gly_dl_w.get(sig_dl, 0.0)
            max_gly_dl = max(gly_dl_w.values(), default=0.0)
            other_dl = max((v for k_s, v in gly_dl_w.items() if k_s != sig_dl), default=0.0)

            w_da = gly_da_w.get(sig_da, 0.0)
            max_gly_da = max(gly_da_w.values(), default=0.0)
            other_da = max((v for k_s, v in gly_da_w.items() if k_s != sig_da), default=0.0)

            if w_full >= 0.60 and w_full == max_gly_full and w_full >= other_full * 1.20 and w_da >= max_gly_da * 0.65:
                bonus -= 0.13 * tier_scale
                has_consensus = True
            elif max_gly_full >= 0.85 and w_full < max_gly_full * 0.50:
                bonus += 0.09 * tier_scale

            if w_dl >= 0.60 and w_dl == max_gly_dl and w_dl >= other_dl * 1.20 and w_da >= max_gly_da * 0.65:
                bonus -= 0.11 * tier_scale
                has_consensus = True
            elif max_gly_dl >= 0.85 and w_dl < max_gly_dl * 0.50:
                bonus += 0.11 * tier_scale

            if w_da >= 0.65 and w_da == max_gly_da and w_da >= other_da * 1.20:
                bonus -= 0.08 * tier_scale
                has_consensus = True
            elif max_gly_da >= 0.90 and w_da < max_gly_da * 0.45:
                bonus += 0.16 * tier_scale

            mech_map = self.cluster_mech.get(c_tier, {})
            tot_mech = sum(mech_map.values())
            if mech > 0 and tot_mech >= 1.2:
                frac_mech = mech_map.get(mech, 0.0) / tot_mech
                if frac_mech >= 0.68:
                    bonus -= 0.03 * tier_scale * frac_mech
                elif frac_mech <= 0.18:
                    bonus += 0.12 * tier_scale

        e_curr = glycan.edges[edge_idx]
        donor_s = glycan.nodes[e_curr.source].sugar

        # Signal 1: (INITIATOR_cluster -> GT_cluster) step-1 anchoring onto the biological reducing end
        if e_curr.target == red_idx and cluster_id > 0:
            step1_edges = [k for k in cand_edges if glycan.edges[k].target == red_idx]
            for c_init in adj_init_clusters:
                for c_tier in (cluster_id, self.fine_to_med.get(cluster_id, 0)):
                    if c_tier > 0 and self.cluster_record_counts.get(c_tier, 0) >= 2:
                        step1_w = self.init_adj_gt_step1.get((c_init, c_tier), {}).get(sig_dl, 0.0)
                        if step1_w >= 0.85:
                            bonus -= 0.07
                            has_consensus = True
                            break
                if self.cluster_record_counts.get(c_init, 0) >= 2 and len(step1_edges) == 1:
                    bonus -= 0.055
                    has_consensus = True
                    break

        # Signal 2: (Precursor_Operon -> GT_cluster) syntenic donor-sugar coupling
        if cluster_id > 0:
            for c_prec in adj_precursor_clusters:
                for c_tier in (cluster_id, self.fine_to_med.get(cluster_id, 0)):
                    if c_tier <= 0 or self.cluster_record_counts.get(c_tier, 0) < 2:
                        continue
                    prec_map = self.precursor_gt_donor.get((c_prec, c_tier), {})
                    if prec_map:
                        max_prec_w = max(prec_map.values(), default=0.0)
                        w_prec = prec_map.get(donor_s, 0.0)
                        if w_prec >= 0.85 and w_prec >= max_prec_w * 0.80:
                            bonus -= 0.06
                            has_consensus = True
                            break

        if adj_precursor_sugars and donor_s in adj_precursor_sugars:
            bonus -= 0.055
            matching_cand_donors = {
                glycan.nodes[glycan.edges[k].source].sugar
                for k in cand_edges
                if glycan.nodes[glycan.edges[k].source].sugar in adj_precursor_sugars
            }
            matching_cand_sigs = {
                _edge_signatures(glycan, k)[1]
                for k in cand_edges
                if glycan.nodes[glycan.edges[k].source].sugar == donor_s
            }
            if len(matching_cand_donors) == 1 and len(matching_cand_sigs) == 1:
                has_consensus = True

        # Signal 3: Syntenic (GT_A -> GT_B) x (Linkage_1 -> Linkage_2) bigrams weighted by Pan-Glycan network
        if cluster_id > 0:
            for c_other, dir_flag in adj_clusters:
                if c_other <= 0:
                    continue
                for c_o_tier, c_s_tier in (
                    (c_other, cluster_id),
                    (self.fine_to_med.get(c_other, 0), self.fine_to_med.get(cluster_id, 0)),
                ):
                    if c_o_tier <= 0 or c_s_tier <= 0:
                        continue
                    pair_key = (c_o_tier, c_s_tier) if dir_flag < 0 else (c_s_tier, c_o_tier)
                    if self.syntenic_pair_record_counts.get(pair_key, 0) < 2:
                        continue
                    pair_map = self.syntenic_pair_dl.get(pair_key, {})
                    if not pair_map:
                        continue
                    for k_other in cand_edges:
                        if k_other == edge_idx:
                            continue
                        e_other = glycan.edges[k_other]
                        if e_curr.target == e_other.source or e_curr.source == e_other.target:
                            _, sig_dl_other, _, _ = _edge_signatures(glycan, k_other)
                            sig_pair = (sig_dl_other, sig_dl) if dir_flag < 0 else (sig_dl, sig_dl_other)
                            sup = pair_map.get(sig_pair, 0.0)
                            if sup >= 0.85:
                                bonus -= 0.06
                                has_consensus = True
                                break

                # Cross-network anchor: if neighbor GT c_other has multi-record consensus on an adjacent edge k_other
                # and (dl_other, dl_curr) is a recurring disaccharide motif in the global Pan-Glycan network
                if not has_consensus and self.cluster_record_counts.get(c_other, 0) >= 2:
                    o_dl_map = self.cluster_sig_dl.get(c_other, {})
                    if o_dl_map:
                        cand_o_w = {
                            _edge_signatures(glycan, k)[1]: o_dl_map.get(_edge_signatures(glycan, k)[1], 0.0)
                            for k in cand_edges
                        }
                        max_o_w = max(cand_o_w.values(), default=0.0)
                        if max_o_w >= 0.80:
                            for k_other in cand_edges:
                                if k_other == edge_idx:
                                    continue
                                _, dl_other, _, _ = _edge_signatures(glycan, k_other)
                                if cand_o_w.get(dl_other, 0.0) < max_o_w:
                                    continue
                                e_other = glycan.edges[k_other]
                                if e_curr.target == e_other.source or e_curr.source == e_other.target:
                                    sig_pair = (dl_other, sig_dl) if dir_flag < 0 else (sig_dl, dl_other)
                                    bg_cnt = self.pan_glycan_bigram_counts.get(sig_pair, 0)
                                    if bg_cnt >= 2:
                                        bonus -= 0.045
                                        adj_sigs = {
                                            _edge_signatures(glycan, k_c)[1]
                                            for k_c in cand_edges
                                            if k_c != k_other
                                            and (
                                                glycan.edges[k_c].target == e_other.source
                                                or glycan.edges[k_c].source == e_other.target
                                            )
                                        }
                                        if len(adj_sigs) == 1:
                                            has_consensus = True
                                        break

        # Hybrid fallback: for singleton/unresolved GT clusters or unclustered genes (cluster_id <= 0),
        # query the factorized 1D-CNN BiosyntheticClassifierSuite (donor, linkage, mechanism)
        if self.classifier_suite is not None and (
            not has_consensus or self.cluster_record_counts.get(cluster_id, 0) < 2
        ):
            clf_b, clf_cons = self.classifier_suite.gt_edge_score(
                cluster_id,
                glycan,
                edge_idx,
                aa_seq=aa_seq,
            )
            scale_clf = 1.0 if self.cluster_record_counts.get(cluster_id, 0) < 2 else 0.55
            bonus += scale_clf * clf_b
            if clf_cons:
                has_consensus = True

        return float(np.clip(bonus, -0.32, 0.45)), has_consensus

    @classmethod
    def fit(
        cls,
        glycans: Sequence[Glycan],
        loci: Loci | Sequence[Locus],
        archetypes: Sequence[LocusArchetype | int] | np.ndarray,
        is_repeating: Sequence[bool] | np.ndarray,
        *,
        assigner: BiosyntheticAssigner | None = None,
        em_iterations: int = 5,
        classifier_suite: Any | None = None,
    ) -> tuple[Self, list[Glycan]]:
        """Fit a `PanGraphConsensusPrior` across all paired records and rotate glycans to biological reducing ends."""
        from collections import Counter, defaultdict

        from scipy.optimize import linear_sum_assignment

        from kbank.biosynthesis.classifier import BiosyntheticClassifierSuite
        from kbank.glycan import Glycans

        prior = cls()
        suite = classifier_suite if classifier_suite is not None else BiosyntheticClassifierSuite.pretrained()
        if suite is not None and suite.is_fitted:
            prior.classifier_suite = suite
        n = len(glycans)
        if n == 0:
            return prior, list(glycans)

        loci_batch = loci if isinstance(loci, Loci) else Loci.from_list(list(loci))
        glycans_batch = Glycans.from_list(list(glycans))

        # 1. Extract pan-locus and pan-glycan consensus transition topologies via .as_pan_graph()
        if len(loci_batch.node_types) > 0:
            role_cluster_mask = int(GeneNode.mask("ROLE_ID", "CLUSTER_ID"))
            pan_l_adj, pan_l_nodes, _ = loci_batch.as_pan_graph(
                node_mask=role_cluster_mask,
                canonicalize_strands=False,
            )
            coo_l = pan_l_adj.tocoo()
            for r_i, c_i, cnt in zip(coo_l.row, coo_l.col, coo_l.data, strict=False):
                c_src = int(GeneNode(int(pan_l_nodes[int(r_i)])).cluster_id)
                c_tgt = int(GeneNode(int(pan_l_nodes[int(c_i)])).cluster_id)
                if c_src > 0 and c_tgt > 0:
                    prior.pan_locus_edge_counts[(c_src, c_tgt)] = prior.pan_locus_edge_counts.get(
                        (c_src, c_tgt), 0
                    ) + int(cnt)

        if len(glycans_batch.node_types) > 0:
            _, pan_g_nodes, pan_g_counts = glycans_batch.as_pan_graph()
            for g_n, g_c in zip(pan_g_nodes, pan_g_counts, strict=False):
                prior.pan_glycan_node_counts[int(g_n)] = int(g_c)

        # Build Medium/Coarse superfamily cluster hierarchy when protein sequences are available (by_role=False)
        cluster_rec_sets: dict[int, set[int]] = defaultdict(set)
        if (
            len(loci_batch.node_types) > 0
            and len(loci_batch.genes) == len(loci_batch.node_types)
            and bool(np.any(loci_batch.cluster_ids > 0))
        ):
            med_loci = loci_batch.cluster_genes(
                mode="mcl",
                min_identity=20.0,
                min_bitscore=35.0,
                inflation=1.6,
                by_role=False,
            )
            for fine_c, med_c in zip(loci_batch.cluster_ids, med_loci.cluster_ids, strict=False):
                fc = int(fine_c)
                mc = int(med_c)
                if fc > 0 and mc > 0:
                    prior.fine_to_med[fc] = mc + 100000

        # 2. Pass 0: Accumulate INITIATOR & Precursor cluster -> sugar specificity & rotate biological reducing ends
        canon_init_names = {
            b"GlcNAc",
            b"GalNAc",
            b"QuiNAc",
            b"FucNAc",
            b"FucN4N",
            b"QuiN4N",
            b"Bac",
            b"Gal",
            b"Glc",
        }
        canon_init_ids = {Sugars.name_to_index[nm] for nm in canon_init_names if nm in Sugars.name_to_index}
        housekeeping_names = {b"Glc", b"Gal", b"GlcNAc"}
        housekeeping_ids = {Sugars.name_to_index[nm] for nm in housekeeping_names if nm in Sugars.name_to_index}
        precursor_roles = {
            GeneRole.OXIDOREDUCTASE,
            GeneRole.EPIMERASE,
            GeneRole.AMINOTRANSFERASE,
            GeneRole.ALDOLASE,
            GeneRole.NUCLEOTIDASE,
            GeneRole.NUCLEOTIDYLTRANSFERASE,
        }
        helper = assigner if assigner is not None else BiosyntheticAssigner()
        init_sugar_acc: dict[int, dict[int, float]] = defaultdict(lambda: defaultdict(float))
        prec_sugar_acc: dict[int, dict[int, float]] = defaultdict(lambda: defaultdict(float))
        mod_spec_acc: dict[int, dict[tuple[int, int, int], float]] = defaultdict(lambda: defaultdict(float))
        mod_roles_set = {
            GeneRole.ACETYLTRANSFERASE,
            GeneRole.PYRUVYLTRANSFERASE,
            GeneRole.METHYLTRANSFERASE,
            GeneRole.PHOSPHOTRANSFERASE,
            GeneRole.FORMYLTRANSFERASE,
            GeneRole.ALANINE_TRANSFERASE,
        }

        raw_locus_list = [loci_batch[i] for i in range(n)]
        for i in range(n):
            loc = raw_locus_list[i]
            gly = glycans[i]
            if len(loc.nodes) == 0 or len(gly.nodes) == 0:
                continue
            for g_node in loc.nodes:
                cid = int(g_node.cluster_id)
                if cid > 0:
                    cluster_rec_sets[cid].add(i)
                    mc = prior.fine_to_med.get(cid, 0)
                    if mc > 0:
                        cluster_rec_sets[mc].add(i)

            n_ac_mods, _ = helper._find_n_acetyl_nodes(gly)
            for g_node in loc.nodes:
                if g_node.role in mod_roles_set and g_node.cluster_id > 0:
                    cid = int(g_node.cluster_id)
                    m_targets = helper.registry.get_rule(g_node.role).target_modifiers
                    for e_mod in gly.edges:
                        if 0 <= e_mod.source < len(gly.nodes) and 0 <= e_mod.target < len(gly.nodes):
                            s_m = gly.nodes[e_mod.source]
                            t_m = gly.nodes[e_mod.target]
                            if (
                                s_m.is_modifier
                                and e_mod.source not in n_ac_mods
                                and s_m.sugar in m_targets
                                and not t_m.is_modifier
                            ):
                                m_key = (
                                    int(s_m.sugar),
                                    int(t_m.sugar),
                                    int(e_mod.edge_type.chemical_int),
                                )
                                for t_cid in (cid, prior.fine_to_med.get(cid, 0)):
                                    if t_cid > 0:
                                        mod_spec_acc[t_cid][m_key] += 1.0

            core_sugars = {nd.sugar for nd in gly.nodes if not nd.is_modifier}
            for g_node in loc.nodes:
                if g_node.role in precursor_roles and g_node.cluster_id > 0:
                    cid = int(g_node.cluster_id)
                    compat = helper.registry.get_rule(g_node.role).compatible_sugars
                    cand_s = (core_sugars & compat) if compat else core_sugars
                    spec_s = (cand_s - housekeeping_ids) or cand_s
                    if spec_s:
                        w_p = 2.0 if len(spec_s) == 1 else (1.0 / len(spec_s))
                        for t_cid in (cid, prior.fine_to_med.get(cid, 0)):
                            if t_cid > 0:
                                for s_idx in spec_s:
                                    prec_sugar_acc[t_cid][int(s_idx)] += w_p

            cyc = gly.cycle_nodes()
            cand_nodes = cyc if cyc else tuple(idx for idx, nd in enumerate(gly.nodes) if not nd.is_modifier)
            cyc_sugars = {gly.nodes[idx].sugar for idx in cand_nodes}
            canon_in_cyc = cyc_sugars & canon_init_ids
            target_sugars = canon_in_cyc if canon_in_cyc else cyc_sugars
            if not target_sugars:
                continue

            w_init = 3.0 if len(target_sugars) == 1 else (1.0 / len(target_sugars))
            for g_node in loc.nodes:
                if g_node.role == GeneRole.INITIATOR and g_node.cluster_id > 0:
                    for s_idx in target_sugars:
                        init_sugar_acc[int(g_node.cluster_id)][int(s_idx)] += w_init

        prior.cluster_record_counts = {c: len(recs) for c, recs in cluster_rec_sets.items()}
        prior.init_cluster_sugar = {c: dict(cnts) for c, cnts in init_sugar_acc.items()}
        prior.precursor_cluster_sugar = {c: dict(cnts) for c, cnts in prec_sugar_acc.items()}
        prior.modifier_cluster_spec = {c: dict(cnts) for c, cnts in mod_spec_acc.items()}

        rotated_glycans: list[Glycan] = []
        locus_list: list[Locus] = []
        bigram_acc: Counter[tuple[tuple[int, ...], tuple[int, ...]]] = Counter()

        for i in range(n):
            gly = glycans[i]
            loc = raw_locus_list[i]
            arch = LocusArchetype(int(archetypes[i]))
            rep = bool(is_repeating[i])
            if len(gly.nodes) > 0 and rep:
                gly = helper.resolve_biological_reducing_end(
                    gly,
                    loc if len(loc.nodes) > 0 else None,
                    archetype=arch,
                    pan_prior=prior,
                )
            if len(loc.nodes) > 0 and len(gly.nodes) > 0:
                loc = helper.rescue_cassette_roles(
                    loc,
                    gly,
                    archetype=arch,
                    is_repeating=rep,
                    pan_cluster_counts=prior.cluster_record_counts,
                )
                for g_node in loc.nodes:
                    cid = int(g_node.cluster_id)
                    if cid > 0:
                        cluster_rec_sets[cid].add(i)
                        mc = prior.fine_to_med.get(cid, 0)
                        if mc > 0:
                            cluster_rec_sets[mc].add(i)

            rotated_glycans.append(gly)
            locus_list.append(loc)

            # Accumulate global Pan-Glycan consecutive linkage bigrams across ALL glycans (paired & unpaired)
            core_e = [
                k
                for k, e in enumerate(gly.edges)
                if not gly.nodes[e.source].is_modifier and not gly.nodes[e.target].is_modifier
            ]
            for k1 in core_e:
                e1 = gly.edges[k1]
                _, dl1, _, _ = _edge_signatures(gly, k1)
                for k2 in core_e:
                    if k1 == k2:
                        continue
                    e2 = gly.edges[k2]
                    if e1.source == e2.target or e1.target == e2.source:
                        _, dl2, _, _ = _edge_signatures(gly, k2)
                        bigram_acc[(dl1, dl2)] += 1

        prior.cluster_record_counts = {c: len(recs) for c, recs in cluster_rec_sets.items()}
        prior.pan_glycan_bigram_counts = dict(bigram_acc)
        prior.rescued_loci = locus_list
        if prior.classifier_suite is not None:
            prior.classifier_suite.precompute_for_loci(locus_list)

        # 3. Iterative EM across all paired records to refine GT cluster_id -> linkage & syntenic pair priors
        paired_indices = [i for i in range(n) if len(locus_list[i].nodes) > 0 and len(rotated_glycans[i].nodes) > 0]
        if not paired_indices:
            return prior, rotated_glycans

        cost_strat = BiosyntheticCostStrategy(rule_registry=helper.registry, pan_prior=prior)

        for _em_round in range(max(1, em_iterations)):
            sig_full_acc: dict[int, dict[tuple[int, ...], float]] = defaultdict(lambda: defaultdict(float))
            sig_dl_acc: dict[int, dict[tuple[int, ...], float]] = defaultdict(lambda: defaultdict(float))
            sig_da_acc: dict[int, dict[tuple[int, ...], float]] = defaultdict(lambda: defaultdict(float))
            mech_acc: dict[int, dict[int, float]] = defaultdict(lambda: defaultdict(float))
            syn_pair_acc: dict[tuple[int, int], dict[tuple[tuple[int, ...], tuple[int, ...]], float]] = defaultdict(
                lambda: defaultdict(float)
            )
            syn_pair_recs: dict[tuple[int, int], set[int]] = defaultdict(set)
            init_step1_acc: dict[tuple[int, int], dict[tuple[int, ...], float]] = defaultdict(
                lambda: defaultdict(float)
            )
            prec_donor_acc: dict[tuple[int, int], dict[int, float]] = defaultdict(lambda: defaultdict(float))

            for i in paired_indices:
                loc = locus_list[i]
                gly = rotated_glycans[i]
                n_res = len(gly.nodes)
                n_edges = len(gly.edges)
                if n_edges == 0:
                    continue

                link_genes = [
                    g_idx
                    for g_idx, g_node in enumerate(loc.nodes)
                    if g_node.role in (GeneRole.GLYCOSYLTRANSFERASE, GeneRole.POLYMERASE)
                ]
                core_edges = [
                    k
                    for k, e in enumerate(gly.edges)
                    if not gly.nodes[e.source].is_modifier and not gly.nodes[e.target].is_modifier
                ]
                if not link_genes or not core_edges:
                    continue

                _, loc_labels = loc.as_graph(as_bipartite=True)
                _, gly_labels = gly.as_graph(as_bipartite=True)
                c_mat = cost_strat.compute_cost(loc_labels, gly_labels, graph_a=loc, graph_b=gly)

                bip_cols = [n_res + k for k in core_edges]
                raw_sub = c_mat[np.ix_(link_genes, bip_cols)]
                sub_c = np.where(np.isfinite(raw_sub) & (raw_sub < 1e4), raw_sub, 10.0)
                r_ind, c_ind = linear_sum_assignment(sub_c)

                red_idx = gly.reducing_end_index
                init_genes_in_loc = [
                    (idx_g, int(nd.cluster_id))
                    for idx_g, nd in enumerate(loc.nodes)
                    if nd.role == GeneRole.INITIATOR and nd.cluster_id > 0
                ]
                prec_genes_in_loc = [
                    (idx_g, int(nd.cluster_id))
                    for idx_g, nd in enumerate(loc.nodes)
                    if nd.role in precursor_roles and nd.cluster_id > 0
                ]

                gt_gene_to_edge: dict[int, tuple[int, float]] = {}
                n_matched = len(r_ind)
                for idx_m in range(n_matched):
                    r_m = int(r_ind[idx_m])
                    c_m = int(c_ind[idx_m])
                    if sub_c[r_m, c_m] >= 9.0:
                        continue
                    g_idx = link_genes[r_m]
                    k_loc = core_edges[c_m]
                    c_id = int(loc.nodes[g_idx].cluster_id)
                    if c_id <= 0:
                        continue

                    c_mm = float(sub_c[r_m, c_m])
                    sig_m = _edge_signatures(gly, k_loc)[1]
                    alt_sum = 0.0
                    for idx_o in range(n_matched):
                        if idx_o == idx_m:
                            continue
                        r_o = int(r_ind[idx_o])
                        c_o = int(c_ind[idx_o])
                        sig_o = _edge_signatures(gly, core_edges[c_o])[1]
                        if sig_m and sig_m == sig_o:
                            continue
                        c_oo = float(sub_c[r_o, c_o])
                        delta_sw = max(0.0, (float(sub_c[r_m, c_o]) + float(sub_c[r_o, c_m])) - (c_mm + c_oo))
                        alt_sum += float(np.exp(-min(delta_sw / 0.12, 30.0)))

                    w_assign = 1.0 / (1.0 + alt_sum)
                    gt_gene_to_edge[g_idx] = (k_loc, w_assign)

                    s_full, s_dl, s_da, mech = _edge_signatures(gly, k_loc)
                    for target_cid in (c_id, prior.fine_to_med.get(c_id, 0)):
                        if target_cid <= 0:
                            continue
                        sig_full_acc[target_cid][s_full] += w_assign
                        sig_dl_acc[target_cid][s_dl] += w_assign
                        sig_da_acc[target_cid][s_da] += w_assign
                        if mech > 0:
                            mech_acc[target_cid][mech] += w_assign

                    e_matched = gly.edges[k_loc]
                    donor_s = gly.nodes[e_matched.source].sugar
                    if e_matched.target == red_idx:
                        for idx_init, c_init in init_genes_in_loc:
                            if abs(g_idx - idx_init) <= 2:
                                for target_cid in (c_id, prior.fine_to_med.get(c_id, 0)):
                                    if target_cid > 0:
                                        init_step1_acc[(c_init, target_cid)][s_dl] += w_assign

                    for idx_prec, c_prec in prec_genes_in_loc:
                        if abs(g_idx - idx_prec) <= 3:
                            for target_cid in (c_id, prior.fine_to_med.get(c_id, 0)):
                                if target_cid > 0:
                                    prec_donor_acc[(c_prec, target_cid)][donor_s] += w_assign

                # Accumulate syntenic transition pairs among GT/Polymerase genes weighted by Pan-Glycan bigrams
                gt_indices = sorted(gt_gene_to_edge.keys())
                for pos_a in range(len(gt_indices) - 1):
                    g_a = gt_indices[pos_a]
                    g_b = gt_indices[pos_a + 1]
                    if g_b - g_a > 3:
                        continue
                    c_a = int(loc.nodes[g_a].cluster_id)
                    c_b = int(loc.nodes[g_b].cluster_id)
                    if c_a <= 0 or c_b <= 0:
                        continue
                    k_a, w_a = gt_gene_to_edge[g_a]
                    k_b, w_b = gt_gene_to_edge[g_b]
                    if k_a == k_b:
                        continue
                    e_a = gly.edges[k_a]
                    e_b = gly.edges[k_b]
                    if e_a.target == e_b.source or e_a.source == e_b.target:
                        _, dl_a, _, _ = _edge_signatures(gly, k_a)
                        _, dl_b, _, _ = _edge_signatures(gly, k_b)
                        bg_cnt = prior.pan_glycan_bigram_counts.get((dl_a, dl_b), 0)
                        bg_mult = 1.0 + 0.18 * float(np.log1p(bg_cnt))
                        for ca_tier, cb_tier in (
                            (c_a, c_b),
                            (prior.fine_to_med.get(c_a, 0), prior.fine_to_med.get(c_b, 0)),
                        ):
                            if ca_tier > 0 and cb_tier > 0:
                                syn_pair_acc[(ca_tier, cb_tier)][(dl_a, dl_b)] += w_a * w_b * bg_mult
                                syn_pair_recs[(ca_tier, cb_tier)].add(i)

            prior.cluster_sig_full = {c: dict(m) for c, m in sig_full_acc.items()}
            prior.cluster_sig_dl = {c: dict(m) for c, m in sig_dl_acc.items()}
            prior.cluster_sig_da = {c: dict(m) for c, m in sig_da_acc.items()}
            prior.cluster_mech = {c: dict(m) for c, m in mech_acc.items()}
            prior.syntenic_pair_dl = {p: dict(m) for p, m in syn_pair_acc.items()}
            prior.syntenic_pair_record_counts = {p: len(rs) for p, rs in syn_pair_recs.items()}
            prior.init_adj_gt_step1 = {p: dict(m) for p, m in init_step1_acc.items()}
            prior.precursor_gt_donor = {p: dict(m) for p, m in prec_donor_acc.items()}

        return prior, rotated_glycans