Haplotype-aware inference of human chromosome abnormalities

Edited by Aravinda Chakravarti, New York University Langone Medical Center, New York, NY, and approved September 30, 2021 (received for review May 25, 2021)
November 12, 2021
118 (46) e2109307118

Significance

Whole-chromosome gains and losses (aneuploidies) are the leading cause of human pregnancy loss and congenital disorders. Recent work has demonstrated that in addition to harmful meiotic aneuploidies, mitotic aneuploidies (which lead to mosaic embryos harboring cells with different numbers of chromosomes) may also be common in preimplantation embryos but potentially compatible with healthy birth. Here we developed and tested a method for distinguishing these forms of aneuploidy using genetic testing data from 8,154 in vitro fertilization (IVF) embryos. We reclassified embryos based on signatures of meiotic and mitotic error, while also revealing lethal forms of chromosome abnormality that were previously hidden. Our method complements standard protocols for preimplantation genetic testing, while offering insight into the biology of early development.

Abstract

Extra or missing chromosomes—a phenomenon termed aneuploidy—frequently arise during human meiosis and embryonic mitosis and are the leading cause of pregnancy loss, including in the context of in vitro fertilization (IVF). While meiotic aneuploidies affect all cells and are deleterious, mitotic errors generate mosaicism, which may be compatible with healthy live birth. Large-scale abnormalities such as triploidy and haploidy also contribute to adverse pregnancy outcomes, but remain hidden from standard sequencing-based approaches to preimplantation genetic testing for aneuploidy (PGT-A). The ability to reliably distinguish meiotic and mitotic aneuploidies, as well as abnormalities in genome-wide ploidy, may thus prove valuable for enhancing IVF outcomes. Here, we describe a statistical method for distinguishing these forms of aneuploidy based on analysis of low-coverage whole-genome sequencing data, which is the current standard in the field. Our approach overcomes the sparse nature of the data by leveraging allele frequencies and linkage disequilibrium (LD) measured in a population reference panel. The method, which we term LD-informed PGT-A (LD-PGTA), retains high accuracy down to coverage as low as 0.05× and at higher coverage can also distinguish between meiosis I and meiosis II errors based on signatures spanning the centromeres. LD-PGTA provides fundamental insight into the origins of human chromosome abnormalities, as well as a practical tool with the potential to improve genetic testing during IVF.
Whole-chromosome gains and losses (aneuploidies) are extremely common in human embryos and are the leading cause of pregnancy loss and congenital disorders, both in the context of in vitro fertilization (IVF) and in natural conception (1, 2). Aneuploidy frequently arises during maternal meiosis due to mechanisms such as classical nondisjunction, premature separation of sister chromatids, and reverse segregation (3). Such meiotic aneuploidies are strongly associated with maternal age, with risk of aneuploid conception increasing exponentially starting around age 35 y. Although not as well characterized, research has also demonstrated that aneuploidy of mitotic origin is prevalent during the initial postzygotic cell divisions, potentially owing to relaxation of cell cycle checkpoints prior to embryonic genome activation (4, 5). Such mitotic errors, which are independent of maternal or paternal age, generate mosaic embryos possessing both normal and aneuploid cells (6). Mechanisms of mitotic aneuploidy include anaphase lag and mitotic nondisjunction (7), but also newly appreciated phenomena such as multipolar mitotic division (8, 9). Such abnormal mitoses are surprisingly common in cleavage-stage embryos (1012) and partially explain the high observed rates of embryonic mortality (∼50%) during preimplantation human development.
In light of these observations, preimplantation genetic testing for aneuploidy (PGT-A) has been developed as an approach to improve IVF outcomes by prioritizing chromosomally normal (i.e., euploid) embryos for transfer, based on the inferred genetic constitution of an embryo biopsy. First introduced in the early 1990s, PGT-A has been the subject of long-standing controversy, with some meta-analyses and clinical trials drawing its benefits into question (13, 14). Meanwhile, technical platforms underlying the test have steadily improved over time, with the current state of the art comprising low-coverage whole-genome sequencing of DNA extracted from 5 to 10 trophectoderm cells of blastocyst-stage embryos from days 5 to 7 postfertilization (15, 16).
The improved sensitivity and resolution of sequencing-based PGT-A have placed a renewed focus on chromosomal mosaicism as a potential confounding factor for diagnosis and interpretation (17). Mosaicism within a PGT-A biopsy may generate a copy number profile that is intermediate between euploid and aneuploid expectations, although it is important to note that technical artifacts can mimic such signatures (18). At the whole-embryo scale, mosaic aneuploidies may be prevalent but systematically underestimated due to the reliance of PGT-A on biopsies of one or few cells, which by chance may fail to sample low-level aneuploid lineages (19). While uniform aneuploidies arising from meiotic errors are unambiguously harmful, mosaic aneuploidies are potentially compatible with healthy live birth (2022). Research using mouse models and human gastruloids has recently offered initial insights into the selective elimination of aneuploid cells from mosaic embryos via lineage-specific mechanisms of apoptosis (2325).
Another underappreciated challenge in the analysis of sequencing-based PGT-A data is the detection of complex abnormalities including errors in genome-wide ploidy (e.g., triploidy and haploidy). Because existing algorithms detect chromosome abnormalities by comparing the normalized counts of aligned sequencing reads (i.e., depth of coverage) across chromosomes within a sample, inferences are compromised when many or all chromosomes are affected. Extreme cases such as haploidy and triploidy may evade detection entirely and be falsely interpreted as normal euploid samples. Such classification errors are particularly concerning given the association of ploidy abnormalities with molar pregnancy and miscarriage (26, 27). Triploidy in particular comprises more than 10% of cytogenetically abnormal miscarriages (28).
The ability to reliably distinguish meiotic- and mitotic-origin aneuploidies, as well as complex and genome-wide errors of ploidy, may thus prove valuable for enhancing IVF outcomes. Notably, trisomy (and triploidy) of meiotic origin is expected to exhibit a unique genetic signature, characterized by the presence of three distinct parental haplotypes (i.e., “both parental homologs” [BPH] from a single parent) in contrast to the mitotic trisomy signature of only two distinct haplotypes chromosome-wide (i.e., a “single parental homolog” [SPH] from each parent; Fig. 1) (29). Conversely, monosomy and haploidy will exhibit genetic signatures of only a single haplotype chromosome- or genome-wide, respectively. To date, few methods have explicitly attempted to use these signatures to distinguish forms of aneuploidy. Exceptions include approaches that require genetic material from both parents and embryo biopsies (3033) as well as targeted sequencing approaches that require alternative methods of library preparation to sequence short amplicons to higher coverage (34).
Fig. 1.
Signatures of various forms of chromosome abnormality with respect to their composition of identical and distinct parental homologs. Normal gametogenesis produces two genetically distinct copies of each chromosome—one copy from each parent—that comprise mosaics of two homologs possessed by each parent. Meiotic-origin trisomies may be diagnosed by the presence of one or more tracts with three distinct parental homologs (i.e., transmission of BPH from a given parent). In contrast, mitotic-origin trisomies are expected to exhibit only two genetically distinct parental homologs chromosome-wide (i.e., duplication of a SPH from a given parent). Triploidy and haploidy will mirror patterns observed for individual meiotic trisomies and monosomies, respectively, but across all 23 chromosome pairs—a pattern that confounds standard coverage-based analysis of PGT-A data.
Here we describe a statistical approach to classify aneuploidies using genotype information encoded in low-coverage whole-genome sequencing data from PGT-A. Inspired by the related challenge of imputation (3540), our method overcomes the sparse nature of the data by leveraging information from a population reference panel. We test our method on simulated and empirical data of varying sequencing depths, meiotic recombination patterns, and patient ancestries, evaluating its strengths and limitations under realistic scenarios. At higher coverage, we further demonstrate that our method permits the mapping of meiotic cross-overs on trisomic chromosomes, as well as the distinction of trisomies originating in meiosis I and meiosis II. Our method reveals biological insight into the fidelity of meiosis and mitosis, while also holding promise for improving preimplantation genetic testing.

Results

Genotypic Signatures of Meiotic and Mitotic Aneuploidy

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 (3540) 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).

LD-PGTA: A Statistical Model to Classify Aneuploidies

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 (104 to 105 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.”
Fig. 2.
(A) Within defined genomic windows, select reads overlapping informative SNPs that tag common haplotype variation. (B) Obtain joint frequencies of corresponding SNPs from a phased panel of ancestry-matched reference haplotypes. (C) Randomly resample 2 to 18 reads and compute probabilities of observed alleles under two competing ploidy hypotheses. (D) Compare the hypotheses by computing a likelihood ratio and estimate the mean and variance by resampling random sets of reads using a bootstrapping approach. (E) Repeat steps AD for consecutive nonoverlapping genomic windows and aggregate the log-likelihood ratios over larger genomic intervals. (F) Use the mean and variance of the combined log-likelihood ratio to estimate a confidence interval and classify the ploidy state of the genomic interval.

Benchmarking LD-PGTA with Simulated Sequences

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.
Fig. 3.
Balanced ROC curves for BPH vs. SPH with matched and random reference panels for nonadmixed embryos, varying depths of coverage. Each balanced ROC curve reflects an average over bins across the genome. We averaged both the BTPR and BFPR for common z scores across bins.
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.
Fig. 4.
Demonstration of the detection of meiotic cross-overs from low-coverage PGT-A data. Trisomies were simulated with varying locations of meiotic cross-overs, as depicted in Top diagrams, and varying depths of coverage (0.01×, 0.05×, and 0.1×). Confidence intervals correspond to a z score of one (confidence level of 68.3%). The size of the genomic windows (GW) varies with the coverage, while the size of the bins is kept constant.

Evaluating Sensitivity to Ancestries of the Target Sample and Reference Panel

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 (4648). 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.

Application to PGT-A Data

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).

Ancestry Inference for Reference Panel Selection

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.
Fig. 5.
Ancestry inference from low-coverage PGT-A data informs the selection of matched reference panels. Principal component axes were defined based on analysis of 1000 Genomes reference samples (A and C) and colored according to superpopulation annotations (African [AFR], Admixed American [AMR], East Asian [EAS], European [EUR], South Asian [SAS]). Low-coverage embryo samples were then projected onto these axes using a Procrustes approach implemented with LASER (v2.0) (50) (B and D) and their ancestries were classified using k-nearest neighbors of the first four principal component axes (k = 10). Projections of first-generation admixed reference samples were approximated as the mean of random samples from each of the component superpopulations. A and B depict principal components 1 and 2, while C and D depict principal components 2 and 3, together capturing continental-scale patterns of global ancestry.

Preliminary Analysis with a Coverage-Based Classifier

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] 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]: β^age= 0.0957, SE = 0.0080, P ¡ 1× 1010), while the rate of mosaic aneuploidy was not (quasibinomial GLM: β^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.
Fig. 6.
Application of LD-PGTA to low-coverage sequencing-based data. (A) LD-PGTA was used to refine classification results from 2,988 chromosomes originally diagnosed as whole-chromosome (i.e., nonsegmental) aneuploid based on standard coverage-based analysis with Bluefuse Multi, including strong validation of putative monosomies, as well as subclassification of BPH and SPH trisomies. (B) Association of BPH and SPH trisomies with maternal age. (C) LD-PGTA classifications of BPH versus SPH trisomy, stratifying by individual chromosome. BPH trisomy is strongly enriched on chromosomes 16 and 22, while signatures of SPH trisomy are more evenly distributed among the various autosomes. Chromosomes were classified using a 50% confidence interval. Dashed line indicates x = y.

Subclassification of Meiotic and Mitotic Trisomies

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: β^age×BPH/SPH= –0.0128, SE = 0.0184, P = 0.486; Fig. 6B) (29). Intriguingly, chromosomes 16 (binomial test: P = 5.3× 106) and 22 (binomial test: P = 2.7× 107) 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).

Hidden Abnormalities in Genome-Wide Ploidy

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.
Fig. 7.
Visualization of results from representative putative diploid (A), triploid (B), and haploid (C) samples. Copy number estimates obtained using a standard coverage-based approach (BlueFuse Multi) are depicted in A–C, Left and are indicative of diploidy. LD-PGTA results are depicted in A–C, Center and Right and suggest genome-wide abnormalities in ploidy for the latter two samples.

Discussion

Aneuploidy affects more than half of human embryos and is the leading cause of pregnancy loss and IVF failure. PGT-A seeks to improve IVF outcomes by prioritizing euploid embryos for transfer. Current low-coverage sequencing-based implementations of PGT-A rely entirely on comparisons of the depth of coverage of reads aligned to each chromosome. In such low-coverage settings, genotype observations are so sparse that they are typically regarded as uninformative.
Here we showed that by leveraging patterns of LD, even sparse genotype data are sufficient to capture signatures of aneuploidy, especially when aggregating over large genomic intervals such as entire chromosomes. In doing so, our method, termed LD-PGTA, reveals additional information that is hidden from standard coverage-based analyses of PGT-A data, including the designation of meiotic and mitotic trisomies, as well as the discovery of errors in genome-wide ploidy (i.e., haploidy and triploidy). Using simulations, we demonstrated the high accuracy of LD-PGTA even at extremely low depths of coverage (0.05×). We then showcased the utility of our method through application to an existing PGT-A dataset composed of 7,046 IVF embryos, stratifying trisomies of putative meiotic and mitotic origin and reclassifying 11 presumed diploid embryos as haploid and an additional 12 as triploid.
The key innovation of LD-PGTA is the use of allele frequencies and patterns of LD as measured in an external reference panel such as the 1000 Genomes Project (57). Because these parameters vary across human populations, the accuracy of our method depends on the correct specification of the reference panel with respect to the ancestry of the target sample. This challenge is analogous to that described in several recent studies investigating the portability of polygenic scores between populations, which are similarly biased by population differences in allele frequencies and patterns of LD (4648). In contrast to polygenic scores, however, our simulations suggest that the practical effects of reference panel misspecification on aneuploidy classification are typically modest. Moreover, our analysis of PGT-A data demonstrates that even at very low depths of coverage, existing extensions of principal components analysis (50) are sufficient to infer sample ancestries and guide the selection of an appropriate matched reference panel. However, this still depends on the availability of those genetic reference panels. While studies such as the 1000 Genomes Project (57) offer a broad sampling of global populations, certain populations remain underrepresented (including many African populations, as well as indigenous populations around the world). We anticipate that the performance of our method will improve as genomic databases grow larger and more diverse.
While LD-PGTA is conceptually related to genotype imputation in the sense that both exploit patterns of LD from a population reference panel, the goals of these two tasks are distinct. Whereas genotype imputation assumes some ploidy state and seeks to infer missing genotypes (e.g., for single-nucleotide polymorphisms [SNPs] that are not assayed by a genotyping microarray) (36, 3840), LD-PGTA uses the haplotype patterns to quantify evidence of different ploidy configurations. By combining evidence over large genomic intervals, LD-PGTA achieves classification power at coverages so low (<0.05×) that accurate genotype imputation is not possible (39). One intriguing area for future development is the use of pedigree information when genotype data from parents and/or siblings are also available, which could bolster the accuracy of LD-PGTA at very low coverages. We also note that while not yet tested for these applications, the LD-PGTA framework should generalize to other low-coverage sequencing scenarios, such as prenatal, postnatal, and miscarriage testing.
The genomic resolution of LD-PGTA is a function of multiple factors, including technical variables such as depth and evenness of sequencing coverage, read lengths, and desired confidence level, as well as biological variables such as the magnitudes of local LD and genetic variation. The low genetic diversity of human populations (π 0.001) is particularly limiting, even with increased sequencing coverage. In practice, an average coverage of 0.01× offers resolution of ∼5 Mbp, which corresponds to 10 nonoverlapping genomic bins on the shortest human chromosome. At higher depths of coverage (∼0.5×), our method permits the mapping of points of transition between different ploidy configurations. Depending on the particular hypotheses under consideration, such transitions could reflect evidence of segmental aneuploidies or—in the case of BPH and SPH trisomy—allow the mapping of meiotic cross-over breakpoints on trisomic chromosomes. Beyond the clinical applications previously discussed, this output of our method will facilitate future research into the factors influencing meiotic recombination, as well as its impacts on aneuploidy formation—topics of longstanding interest in basic reproductive biology (58).
Standard sequencing-based implementations of PGT-A infer the copy number of each chromosome based on variation in the local depth of coverage and typically support coverages as low as 0.005×. In contrast, LD-PGTA requires mean coverage of approximately 0.05× to yield high accuracy for any individual embryo, as demonstrated by simulations across a wide range of coverage, read length, and ancestry and admixture scenarios. Nevertheless, when applied to large datasets such as that investigated in our study, global patterns emerge at even lower coverage that offer rough stratification of meiotic and mitotic errors and provide insight into the biological origins of aneuploidy, beyond the binary classification of aneuploid and euploid. As costs of sequencing continue to plummet, application of the method to higher-coverage datasets and additional stages of development will further unlock its potential for both research and diagnostic aims.
We envision our method complementing rather than supplanting current coverage-based approaches to PGT-A, whose performance remains superior for tasks such as classification of monosomies and trisomies of few chromosomes. Meanwhile, the application of LD-PGTA can be used to flag haploidies or triploidies that remain invisible to current coverage-based approaches. Additionally, the subclassification of meiotic and mitotic trisomies may prove valuable as orthogonal evidence to distinguish potentially viable mosaic embryos from those possessing lethal meiotic errors. This knowledge is particularly relevant to the many patients with no euploid-testing embryos available for transfer (59). Together, our method offers an approach for extracting useful hidden information from standard preimplantation and prenatal genetic testing data, toward the goal of improving pregnancy outcomes.

Materials and Methods

Comparing Hypotheses with the Likelihood Ratio

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:
1)
A single copy of a chromosome, namely monosomy, which may arise by meiotic mechanisms such as nondisjunction, premature separation of sister chromatids, and reverse segregation or mitotic mechanisms such as mitotic nondisjunction and anaphase lag.
2)
Two distinct homologous chromosomes, namely disomy, the outcome of normal meiosis, fertilization, and embryonic mitosis.
3)
Two identical homologs with a third distinct homolog, denoted as SPH. SPH may originate from mitotic error (Fig. 1) or rare meiotic errors without recombination (29).
4)
Three distinct homologous chromosomes, denoted as BPH. BPH observed in any portion of a chromosome is a clear indication of meiotic error (Fig. 1) (29).
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:
Pmonosomy(AB)=f(AB),
[1]
Pdisomy(AB)=12f(AB)+12f(A)f(B),
[2]
PSPH(AB)=59f(AB)+49f(A)f(B),
[3]
PBPH(AB)=13f(AB)+23f(A)f(B),
[4]
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:
γ(A,B)=logPhypothesis 1(AB)Phypothesis 0(AB).
[5]
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).

Quantifying Uncertainty by Bootstrapping

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,
γ¯(w)=1mswγs,
[6]
Var(γ(w))=1m1sw(γsγ¯(w))2.
[7]
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.

Aggregating Log-Likelihood Ratios across Consecutive Windows

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,
Γbin=wbinγ¯(w),
[8]
where γ¯(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,
Var(Γbin)=wbinVar(γ(w)),
[9]
where Var(γ(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 Γbin±zVar(Γbin), where z=Φ1(1α) is the z score, Φ is the cumulative distribution function of the standard normal distribution, and C=100(12α)% 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, γ¯bin=Γbin/g. Thus, the variance of the mean LLR per window is Var(γ¯bin)=Var(Γbin)/g2. 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.

Evaluating Model Performance on Simulated Data

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 γ¯bin±zVar(γ¯bin), and z is referred to as the z score. Thus, we classify a bin as exhibiting support for hypothesis 1 when
γ¯binzVar(γ¯bin)>0,
[10]
and for hypothesis 0 when
γ¯bin+zVar(γ¯bin)<0,
[11]
where 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 as
BTPR=12(classifiedas H0H0 instancesH0 instances+classifiedas H1H1 instancesH1 instances),
[12]
BFPR=12(classifiedas H1H0 instancesH0 instances+classifiedas H0H1 instancesH1 instances),
[13]
where an Hi 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: H0, H1, 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 H0 and H1 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.

Detecting Haploidy and Triploidy

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,
γ¯bin+zVar(γ¯bin)>0.
[14]
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 0.05× with a correctly specified reference panel.

Ethics Approval and Consent to Participate

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

Data Availability

Software for implementing our method and reproducing all analysis, Bluefuse Multi aneuploidy data, LD-PGTA results, and aneuploidy calls (i.e., PGT-A results for each chromosome of each sample) for 8,154 deidentified embryos have been deposited in GitHub, https://github.com/mccoy-lab/LD-PGTA/tree/master/data. The repository is also archived on Zenodo, https://doi.org/10.5281/zenodo.5610347.

Acknowledgments

We thank the Origins of Aneuploidy Research Consortium, Shai Carmi, and Alexander Zaranek, as well as members of the McCoy laboratory for constructive feedback. We also thank the staff of the Maryland Advanced Research Computing Center for computing support. Research reported in this publication was supported by the National Institute of General Medical Sciences of the National Institutes of Health under Award R35GM133747. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Supporting Information

Appendix 01 (PDF)

References

1
T. Hassold, P. Hunt, To err (meiotically) is human: The genesis of human aneuploidy. Nat. Rev. Genet. 2, 280–291 (2001).
2
S. I. Nagaoka, T. J. Hassold, P. A. Hunt, Human aneuploidy: Mechanisms and new insights into an age-old problem. Nat. Rev. Genet. 13, 493–504 (2012).
3
J. R. Gruhn et al., Chromosome errors in human eggs shape natural fertility over reproductive life span. Science 365, 1466–1469 (2019).
4
R. H. Harrison, H. C. Kuo, P. N. Scriven, A. H. Handyside, C. M. Ogilvie, Lack of cell cycle checkpoints in human cleavage stage embryos revealed by a clonal pattern of chromosomal mosaicism analysed by sequential multicolour FISH. Zygote 8, 217–224 (2000).
5
A. A. Kiessling et al., Evidence that human blastomere cleavage is under unique cell cycle control. J. Assist. Reprod. Genet. 26, 187–195 (2009).
6
R. C. McCoy, Mosaicism in preimplantation human embryos: When chromosomal abnormalities are the norm. Trends Genet. 33, 448–463 (2017).
7
E. Mantikou, K. M. Wong, S. Repping, S. Mastenbroek, Molecular origin of mitotic aneuploidies in preimplantation embryos. Biochim. Biophys. Acta 1822, 1921–1930 (2012).
8
C. S. Ottolini et al., Tripolar mitosis and partitioning of the genome arrests human preimplantation development in vitro. Sci. Rep. 7, 9744 (2017).
9
R. C. McCoy et al., Tripolar chromosome segregation drives the association between maternal genotype at variants spanning PLK4 and aneuploidy in human preimplantation embryos. Hum. Mol. Genet. 27, 2573–2585 (2018).
10
Q. Zhan, Z. Ye, R. Clarke, Z. Rosenwaks, N. Zaninovic, Direct unequal cleavages: Embryo developmental competence, genetic constitution and clinical outcome. PLoS One 11, e0166398 (2016).
11
T. Cavazza et al., Parental genome unification is highly error-prone in mammalian embryos. Cell 184, 2860–2877.e22 (2021).
12
E. Ford et al., The first mitotic division of the human embryo is highly error-prone. bioRxiv [Preprint] (2020). https://doi.org/10.1101/2020.07.17.208744 (Accessed 1 December 2020).
13
S. Mastenbroek et al., In vitro fertilization with preimplantation genetic screening. N. Engl. J. Med. 357, 9–17 (2007).
14
S. Munné et al., Preimplantation genetic testing for aneuploidy versus morphology as selection criteria for single frozen-thawed embryo transfer in good-prognosis patients: A multicenter randomized clinical trial. Fertil. Steril. 112, 1071–1079.e7 (2019).
15
J. R. Vermeesch, T. Voet, K. Devriendt, Prenatal and pre-implantation genetic diagnosis. Nat. Rev. Genet. 17, 643–656 (2016).
16
M. Viotti, Preimplantation genetic testing for chromosomal abnormalities: Aneuploidy, mosaicism, and structural rearrangements. Genes (Basel) 11, 602 (2020).
17
M. Popovic, L. Dhaenens, A. Boel, B. Menten, B. Heindryckx, Chromosomal mosaicism in human blastocysts: The ultimate diagnostic dilemma. Hum. Reprod. Update 26, 313–334 (2020).
18
A. Capalbo, F. M. Ubaldi, L. Rienzi, R. Scott, N. Treff, Detecting mosaicism in trophectoderm biopsies: Current challenges and future possibilities. Hum. Reprod. 32, 492–498 (2017).
19
M. R. Starostik, O. A. Sosina, R. C. McCoy, Single-cell analysis of human embryos reveals diverse patterns of aneuploidy and mosaicism. Genome Res. 30, 814–825 (2020).
20
E. Greco, M. G. Minasi, F. Fiorentino, Healthy babies after intrauterine transfer of mosaic aneuploid blastocysts. N. Engl. J. Med. 373, 2089–2090 (2015).
21
A. R. Victor et al., One hundred mosaic embryos transferred prospectively in a single clinic: Exploring when and why they result in healthy pregnancies. Fertil. Steril. 111, 280–293 (2019).
22
M. Viotti et al., Using outcome data from one thousand mosaic embryo transfers to formulate an embryo ranking system for clinical use. Fertil. Steril. 115, 1212–1224 (2021).
23
H. Bolton et al., Mouse model of chromosome mosaicism reveals lineage-specific depletion of aneuploid cells and normal developmental potential. Nat. Commun. 7, 11165 (2016).
24
S. Singla, L. K. Iwamoto-Stohl, M. Zhu, M. Zernicka-Goetz, Autophagy-mediated apoptosis eliminates aneuploid cells in a mouse model of chromosome mosaicism. Nat. Commun. 11, 2958 (2020).
25
M. Yang et al., Depletion of aneuploid cells in human embryos and gastruloids. Nat. Cell Biol. 23, 314–321 (2021).
26
S. D. Lawler, S. Povey, R. A. Fisher, V. J. Pickthall, Genetic studies on hydatidiform moles. II. The origin of complete moles. Ann. Hum. Genet. 46, 209–222 (1982).
27
P. A. Jacobs, A. E. Szulman, J. Funkhouser, J. S. Matsuura, C. C. Wilson, Human triploidy: Relationship between parental origin of the additional haploid complement and development of partial hydatidiform mole. Ann. Hum. Genet. 46, 223–231 (1982).
28
B. Levy et al., Genomic imbalance in products of conception: Single-nucleotide polymorphism chromosomal microarray analysis. Obstet. Gynecol. 124, 202–209 (2014).
29
R. C. McCoy et al., Evidence of selection against complex mitotic-origin aneuploidy during preimplantation development. PLoS Genet. 11, e1005601 (2015).
30
A. H. Handyside et al., Karyomapping: A universal method for genome wide analysis of genetic disease based on mapping crossovers between parental haplotypes. J. Med. Genet. 47, 651–658 (2010).
31
D. S. Johnson et al., Preclinical validation of a microarray method for full molecular karyotyping of blastomeres in a 24-h protocol. Hum. Reprod. 25, 1066–1075 (2010).
32
M. Rabinowitz et al., Origins and rates of aneuploidy in human blastomeres. Fertil. Steril. 97, 395–401 (2012).
33
M. Zamani Esteki et al., Concurrent whole-genome haplotyping and copy-number profiling of single cells. Am. J. Hum. Genet. 96, 894–912 (2015).
34
D. Marin et al., Validation of a targeted next generation sequencing-based comprehensive chromosome screening platform for detection of triploidy in human blastocysts. Reprod. Biomed. Online 36, 388–395 (2018).
35
J. Marchini, B. Howie, Genotype imputation for genome-wide association studies. Nat. Rev. Genet. 11, 499–511 (2010).
36
C. Fuchsberger, G. R. Abecasis, D. A. Hinds, minimac2: Faster genotype imputation. Bioinformatics 31, 782–784 (2015).
37
S. Das, G. R. Abecasis, B. L. Browning, Genotype imputation from large reference panels. Annu. Rev. Genomics Hum. Genet. 19, 73–96 (2018).
38
B. L. Browning, Y. Zhou, S. R. Browning, A one-penny imputed genome from next-generation reference panels. Am. J. Hum. Genet. 103, 338–348 (2018).
39
S. Rubinacci, D. M. Ribeiro, R. J. Hofmeister, O. Delaneau, Efficient phasing and imputation of low-coverage sequencing data using large reference panels. Nat. Genet. 53, 120–126 (2021).
40
R. W. Davies et al., Rapid genotype imputation from sequence with reference panels. Nat. Genet. 53, 1104–1111 (2021).
41
S. H. Vohr, C. F. Buen Abad Najar, B. Shapiro, R. E. Green, A method for positive forensic identification of samples from extremely low-coverage sequence data. BMC Genomics 16, 1034 (2015).
42
N. Li, M. Stephens, Modeling linkage disequilibrium and identifying recombination hotspots using single-nucleotide polymorphism data. Genetics 165, 2213–2233 (2003).
43
N. E. Lamb et al., Susceptible chiasmate configurations of chromosome 21 predispose to non-disjunction in both maternal meiosis I and meiosis II. Nat. Genet. 14, 400–405 (1996).
44
M. R. Chernick, Bootstrap Methods: A Guide for Practitioners and Researchers (John Wiley & Sons, Hoboken, NJ, ed. 2, 2007), vol. 619.
45
W. J. Krzanowski, D. J. Hand, ROC Curves for Continuous Data (CRC Press, 2009).
46
A. R. Martin et al., Human demographic history impacts genetic risk prediction across diverse populations. Am. J. Hum. Genet. 100, 635–649 (2017).
47
A. R. Martin et al., Clinical use of polygenic scores will risk exacerbating health disparities. Nat. Genet. 51, 584–591 (2019).
48
Y. Wang et al., Theoretical and empirical quantification of the accuracy of polygenic scores in ancestry divergent populations. Nat. Commun. 11, 3865 (2020).
49
C. Wang et al., Ancestry estimation and control of population stratification for sequence-based association studies. Nat. Genet. 46, 409–415 (2014).
50
C. Wang, X. Zhan, L. Liang, G. R. Abecasis, X. Lin, Improved ancestry estimation for both genotyping and sequencing data using projection Procrustes analysis and genotype imputation. Am. J. Hum. Genet. 96, 926–937 (2015).
51
A. Abyzov, A. E. Urban, M. Snyder, M. Gerstein, CNVnator: An approach to discover, genotype, and characterize typical and atypical CNVs from family and population genome sequencing. Genome Res. 21, 974–984 (2011).
52
J. M. Franasiak et al., Aneuploidy across individual chromosomes at the embryonic level in trophectoderm biopsies: Changes with patient age and chromosome structure. J. Assist. Reprod. Genet. 31, 1501–1509 (2014).
53
K. Mutia et al., The frequency of chromosomal euploidy among 3PN embryos. J. Reprod. Infertil. 20, 127–131 (2019).
54
C. Staessen, A. C. Van Steirteghem, The chromosomal constitution of embryos developing from abnormally fertilized oocytes after intracytoplasmic sperm injection and conventional in-vitro fertilization. Hum. Reprod. 12, 321–327 (1997).
55
B. Rosenbusch, The chromosomal constitution of embryos arising from monopronuclear oocytes in programmes of assisted reproduction. Int. J. Reprod. Med. 2014, 418198 (2014).
56
A. S. T. Lim, V. H. H. Goh, C. L. Su, S. L. Yu, Microscopic assessment of pronuclear embryos is not definitive. Hum. Genet. 107, 62–68 (2000).
57
The 1000 Genomes Project Consortium, A global reference for human genetic variation. Nature 526, 68–74 (2015).
58
T. J. Hassold, P. A. Hunt, Missed connections: Recombination and human aneuploidy. Prenat. Diagn. 41, 584–590 (2021).
59
Z. P. Demko, A. L. Simon, R. C. McCoy, D. A. Petrov, M. Rabinowitz, Effects of maternal age on euploidy rates in a large cohort of embryos analyzed with 24-chromosome single-nucleotide polymorphism-based preimplantation genetic screening. Fertil. Steril. 105, 1307–1313 (2016).
60
D. Ariad et al., LD-PGTA code and data. Zenodo. https://doi.org/10.5281/zenodo.5610347. Deposited 29 October 2021.

Information & Authors

Information

Published in

Go to Proceedings of the National Academy of Sciences
Go to Proceedings of the National Academy of Sciences
Proceedings of the National Academy of Sciences
Vol. 118 | No. 46
November 16, 2021
PubMed: 34772814

Classifications

Data Availability

Software for implementing our method and reproducing all analysis, Bluefuse Multi aneuploidy data, LD-PGTA results, and aneuploidy calls (i.e., PGT-A results for each chromosome of each sample) for 8,154 deidentified embryos have been deposited in GitHub, https://github.com/mccoy-lab/LD-PGTA/tree/master/data. The repository is also archived on Zenodo, https://doi.org/10.5281/zenodo.5610347.

Submission history

Accepted: September 16, 2021
Published online: November 12, 2021
Published in issue: November 16, 2021

Keywords

  1. trisomy
  2. triploidy
  3. haploidy
  4. mosaicism
  5. in vitro fertilization

Acknowledgments

We thank the Origins of Aneuploidy Research Consortium, Shai Carmi, and Alexander Zaranek, as well as members of the McCoy laboratory for constructive feedback. We also thank the staff of the Maryland Advanced Research Computing Center for computing support. Research reported in this publication was supported by the National Institute of General Medical Sciences of the National Institutes of Health under Award R35GM133747. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Notes

Competing interest statement: D.A., M.V., and R.C.M. are coinventors of the method described herein, which is the subject of a provisional patent application by Johns Hopkins University.
This article is a PNAS Direct Submission.

Authors

Affiliations

Department of Biology, Johns Hopkins University, Baltimore, MD 21218;
Department of Biology, Johns Hopkins University, Baltimore, MD 21218;
Andrea R. Victor
Zouves Fertility Center, Foster City, CA 94131;
Zouves Fertility Center, Foster City, CA 94131;
Christo G. Zouves
Zouves Fertility Center, Foster City, CA 94131;
Zouves Foundation for Reproductive Medicine, Foster City, CA 94131
Department of Biology, Johns Hopkins University, Baltimore, MD 21218;

Notes

1
To whom correspondence may be addressed. Email: [email protected] or [email protected].
Author contributions: D.A., M.V., and R.C.M. designed research; D.A. and R.C.M. performed research; A.R.V., F.L.B., and C.G.Z. contributed new reagents/analytic tools; D.A. and R.C.M. analyzed data; D.A., S.M.Y., M.V., and R.C.M. wrote the paper; and A.R.V., F.L.B., C.G.Z., and M.V. performed data collection and data curation.

Metrics & Citations

Metrics

Note: The article usage is presented with a three- to four-day delay and will update daily once available. Due to ths delay, usage data will not appear immediately following publication. Citation information is sourced from Crossref Cited-by service.


Citation statements




Altmetrics

Citations

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

Cited by

    Loading...

    View Options

    View options

    PDF format

    Download this article as a PDF file

    DOWNLOAD PDF

    Get Access

    Login options

    Check if you have access through your login credentials or your institution to get full access on this article.

    Personal login Institutional Login

    Recommend to a librarian

    Recommend PNAS to a Librarian

    Purchase options

    Purchase this article to get full access to it.

    Single Article Purchase

    Haplotype-aware inference of human chromosome abnormalities
    Proceedings of the National Academy of Sciences
    • Vol. 118
    • No. 46

    Media

    Figures

    Tables

    Other

    Share

    Share

    Share article link

    Share on social media