From 60761c77e92f5e3b45447897240c6e6eb63f978c Mon Sep 17 00:00:00 2001 From: Taeyoon Kim Date: Thu, 20 Aug 2026 09:35:10 +0900 Subject: [PATCH 1/3] refactor: classify isotype before EU-numbered allotype detection --- src/abcat/cdomain.py | 287 ++++++++++++++++++++++++++++--------------- 1 file changed, 189 insertions(+), 98 deletions(-) diff --git a/src/abcat/cdomain.py b/src/abcat/cdomain.py index aa18237..7a1fce0 100644 --- a/src/abcat/cdomain.py +++ b/src/abcat/cdomain.py @@ -10,7 +10,6 @@ logger = logging.getLogger(__name__) -# Data directory path DATA_DIR = Path(__file__).parent.parent.parent / "data" _C_GENE_DB: dict[str, Any] | None = None @@ -38,63 +37,159 @@ def _load_databases() -> tuple[dict[str, Any], dict[str, Any]]: return _C_GENE_DB, _ALLOTYPE_DB +def _make_aligner() -> PairwiseAligner: + aligner = PairwiseAligner() + aligner.mode = "global" + aligner.match_score = 2 + aligner.mismatch_score = -1 + aligner.open_gap_score = -5 + aligner.extend_gap_score = -1 + return aligner + + +def _isotype_candidates( + chain_category: str, + c_gene_db: dict[str, Any], +) -> dict[str, dict[str, tuple[str, dict[str, Any]]]]: + """Group constant-gene references by isotype before subclass matching.""" + grouped: dict[str, dict[str, tuple[str, dict[str, Any]]]] = {} + for gene_name, gene_info in c_gene_db.get(chain_category, {}).items(): + isotype = gene_info.get("isotype", "Unknown") + grouped.setdefault(isotype, {})[gene_name] = (gene_name, gene_info) + return grouped + + +def _classify_isotype( + c_seq: str, + chain_category: str, + c_gene_db: dict[str, Any], + aligner: PairwiseAligner, +) -> tuple[str, float]: + """Determine isotype first, using the N-terminal constant-region signature.""" + grouped = _isotype_candidates(chain_category, c_gene_db) + signature_len = min(70, len(c_seq)) + signature = c_seq[:signature_len] + + best_isotype = "Unknown" + best_score = float("-inf") + for isotype, genes in grouped.items(): + iso_score = float("-inf") + for _, gene_info in genes.values(): + ref_signature = gene_info["sequence"][:signature_len] + alignments = aligner.align(signature, ref_signature) + if alignments: + iso_score = max(iso_score, alignments[0].score) + if iso_score > best_score: + best_score = iso_score + best_isotype = isotype + + return best_isotype, best_score + + +def _select_subclass( + c_seq: str, + chain_category: str, + isotype: str, + c_gene_db: dict[str, Any], + aligner: PairwiseAligner, +) -> tuple[str | None, dict[str, Any] | None, float, Any | None]: + """Match the constant sequence only against genes belonging to the chosen isotype.""" + best_gene: str | None = None + best_info: dict[str, Any] | None = None + best_score = float("-inf") + best_alignment: Any | None = None + + for gene_name, gene_info in c_gene_db.get(chain_category, {}).items(): + if gene_info.get("isotype") != isotype: + continue + alignments = aligner.align(c_seq, gene_info["sequence"]) + if alignments and alignments[0].score > best_score: + best_score = alignments[0].score + best_gene = gene_name + best_info = gene_info + best_alignment = alignments[0] + + return best_gene, best_info, best_score, best_alignment + + +def _infer_eu_start(reference_c_seq: str, markers: list[dict[str, Any]]) -> int | None: + """Infer the EU numbering origin for the reference constant sequence.""" + del reference_c_seq # Reserved for future explicit EU-numbered reference tables. + starts = { + int(m["subclass_offset"]) - int(m["eu_position"]) + 1 + for rule in markers + for m in rule.get("markers", []) + if "subclass_offset" in m and "eu_position" in m + } + if len(starts) == 1: + return starts.pop() + if not starts: + logger.warning("No EU origin metadata available for allotype markers") + else: + logger.warning("Inconsistent EU origins in allotype markers: %s", sorted(starts)) + return None + + +def _get_query_index_for_eu_position( + alignment: Any, + eu_position: int, + reference_eu_start: int | None, +) -> int | None: + """Map a canonical EU position through reference and query pairwise alignment.""" + if reference_eu_start is None: + return None + ref_idx = eu_position - reference_eu_start + if ref_idx < 0: + return None + return _get_query_index_for_ref_offset(alignment, ref_idx, 0) + + +def _get_query_index_for_ref_offset( + alignment: Any, ref_idx: int, default_len: int +) -> int | None: + """Map a 0-indexed reference position to query position using aligned blocks.""" + if hasattr(alignment, "aligned") and len(alignment.aligned) >= 2: + query_blocks, ref_blocks = alignment.aligned[0], alignment.aligned[1] + for (q_start, q_end), (r_start, r_end) in zip(query_blocks, ref_blocks): + if r_start <= ref_idx < r_end: + return q_start + (ref_idx - r_start) + return None + if 0 <= ref_idx < default_len: + return ref_idx + return None + + def analyze_cdomain( sequence: str, v_domain_len: int = 0, chain_hint: ChainType | None = None, ) -> CAnalysisResult: - """ - Analyzes the Constant Domain of an antibody sequence. - Determines Isotype, Subclass, Isoallotype, and Allotype markers. - """ + """Analyze the constant domain in order: isotype -> subclass -> EU allotype.""" cleaned = clean_sequence(sequence) - - # Extract constant domain sequence if full sequence is provided - if v_domain_len > 0 and len(cleaned) > v_domain_len: - c_seq = cleaned[v_domain_len:] - else: - c_seq = cleaned + c_seq = cleaned[v_domain_len:] if v_domain_len > 0 and len(cleaned) > v_domain_len else cleaned c_gene_db, allotype_db = _load_databases() - aligner = PairwiseAligner() - aligner.mode = "global" - aligner.match_score = 2 - aligner.mismatch_score = -1 - aligner.open_gap_score = -5 - aligner.extend_gap_score = -1 - - best_match_key = None - best_c_info = None - best_score = -1e9 - best_alignment = None - chain_category = "heavy" + aligner = _make_aligner() - # Filter targets based on chain hint if provided - search_groups = [] if chain_hint == ChainType.HEAVY: - search_groups = [("heavy", c_gene_db.get("heavy", {}))] + search_groups = ["heavy"] elif chain_hint in (ChainType.KAPPA, ChainType.LAMBDA): - search_groups = [("light", c_gene_db.get("light", {}))] + search_groups = ["light"] else: - search_groups = [ - ("heavy", c_gene_db.get("heavy", {})), - ("light", c_gene_db.get("light", {})), - ] - - for cat_name, genes in search_groups: - for gene_name, gene_info in genes.items(): - ref_seq = gene_info["sequence"] - alignments = aligner.align(c_seq, ref_seq) - if alignments: - top = alignments[0] - if top.score > best_score: - best_score = top.score - best_match_key = gene_name - best_c_info = gene_info - best_alignment = top - chain_category = cat_name - - if not best_match_key or not best_c_info or not best_alignment: + search_groups = ["heavy", "light"] + + # Stage 1: determine isotype without letting allotype/subclass differences influence the class call. + best_chain: str | None = None + best_isotype = "Unknown" + best_iso_score = float("-inf") + for category in search_groups: + isotype, score = _classify_isotype(c_seq, category, c_gene_db, aligner) + if score > best_iso_score: + best_iso_score = score + best_isotype = isotype + best_chain = category + + if best_chain is None or best_isotype == "Unknown": return CAnalysisResult( c_region_sequence=c_seq, isotype="Unknown", @@ -105,19 +200,33 @@ def analyze_cdomain( isoallotypes=[], ) + # Stage 2: subclass is matched only inside the selected isotype. + best_match_key, best_c_info, best_score, best_alignment = _select_subclass( + c_seq, best_chain, best_isotype, c_gene_db, aligner + ) + if not best_match_key or not best_c_info or best_alignment is None: + return CAnalysisResult( + c_region_sequence=c_seq, + isotype=best_isotype, + subclass="Unknown", + matched_c_gene="None", + alignment_identity=0.0, + allotypes=[], + isoallotypes=[], + ) + ref_len = len(best_c_info["sequence"]) identity = round(min(1.0, max(0.0, best_score / (2.0 * ref_len))), 4) - - isotype = best_c_info["isotype"] subclass = best_c_info["subclass"] - # Call Allotypes and Isoallotypes + # Stage 3: determine allotypes from canonical EU-numbered positions. allotype_calls, isoallotype_calls = _call_allotypes( query_c_seq=c_seq, matched_subclass=subclass, - chain_category=chain_category, + chain_category=best_chain, allotype_db=allotype_db, alignment=best_alignment, + reference_c_seq=best_c_info["sequence"], ) present_allotypes = [a.allotype for a in allotype_calls if a.status == "present"] @@ -125,7 +234,7 @@ def analyze_cdomain( return CAnalysisResult( c_region_sequence=c_seq, - isotype=isotype, + isotype=best_isotype, subclass=subclass, matched_c_gene=best_match_key, alignment_identity=identity, @@ -147,13 +256,9 @@ def format_allotype_summary(present_allotypes: list[str]) -> str: match = re.match(r"^([A-Za-z]+[0-9]*m|\w+[\*\_]?)(.*)$", allo) if match: prefix, suffix = match.groups() - if prefix not in groups: - groups[prefix] = [] - groups[prefix].append(suffix) + groups.setdefault(prefix, []).append(suffix) else: - if "other" not in groups: - groups["other"] = [] - groups["other"].append(allo) + groups.setdefault("other", []).append(allo) formatted_parts = [] for prefix, suffixes in groups.items(): @@ -162,38 +267,20 @@ def format_allotype_summary(present_allotypes: list[str]) -> str: elif len(suffixes) == 1: formatted_parts.append(f"{prefix}{suffixes[0]}") else: - joined_suffixes = ",".join(suffixes) - formatted_parts.append(f"{prefix}{joined_suffixes}") + formatted_parts.append(f"{prefix}{','.join(suffixes)}") return ", ".join(formatted_parts) -def _get_query_index_for_ref_offset( - alignment: Any, ref_idx: int, default_len: int -) -> int | None: - """ - Maps 0-indexed reference C-gene position (ref_idx) to 0-indexed position in query_c_seq using PairwiseAligner alignment. - alignment.aligned[0] is query blocks, alignment.aligned[1] is target (ref) blocks. - """ - if hasattr(alignment, "aligned") and len(alignment.aligned) >= 2: - query_blocks, ref_blocks = alignment.aligned[0], alignment.aligned[1] - for (q_start, q_end), (r_start, r_end) in zip(query_blocks, ref_blocks): - if r_start <= ref_idx < r_end: - return q_start + (ref_idx - r_start) - return None - if 0 <= ref_idx < default_len: - return ref_idx - return None - - def _call_allotypes( query_c_seq: str, matched_subclass: str, chain_category: str, allotype_db: dict[str, Any], alignment: Any, + reference_c_seq: str, ) -> tuple[list[AllotypeMarkerCall], list[AllotypeMarkerCall]]: - """Scans key polymorphic position markers to assign Allotypes and Isoallotypes.""" + """Call allotype/isoallotype markers using canonical EU positions.""" allotypes: list[AllotypeMarkerCall] = [] isoallotypes: list[AllotypeMarkerCall] = [] @@ -201,7 +288,10 @@ def _call_allotypes( if not subclass_markers: return allotypes, isoallotypes - # Align query C-region to reference + eu_start = _infer_eu_start(reference_c_seq, subclass_markers) + if eu_start is None: + return allotypes, isoallotypes + for rule in subclass_markers: allotype_name = rule["allotype"] opposing = rule.get("opposing_allotype") @@ -209,29 +299,30 @@ def _call_allotypes( category = rule.get("type", "allotype") markers = rule.get("markers", []) - matches = {} + matches: dict[str, str] = {} all_matched = True + any_observed = False + + for marker in markers: + eu_pos = int(marker["eu_position"]) + expected_aa = marker["amino_acid"] + query_idx = _get_query_index_for_eu_position(best_alignment := alignment, eu_pos, eu_start) + if query_idx is None or not (0 <= query_idx < len(query_c_seq)): + all_matched = False + continue - for m in markers: - ref_offset = ( - m["subclass_offset"] - 1 - ) # 0-indexed offset in reference C-gene - expected_aa = m["amino_acid"] - eu_pos = m["eu_position"] - - query_idx = _get_query_index_for_ref_offset( - alignment, ref_offset, len(query_c_seq) - ) - - if query_idx is not None and 0 <= query_idx < len(query_c_seq): - actual_aa = query_c_seq[query_idx] - matches[f"EU_{eu_pos}"] = actual_aa - if actual_aa != expected_aa: - all_matched = False - else: + actual_aa = query_c_seq[query_idx] + matches[f"EU_{eu_pos}"] = actual_aa + any_observed = True + if actual_aa != expected_aa: all_matched = False - status = "present" if all_matched else "absent" + if all_matched and any_observed: + status = "present" + elif any_observed: + status = "absent" + else: + status = "inconclusive" call = AllotypeMarkerCall( allotype=allotype_name, From 16bb22f2a6ebbca4dec04638db46a0301a8552db Mon Sep 17 00:00:00 2001 From: Taeyoon Kim Date: Thu, 20 Aug 2026 09:35:47 +0900 Subject: [PATCH 2/3] test: verify EU-position based allotype mapping --- tests/test_cdomain_eu.py | 49 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) create mode 100644 tests/test_cdomain_eu.py diff --git a/tests/test_cdomain_eu.py b/tests/test_cdomain_eu.py new file mode 100644 index 0000000..50a8607 --- /dev/null +++ b/tests/test_cdomain_eu.py @@ -0,0 +1,49 @@ +from copy import deepcopy + +from abcat.cdomain import ( + _call_allotypes, + _load_databases, + _make_aligner, + analyze_cdomain, +) +from abcat.schemas import ChainType + + +def test_isotype_is_selected_before_subclass(): + c_gene_db, _ = _load_databases() + igg1 = c_gene_db["heavy"]["IGHG1"]["sequence"] + + result = analyze_cdomain(igg1, chain_hint=ChainType.HEAVY) + + assert result.isotype == "IgG" + assert result.subclass == "IgG1" + assert result.matched_c_gene == "IGHG1" + + +def test_allotype_matching_uses_eu_position_not_legacy_offset(): + c_gene_db, allotype_db = _load_databases() + reference = c_gene_db["heavy"]["IGHG1"]["sequence"] + markers = deepcopy(allotype_db["heavy"]["IgG1"]) + + # Corrupt the legacy subclass-relative offset for the EU 214 marker. The + # canonical EU coordinate must still resolve G1m17 correctly. + for rule in markers: + for marker in rule.get("markers", []): + if marker["eu_position"] == 214: + marker["subclass_offset"] = 1 + + alignment = _make_aligner().align(reference, reference)[0] + allotypes, isoallotypes = _call_allotypes( + query_c_seq=reference, + matched_subclass="IgG1", + chain_category="heavy", + allotype_db={"heavy": {"IgG1": markers}}, + alignment=alignment, + reference_c_seq=reference, + ) + + present = {a.allotype for a in allotypes if a.status == "present"} + assert "G1m17" in present + g1m17 = next(a for a in allotypes if a.allotype == "G1m17") + assert g1m17.matched_residues["EU_214"] == "K" + assert not any(a.allotype == "G1m17" for a in isoallotypes) From 8fd7dadee69898f858b8b2414cf594f32d9b3a3d Mon Sep 17 00:00:00 2001 From: Taeyoon Kim Date: Thu, 20 Aug 2026 09:36:35 +0900 Subject: [PATCH 3/3] fix: use explicit EU numbering for Fc allotype calls --- src/abcat/cdomain.py | 70 ++++++++++++++++++++-------------------- tests/test_cdomain_eu.py | 6 ++-- 2 files changed, 37 insertions(+), 39 deletions(-) diff --git a/src/abcat/cdomain.py b/src/abcat/cdomain.py index 7a1fce0..94e6e9d 100644 --- a/src/abcat/cdomain.py +++ b/src/abcat/cdomain.py @@ -15,6 +15,18 @@ _C_GENE_DB: dict[str, Any] | None = None _ALLOTYPE_DB: dict[str, Any] | None = None +# EU numbering origins for the constant-domain reference sequences used by the +# allotype database. These are explicit so allotype detection never depends on +# the legacy subclass-relative offsets stored in historical marker records. +_EU_REFERENCE_STARTS: dict[tuple[str, str], int] = { + ("heavy", "IgG1"): 118, + ("heavy", "IgG2"): 118, + ("heavy", "IgG3"): 118, + ("heavy", "IgG4"): 118, + ("heavy", "IgE"): 99, + ("light", "IGKC"): 109, +} + def _load_databases() -> tuple[dict[str, Any], dict[str, Any]]: global _C_GENE_DB, _ALLOTYPE_DB @@ -112,36 +124,9 @@ def _select_subclass( return best_gene, best_info, best_score, best_alignment -def _infer_eu_start(reference_c_seq: str, markers: list[dict[str, Any]]) -> int | None: - """Infer the EU numbering origin for the reference constant sequence.""" - del reference_c_seq # Reserved for future explicit EU-numbered reference tables. - starts = { - int(m["subclass_offset"]) - int(m["eu_position"]) + 1 - for rule in markers - for m in rule.get("markers", []) - if "subclass_offset" in m and "eu_position" in m - } - if len(starts) == 1: - return starts.pop() - if not starts: - logger.warning("No EU origin metadata available for allotype markers") - else: - logger.warning("Inconsistent EU origins in allotype markers: %s", sorted(starts)) - return None - - -def _get_query_index_for_eu_position( - alignment: Any, - eu_position: int, - reference_eu_start: int | None, -) -> int | None: - """Map a canonical EU position through reference and query pairwise alignment.""" - if reference_eu_start is None: - return None - ref_idx = eu_position - reference_eu_start - if ref_idx < 0: - return None - return _get_query_index_for_ref_offset(alignment, ref_idx, 0) +def _get_eu_reference_start(chain_category: str, matched_subclass: str) -> int | None: + """Return the explicit EU-numbering origin for the selected reference.""" + return _EU_REFERENCE_STARTS.get((chain_category, matched_subclass)) def _get_query_index_for_ref_offset( @@ -159,6 +144,20 @@ def _get_query_index_for_ref_offset( return None +def _get_query_index_for_eu_position( + alignment: Any, + eu_position: int, + reference_eu_start: int | None, +) -> int | None: + """Map a canonical EU position through reference and query pairwise alignment.""" + if reference_eu_start is None: + return None + ref_idx = eu_position - reference_eu_start + if ref_idx < 0: + return None + return _get_query_index_for_ref_offset(alignment, ref_idx, 0) + + def analyze_cdomain( sequence: str, v_domain_len: int = 0, @@ -178,7 +177,7 @@ def analyze_cdomain( else: search_groups = ["heavy", "light"] - # Stage 1: determine isotype without letting allotype/subclass differences influence the class call. + # Stage 1: determine isotype without letting subclass/allotype differences influence the class call. best_chain: str | None = None best_isotype = "Unknown" best_iso_score = float("-inf") @@ -219,7 +218,7 @@ def analyze_cdomain( identity = round(min(1.0, max(0.0, best_score / (2.0 * ref_len))), 4) subclass = best_c_info["subclass"] - # Stage 3: determine allotypes from canonical EU-numbered positions. + # Stage 3: determine allotypes from canonical EU-numbered constant/Fc positions. allotype_calls, isoallotype_calls = _call_allotypes( query_c_seq=c_seq, matched_subclass=subclass, @@ -288,8 +287,9 @@ def _call_allotypes( if not subclass_markers: return allotypes, isoallotypes - eu_start = _infer_eu_start(reference_c_seq, subclass_markers) - if eu_start is None: + reference_eu_start = _get_eu_reference_start(chain_category, matched_subclass) + if reference_eu_start is None: + logger.warning("No EU reference origin for %s/%s", chain_category, matched_subclass) return allotypes, isoallotypes for rule in subclass_markers: @@ -306,7 +306,7 @@ def _call_allotypes( for marker in markers: eu_pos = int(marker["eu_position"]) expected_aa = marker["amino_acid"] - query_idx = _get_query_index_for_eu_position(best_alignment := alignment, eu_pos, eu_start) + query_idx = _get_query_index_for_eu_position(alignment, eu_pos, reference_eu_start) if query_idx is None or not (0 <= query_idx < len(query_c_seq)): all_matched = False continue diff --git a/tests/test_cdomain_eu.py b/tests/test_cdomain_eu.py index 50a8607..3e27b05 100644 --- a/tests/test_cdomain_eu.py +++ b/tests/test_cdomain_eu.py @@ -25,12 +25,10 @@ def test_allotype_matching_uses_eu_position_not_legacy_offset(): reference = c_gene_db["heavy"]["IGHG1"]["sequence"] markers = deepcopy(allotype_db["heavy"]["IgG1"]) - # Corrupt the legacy subclass-relative offset for the EU 214 marker. The - # canonical EU coordinate must still resolve G1m17 correctly. + # EU position is canonical; legacy subclass offsets should not be required. for rule in markers: for marker in rule.get("markers", []): - if marker["eu_position"] == 214: - marker["subclass_offset"] = 1 + marker.pop("subclass_offset", None) alignment = _make_aligner().align(reference, reference)[0] allotypes, isoallotypes = _call_allotypes(