|
| 1 | +from typing import NamedTuple |
| 2 | + |
| 3 | +import requests |
| 4 | + |
| 5 | +from models import Collection, Variant |
| 6 | +from sources import Source |
| 7 | + |
| 8 | +RSV_A_TREE_URL = ( |
| 9 | + "https://raw.githubusercontent.com/nextstrain/nextclade_data" |
| 10 | + "/refs/heads/master/data_output/nextstrain/rsv/a" |
| 11 | + "/EPI_ISL_412866/2026-04-14--11-55-23Z/tree.json" |
| 12 | +) |
| 13 | + |
| 14 | +RSV_B_TREE_URL = ( |
| 15 | + "https://raw.githubusercontent.com/nextstrain/nextclade_data" |
| 16 | + "/refs/heads/master/data_output/nextstrain/rsv/b" |
| 17 | + "/EPI_ISL_1653999/2026-04-14--11-55-23Z/tree.json" |
| 18 | +) |
| 19 | + |
| 20 | + |
| 21 | +class _RsvNextcladeLineagesBase(Source): |
| 22 | + owned_tag = "#nextclade-lineage" |
| 23 | + _tree_url: str |
| 24 | + _organism_label: str |
| 25 | + |
| 26 | + def get_collections(self) -> list[Collection]: |
| 27 | + print(f"Fetching Nextclade tree from {self._tree_url} ...") |
| 28 | + response = requests.get(self._tree_url, timeout=60) |
| 29 | + response.raise_for_status() |
| 30 | + tree_json = response.json() |
| 31 | + collections = _build_collections(tree_json, self.organism, self._organism_label, self.owned_tag) |
| 32 | + print(f" Loaded {len(collections)} clade(s).") |
| 33 | + return collections |
| 34 | + |
| 35 | + |
| 36 | +class RsvANextcladeLineagesSource(_RsvNextcladeLineagesBase): |
| 37 | + name = "rsv-a-nextclade-lineages" |
| 38 | + organism = "rsvA" |
| 39 | + _tree_url = RSV_A_TREE_URL |
| 40 | + _organism_label = "RSV-A" |
| 41 | + |
| 42 | + |
| 43 | +class RsvBNextcladeLineagesSource(_RsvNextcladeLineagesBase): |
| 44 | + name = "rsv-b-nextclade-lineages" |
| 45 | + organism = "rsvB" |
| 46 | + _tree_url = RSV_B_TREE_URL |
| 47 | + _organism_label = "RSV-B" |
| 48 | + |
| 49 | + |
| 50 | +class _CladeInfo(NamedTuple): |
| 51 | + clade_name: str |
| 52 | + parent_clade: str | None |
| 53 | + branch_nuc: list[str] |
| 54 | + branch_aa: list[str] |
| 55 | + full_nuc: list[str] |
| 56 | + full_aa: list[str] |
| 57 | + |
| 58 | + |
| 59 | +def _build_collections(tree_json: dict, organism: str, organism_label: str, owned_tag: str) -> list[Collection]: |
| 60 | + """Build one Collection per clade in the Nextclade reference tree. |
| 61 | +
|
| 62 | + Each collection gets four variants, mirroring the shape used for COVID Pango lineages: |
| 63 | +
|
| 64 | + - "Nucleotide substitutions" — full set from reference root (for searching |
| 65 | + sequences that carry all defining substitutions |
| 66 | + of this clade) |
| 67 | + - "Amino acid substitutions" — full AA set from reference root |
| 68 | + - "New nucleotide substitutions" — branch-only mutations (what is newly introduced |
| 69 | + in this clade step relative to its parent clade) |
| 70 | + - "New amino acid substitutions" — branch-only AA mutations |
| 71 | + """ |
| 72 | + collections = [] |
| 73 | + for clade in _extract_clades(tree_json["tree"]): |
| 74 | + parent_str = clade.parent_clade or "—" |
| 75 | + variants: list[Variant] = [ |
| 76 | + { |
| 77 | + "type": "filterObject", |
| 78 | + "name": "Nucleotide substitutions", |
| 79 | + "filterObject": {"nucleotideMutations": clade.full_nuc}, |
| 80 | + }, |
| 81 | + { |
| 82 | + "type": "filterObject", |
| 83 | + "name": "Amino acid substitutions", |
| 84 | + "filterObject": {"aminoAcidMutations": clade.full_aa}, |
| 85 | + }, |
| 86 | + { |
| 87 | + "type": "filterObject", |
| 88 | + "name": "New nucleotide substitutions", |
| 89 | + "filterObject": {"nucleotideMutations": clade.branch_nuc}, |
| 90 | + }, |
| 91 | + { |
| 92 | + "type": "filterObject", |
| 93 | + "name": "New amino acid substitutions", |
| 94 | + "filterObject": {"aminoAcidMutations": clade.branch_aa}, |
| 95 | + }, |
| 96 | + ] |
| 97 | + description = ( |
| 98 | + f"{organism_label} Nextclade clade {clade.clade_name}. " |
| 99 | + f"Parent clade: {parent_str}. " |
| 100 | + f"{owned_tag}" |
| 101 | + ) |
| 102 | + collections.append({ |
| 103 | + "name": clade.clade_name, |
| 104 | + "organism": organism, |
| 105 | + "description": description, |
| 106 | + "variants": variants, |
| 107 | + }) |
| 108 | + return collections |
| 109 | + |
| 110 | + |
| 111 | +def _extract_clades(node, parent_clade=None, accum_nuc=None, accum_aa=None): |
| 112 | + """Walk the Nextclade reference tree recursively, yielding one _CladeInfo per introduced clade. |
| 113 | +
|
| 114 | + Fields: |
| 115 | + - clade_name: the clade label from branch_attrs.labels.clade, e.g. "A.D.3.5" |
| 116 | + - parent_clade: node_attrs.clade_membership of the parent node (None for the root clade) |
| 117 | + - branch_nuc: nucleotide mutations on this branch only (relative to parent node), |
| 118 | + e.g. ["C982T", "G108A"] |
| 119 | + - branch_aa: AA mutations on this branch only, formatted as GENE:MUT, |
| 120 | + e.g. ["F:T8A", "G:T4N"] |
| 121 | + - full_nuc: all nucleotide mutations from the reference root down to and including |
| 122 | + this clade's branch, expressed relative to the root reference sequence, |
| 123 | + e.g. ["A982T", "G108A", "C241T"] |
| 124 | + - full_aa: all AA mutations from root to this clade, expressed relative to root, |
| 125 | + e.g. ["F:T8G", "G:T4N"] |
| 126 | +
|
| 127 | + A clade is introduced at any node where branch_attrs.labels.clade is set. |
| 128 | +
|
| 129 | + The root reference sequence (e.g. EPI_ISL_412866 for RSV-A) defines position zero |
| 130 | + for all accumulated mutations — i.e. full_nuc and full_aa give the complete set of |
| 131 | + substitutions needed to go from the reference to this clade. |
| 132 | +
|
| 133 | + accum_nuc / accum_aa carry the accumulated state downward. Because _apply_*_mutations |
| 134 | + always returns a new dict rather than mutating in place, each recursive call receives |
| 135 | + its own independent copy of the state — sibling subtrees cannot interfere with each other. |
| 136 | + """ |
| 137 | + if accum_nuc is None: |
| 138 | + accum_nuc = {} |
| 139 | + if accum_aa is None: |
| 140 | + accum_aa = {} |
| 141 | + |
| 142 | + labels = node.get("branch_attrs", {}).get("labels", {}) |
| 143 | + muts = node.get("branch_attrs", {}).get("mutations", {}) |
| 144 | + node_clade = node.get("node_attrs", {}).get("clade_membership", {}).get("value") |
| 145 | + |
| 146 | + # Separate nucleotide mutations (key "nuc") from amino acid mutations (all other keys |
| 147 | + # are gene names such as "F", "G", "L", etc.). |
| 148 | + branch_nuc = muts.get("nuc", []) |
| 149 | + branch_aa_by_gene = {k: v for k, v in muts.items() if k != "nuc"} |
| 150 | + |
| 151 | + # Fold this branch's mutations into the running accumulated state. |
| 152 | + # The returned dicts are new objects — the parent's dicts are untouched. |
| 153 | + accum_nuc = _apply_nuc_mutations(branch_nuc, accum_nuc) |
| 154 | + accum_aa = _apply_aa_mutations(branch_aa_by_gene, accum_aa) |
| 155 | + |
| 156 | + if "clade" in labels: |
| 157 | + # Flatten branch-level AA mutations into GENE:MUT strings for the "new" variant. |
| 158 | + branch_aa_flat = [ |
| 159 | + f"{gene}:{m}" |
| 160 | + for gene, gene_muts in branch_aa_by_gene.items() |
| 161 | + for m in gene_muts |
| 162 | + ] |
| 163 | + yield _CladeInfo( |
| 164 | + clade_name=labels["clade"], |
| 165 | + parent_clade=parent_clade, |
| 166 | + branch_nuc=branch_nuc, |
| 167 | + branch_aa=branch_aa_flat, |
| 168 | + full_nuc=_format_accum_nuc(accum_nuc), |
| 169 | + full_aa=_format_accum_aa(accum_aa), |
| 170 | + ) |
| 171 | + |
| 172 | + # Pass the current node's clade_membership as the parent context for children. |
| 173 | + # If a node has no clade_membership (e.g. the synthetic root), fall back to whatever |
| 174 | + # was passed in from above. |
| 175 | + next_parent = node_clade or parent_clade |
| 176 | + for child in node.get("children", []): |
| 177 | + yield from _extract_clades(child, next_parent, accum_nuc, accum_aa) |
| 178 | + |
| 179 | + |
| 180 | +def _apply_nuc_mutations( |
| 181 | + branch_muts: list[str], |
| 182 | + accum: dict[str, tuple[str, str]], |
| 183 | +) -> dict[str, tuple[str, str]]: |
| 184 | + """Apply a list of branch-level nucleotide mutations on top of an accumulated mutation dict. |
| 185 | +
|
| 186 | + Each mutation string has the form REF_BASE + POSITION + NEW_BASE, e.g. "C982T", where: |
| 187 | + - REF_BASE is the base in the *parent* node at this position. |
| 188 | + - POSITION is the 1-based coordinate in the genome. |
| 189 | + - NEW_BASE is the base introduced by this branch. |
| 190 | +
|
| 191 | + The accumulated dict maps POSITION (str) → (ORIG_REF_BASE, CURRENT_BASE), where: |
| 192 | + - ORIG_REF_BASE is the base in the root reference genome (not the parent). |
| 193 | + - CURRENT_BASE is the base after all mutations from root to here. |
| 194 | +
|
| 195 | + Three cases when applying a branch mutation at a position: |
| 196 | +
|
| 197 | + 1. Position not yet in accum (first time this position is mutated on this path): |
| 198 | + The parent base IS the reference base, so we record (ref_base_from_mut, new_base). |
| 199 | +
|
| 200 | + 2. Position already in accum (mutated earlier on the path from root): |
| 201 | + The parent base is no longer the reference base. We keep the original reference |
| 202 | + base from the existing entry and update only the current base. |
| 203 | +
|
| 204 | + 3. New base equals original reference base (reversion to reference): |
| 205 | + Net effect from root is zero — remove the position from the dict entirely. |
| 206 | +
|
| 207 | + Returns a new dict; the input is not modified (so sibling branches are unaffected). |
| 208 | + """ |
| 209 | + result = dict(accum) |
| 210 | + for mut in branch_muts: |
| 211 | + ref_base = mut[0] # base in the parent (= reference if first mutation here) |
| 212 | + pos = mut[1:-1] # genome position as string, e.g. "982" |
| 213 | + new_base = mut[-1] # base introduced by this branch |
| 214 | + existing = result.get(pos) |
| 215 | + # If this position was already mutated earlier on the path, keep the original |
| 216 | + # reference base; otherwise the parent base is the reference base. |
| 217 | + orig_ref = existing[0] if existing else ref_base |
| 218 | + if new_base == orig_ref: |
| 219 | + # Reverted back to the reference state — remove from accumulated set. |
| 220 | + result.pop(pos, None) |
| 221 | + else: |
| 222 | + result[pos] = (orig_ref, new_base) |
| 223 | + return result |
| 224 | + |
| 225 | + |
| 226 | +def _apply_aa_mutations( |
| 227 | + branch_muts_by_gene: dict[str, list[str]], |
| 228 | + accum: dict[tuple[str, str], tuple[str, str]], |
| 229 | +) -> dict[tuple[str, str], tuple[str, str]]: |
| 230 | + """Apply branch-level amino acid mutations on top of an accumulated AA mutation dict. |
| 231 | +
|
| 232 | + branch_muts_by_gene maps gene name → list of mutation strings, e.g. |
| 233 | + {"F": ["T8A", "L20F"], "G": ["T4N"]}. |
| 234 | +
|
| 235 | + Each mutation string has the form REF_AA + POSITION + NEW_AA, e.g. "T8A", where: |
| 236 | + - REF_AA is the amino acid in the *parent* node at this codon. |
| 237 | + - POSITION is the 1-based codon position within the gene. |
| 238 | + - NEW_AA is the amino acid introduced by this branch. |
| 239 | +
|
| 240 | + The accumulated dict maps (GENE, POSITION) → (ORIG_REF_AA, CURRENT_AA), where: |
| 241 | + - ORIG_REF_AA is the amino acid in the root reference genome. |
| 242 | + - CURRENT_AA is the amino acid after all mutations from root to here. |
| 243 | +
|
| 244 | + The same three accumulation cases apply as for nucleotide mutations. |
| 245 | +
|
| 246 | + Returns a new dict; the input is not modified. |
| 247 | + """ |
| 248 | + result = dict(accum) |
| 249 | + for gene, muts in branch_muts_by_gene.items(): |
| 250 | + for mut in muts: |
| 251 | + ref_aa = mut[0] # amino acid in the parent |
| 252 | + pos = mut[1:-1] # codon position within the gene |
| 253 | + new_aa = mut[-1] # amino acid introduced by this branch |
| 254 | + key = (gene, pos) |
| 255 | + existing = result.get(key) |
| 256 | + orig_ref = existing[0] if existing else ref_aa |
| 257 | + if new_aa == orig_ref: |
| 258 | + result.pop(key, None) |
| 259 | + else: |
| 260 | + result[key] = (orig_ref, new_aa) |
| 261 | + return result |
| 262 | + |
| 263 | + |
| 264 | +def _format_accum_nuc(accum: dict[str, tuple[str, str]]) -> list[str]: |
| 265 | + """Render accumulated nucleotide mutations as ORIG_REF_BASE + POSITION + CURRENT_BASE strings. |
| 266 | +
|
| 267 | + E.g. {"982": ("A", "T")} → ["A982T"], meaning: at genome position 982 the reference |
| 268 | + has A and this clade (from root) has T. |
| 269 | + """ |
| 270 | + return [f"{ref}{pos}{cur}" for pos, (ref, cur) in accum.items()] |
| 271 | + |
| 272 | + |
| 273 | +def _format_accum_aa(accum: dict[tuple[str, str], tuple[str, str]]) -> list[str]: |
| 274 | + """Render accumulated AA mutations as GENE:ORIG_REF_AA + POSITION + CURRENT_AA strings. |
| 275 | +
|
| 276 | + E.g. {("F", "8"): ("T", "G")} → ["F:T8G"], meaning: in gene F at codon 8 the reference |
| 277 | + has T and this clade (from root) has G. |
| 278 | + """ |
| 279 | + return [f"{gene}:{ref}{pos}{cur}" for (gene, pos), (ref, cur) in accum.items()] |
0 commit comments