Multi-nucleotide variants not split correctly – potential normalization issue?

[longfei-stat]

Hi team,

Thanks for developing and maintaining this great tool.

I’ve noticed that some nearby variants are not being split into individual SNVs as expected. For example:

15850TA>CG
This appears to be a multi-nucleotide variant (MNV), but I would expect it to be normalized into two SNVs:

• 15850T>C

• 15851A>G

185_GCGAA_ACGAG
This looks like it should correspond to two SNVs:

• 185G>A

• 189A>G

This behavior makes downstream analysis and annotation more difficult. I suspect this might be related to how variant normalization is being handled in the pipeline.

Could you clarify whether these cases are expected behavior? And if not, is there a way to enable proper decomposition of MNVs into SNVs?

Thanks again for your support!

Kind regards,
Fei

[trentzz]

Hi @longfei-stat,

Our normalise module is built around bcftools norm bcftools norm website.

From our testing we found the -m-both option to perform the best.

More info from the website:

-m, --multiallelics -|+[snps|indels|both|any]

split multiallelic sites into biallelic records (-) or join biallelic sites into multiallelic records (+). An optional type string can follow which controls variant types which should be split or merged together: If only SNP records should be split or merged, specify snps; if both SNPs and indels should be merged separately into two records, specify both; if SNPs and indels should be merged into a single record, specify any.

Note that multiallelic sites with both SNPs and indels will be split into biallelic sites with both -m -snps and -m -indels.

Would you be able to run mity call without the --normalise option and then try bcftools norm -f <fasta> -m-both <vcf> with the other options, e.g. -m-any, -m-snps, -m-any instead of -m-both.

If you find that any one of those gets you the results you're looking for, I can make a change to allow users to choose which -m option they want to use for normalisation.

Thank you!
Trent Zeng

[longfei-stat]

Hi @trentzz,

Thanks for the information. I tried all the -m options, but unfortunately none of them produced the desired result. It seems that bcftools norm works well for cases like:

100  CC  C,GG

which gets normalised to:

100	 C	 G 
100	 CC	 C
101	 C	 G

However, it doesn't appear to handle variant decomposition—such as breaking down complex or compound variants into simpler constituent parts. Could this be a limitation of bcftools norm, or possibly an issue related to how FreeBayes represents variants? Perhaps need to explore using a variant decomposition tool?

HG00405.mity.call.vcf.gz
Here I attached the data I tested with the 15850TA>CG.

Thanks,
Fei

[trentzz]

Hi @longfei-stat,

Thanks, that's good information for us to work off of. I think we may need some time to get back to you to investigate bcftools norm, the output of freebayes and decomp tools more.

Thanks!
Trent Zeng

[anjalijain22]

Hi @trentzz and @longfei-stat!
I was also encountering this same issue with mity's normalise command. After doing a bit of digging, I found that Freebayes can generate multinucleotide polymorphisms (MNPs). My understanding is that bcftools norm decomposes multiallelic variants but not MNPs. In the earlier version of mity (v0.4.0) this wasn't an issue since it was using a custom script which would decompose these MNPs. To solve this issue, I added an extra step in our WGS pipeline which uses vt decompose_blocksub on mity's normalized VCF to decompose biallelic block substitutions into its constituent SNPs. For example:

MNP:
MT-51-CT >TC

Constituent SNPs:
MT 151 C T
MT 152 T C

I thought I'd share this as a possible solution while you look into other options. Happy to provide more details if needed!