Reproducing norovirus VP1 alignment/reference system for barcode generation

[adiama0]

Hello,

We are interested in trying to generate a custom barcode for norovirus using Freyja’s BarcodeForge workflow and the published norovirus VP1 dataset as the starting point. To do this, we are retrieving the sequence data from the Nextstrain workflow files and the tree/metadata from the published norovirus Auspice JSON. However, BarcodeForge fails because some positions become multiallelic across lineages during barcode construction.

We are trying to understand whether this is an expected property of the published norovirus VP1 build, or whether it is something introduced by our local VP1 extraction/alignment workflow. In our local workflow, we isolate the VP1 region from the retrieved sequences and align those VP1 sequences with MAFFT before running BarcodeForge. Could you clarify what alignment procedure was used upstream for the published VP1 dataset?

I am happy to share additional details about our workflow if that would be helpful for diagnosing the issue.

Thank you so much for your help!

[joverlee521]

Hi @adiama0,

I'm assuming you are retrieving the following as the starting files:

1. sequences from https://data.nextstrain.org/files/workflows/norovirus/sequences.fasta.zst
2. tree/metadata from https://nextstrain.org/norovirus/all/VP1

The sequences.fasta file is an output for our ingest workflow and these sequences have not been aligned.

Our phylogenetic workflow uses the sequences.fasta file as input to filter and align sequences to different references for each genogroup and the aligned sequences are output in an intermediate file before the final Auspice JSON. Unfortunately, we do not save these intermediate aligned files in a permanent location. If you need to start from the aligned file, I suggest running the phylogenetic workflow locally and use the intermediate files that were produced.

[huddlej]

@adiama0, just to add to @joverlee521's note above, we align norovirus sequences with augur align as shown in the workflow rule below:

norovirus/phylogenetic/rules/prepare_sequences.smk

Lines 99 to 125 in a099a6c

	rule align:
	"""
	Aligning sequences to {input.reference}
	- filling gaps with N
	"""
	input:
	sequences = "results/{group}/{gene}/filtered.fasta",
	reference = "results/defaults/norovirus_reference_{group}_{gene}.gb"
	output:
	alignment = "results/{group}/{gene}/alignment.fasta"
	benchmark:
	"benchmarks/{group}/{gene}/align.txt",
	log:
	"logs/{group}/{gene}/align.txt",
	threads: 1
	shell:
	r"""
	exec &> >(tee {log:q})
	augur align \
	--sequences {input.sequences:q} \
	--reference-sequence {input.reference:q} \
	--output {output.alignment:q} \
	--fill-gaps \
	--nthreads {threads} \
	--remove-reference
	"""

The augur align command is a light wrapper around MAFFT that adds the reference sequence to the input sequences prior to alignment, fills any gaps relative to the reference with Ns, and then removes the reference sequence prior to saving the alignment to disk.

The gene-specific trees like VP1 use a reference GenBank record parsed from the group-specific whole genome GenBank record:

norovirus/phylogenetic/rules/prepare_sequences.smk

Lines 77 to 97 in a099a6c

	rule parse_reference:
	"""
	parsing gene from reference genbank
	"""
	input:
	reference = config['reference']
	output:
	output = "results/defaults/norovirus_reference_{group}_{gene}.gb"
	benchmark:
	"benchmarks/{group}/{gene}/parse_reference.txt",
	log:
	"logs/{group}/{gene}/parse_reference.txt",
	shell:
	r"""
	exec &> >(tee {log:q})
	python scripts/reference_parsing.py \
	--reference {input.reference} \
	--gene {wildcards.gene} \
	--output {output.output}
	"""

The whole genome GenBank record for the "all" group lives here:

https://github.com/nextstrain/norovirus/blob/main/phylogenetic/defaults/all/reference.gb

[adiama0]

Hi @joverlee521 and @huddlej,

Thank you both again for your help on this. We were able to resolve the earlier issue and your clarification on the alignment procedure was very helpful.

At this point we seem to be running into a different problem. The coverage for our samples is consistently low, even after updating our workflow to better match your approach.

I am working on this with @cimendes and trying to understand whether the issue is related to the reference choice, our rooting/ancestral-sequence strategy, or something about the VP1 region of the reads. Have you seen low coverage when working with the norovirus VP1 dataset, or do you have any recommendations on how you would troubleshoot it from here?

If it would be helpful, we would be happy to share more details or set up a brief conversation to discuss the workflow more directly.

[joverlee521]

Prior to aligning, we filter for sequences that have a minimum coverage of 0.8. The coverage comes from a makeshift Nextclade run with the NC_039477_REF GFF file. These coverage values are available in the metadata file at https://data.nextstrain.org/files/workflows/norovirus/metadata.tsv.zst. Particularly, the VP1_coverage column shows only roughly 16% of all norovirus records have coverage >= 0.8.

[huddlej]

@adiama0 @cimendes Does @joverlee521's response above help clarify the issue with coverage? We are both available to meet to discuss your analysis, if you still have questions. You can email hello@nextstrain.org to coordinate on a time.

[adiama0]

Hello, thank you for your suggestions! I am currently running the workflow with the added filter to see if our results improve. I will reach out if we run into the same issue with the modification.