Haplotype-aware inference of human chromosome abnormalities

Most contemporary implementations of PGT-A are based on low-coverage (<0.05×) whole-genome sequencing of 5 to 10 trophectoderm cells biopsied from blastocyst-stage embryos at day 5, 6, or 7 postfertilization (15, 16). Given the level of coverage and the paucity of heterozygous sites within human genomes, standard approaches for analyzing such data typically ignore genotype information and instead infer aneuploidies based on deviations in relative coverage across chromosomes within a sample. Inspired by genotype imputation (35–40) and related forensic methods (41), we hypothesized that even low-coverage data could provide orthogonal evidence of chromosome abnormalities, complementing and refining the inferences obtained from coverage-based methods. As in imputation, the information content of the genotype data is greater than first appears by virtue of linkage disequilibrium (LD)—the population genetic correlation of alleles at two or more loci in genomic proximity, which together compose a haplotype (42). We developed an approach to quantify the probability that two sequencing reads originated from the same chromosome, based on known patterns of LD among alleles observed on those reads. Hereafter, we refer to our method as LD-informed PGT-A (LD-PGTA).

Specifically, in the case of trisomy, we sought to identify a signature of meiotic error wherein portions of the trisomic chromosome are composed of two distinct homologs from one parent and a third distinct homolog from the other parent (BPH; Fig. 1). In contrast, trisomies arising via postzygotic mitotic errors are composed of two identical copies of one parental homolog and a second distinct homolog from the other parent (SPH). Although this chromosome-wide SPH signature may also capture rare meiotic errors with no recombination (Fig. 1), BPH serves as an unambiguous signature of deleterious meiotic aneuploidy. In the presence of recombination, meiotic trisomies will comprise a mixture of BPH and SPH tracts. BPH tracts that span the centromere are consistent with missegregation of homologous chromosomes during meiosis I, while BPH tracts that lie distal to the centromere are consistent with missegregation of sister chromatids during meiosis II (43).

Distinguishing these scenarios in low-coverage data (¡1×) is challenging due to the fact that reads rarely overlap one another at sites that would distinguish the different homologs. Our classifier therefore uses patterns of allele frequencies and LD from large population reference panels to account for potential relationships among the sparse reads. Specifically, the length of the trisomic chromosome is divided into genomic windows, on a scale consistent with the length of typical human haplotypes (10${}^4$ to 10${}^5$ bp). For each genomic window, several reads are sampled, and the probabilities of the observed alleles under the BPH trisomy and SPH trisomy hypotheses (i.e., likelihoods) are computed.

Our likelihood functions are based on the premise that the probability of drawing reads from the same haplotype differs under different ploidy hypotheses, which are defined by various configurations of genetically distinct or identical homologs (Fig. 1). If a pair of reads originates from two identical homologs, the probability of observing the alleles on these reads is given by the joint frequency of the linked alleles. On the other hand, if a pair of reads originates from two different homologs, the probability of observing their associated alleles is simply equal to the product of the allele frequencies, since the alleles are unlinked. Because the BPH and SPH hypotheses are defined by different ratios of identical and distinct homologs, the probability of sampling reads from the same homolog also differs under each hypothesis (1/3 and 5/9, respectively; Fig. 1).

Allele frequencies and joint allele frequencies (i.e., haplotype frequencies) are in turn estimated through examination of a phased population genetic reference panel, ideally matched to the ancestry of the target sample. Such estimates have the virtue that they do not depend on theoretical assumptions, but simply on the sample having been randomly drawn from the genetically similar populations. A drawback, which we take into account, is that reliable estimates of the frequencies of rare alleles or haplotypes require large reference panels. In practice, it is sufficient to construct the reference panels using statistically phased genotypes from large surveys of human genetic diversity, such as the 1000 Genomes Project (n = 2,504 individuals).

Given the importance of admixture in contemporary populations, we also generalized our models to allow for scenarios involving admixture in the parental generation, as well as more distant and complex admixture scenarios (SI Appendix). These admixture-aware models require only knowledge of the ancestry of the target sample (i.e., embryo), which has a practical advantage for PGT-A where parental genotype data are typically unavailable.

Within each genomic window, we compare the likelihoods under the BPH and SPH hypotheses by computing a log-likelihood ratio (LLR) and use the m out of n bootstrap procedure (44) to assess uncertainty (Materials and Methods and Fig. 2). LLRs are then aggregated across consecutive genomic windows comprising larger intervals (i.e., “bins”). The mean and variance of the combined LLR are used to compute a confidence interval for each bin, which can then be classified as supporting the BPH or SPH hypothesis, depending on whether the bounds of the confidence interval are positive or negative, respectively. Confidence intervals that span zero remain inconclusive and are classified as “ambiguous.”

To evaluate our method, we simulated sequencing reads from BPH and SPH trisomies according to their defining haplotype configurations (Fig. 1). Our simulations assumed uniform depth of coverage, random mating (by randomly drawing haplotypes from the 1000 Genomes Project), and equal probability of drawing a read from any of the homologs. We varied the mean depths of coverage (0.01×, 0.05×, and 0.1×) and read lengths (36 to 250 bp), testing model assumptions and performance over a range of plausible scenarios.

To benchmark and optimize LD-PGTA, we developed a generalization of the receiver operating characteristic (ROC) curve (45) for scenarios that include an ambiguous class (i.e., bins with a confidence interval spanning zero), which we hereafter denote as the “balanced” ROC curve. For a given discrimination threshold, a balanced true (false) positive rate, denoted as BTPR (BFPR), is defined as the average of the true (false) positive rate of predicting BPH and SPH trisomy (SI Appendix, Fig. S1 A and B). The balanced ROC curve thus depicts the relationship between the BTPR and BFPR at various confidence thresholds.

We generated balanced ROC curves for bins within simulated trisomic chromosomes (BPH and SPH) for samples of nonadmixed ancestry and a correctly specified reference panel (Fig. 3). At depths of 0.01×, 0.05×, and 0.1×, LD-PGTA achieved a mean BTPR across ancestry groups of 33.4%, 88.6%, and 97.5%, respectively, with a BFPR of 10%. Expanding our view across all classification thresholds, we found that the area under the balanced ROC curves was 0.72, 0.95, and 0.99 at depths of 0.01×, 0.05×, and 0.1×, respectively. Our results thus demonstrate that as the depth of coverage increases from 0.01× to 0.1×, the balanced ROC curve approaches a step function, nearing ideal classification performance.

We also observed that short read lengths performed better than longer reads at low coverage (SI Appendix, Fig. S2). While potentially counterintuitive, this relationship arises due to the fact that coverage is a linear function of the read length as well as the number of reads. Thus, holding coverage constant, a large number of shorter reads will achieve more uniform coverage than a small number of longer reads.

For meiotic-origin trisomies, meiotic cross-overs should manifest as switches between tracts of BPH and SPH trisomy. To test the utility of our method for revealing such recombination events, we simulated trisomic chromosomes with random mixtures of BPH and SPH tracts, while varying the sequencing depth, as well as the chromosome length and size of the corresponding genomic windows (Fig. 4). Applying LD-PGTA to the simulated data, we examined the resulting changes in the LLRs of consecutive bins and their relation to simulated meiotic cross-overs. While muted signatures could be discerned at the lowest coverage of 0.01×, the signatures were difficult to distinguish from background noise and spatial resolution was poor. In contrast, cross-overs were more pronounced at higher sequencing depths of 0.05× and 0.1× and closely corresponded to simulated cross-over breakpoints. For certain chromosomes (e.g., chromosome [Chr]11, Chr13, and Chr18) and at the higher sequencing depths, the LLRs observed in centromere-flanking bins provided an indication of the originating stage of the trisomy (here simulated as meiosis II errors). We note, however, that the bin encompassing the centromere itself was rarely informative, as such genomic regions are typically composed of repetitive heterochromatin that is inaccessible to standard short-read mapping and genotyping methods.

We next extended our simulations to test the sensitivity of our classifier to specification of a reference panel that is appropriately matched to the ancestry of the target sample. Specifically, we simulated samples by drawing haplotypes from a reference panel composed of one of the five superpopulations of the 1000 Genomes Project and then applied LD-PGTA while specifying a reference panel composed of haplotypes either from the same superpopulation (i.e., “matched”) or from a mixture of all superpopulations (i.e., “mismatched”). Notably, any phased genomic dataset of matched ancestry could be used in place of the 1000 Genomes Project dataset, with performance maximized in large datasets with closely matched ancestries.

While classification performance was high across all continental ancestry groups under the matched-ancestry scenario, performance moderately declined in all scenarios where the reference panel was misspecified (Fig. 3 and SI Appendix, Figs. S3–S5). Due to such misspecification, the mean BTPR across ancestry groups declined by 8.5%, 28.3%, and 26.4%, at a BFPR of 10% and depths of 0.01×, 0.05×, and 0.1×, respectively. Moreover, the area under the balanced ROC curves decreased by 0.06, 0.10, and 0.10 at depths of 0.01×, 0.05×, and 0.1×, respectively. This performance decline mimics that of related methods such as polygenic scores and arises by consequence of systematic differences in allele frequencies and LD structure between the respective populations (46–48). As an example of such effects, we observed that using an African reference panel for a European target sample produced a bias in favor of SPH trisomy (SI Appendix, Fig. S6). Conversely, using a European reference panel for an African target sample produced a bias in favor of BPH trisomy.

Seeking to test classification performance on individuals of admixed ancestries, we extended our simulation procedure to generate test samples composed of chromosomes drawn from distinct superpopulations. Following our previous analysis, we tested our method with and without correct specification of the component ancestries contributing to the admixture. Specifically, we implemented the latter scenario by specifying a single reference panel that matched the ancestry of either the maternal or the paternal haplotype (i.e., “partially matched”; SI Appendix, Fig. S7). Our results demonstrated that the performance of LD-PGTA for correctly specified admixture scenarios was comparable to that observed for nonadmixed scenarios. The ratio of mean area under the curve (AUC) across ancestries and coverages for nonadmixed samples versus admixed samples was 1.0. The impacts of reference panel misspecification were again greatest for admixed individuals with African ancestry, likely reflecting differences in structure of LD (shorter haplotypes) among African populations. For the African–European admixture scenario the BTPR decreased by 5.1%, 30.9%, and 28.7% at a BFPR of 10% and depths of 0.01×, 0.05×, and 0.1× when admixture was misspecified.

We next proceeded to apply our method to existing data from PGT-A, thereby further evaluating its performance and potential utility under realistic clinical scenarios. Data were obtained from IVF cases occurring between 2015 and 2020 at the Zouves Fertility Center. Embryos underwent trophectoderm biopsy at day 5, 6, or 7 postfertilization, followed by PGT-A using the Illumina VeriSeq PGS kit and protocol, which entails sequencing on the Illumina MiSeq platform (36-bp single-end reads). The dataset comprised a total of 8,154 embryo biopsies from 1,640 IVF cases, with maternal age ranging from 22 to 56 y (median = 38). Embryos were sequenced to a median coverage 0.0083× per sample, which we note lies near the lower limit of LD-PGTA and corresponds to an expected false positive rate of ∼0.02 at a 50% confidence threshold for classification of chromosome-scale patterns of SPH and BPH trisomy (SI Appendix, Fig. S8).

Given the aforementioned importance of the ancestry of the reference panel, we used LASER (49, 50) to perform automated ancestry inference of each embryo sample from the low-coverage sequencing data. LASER applies principal components analysis (PCA) to genotypes of reference individuals with known population affiliations. It then projects target samples onto the reference PCA space, using a Procrustes analysis to overcome the sparse nature of the data. Ancestry of each target sample is then deduced using a k-nearest neighbors approach.

Our analysis revealed a diverse patient cohort, consistent with the demographic composition of the local population (Fig. 5). Specifically, we inferred a total of 2,037 (25.0%) embryos of predominantly East Asian ancestries, 900 (11.0%) of South Asian ancestries, 2,554 (31.3%) of European ancestries, 669 (8.2%) of admixed American ancestries, and 44 (0.5%) of African ancestries, according to the reference panels defined by the 1000 Genomes Project. Interestingly, we also observed 1,936 (23.7%) embryos with principal component scores indicating admixture between parents of ancestries from one of the aforementioned superpopulations, which we thus evaluated using the generalization of PGT-A for admixed samples. More specifically, we inferred a total of 1,264 (15.5%) embryos to possess admixture of East Asian and European ancestries, 447 (5.5%) of South Asian and European ancestries, 129 (1.6%) of East and South Asian ancestries, and 96 (1.2%) of African and European ancestries. The ancestry of 14 (0.2%) embryos remained undetermined. To further test the robustness of these ancestry inferences with low amounts of input data, we separately analyzed chromosome 1 and compared the result to inferences for the rest of the genome. Even when restricting to this small proportion of the genome (¡10%), we observed concordances of 87.1% and 81.9% for putative nonadmixed and admixed individuals, respectively.

Ploidy status of each chromosome from each embryo biopsy was first inferred using the Illumina BlueFuse Multi Software suite in accordance with the VeriSeq protocol. Similar in principle to several open-source tools (51), BlueFuse Multi detects aneuploidies based on the coverage of reads aligned within 2,500 bins that are distributed along the genome, normalizing and adjusting for local variability in GC content and other potential biases. It then reports the copy number estimates for each bin (and interpolates between bins), as well as the copy number classification of each full chromosome (as a “gain,” a “loss,” or normal copy number) and numerous auxiliary metrics to assist with quality control.

Of the original 8,154 embryos, we excluded 138 embryos with low signal-to-noise ratios (derivative log ratio [DLR] $\ge$ 0.4), as well as an additional 970 embryos with low coverages that were unsuitable for analysis with LD-PGTA (mean coverage ¡0.005×). This resulted in a total of 7,046 embryos used for all downstream analyses. We focused our analysis on aneuploidies affecting entire chromosomes (i.e., nonsegmental aneuploidies; Materials and Methods), as these are the hypotheses that LD-PGTA is designed to test at an extremely low coverage (Fig. 1). Aggregating the BlueFuse Multi results from these embryos, we identified 2,006 embryos that possessed one or more nonmosaic whole-chromosome aneuploidies (28.5%), including 869 (12.3%) with one or more chromosome gains and 1,186 (16.8%) with one or more chromosome losses. Of these, a total of 494 embryos (7.0%) possessed aneuploidies involving two or more chromosomes. A total of 307 (4.4%) embryos possessed one or more chromosomes in the mosaic range, including 94 (1.3%) embryos that also harbored nonmosaic aneuploidies. The rate of nonmosaic aneuploidy was significantly associated with maternal age (quasibinomial generalized linear model [GLM]: ${\widehat{\beta }}_{\text{age}}=$ 0.0957, SE = 0.0080, P ¡ 1× 10${}^{-10}$), while the rate of mosaic aneuploidy was not (quasibinomial GLM: ${\widehat{\beta }}_{\text{age}}=$ 0.0278, SE = 0.0143, P = 0.0517), replicating well-described patterns from the literature (1). On an individual chromosome basis, the above results correspond to a total of 1,657 whole-chromosome losses (1,456 complete and 201 in the mosaic range) and 1,331 whole-chromosome gains (1,063 complete and 268 in the mosaic range) (Fig. 6A). These conventional coverage-based results served as the starting point for refinement and subclassification with LD-PGTA.

LD-PGTA is intended to complement coverage-based methods, using orthogonal evidence from genotype data to support, refute, or refine initial ploidy classifications. The relevant hypotheses for LD-PGTA to test are therefore based on the initial coverage-based classification. Applying LD-PGTA to the 1,456 putative nonmosaic whole-chromosome losses, we confirmed the monosomy diagnosis for 1,131 (77.7%) chromosomes, designated 290 as ambiguous (i.e., 50% confidence interval spanning 0), and obtained conflicting evidence favoring disomy for only 35 chromosomes (2.4%). We note that the latter group is within the expected margin of error (FPR = 5.1%) if all 1,456 chromosomes were in fact monosomic, given the selected confidence threshold of 50%.

Seeking to better understand the molecular origins of trisomy, we next applied LD-PGTA to the 1,331 originally diagnosed whole-chromosome gains (1,063 complete and 268 in the mosaic range). Of the nonmosaic gains, a total of 420 (39.5%) chromosomes exhibited evidence of BPH trisomy (i.e., meiotic origin), 217 (20.4%) chromosomes exhibited evidence of full SPH trisomy (i.e., mitotic origin or meiotic origin lacking recombination), and 426 (40.1%) chromosomes exhibited evidence of a mixture of BPH and SPH tracts (i.e., meiotic origin with recombination) or uncertain results. In contrast, the vast majority of the 268 putative mosaic chromosome gains exhibited evidence of SPH trisomy (152 chromosomes [56.7%]), while only 28 (10.4%) exhibited evidence of BPH trisomy and 88 (32.8%) remained ambiguous. We observed that BPH trisomies exhibited a nominally stronger association with maternal age than SPH trisomies, consistent with previous literature demonstrating a maternal age effect on meiotic but not mitotic aneuploidy, although the difference was not statistically significant (quasibinomial GLM: ${\widehat{\beta }}_{\text{age}\times \text{BPH/SPH}}=$ –0.0128, SE = 0.0184, P = 0.486; Fig. 6B) (29). Intriguingly, chromosomes 16 (binomial test: P = 5.3× 10${}^{-6}$) and 22 (binomial test: P = 2.7× 10${}^{-7}$) were strongly enriched for BPH versus SPH trisomies, consistent with a known predisposition of these two chromosomes to missegregation during maternal meiosis, including based on their strong maternal age effect (Fig. 6C) (29, 52).

Scanning along consecutive bins across a trisomic chromosome, switches between tracts of BPH and SPH trisomy indicate the occurrence of recombination, offering a method for mapping meiotic cross-overs and studying their role in the genesis of aneuploidy. To detect these cross-overs in practice, we omitted all ambiguous bins and then approximated each cross-over breakpoint as the midpoint between the centers of nearest-neighbor bins, where one was unambiguously classified as BPH and the other as SPH. Although the low-coverage nature of this particular dataset offered very coarse resolution (∼5 Mbp; Discussion), we nevertheless identified 1,614 putative meiotic cross-overs and mapped their heterogeneous distribution across 2,149 trisomic chromosomes (conditioning on coverage >0.01× and based on a ±1 SD confidence threshold; SI Appendix, Fig. S9).

Because of their reliance on differences in normalized coverage of reads aligned to each chromosome, sequencing-based implementations of PGT-A are typically blind to haploidy, triploidy, and other potential genome-wide aberrations. While differences in coverage between the X and Y chromosome may offer a clue about certain forms of triploidy (53), this signature does not apply to haploidy or to forms of triploidy where the Y chromosome is absent. Moreover, even when present, the short length of the Y chromosome diminishes the ratio of signal to noise and limits its diagnostic utility. As such, haploid and triploid embryos (especially those lacking Y chromosomes) are routinely misclassified as diploid by coverage-based methods for PGT-A analysis (34).

To overcome these challenges, we extended LD-PGTA to detect abnormalities in genome-wide ploidy, effectively combining evidence for or against specific chromosome abnormalities across the entire genome. In the case of haploidy, this entailed aggregating the LLRs of the comparison between disomy and monosomy across all genomic bins, while in the case of triploidy we aggregated the LLRs of the comparison between disomy and BPH trisomy across all genomic bins. Because the chromosomes of triploid samples will possess tracts of both BPH and SPH trisomy, our power for detection is lower than for haploidy, and it thus requires a less stringent threshold (Materials and Methods and SI Appendix, Fig. S10).

We started our analysis with the 4,495 embryos that were initially classified as euploid based on conventional coverage-based tests (i.e., Bluefuse Multi). To take a more conservative approach, we applied additional filtering to consider only embryos with >0.06× and >0.03× with at least 1,000 and 250 informative genomic windows for the classification of triploidy and haploidy, respectively. In the case of triploidy, this confined the initial pool of embryos to 3,821, among which we identified 12 (0.31%) as likely triploid (Materials and Methods, Eq. 14, and Fig. 7). In the case of haploidy, our stringent filtering criteria restricted our analysis to an initial pool of 1,320 embryos, among which we identified 11 (0.83%) as likely haploid (Materials and Methods, Eq. 10, and Fig. 7). Of the 12 putative triploid embryos, 5 (42%) exhibited an XXX composition of sex chromosomes. In contrast, all of the haploid embryos possessed an X chromosome and no Y chromosome, consistent with previous studies of IVF embryos that suggested predominantly gynogenetic origins of haploidy, in the context of both intracytoplasmic sperm injection (ICSI) and conventional IVF (54, 55). Among the 12 putative triploid embryos, we observed 5 embryos with coverage-based signatures of XXY sex chromosome complements. Moreover, 6 of the 12 embryos that we classified as haploid were independently noted as possessing only a single pronucleus (i.e., monopronuclear [1PN]) at the zygote stage upon retrospective review of notes from the embryology laboratory. Importantly, however, microscopic assessment of pronuclei is known to be prone to error and is imperfectly associated with ploidy status and developmental potential (56). Nevertheless, these results offer independent validation of our statistical approach, which may be used to flag diploid-testing embryos as harboring potential abnormalities in genome-wide ploidy.

By virtue of LD, observations of a set of alleles from one read may provide information about the probabilities of allelic states in another read that originated from the same DNA molecule (i.e., chromosome). In contrast, when comparing reads originating from distinct homologous chromosomes, allelic states observed in one read will be uninformative of allelic states observed in the other read. As different ploidy configurations are defined by different counts of identical and distinct homologous chromosomes, sparse genotype observations may provide information about the underlying ploidy status, especially when aggregated over large chromosomal regions. For a set of reads aligned to a defined genomic region, we compare the likelihoods of the observed alleles under four competing hypotheses:

Our statistical models are based on the premise that the odds of two reads being drawn from the same haplotype differ under the different scenarios. Specifically, for disomy, the odds are $1:1$; for monosomy, the odds are $1:0$; for BPH trisomy, the odds are $1:2$; and for SPH trisomy, the odds are $5:4$ (Fig. 2). If a pair of reads is drawn from identical homologs, the probability of observing the two alleles is given by the joint frequency of these two alleles (i.e., the frequency of the haplotype that they define) in the reference panel. In contrast, if a pair of reads is drawn from distinct homologs, then the probability of observing the two alleles is simply the product of the frequencies of the two alleles in the reference panel:where f(A) and f(B) are the frequencies of alleles A and B in the population. Likelihoods of two of the above hypotheses are compared by computing a log-likelihood ratio:

When a read overlaps with multiple SNPs, f(A) should be interpreted as the joint frequency of all SNP alleles that occur in read A (i.e., the frequency of haplotype A). Similarly, f(AB) would denote the joint frequency of all SNP alleles occurring in reads A and B. The equations above were extended to consider up 18 reads per window, as described in SI Appendix, Generalization to arbitrary ploidy hypotheses. Estimates of allele and haplotype frequencies from a reference panel do not depend on theoretical assumptions, but rely on the idea that the sample is randomly drawn from a population with similar ancestry. One limitation, which we consider, is that reliable estimates of probabilities near zero or one require large reference panels, such as the 1000 Genomes Project (57).

To quantify uncertainty in our likelihood estimates, we performed m over n bootstrapping by iteratively resampling reads within each window (44). Resampling was performed without replacement to comply with the assumptions of the statistical models about the odds of drawing two reads from the same haplotype. Thus, in each iteration, only subsets of the available reads can be resampled. Specifically, within each genomic window, up to 18 reads with a priority score exceeding a defined threshold are randomly sampled with equal probabilities. The likelihood of the observed combination of SNP alleles under each competing hypothesis is then calculated, and the hypotheses are compared by computing a LLR. The sample mean and the unbiased sample variance (i.e., with Bessel’s correction) of the LLR in each window are calculated by repeating this process using a bootstrapping approach,where γ_s is the log-likelihood ratio for the sth subsample of reads from the wth genomic window and m is the number of subsamples. Because the number of terms in the statistical models grows exponentially, we subsample at most 18 reads per window. Moreover, accurate estimates of joint frequencies of many alleles require a very large reference panel. Given the rate of heterozygosity in human populations and the size of the 1000 Genomes dataset, 18 reads are generally sufficient to capture one or more heterozygous SNPs that would inform our comparison of hypotheses.

Even when sequences are generated according to one of the hypotheses, a fraction of genomic windows will emit alleles that do not support that hypothesis and may even provide modest support for an alternative hypothesis. This phenomenon is largely driven by the sparsity of the data as well as the low rates of heterozygosity in human genomes, which together contribute to random noise. Another possible source of error is a local mismatch between the ancestry of the reference panel and the tested sequence. Moreover, technical errors such as spurious alignment and genotyping could also contribute to poor results within certain genomic regions (e.g., near the centromeres). To overcome this noise, we binned LLRs across consecutive genomic windows, thereby reducing biases and increasing the classification accuracy at the cost of a lower resolution. Specifically, we aggregated the mean LLRs of genomic windows within a larger bin,where ${\overline{\gamma }}^{(w)}$ is the mean of the LLRs associated with the wth genomic window. In addition, using the Bienaymé formula, we calculated the variance of the aggregated LLRs,where $\text{Var}({\gamma }^{(w)})$ is the variance of the LLRs associated with the wth window. For a sufficiently large bin, the confidence interval for the aggregated LLR is ${\mathrm{\Gamma }}_{\text{bin}}\pm z\sqrt{\text{Var}({\mathrm{\Gamma }}_{\text{bin}})}$, where $z={\mathrm{\Phi }}^{-1}\left(1-\alpha \right)$ is the z score, $\mathrm{\Phi }$ is the cumulative distribution function of the standard normal distribution, and $C=100(1-2\alpha )\%$ is the confidence level. The confidence level is chosen based on the desired sensitivity vs. specificity. We normalized the aggregated LLRs by the number of genomic windows that compose each bin, ${\overline{\gamma }}_{\text{bin}}={\mathrm{\Gamma }}_{\text{bin}}/g$. Thus, the variance of the mean LLR per window is $\text{Var}({\overline{\gamma }}_{\text{bin}})=\text{Var}({\mathrm{\Gamma }}_{\text{bin}})/g^2$. These normalized quantities can be compared across different regions of the genome, as long as the size of the genomic window is the same on average.

We developed a classification scheme to infer the ploidy status of each genomic bin. Each class is associated with a hypothesis about the number of distinct homologs and their degeneracy (i.e., copy number of nonunique homologs). To this end, we compare pairs of hypotheses by computing LLRs of competing statistical models. In general, the specific models that we compare are informed by orthogonal evidence obtained using standard coverage-based approaches to aneuploidy detection (i.e., tag counting).

The confidence interval for the mean LLR is ${\overline{\gamma }}_{\text{bin}}\pm z\sqrt{\text{Var}({\overline{\gamma }}_{\text{bin}})}$, and z is referred to as the z score. Thus, we classify a bin as exhibiting support for hypothesis 1 when

and for hypothesis 0 whenwhere the first (second) criterion is equivalent to requiring that the bounds of the confidence interval lie on the positive (negative) side of the number line.

For a given depth of coverage and read length, we simulate an equal number of sequences generated according to hypotheses 0 and 1, as explained in the previous section. Based on these simulations, we generate two ROC curves for each bin. The first ROC curve is produced by defining true positives as simulations where sequences generated under hypothesis 0 are correctly classified. For the second ROC curve, true positives are defined as simulations where sequences generated under hypothesis 1 are correctly classified. These two ROC curves can be combined into a single curve, which we term the “balanced ROC curve.” The balanced true and false positive rates for a bin are defined aswhere an H_i instance with i = 1, 2 is a sequence that was generated under hypothesis i.

The balanced ROC curve is particularly suited for classification tasks with three possible classes: H₀, H₁, and ambiguous. The ambiguous class contains all instances that do not fulfill the criteria in Eqs. 10 and 11, i.e., where the boundaries of the confidence interval span zero. This classification scheme allows us to optimize the classification of both H₀ and H₁ instances, at the expense of leaving out ambiguous instances. The advantage of this optimization is a reduction in the rate of spurious classification. To generate each curve, we varied the confidence level and the number of bins.

A key aim of our method is the reclassification of sequences that were initially identified as euploid by standard coverage-based approaches as exhibiting potential abnormalities of genome-wide ploidy. Specifically, coverage-based methods are based on the correlation between the depth of coverage of reads aligned to a genomic region and the copy number of that region. In the cases of (near) haploidy and triploidy, no such correlation is evident because (nearly) all chromosomes exhibit the same abnormality, resulting in erroneous classification.

To this end, following the binning procedure that was introduced in Eqs. 8 and 9, we aggregate the LLRs of each genomic window along the entire genome. To identify haploids, we test for cases where the confidence interval lies on the positive side of the number line, as formulated in Eq. 10.

In the idealized case where triploids are composed of entirely BPH sequence, triploids are easily distinguished from diploids. However, true cases of triploidy generally originate from retention of the polar body, leading to a more realistic scenario in which regions of BPH are relegated to the ends of chromosomes, while the rest of the genome exhibits SPH. The remedy is to introduce a binary classifier that assesses whether there is enough evidence of genome-wide BPH to question the hypothesis of genome-wide diploidy,

In other words, it is sufficient to show that the confidence interval spans the zero to classify the case as a triploid. We evaluated this criterion by applying our method to simulations of triploid embryos with realistic tracts of BPH and SPH trisomy. ROC curves constructed from these simulations demonstrate near-ideal performance at coverage $\ge$0.05× with a correctly specified reference panel.

The Homewood Institutional Review Board of Johns Hopkins University determined that this work does not qualify as federally regulated human subjects research (HIRB00011431).