H3N2 sequences classified as J.2.4 in Nextclade appear only as J.2 in my Nextstrain tree

[RasmusKoRiis]

#Description / Problem Statement
When I analyze certain H3N2 sequences with Nextclade, they are classified into subclade J.2.4.
However, when I build a tree with Nextstrain, those same sequences are shown only as J.2 (the parent clade).

This looks like a mismatch between the clade definitions used in Nextclade vs those used in the Nextstrain seasonal influenza build.

#Steps to Reproduce / Workflow

I have used the both nextclade in CLI and online version to get J.2.4 classification of my H3N2 seqeunces.

I then run the same H3N2 seqeunces with references using the GISAID-workflow in your nextstrain (seasonl flu) pipeline and get the seqeunces classified as J.2.

#Environment / Versions

Nextstrain CLI: 10.2.1.post1

OS: Linux (Ubuntu-based)

Expected Behavior
I expected the Nextstrain tree build to classify these sequences as J.2.4, consistent with Nextclade.

Actual Behavior
In the Nextstrain tree, they appear only as J.2.

[huddlej]

Thank you for reporting this issue, @RasmusKoRiis! The seasonal-flu workflow downloads both clade and subclade annotations from the influenza-clade-nomenclature repository for each subtype and gene segment combination (for example: this is the H3N2 HA repo) and then runs augur clades with those downloaded files. This implementation means that the downloaded copies of the clade definitions can get stale without Snakemake realizing it. There are other possible reasons we could get mismatched results between the seasonal-flu workflow and Nextclade, but this seems like the simplest explanation to check first.

Can you share the contents of the last 5 lines of your local copy of config/h3n2/ha/subclades.tsv? The lines should look like this:

J.3	HA2	176	I
J.4	clade	J
J.4	HA1	173	R
J.4	HA1	276	E
J.4	SigPep	3	T

If you don't see J.4 at the end of the file, you can force the workflow to redownload the clade definitions by adding the --forceall Snakemake flag to the end of the command from the quickstart guide like so:

nextstrain build . --configfile profiles/gisaid/builds.yaml --forceall

[huddlej]

As a separate note for the rest of the Nextstrain team, @joverlee521 pointed out that we should be able to fix this issue of clades files getting stale by defining an input for the rules that download clades that uses the shared path_or_url function. If we choose that solution, we'll need to setup vendoring of the shared repo in this seasonal-flu repo.

[RasmusKoRiis]

Thank you for your response @huddlej , my config/h3n2/ha/subclades.tsv file seems okay with this at the end of the file:

J.3 clade J
J.3 HA2 176 I
J.4 clade J
J.4 HA1 173 R
J.4 HA1 276 E
J.4 SigPep 3 T

[huddlej]

Ah, I'm sorry, @RasmusKoRiis, I should have provided another test than checking the end of the file; the J.4 clade is older than J.2.4 which appears in the earlier in that file. Do you also see the following lines in your subclades file?

J.2.4	clade	J.2
J.2.4	HA1	189	R
J.2.4	HA1	135	K
J.2.4	nuc	469	A
J.2.4	nuc	1568	A

If you do, my next best guess is that your input sequences don't reflect enough of the genetic diversity required to label the J.2.4 clade. The following screenshot shows an example of how this can happen with a subset of 19 recent H3N2 HA sequences colored by the subclade assignment and with branch labels showing amino acid substitutions. The clade at the top of the tree has J.4's 189R and 135K substitutions and also the 469A and 1568A nucleotide substitutions, but it doesn't get labeled as J.2.4:

This happens because augur clades assigns the J.2.4 nodes to a more ancestral clade (in this example, A.2) and assigns the more ancestral nodes to J.2. The algorithm notes that A.2 is not nested inside J.2, so it overwrites the A.2 assignments with J.2. This bug appears because augur clades relies on the inferred ancestral sequences in the tree to assign clades. When there isn't enough of the ancestral diversity associated with all clade assignments in the input data, the algorithm assigns spurious labels.

Nextclade does not rely on inferred ancestral sequences, so it isn't affected by the amount of genetic diversity in your input sequences. Instead, it takes an existing tree structure whose nodes have already been labeled with clades, maps each input sequence to the closest sequence in the tree, and assigns the input sequence the clade label of that closest sequence.

The long-term solution to this issue is to use Nextclade to assign clade labels in the phylogenetic workflow instead of augur clades.

The short-term solution for your analysis could be to include ~20-30 more sequences from the last six months that represent the diversity of currently circulating clades. Adding these additional contextual sequences to your analysis should allow the workflow to properly infer ancestral sequences that are needed to label recent clades.

Here is what I get when I rerun the analysis above with 82 genomes instead of 19:

[RasmusKoRiis]

Thanks @huddlej for the detailed explanation.

Yes — my subclades.tsv does include the lines you mention:

J.2.4 clade J.2 J.2.4 HA1 189 R J.2.4 HA1 135 K J.2.4 nuc 469 A J.2.4 nuc 1568 A

In my tree, one sample does get J.2.4, but most of the recent sequences collapse to J.2, while those same sequences are called J.2.4 by Nextclade. I included some older strains and a set of HA sequences from the last year from Norway to try to capture the diversity, but I still see the collapse.

Given that context, do you have a sense of why the subset might still be assigned to J.2 rather than J.2.4?

I’ve attached a screenshot of my current tree for reference.

[huddlej]

Thank you for these screenshot examples, @RasmusKoRiis. Based on that second figure, I suspect the issue is still that there aren't enough J.2.4 samples to get ancestral sequences that augur clades can properly label the nodes you'd expect. I can also see how this type of problem is harder to work around by the "add more samples" approach when you want to focus on a specific geographic location like Norway where the diversity might be limited compared to global diversity.

The workflow already generates the Nextclade subclade annotations and stores them in the metadata TSV that we pass to augur export to generate the final Auspice JSON. The only reason those Nextclade subclade annotations don't appear in the Auspice JSONs by default is that the augur clades node data JSON defines a subclade attribute that takes precedent over the metadata field of the same name.

I just pushed a feature branch to this repository called export-nextclade-subclades which makes a copy of the Nextclade subclade column (named subclade_nextclade) and adds that new column to the Auspice config JSON for H3N2 HA. If you run the GISAID quickstart guide from this branch with the 19 recent sequences I showed above, you get annotations from both augur clades and from Nextclade. The screenshot below shows the scatterplot of augur clades output on the x-axis and Nextclade outputs on the y-axis.

Would you be willing to fetch this branch and try running your workflow to see if you get the desired clade annotations with it and your data? If this seems to work well enough for you (given the caveats that you don't get internal node colorings or branch labels), I could formalize this approach in the workflow and that will at least resolve this issue in the short term.

[RasmusKoRiis]

@huddlej Thank you for this! I have tried out the new branch with my data (see images) and as far as I can see the strains are now annotated with the correct subclade under "Subclade "Nextclade".