Charles R. ,

Thanks for your post, which relates to recent, key findings in genomics. First, to be thorough, I have to acknowledge the possibility that a variety of technical issues could contribute to findings like what you have described, e.g. an early cycle PCR amp error could give you reads precisely half of which have a C the other an A. However, for the purposes of this post, I will assume that the quality metrics like those you provide are accurate, and that the read such stats (e.g. quality score) in fact perfectly reflects the primary sequence of the gDNA.

The problem is that, even if those statistics are entirely accurate, the data may still prove very hard to interpret correctly due to specific sequence characteristics of ChrY (as well as other chromosomes). To explain this, will need a bit of background:

Both the X and the Y chromosomes have unique sequence properties that have complicated their analysis - and thus also our understanding of the variation in them, until very recently. The completion of the T2T project and creation of the T2T v2 assembly concluded with the inclusion of chrY to T2T CHM13 v1.3 which with its addition then became T2T CHM13 v2, the best assembly ever until T2T - merged with the Human Pangenome Reference Consortium (HPRC). At last count they have several hundred gapless, phased human haplotypes, and as such are well on their way to the first human pangenome reference.

ChrY in particular displays overlapping as well as distinct types of genomic fluidity than autosomes as the forces imposed by homologous recombination do not act on autosomes and sex chromosomes in precisely the same way. In fact, ChrY has several distinct types of variation that are very hard to deal with using short read (next-gen) sequencing data. Thus, while GRCh38 was >90% complete in other chromosomes, about half of the sequence for ChrY was missing from even the most recent references based on short-read sequencing data, e.g. GRCh38.p13 (!!).

Specifically, complex genomic structures including satellite regions, long palindromic sequences, tandem repeats, and segmental duplications remained opaque to NGS for the entire ~20 year period between completion of HGP and T2T-CHM13-v2.0. To understand why NGS fails in these regions, let's take just one example: segmental duplications, or SDs. SDs are repeat regions of 10-200kB that bear very high homology to one another (95-99.9%), and arise due to specific genomic forces, e.g. the presence of inversion events (porubsky et al 2022).

This ultra high homology alone is a formidable a challenge for short read sequencing technologies, however, additional phenomena, such as gene conversion events, further influence polymorphism at these loci and have greatly complicated their analysis. While with the advent of 3rd generation sequencing techniques (SMRT, nanopore) we now have the ability to obtain phased, gapless assemblies of these regions, in reality, the praxis of nearly all of genomic science lags behind this. What I mean is, insofar as the vast majority of sequencing data are still generated using NGS (short-read, 2nd generation) sequencing, most datasets are still unable to resolve these regions.

We are now at last in a position to address your question.

There are two broad classes of structural variants, copy number neutral SV and copy number altering SV. Segmental Duplications are the later type - they increase copy number at a locus. So, while a given human individual may have only one chr Y, this in no way implies that they have only one copy of a given sequence. Rather, consider these two statements:

1) SDs increase copy number at a locus/loci 2) Polymorphisms account for the 0.1-5% their sequence, but they are otherwise identical.

These two conditions, together, strongly suggest that NGS-based sequencing of two homologous SDs will produce reads that are identical except for one or a few NTs. In this case, even in high-quality sequencing data, this could potentially result in 2 different base calls at a locus that, in the reference being used, only has a copy number of 1.

Actually, though, we now know that having just two copies of an SD is not all that surprising. Recent evidence spanning several fields suggests that SDs evolved recently in the speciation of humans from great apes in order to fuel the evolutionary thirst for greater and greater intelligence (as well as other reasons). In this case, SDs have increased copy number of certain genes/gene families from ~2 to ~30 in just 200,000 (200kYa). In loci like these, a person might not only have two alleles, {she, he} might have MANY. This could show up in a .vcf file in a variety of ways, but detailed studies of alignment of such reads shows that prior assembly structure was a contributor to inaccurate read mapping in cases when a SD or other SV had not yet been ascertained.

The question to consider is, "but if that were true, how would it show up in my data?" Simply put, the case you mention is a possible example. In other words, it is possible that your data contain two reads that map to the same portion of two SDs, even with these characteristics:

Allele quality of a SNP is over 400. The Allele Freq. 0.5, Depth 30, read strands even balanced, Genotype quality 99.

Just one last thing to mention. The reason why these results may be counterintuitive is because we are used to thinking of ourselves as diploid, and because we are not accustomed to thinking about what variant calling accuracy metrics mean in the context of a reference genome that is itself flawed. Metrics like quality score, depth etc. provide information concerning the level of accuracy of the nucleotides on a read and the number of reads generated. They aren't, however, useful for predicting whether the reference genome being used is itself accurate in such an area.

To use different words, what I mean is that having an accurate variant call in such a locus tells you nothing about whether that read maps to one SD, another SD, etc. Thus, you have to be really careful about what you conclude from reads in these regions, irrespective of their quality scores. Why? Because when such reads map to structurally heteromorphic loci (like SDs) their quality scores lose meaning insofar as they neglect higher-level tasks like genotype phasing, haplotype resolution, or read mapping, which ultimately influence most applications clinicians and scientists have for such information...