Meaning of warning

[josephhughes]

Hi,

What does the following warning message mean:

WARNING: Could not assign a node to clade 6 with 86 samples since placement node(s) did not overlap with enough clade samples!

Thanks in advance

[josephhughes]

I'm trying to build a backbone phylogeny of HCV with genotype and subgentotype annotations. I don't understand why I am getting the following error:
WARNING: Could not assign a node to clade 1a with 6 samples since placement node(s) did not overlap with enough clade samples!

Especially as it is clear from the tree that the sub genotype is monophyletic.

I faced similar problems when trying to annotate at the genotype level.

Any advice will be most welcome.

[AngieHinrichs]

Hi Joseph, the warning means that matUtils annotate's algorithm for finding the best node to match your representative samples does not find a node that passes its thresholds for --allele-frequency (by default, a frequency of at least 0.8 in representative samples means that a mutation is required for candidate nodes) and --set-overlap (by default, at least 0.6 * number of representative samples should be descendants of candidate nodes), influenced by --mask-frequency and --clip-sample-frequency. The 1a example does look like it has an easily identified best node, provided there aren't some stray 1a's elsewhere in the tree? (Or perhaps that node includes a mutation that is reverted in a significant proportion of samples?)

I'd need more information to tell what's really going on. If you can share the exact matUtils annotate command that you ran, and the input files, then I can take a look after I return from vacation mid next week. Alternatively, if you want to drill down into the algorithm and look at the parameter values measured on each candidate node then you can run with --write-details details.tsv and try to relate those values to the other output (such as which nodes were identified as candidates) and to the threshold values.

[josephhughes]

Hi Angie,
I've tried to understand what is going on with no luck. I'm getting the same issue with a different virus alignment and tree which I am sharing with you.
ref_padded_aln.fa is a FASTA alignment of the reference sequences (not all are complete genomes)
ref_padded_aln_midpoint.treefile is a Newick tree of the reference sequences
major_clade_and_accession_noNULL.tsv is a TSV with the major clade assignments
minor_clade_and_accession_noNULL.tsv is a TSV of the minor clade assignments
query_ids_and_ref_padded_aln.fa is a FASTA file with my query sequences to assign clades to.
query_clades.txt are the clade assignments I am expecting

This is the code I am using where they all get assigned to Cosmopolitan major clade and undefined minor clade:

faToVcf -ref=NC_001542 ref_padded_aln.fa ref_padded_aln.vcf

usher --vcf ref_padded_aln.vcf --tree ref_padded_aln_midpoint.treefile -o RABV_backbone.pb

matUtils annotate -i RABV_backbone.pb -c major_clade_and_accession_noNULL.tsv -o RABV_backbone_major.pb --write-details major_details.tsv
 
matUtils annotate -i RABV_backbone_major.pb -c minor_clade_and_accession_noNULL.tsv -o RABV_backbone_major_and_minor.pb --write-details minor_details.tsv
 
#converting the query sequences to vcf

faToVcf -ref=NC_001542 query_ids_and_ref_padded_aln.fa query_aln_to_NC_001542.vcf

#getting the clade assignments for the query
usher -i RABV_backbone_major_and_minor.pb -v query_aln_to_NC_001542.vcf -u -d query_output/

I'd be very grateful for any advice.

query_clades.txt

query_ids_and_ref_padded_aln.fa.txt

minor_clade_and_accession_noNULL.tsv.txt

major_clade_and_accession_noNULL.tsv.txt

ref_padded_aln_midpoint.treefile.txt

ref_padded_aln.fa.txt

[josephhughes]

Here is the data for the HCV genotypes:

genotype_Accession.txt
HCV_aln_to_master.fa.txt
query_aln_to_AF009606.fa.txt
subgenotype_Accession.txt
HCV_ref_midpoint.treefile.txt

The following is the code to run the genotype assignment:

#get variants of all sequences relative to my master ref
faToVcf -ref=AF009606 HCV_aln_to_master.fa HCV_aln_to_master.vcf

#building the protobuf
usher --vcf HCV_aln_to_master.vcf --tree HCV_ref_midpoint.treefile -o HCV_ref_backbone.pb

#creating a protobuf with the genotype information
matUtils annotate -i HCV_ref_backbone.pb -c genotype_Accession.txt -o HCV_ref_backbone_genotype.pb
 
#creating a protobuf with the subgenotype information
matUtils annotate -i HCV_ref_backbone_genotype.pb -c subgenotype_Accession.txt -o HCV_ref_backbone_genotype_subgenotype.pb

#converting the query sequences to vcf
faToVcf -ref=AF009606 query_aln_to_AF009606.fa query_aln_to_AF009606.vcf

#getting the clade assignments for the query
usher -i HCV_ref_backbone_genotype_subgenotype.pb -v query_aln_to_AF009606.vcf -u -d query_output/

[AngieHinrichs]

I've worked through your RABV example and found that all of the major clades and most of the minor clades can be annotated if you add these options to both matUtils annotate commands: -m 0.01 -f 0.95 (more about that below). The minor clades YUGFOX, YUGCOW, Vac, LX, LS, LB1 can't be annotated because there is only one sequence for each of them in the tree, and the matUtils annotate algorithm for -c requires a node with multiple descendant sequences. Adding more sequences from those minor clades to the initial tree should help.

BTW KX148261 is associated with both SEA3 and SEA4 in minor_clade_and_accession_noNULL.tsv -- matUtils annotate didn't complain, but it caused trouble for the tree viewer Taxonium so I omitted KX148261 when annotating minor clades. (KX148261 is on the branch annotated SEA4, so that's probably the correct assignment.)

Both major_clade_and_accession_noNULL.tsv and minor_clade_and_accession_noNULL.tsv have many duplicated lines, which caused strange numbers in some of the statistics about how many sequences were matched, so I also removed duplicate lines.

When sequences in query_aln_to_AF009606.vcf were added to the tree, they were placed in a cluster around NC_001542. Here's an excerpt from the file query_output/uncondensed-final-tree.nh:

((NC_001542:33,
  (AB246203:3,
   (AB246204:0,
    (AB010496:26,
     (AB201802:21,AB201803:12)node_402:156)node_400:15)node_399:5)node_398:22)node_397:3,
 (AB010495:27,
  (AB206408:0,AB206407:0)node_396:74)node_401:7)node_395:447

I don't know why they were all placed around NC_001542 instead of dispersed around the tree as one would expect from your query_clades.txt file. (Some of those nodes have a large number of mutations annotated on them, too.) Was the alignment/padding method that generated query_ids_and_ref_padded_aln.fa exactly the same as the method that generated ref_padded_aln.fa?

However, since they were all placed around NC_001542, I can explain why query_output/clades.txt has Cosmopolitan for their major clade and "UNDEFINED" for their minor clade. NC_001542 belongs to major clade Cosmopolitan and minor clade Vac -- but Vac is one of the minor clades that could not be annotated (because there is only one sequence in the original tree for Vac). So while the tree has the Cosmopolitan annotation for NC_001542 and the new sequences placed around it, it does not have any minor clade annotation.

Now, how did I get to -m 0.01 -f 0.95? Ultimately, trial and error, but here's a little background. The default parameter values (-m 0.2 -f 0.8) were tuned for SARS-CoV-2 Pango lineages and the full SARS-CoV-2 UShER tree in early 2021. The Pango lineage designated sequences were selected using a tree built by a different algorithm and didn't match up with branches of the UShER tree as neatly as your clade representatives match up with branches of your tree. Therefore the matUtils annotate algorithm and parameters had to tolerate a fair amount of disagreement. Also, the annotate algorithm uses the core usher placement algorithm to identify candidate nodes; in order to find the right candidate node, the algorithm needs to construct the right set of representative mutations from all of the representative sequences. So the -m and -f parameters, which are used to determine which positions to ignore and which mutations to require, have defaults that allow a fair amount of 'slosh'. (And then there's the relatively low divergence of SARS-CoV-2 sequences compared to most pathogens that have been transmitting for a long time, and the fact that the SARS-CoV-2 RefSeq is so close to ancestral that it can be used as the root of the tree.) But when applying the same algorithm to a tree with fewer sequences, many more mutations, and closer agreement of representative sequences and tree structure, that 'slosh' is counterproductive, and the parameters need to be made more stringent in order to find the correct node.

So, for the trial and error part, I started with the major clade annotation step. I first tried increasing -f from the default of 0.8 to 0.9 and 0.95 -- but by itself, that caused no improvement. So I tried decreasing -m to 0.1, which rescued Bats but not Asian. I kept decreasing -m, and the node found for Asian got a little closer to what it should have been, but was still too deep in the branch. So I kept the low -m value and tried increasing -f, and finally with -m 0.01 -f 0.95 it got all major clades. Then I tried using the same values for the minor clade step, and it got all of the clades except for the ones with only one sequence each.

I'm still not sure what's going on with the placement around NC_001542 instead of across the tree.

[AngieHinrichs]

I'm still not sure what's going on with the placement around NC_001542 instead of across the tree.

As a simple experiment, I tried generating a VCF that renames all of the original sequences, to get the effect of adding duplicate sequences:

sed -re '/^#CHROM/  s/([A-Z]{2}[0-9]{6})/\1_dup/g' ref_padded_aln.vcf > ref_padded_aln_dup.vcf

Then I ran usher with ref_padded_aln_dup.vcf to see where the renamed/duplicate sequences were placed in the uncondensed-final-tree.nh: would they be placed in a new tree around NC_001542, or would each duplicate be placed with its original?

usher -i RABV_backbone_major_and_minor.pb -v ref_padded_aln_dup.vcf -u -d dup_output/

Looking at dup_output/uncondensed-final-tree.nh, it seems that each duplicate is placed with its original (not in a new clump around NC_001542) which is good. So I think that implies there's an issue with query_ids_and_ref_padded_aln.fa (though I still can't identify where it's going wrong).

[AngieHinrichs]

-m 0.01 -f 0.95 seems to help HCV too, except that it has many subgenotypes with only one sequence in the tree, so the tree definitely will need more sequences to capture all the subgenotypes.

[josephhughes]

HCV genotyping is working well now and will work even better once I have added some additional subgenotype sequences.
For RABV, I need to investigate the alignment further. We are building this in a non-standard way as we have selected a range of phylogenetically informative references as a starting point and then query sequences are aligned using nextalign to the reference chosen by the top blast hit. Each sub alignment is then collated and padded according to the original alignment of all the references. It is possible that some informative sites are being lost in subsequent steps as a result.
Thanks for all you help with this.

[AngieHinrichs]

Glad to hear it's working for HCV! Using multiple references for RABV sounds challenging but doable. An UShER tree needs a single reference, which I guess would be the consensus of the orginal alignment of all references? (Are the sequences not alignable to that constructed reference?)