Skip to content

[Rabies] Initialize Lyssavirus rabies all-clades community dataset #333

Open
xonq wants to merge 7 commits into
nextstrain:masterfrom
xonq:master
Open

[Rabies] Initialize Lyssavirus rabies all-clades community dataset #333
xonq wants to merge 7 commits into
nextstrain:masterfrom
xonq:master

Conversation

@xonq

@xonq xonq commented Aug 6, 2025

Copy link
Copy Markdown
Contributor

This pull request initializes a Lyssavirus rabies (rabies) Nextclade dataset with clade-subclade resolution. Created in collaboration with @kimandrews and with subject matter expertise/user input from Massachusetts Department of Public Health. Please review the README.md for information on dataset creation and citations.

@ivan-aksamentov

Copy link
Copy Markdown
Member

Thanks! Seems to be working

https://master.clades.nextstrain.org/?dataset-server=gh:xonq/nextclade_data@master@/data_output&dataset-name=community/theiagen/rabies/all-clades

As a dev I can only review the technical side. And I will let our scientists to check the sciency bits :)

The virus is quite diverse it seems - lots of mutations. But this is probably expected.

If you have an public repo where you prepare trees and other data for the dataset, it would be a great help to the users of your dataset if you add it to the readme. We typically use a boilerplate like this in Nextstrain datasets:

| Key | Value |
| ----------------- | ---------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| authors | [Cornelius Roemer](https://neherlab.org), [Richard Neher](https://neherlab.org), [Nextstrain](https://nextstrain.org) |
| reference | `Wuhan-Hu-1/2019` |
| workflow | https://github.com/neherlab/nextclade_data_workflows/tree/v3-sc2/sars-cov-2 |
| path | `nextstrain/sars-cov-2/orfs` |
| clade definitions | [Nextstrain clades](https://nextstrain.org/blog/2022-04-29-SARS-CoV-2-clade-naming-2022) and [Pango lineages](https://www.nature.com/articles/s41564-020-0770-5) |

But that's not mandatory.

All looks good to me. Smooth work!

@xonq

xonq commented Aug 7, 2025

Copy link
Copy Markdown
Contributor Author

Thank you. Rabies is indeed very diverse - I contemplated creating independent datasets for each clade, but the genotyping of this "all-clades" dataset has been sufficient for our SME partners. Additionally, there were issues with sub-clade metadata quality that limit the improvements more refined datasets may provide.

I do not have a repository for tree building - I built the Nextclade dataset from the Nextstrain rabies build as a template, though I ended up deviating with the tree building methodology and metadata acquisition. The methodology is hopefully adequately documented for users in this PR's README.

@rneher

rneher commented Aug 11, 2025

Copy link
Copy Markdown
Member

Thanks a lot for contributing this dataset! Overall, this looks very good. But I have a few suggestions to make it better.

  • the amino acid annotation is inconsistent with the annotation in Nextclade. On the tree, CDS are called RABVgp1_N etc, in the Nextclade annotation you have things like NP_062343. Harmonizing these annotation would be important.
  • Given the diversity of the virus, I would recommend different alignment parameters. See here for an example.
  • I am not sure how exactly you do things in your workflow, but we recommend aligning and translating the sequences in the tree exactly like Nextclade will align query sequences later. see here for an example. The translation can be later used in ancestral reconstruction.

@xonq

xonq commented Aug 20, 2025

Copy link
Copy Markdown
Contributor Author

hey @rneher, just wanted to reply and inform you that I cannot return to this to address your concerns until a later date. not sure when, but hopefully within the next several weeks. Thanks for your suggestions and my apologies for my ignorance to some of the standardized procedures.

RE: alignment parameters: I'm not really certain how to systematically adjust these parameters - do you have specific recommendations/procedures to determine what parameters are more ideal, or do you suggest dragging and dropping the linked pathogen.json you sent?

RE: Apart from the tree-building, the workflow was performed with AUGUR. With this in mind, do these steps deviate from Nextclade like you're suggesting?:

Alignment:

# align genomes relative to the reference
augur align \
  -s <GENOMES.fa> \
  --nthreads <CPUS> \
  --output aligned.fa \
  --reference-name <REFERENCE_ACCESSION> \
  --debug

Tree building:

performed as discussed in the README

Refinement:

augur refine \
  --tree <IQTREE.contree> \
  --alignment <ALIGNMENT.fa> \
  --metadata <METADATA.tsv> \
  --output-tree refined.nwk \
  --output-node-data branch_lengths.json \
  --keep-root \
  --metadata-id-columns <TIP_NAMES_COLUMN>

Trait application:

augur traits \
	--tree refined.nwk \
	--metadata <METADATA.tsv> \
	--output-node-data traits.json \
	--columns <COLUMN1> <COLUMN2> .. <COLUMNn> \
	--confidence \
	--metadata-id-columns <TIP_NAMES_COLUMN>

Nucleotide mutation calling:

augur ancestral \
  --tree refined.nwk \
  --alignment <ALIGNMENT.fa> \
  --output-node-data nt_muts.json \
  --root-sequence <REFERENCE.fasta> \
  --inference joint

Translation:

augur translate \
  --tree refined.nwk \
  --ancestral-sequences nt_muts.json \
  --reference-sequence <REFERENCE.gbk> \
  --output-node-data aa_muts.json

Clade mutation extraction (non-AUGUR):

extract_nextclades.py \
  -t refined.nwk \
  -m <METADATA.tsv> \
  -cc <CLADE_COLUMN1> <CLADE_COLUMN2> .. <CLADE_COLUMNn> \
  -tc <TIP_NAMES_COLUMN> \
  -aa aa_muts.json \
  -nt nt_muts.json \
  -n

Clade mutation application:

augur clades \
 --tree refined.nwk \
 --mutations nt_muts.json aa_muts.json \
 --clades clades.tsv \
 --output-node-data clades.json

Export:

augur export v2 \
 --tree refined.nwk \
 --metadata <METADATA.tsv> \
 --node-data branch_lengths.json \
   traits.json \
   nt_muts.json \
   aa_muts.json \
   clades.json \
 --lat-longs <LATITUDE_LONGITUDES.tsv> \
 --colors <COLORS.tsv> \
 --auspice-config <AUSPICE_CONFIG>.json \
 --output <FINAL_AUGUR_OUTPUT.json> \
 --metadata-id-columns <TIP_NAMES_COLUMN>

@xonq

xonq commented Jul 21, 2026

Copy link
Copy Markdown
Contributor Author

Hey all, just chiming in to see if I could get some guidance on a couple points to move this forward:

RE: alignment parameters: I'm not really certain how to systematically adjust these parameters - do you have specific recommendations/procedures to determine what parameters are more ideal, or do you suggest dragging and dropping the linked pathogen.json you sent?

RE: Apart from the tree-building, the workflow was performed with AUGUR. With this in mind, do these steps deviate from Nextclade like you're suggesting?:

@rneher

rneher commented Jul 22, 2026

Copy link
Copy Markdown
Member

Hi @xonq

sorry for not responding earlier. I had read the first part of your message "... later day." and mentally shelved it under 'later'.

Re alignment parameters: you can add these to the pathogen json. There is also a pre-set high-diversity.

related, we recommend doing alignment and translation for your tree using nextclade rather than augur align or translate. This is mostly to be sure that what nextclade produces for query sequences are consistent with what is in the tree. for the measles build linked about this is happening here and here.

The latter is probably also responsible for mismatch of the amino acid mutations in the alignment/output and the reference tree. In the reference tree, CDS are named like RABVgp5_L but in the nextclade alignment as NP_056795.1

@ivan-aksamentov

ivan-aksamentov commented Jul 22, 2026

Copy link
Copy Markdown
Member

Posting robotic review below.

It mostly parrots what Richard already said, plus some low prio typos/citation oopsies (to be confirmed!). Nothing much to stress about


⚠️ AI-generated content below. Verify all claims.

First review of data/community/theiagen/rabies/all-clades/ at head f6d98750.

Testing

Observed

Whole-genome Lyssavirus rabies genotyping dataset -- 3,296 tips, 48 clades, reference NC_001542.1. Diff matches the PR description. One undocumented item: a local build path in the tree's seqid fields (F5).

Blocking issues

🔴 F1. Amino-acid CDS names in the tree do not match the genome annotation [click to expand]

The reference tree and the genome annotation name the five CDS differently:

CDS tree.json (meta.genome_annotations and branch_attrs.mutations keys) genome_annotation.gff3 (Name= / protein_id)
N RABVgp1_N NP_056793.1
M1 RABVgp2_M1 NP_056794.1
M2 RABVgp3_M2 NP_056795.1
G RABVgp4_G NP_056796.1
L RABVgp5_L NP_056797.1

Coordinates match (L = 5418-11846 in both), only the names differ. Nextclade derives CDS names from the GFF3, so query AA mutations show up as NP_056797.1:... while the tree carries RABVgp5_L:.... The two never reconcile. Nuc-based clade assignment still works; the AA mutation panel does not.

Same mismatch @rneher raised twice (2025-08-11, 2026-07-22).

Suggestions:

  • Align and translate the tree sequences with Nextclade itself, so both sides use the same GFF3 names -- see the measles build: nextstrain/measles Snakefile#L67-L88
  • Or rename the GFF3 CDS to match the tree (Name=RABVgp1_N, etc.)

Non-blocking issues

🟡 F2. README cites an unrelated paper (wrong Campbell 2022) [click to expand]

The README.md References section lists:

Campbell F, Blank L, Cantrell A, Baxter S, Blackmore C, Dixon J, Goyder E. Factors that influence mental health of university and college students in the UK: a systematic review. BMC Public Health. 2022 ... doi: 10.1186/s12889-022-13943-x.

That DOI resolves to a real but entirely unrelated paper on mental health of UK university students. The dataset's subclade labels (AF1a, AF1b, AM1, SEA1b, ME1a, CE, WE, ...) are the Pango-style RABV lineage nomenclature from Campbell et al., PLoS Pathogens 2022, "Making genomic surveillance deliver: A lineage classification and nomenclature system to inform rabies elimination" (doi:10.1371/journal.ppat.1010023), which is integrated with RABV-GLUE -- the resource the README already cites for clade metadata. Users looking up the nomenclature source land on a mental-health paper instead of the rabies lineage system.

Suggestions:

🟡 F3. No alignment tuning for a divergent virus [click to expand]

pathogen.json sets no alignmentParams. Running the 21 example sequences through Nextclade confirms the problem: 20/21 produce frameshifts (23 in the L gene, 3 in G), and 12/21 score "bad" QC overall. The reference itself passes clean (0 mutations, "good"). This pattern -- systematic frameshifts concentrated in the longest gene across nearly all examples -- is a textbook sign that the aligner needs higher gap penalties.

@rneher recommended tuning alignment or using the high-diversity preset. @ivan-aksamentov noted "lots of mutations" when testing.

Suggestions:

  • Quick fix: add "alignmentPreset": "high-diversity" to pathogen.json
  • Custom tuning: add an "alignmentParams" block. The measles dataset is a good starting point for RABV:
    "alignmentParams": {
      "penaltyGapOpen": 10,
      "penaltyGapOpenInFrame": 15,
      "penaltyGapOpenOutOfFrame": 17,
      "kmerLength": 6,
      "minSeedCover": 0.1
    }
  • See the Nextclade alignment docs for all available parameters
  • Rerun examples after tuning -- the L-gene frameshifts should largely disappear
🟡 F4. QC configuration is minimal compared with sibling datasets [click to expand]

Only missingData and mixedSites are configured. The official measles dataset (which @rneher linked as a reference) configures all six QC rules: privateMutations, frameShifts, stopCodons, snpClusters, missingData, mixedSites. Without the missing four, broken ORFs and divergent sequences pass without a flag.

Suggestions:

  • Add privateMutations, frameShifts, stopCodons, snpClusters rules, scaled to genome size
  • privateMutations needs a generous cutoff given RABV diversity
  • Example: measles pathogen.json QC block (the dataset @rneher pointed to, with all six QC rules configured)
  • See Nextclade QC docs
🔵 F5. Local build path leaked into tree annotations [click to expand]

Every entry in tree.json meta.genome_annotations has "seqid": "../independent_attempt/manual_iqtree/archive/ref_genes.gbk". Nextclade ignores this field, but the path is visible in Auspice and exposes build internals.

Suggestions:

  • Set seqid to the reference accession (NC_001542.1)
  • Or regenerate the tree annotations from the dataset GFF3 (also addresses F1)
🔵 F6. Boilerplate/description mismatch in README and dataset name [click to expand]

README.md opens "This is a minimal example dataset for Lyssavirus rabies" -- leftover template text that contradicts a 3,296-tip full-genome genotyping dataset. attributes.name is "Lyssavirus rabies GenBank WGS" -- the "GenBank WGS" suffix reads as build jargon. Both show up in the Nextclade dataset picker, so users see a confusing name and description.

Suggestions:

  • Rewrite the README opening to describe the actual dataset
  • Consider a cleaner display name, e.g. "Lyssavirus rabies (all clades)"
🔵 F7. Citation formatting errors [click to expand]

In README.md References: the Baynard & Fooks entry has a malformed DOI doi:doi.org/10.1016/B978-0-12-809633-8.20936-9 (double doi); the correct DOI resolves. The Troupin entry date is garbled ("PLoS Pathog 2016 Dec, 15;12(12)"); the paper is PLoS Pathog 2016 Dec 15;12(12):e1006041. The malformed DOI won't resolve if someone clicks it.

Suggestions:

  • Remove the duplicate doi: prefix from the Baynard entry
  • Fix the Troupin volume/date formatting
🔵 F8. Missing trailing newlines [click to expand]

README.md and tree.json lack a final newline.

Suggestions:

  • Add trailing newlines to both files

Clade distribution

48 clades, 3,296 tips [click to expand]

Coverage spans all major RABV clades. Dominant clades reflect sequencing effort (North American RAC-SK, Arctic, and SE-Asian lineages are heavily represented). Several README-included subclades are present but sparsely sampled (<3 tips), which limits confident assignment for those lineages.

Clade Tips %
RAC-SK 740 22.5
Arctic A 327 9.9
Asian SEA4 313 9.5
Asian SEA1b 238 7.2
Cosmopolitan Vac 205 6.2
Arctic AL1a 201 6.1
Cosmopolitan AF1b 184 5.6
Asian 109 3.3
Cosmopolitan NEE 107 3.2
(29 mid clades) 3--78 ea.
Bats TB2 2 0.1
Bats PS 2 0.1
Bats LI 2 0.1
Bats MYsp 2 0.1
Cosmopolitan AF4 2 0.1
Cosmopolitan AM2b 2 0.1
Cosmopolitan ME1b 2 0.1
Total 3,296

Validation summary

Validation checks [click to expand]

GFF3 annotation

  • 5/5 CDS divisible by 3, all ATG starts, all valid stops (TAA/TGA)
  • Coordinates match ##sequence-region 1 11932
  • Fail: CDS names differ between tree (RABVgpN_X) and GFF3 (NP_0567xx.1) -- see F1

Reference

  • 11,932 nt, 0 ambiguous bases, header matches accession NC_001542.1

Tree

  • Reference is a leaf (among 3,296 tips)
  • genome_annotations coordinates match GFF3
  • No root_sequence field (Nextclade uses reference.fasta)

Generated output

  • Source and data_output/ text files identical; generated pathogen.json differs only by version.tag

Registration

  • Present in data/community/collection.json and data_output/index.json
  • Path follows convention (lowercase/dashes, community/<org>/<pathogen>/<features>)

Nextclade CLI run

  • Reference: 0 mutations, QC "good"
  • Example sequences (21 total):
    • 1 "good"
    • 8 "mediocre"
    • 12 "bad"
    • 20/21 have frameshifts in the L gene (see F3)

Notes

Click to expand
  • N1. Reference is in the tree as a leaf -- no spurious private mutations when users analyze the reference itself.
  • N2. All five CDS structurally clean (div-by-3, ATG start, valid stop). Only the naming is wrong (F1).
  • N3. Clade coverage matches published RABV scheme. The README's exclusion-criteria table (non-monophyletic subclades struck out with rationale) is nice.
  • N4. No public workflow repo (author confirmed). README documents the commands. build_url (github.com/theiagen/rabies) resolves.
  • N5. pathogen.json omits $schema, compatibility, defaultCds. Optional, but defaultCds set to N or G would give a better default protein view.
  • N6. F1 and F3 are still open from maintainer feedback. Consider "experimental": true until those are addressed.
  • N7. Path community/theiagen/rabies/all-clades is immutable after release -- all-clades leaves room for clade-specific datasets later.

Background

Taxonomy, genome, and clade system [click to expand]

Taxonomy and reference strain

Rabies virus is the type species of genus Lyssavirus, family Rhabdoviridae, order Mononegavirales; non-segmented negative-sense ssRNA (Baltimore group V). ICTV binomial: Lyssavirus rabies (formerly "Rabies lyssavirus"); NCBI taxid 11292. attributes.name uses the current binomial.

The reference NC_001542.1 is NCBI RefSeq (11,932 bp), derived from M13215 -- the Pasteur virus (PV) lab strain (Tordo et al. 1986; 1988). Old lab strain, so circulating clades diverge from it (relevant to F3).

Genome organization

Five genes in conserved order 3'-N-P-M-G-L-5' (Wikipedia). The GFF3 labels P as "M1" and M as "M2" -- that's the NC_001542.1 RefSeq's historical product names ("phosphoprotein M1" / "M2 protein"), not modern convention. Standard typing targets: N and G genes (Johnson et al. 2002; Wu et al. 2007); this whole-genome dataset covers both.

Clade / lineage classification

RABV splits into bat-related and dog-related groups with different evolutionary dynamics (Troupin et al. 2016). Dog-related major clades: Cosmopolitan, Arctic, Asian, Africa-2, Africa-3, Indian subcontinent. Bat-related: bat and RAC-SK (raccoon-skunk). The finer subclade labels (AF1a, AM1, SEA1b, ME1a, CE, WE, Vac, ...) are the Pango-style nomenclature from Campbell et al. 2022, integrated with RABV-GLUE. Clade exclusion rationale (non-monophyletic / non-unique subclades) is documented in the README.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants