InvalidReferenceBases error for sequences present in GRCh38

[yuliamostovoy]

I'm getting a few of these errors:

 2024-05-23T15:58:27.206Z ERROR kanpig::kplib::vcfreader > skipping invalid VCF entry Custom { kind: InvalidData, error: InvalidReferenceBases(InvalidBase(Invalid('R'))) }
 2024-05-23T15:58:35.842Z ERROR kanpig::kplib::vcfreader > skipping invalid VCF entry Custom { kind: InvalidData, error: InvalidReferenceBases(InvalidBase(Invalid('Y'))) }
 2024-05-23T15:58:35.843Z ERROR kanpig::kplib::vcfreader > skipping invalid VCF entry Custom { kind: InvalidData, error: InvalidReferenceBases(InvalidBase(Invalid('Y'))) }

I've checked and these Rs and Ys are just taken from GRCh38, so it's not a VCF issue. It might be good to add support for these ambiguous bases in future updates.

[ACEnglish]

This error is thrown by noodles-vcf. That package has been updated recently so may not throw the error anymore. However, the updates completely changed the API and some features kanpig uses are not longer available, so it's a huge refactor to pull it in.

I will put this on the development roadmap.

[ACEnglish]

The new version of noodles-vcf did indeed fix this issue. This change is in develop and staged for v0.2.1

Note that the behavior for kmer-featurizing non-ATCG characters is to assume the nucleotide is an 'A'. I'll need to document this somewhere in-case it becomes relevant to some users' debugging of their formerly invalid reference bases.

[ACEnglish]

So, here's a fun update.... tl;dr; is the problem isn't fixed and these invalid reference bases will now cause a corrupted, unusable VCF output.

noodles-vcf can now read IUPAC codes in the REF, but after further testing, I found that it cannot write them. Looking into why noodles can't handle this I found it validates reference bases according to the vcf spec § 1.6.1.4 "Each base must be one of A,C,G,T,N (case insensitive).".

This means it technically is a VCF issue and these lines can be considered invalid. The update to the new noodles-vcf means that instead of throwing out the invalid line, kanpig is asking noodles to write an invalid line, noodles is raising an error, and then writing a corrupted vcf. Not just corrupted in that the IUPAC code is in REF, but noodles stops halfway through writing the line and produces e.g.

chr20	24408074	.	GTATTTCCAGCCTTGGCCATAGchr20	24408077	.	TTTCCAGCCTTGGCC

Somehow fixing the problem just made it worse. Now I have to figure out what to do next. The options are:

• go back to checking reference bases during reading and throw out 'invalid' entries
• correct reference bases during reading/writing according to the vcf spec (pasted below)

Option number two is likely more seamless for the user-experience, but I don't know if data should be "sliently" changing since I can imagine a user asking "Where is my Y?". It also creates a little extra overhead checking every base.

Once the option is picked, I need to figure out if kanpig should be doing it or I put in a noodles-vcf pull request.

vcf spec for how to handle IUPAC in REF

If the reference sequence contains IUPAC ambiguity codes not allowed by this specification (such as R = A/G),
the ambiguous reference base must be reduced to a concrete base by using the one that is first alphabetically
(thus R as a reference base is converted to A in VCF.)

[ACEnglish]

These vcf entries will now be parsed and fixed according to vcf spec v4.4 in the output. A warning is printed to the user if any iupac bases were changed.