Some packages contain incongruent SNPs compared to the rest of the archive

[nevrome]

The new improved detection of incongruent SNPs and strand flips in trident (poseidon-framework/poseidon-hs#375) indicates that some packages have different alternative or reference bases than the others for individual SNPs. I detected this with trident validate -d . --forgeTest --strandCheck, which I ran incrementally, deleting the offending package after each run, to explore the whole SNP space.

With this I detected incongruent SNPs in the following packages:

• 2024_Ravasini_Picene (added by @stschiff in Adding 2024_Ravasini_Picene #238)
• 2021_SaupeCurrBiol (added by @sc13-bioinf in add 2021_SaupeCurrBiol #82)
• 2024_Allentoft_PostGlacial_WestEurasia (added by @smpeltola in Add package for Allentoft at al. 2024 #202)

It seems these all include author-provided genotype data, correct? And they differ to all the other packages, and to each other. This is not the result of strand flips, as these are automatically addressed with --strandCheck.

In the future it will only be possible to trident forge these packages with --skipIncongruentSNPs, so that the relevant SNPs are dropped entirely from the output dataset. Currently this is already happening silently.

Is there anything we can and should do about this? How does this even happen?

[stschiff]

Yes, very annoying. I think one simple check would be to directly compare the BIM files of these packages with a bona-fide package and check how many of these conflicting SNPs there are. If it's just a handful, I would suggest to simply remove them from the offending packages, which would then result in missing data at those SNPs upon forging.

One possible explanation here are triallelic SNPs, so SNPs with a third allele at a very low minor allele frequency (which might also be an error due to damage or so).

[nevrome]

I used 2025_Gretzinger_SlavicPeriod as the suggested bona-fide package and classified SNPs in the three problemantic packages based on how they behave in relation to it:

category	dataset	nr_snps
not present in gretzinger	allentoft	1720
not present in gretzinger	ravaisni	0
not present in gretzinger	saupe	525
identical	allentoft	0
identical	ravaisni	257018
identical	saupe	176087
switched	allentoft	1075557
switched	ravaisni	544916
switched	saupe	307054
probably identical	allentoft	0
probably identical	ravaisni	39165
probably identical	saupe	115600
probably switched	allentoft	0
probably switched	ravaisni	356169
probably switched	saupe	576631
ambiguous palindromic	allentoft	0
ambiguous palindromic	ravaisni	0
ambiguous palindromic	saupe	0
fixable strand flip	allentoft	0
fixable strand flip	ravaisni	0
fixable strand flip	saupe	1
fixable strand flip + switched	allentoft	0
fixable strand flip + switched	ravaisni	1
fixable strand flip + switched	saupe	1
unresolved allele mismatch	allentoft	169
unresolved allele mismatch	ravaisni	29050
unresolved allele mismatch	saupe	10683

• identical: Same two alleles in the same order.
• switched: As identical, but in reversed order.
• probably identical/switched: As identical/switched, but N acts as a wildcard.
• ambiguous palindromic: SNP is A/T or C/G palindromic and strand orientation cannot be determined.
• fixable strand flip: SNP becomes identical after applying the strand complement.
• fixable strand flip + switched: As fixable strand flip, but in reversed order.
• unresolved allele mismatch: SNP can not be reconciled with any of the steps above.

For allentoft we could do what you suggest and just kick out the 169 unfixable SNPs. But for ravasini and (the really odd) saupe that's probably not the best approach.

Edit 2026-06-12: Fixed the count for the saupe paper, I had not considered 0 for N, as pointed out by @TCLamnidis.

[stschiff]

Hmm, very interesting! We should look at this in the next Poseidon meeting, ideally together with Thiseas.

[nevrome]

I had a look into the Ravisini paper with Thiseas. The authors state the following:

Variant calling for autosomes was performed with ANGSD-0.917 [72]. We called haploid genotypes sampling a random base (-doHaploCall 1) for each position that is present in the 1240K SNPs panel using the -sites options. In addition, we specified the major and minor alleles as they are indicated in the 1240K SNPs panel with the option -doMajorMinor 3. The function haploToPlink was used to convert the .haplo file, resulting after the SNP calling, into Plink [82] format files.

It's unlikely that the resulting data would deviate that strongly from the "norm", if we understand angsd -doMajorMinor 3 correctly. How did you obtain the genotype data, @stschiff? Did you do any data transformation that could explain this?

[stschiff]

I did not do any data transformation. The authors sent me the plink data and I included it for this package. I would like to understand what these unresolvable SNP pairs are. Will investigate a bit now.

[stschiff]

OK, so I can confirm that these are really the wrong alleles in some positions. I don't know why. I looked back, and my email trail to the authors just includes the email where they sent me the files, but they were sent through a download link in December 2024, and that link has of course expired. But I definitely haven't touched the file, so I think the erroneous alleles must have entered through them.

Not sure what to do about it...

[nevrome]

@TCLamnidis was so kind to do some research on 2021_SaupeCurrBiol.

The data can be downloaded here: https://evolbio.ut.ee/Saupe_2021. I assume that @sc13-bioinf got it from there, but maybe it's worth checking if switching to this would resolve the issue right away. (Edit: I did, and .bim and .bed from that page are identical to the ones we have in the archive. So no easy solution here.)

There is some documentation on how the data was processed: https://evolbio.ut.ee/master_file_creation_log.txt
Thiseas pointed out this interesting detail:

A few additional strand flips were necessary at each step. Importantly, because PLINK merges monomorphic markers from different strands incorrectly (0-C 0-G gets to be merged as G-C and no strand conflict is reported), we flipped the strand of such (with strand conflict checked for outside of PLINK) monomorphic markers in the dataset to be merged prior to merger.

This sounds like another way how these allele combinations can emerge in the first place, besides triallelic SNPs.

[stschiff]

Oh dear... yes indeed, that's scary, actually... I still don't know what exactly to do about it.