Predicting embryonic aneuploidy rate in IVF patients using whole-exome sequencing.

Patient DNA samples were obtained from Reproductive Medicine Associates of New Jersey (RMANJ) DNA Bank and was approved by the IRB #RMA1-09-165 at Copernicus Group IRB and IRB #Pro2018000106 at Rutgers University. Two published whole-exome sequencing sets were used in this study ( Tyc et al. 2020a , b ; Tyc et al. 2021 ). Details of the sequencing and variant calling were described in previous publications. Briefly, the dataset from ( Tyc et al. 2020a , b ), referred to as the “IonTorrent” in the following text, was generated on an Ion Proton instrument (Thermo Fisher Scientific, Waltham, MA, USA) using an Ion AmpliSeq Exome kit (Thermo Fisher Scientific, Waltham, MA, USA). The variant calling was performed using the Ion Torrent Suite software v4.4 (Thermo Fisher Scientific, Waltham, MA, USA) with default parameters ( Tyc et al. 2020a , b ). The dataset from ( Tyc et al. 2021 ), referred to as the “Illumina” in the following text, was sequenced on the Illumina sequencing platform (Illumina, San Diego, CA, USA) using the Agilent SureSelect Human All Exon V6 kit (Agilent Technologies, CA, USA). The variant calling was performed using the GATK v3.8 pipeline following the GATK best practices ( Tyc et al. 2021 ). An additional 62 unpublished samples that were sequenced and from which variants were identified using the same procedure as the IonTorrent data set were included in this study. These samples will be referred to as “Intermediate” in the following text. The “Intermediate” dataset is defined by the aneuploidy phenotype of the patients in the group (see Individual group assignment section below for more detail). AVA,Dx pipeline [prototype described in ( Wang et al. 2019 ); pipeline manuscript in preparation; proprietary pipeline code available upon request] was used to conduct individual sample and variant quality control (QC), variant score assignment, gene score calculation, feature selection, and model building and testing (described below). For the IonTorrent dataset, the genotype QC was originally performed as a part of the Ion Torrent variant calling pipeline. In the current study, variants were further filtered based on site-wise qualities (quality (QUAE) > 30, mean Read Depth (DP) > = 6, and mean DP = 20%) were also removed. The same QC procedures were applied to the Intermediate dataset for downstream analyses. For the Illumina dataset, variants that failed the GATK’s Variant Quality Score Recalibration (VQSR) standard X( i.e. , FILTER ≠ “PASS” in the VCF file) were removed. Each variant was then filtered based on site-wise quality (QUAL > 30, mean DP > = 6, and mean DP 0.3 and AB 15, and DP > 4). Genotypes that failed genotype quality filters were converted into missing (“./.”) and variant sites with a high missing rate (MR > = 20%) were removed. In all datasets, mitochondrial and sex chromosome variants were removed. After variant QC, an individual’s ethnicity was inferred with EthSEQ ( Romanel et al. 2017 ) and individual relatedness within a cohort was inferred with SNPRelate ( Zheng et al. 2012 ), as parts of the AVA,Dx pipeline. EthSEQ infers an individual’s ancestry based on reference samples from the 1000 Genomes project ( Lowy-Gallego et al. 2019 ). Individuals that are labeled as “inside” of the 1000 Genomes European ancestry cluster (EUR) were selected. The initial IonTorrent dataset consisted of 166 individuals ( Tyc et al. 2020a , b ) and 139 individuals were selected following EthSEQ. The initial Illumina dataset consisted of 160 individuals ( Tyc et al. 2021 ), of which 142 individuals were selected. Among the 62 patients in the initial Intermediate dataset, 50 were inferred as EUR and selected for analysis. For SNPRelate analysis, individuals that are expected to be first-degree relatives (kinship coefficient > = 0.25) were considered related. In each cohort, no individual pair had kinship coefficient > = 0.25 so no individual was excluded. Detailed information of individuals that passed QC (age, aneuploidy rate, etc.) is provided in Supplemental Table 1 . After individual QC, variant QC were performed again to select the final set of variants. After individual QC and filtering, individuals were grouped as either producing low or high proportions of aneuploid blastocysts and referred to as low rate group (LRG) and high rate group (HRG), respectively. For both IonTorrent and Illumina datasets, individuals with ≤ 30% aneuploid blastocysts were defined as LRG, and individuals with ≥ 50% aneuploid blastocysts as HRG. All individuals included in both datasets had at least 4 embryos tested and used for aneuploidy rate calculation ( Tyc et al. 2020a , b ; Tyc et al. 2021 ). The IonTorrent dataset included 82 LRG and 57 HRG individuals, while the Illumina dataset included 68 LRG and 74 HRG individuals. The intermediate dataset contains 50 individuals with intermediate aneuploidy rates falling between 30 and 50%. Variant scores were standardized to fit a 0 to 1 range as described previously ( Wang et al. 2019 ). A synonymous single nucleotide variant (sSNV) was assigned a score of 0.05, indicating that the variant is likely to have a weak functional effect. For non-synonymous single nucleotide variants (nsSNVs), SNAP (screening for non-acceptable polymorphisms) scores ( Bromberg and Rost 2007 ) were used. SNAP uses a neural network-based method to predict the functional impact of nsSNVs and SNAP scores can range between − 1 and 1 ( Bromberg and Rost 2007 ). Because sSNV was assigned a score of 0.05, SNAP scores for nsSNVs were transformed to range between 0.055 and 1 as follows: 0.055 for neutral variants (SNAP score ≤ 0) and 0.06 + SNAP score*0.94 for non-neutral ones (SNAP score > 0). For stop-loss/gain SNVs and insertion/deletions (INDELs), variant scores were assigned as stop-loss/gain SNVs = 1; frameshift INDELs = 1; and non-frameshift INDELs = 0.5. Although non-coding variants could be functionally important and have clinical significance, the functional prediction for non-coding variants is still difficult ( Novikova et al. 2021 ; French and Edwards 2020 ). Therefore, non-coding variants were not included in the current analysis. To score genes for model building, the AVA,Dx pipeline assigns the zygosity variable to 0.25 for heterozygous variants and to 1 for homozygous ones, as optimized in previous work ( Wang et al. 2019 ). AVA,Dx offers two options for scoring genes on the basis of the zygosity -adjusted variant scores. For one individual at a time, for a given gene, all variant scores are either (1) summed (referred to as “Sum”) to quantify how much the protein encoded by the gene is functionally affected based on all variants the gene has; or (2) inverted to indicate remaining functionality and calculates the product of all kept functions within the gene (i.e., for each variant the transformation 1—zygosity*variant score was applied, and all transformed scores for the gene were multiplied, referred to as “Product”). Both options in scoring were evaluated in our cohorts. With the calculated gene scores, leave-one-out cross-validation (LOOCV) method was used to assign different patients into training and testing sets for feature selection and model building. For each iteration of the model training, the Kolmogorov–Smirnov (K–S) test was used to rank genes based on the score differences between patients with high and low aneuploidy rates, and top genes with p value < 0.2 (200 maximum) were selected for model building ( Wang et al. 2019 ). While this is a very high threshold for significance, our aim was to maximize gene set size while removing most uninformative genes. Random Forest (RF) and Support Vector Machine (SVM) models were then built using different numbers of selected genes from K–S test as features (top 5 to 200 genes in increments of 5). For both RF and SVM, HRG and LRG class weights (wLRG, wHRG) were calculated based on the sample size of each group (nLRG, nHRG) by the equations: wLRG = (1/nLRG)/(1/nHRG + 1/nLRG); wHRG = (1/nHRG)/(1/nHRG + 1/nLRG). The calculated weights were assigned by applying the parameter values class_weight = {0: 0.41, 1: 0.59} in IonTorrent and {0: 0.52, 1: 0.48} in Illumina, respectively. For RF, the number of trees in the forest was set to 207 in IonTorrent and 211 in Illumina as determined by the parameter n_estimators = 1.5 * the number of patients. The whole dataset is used to build each tree (bootstrap=False), and the maximum depth of the tree is 1 (max_depth = 1). For SVM probability estimates were enabled by probability = True. Other settings remained default as in the AVA,Dx package documentation. For each model architecture (e.g., top-5-gene, top-10-gene, etc.), individual predictions obtained from each iteration of the model were used to evaluate the model performance. The model architecture with the best ROC-AUC (Area Under the Receiver Operator Characteristic curve) was selected and genes from each iteration were included in the final list of candidate genes for the gene analysis. The Intermediate dataset was used to test the IonTorrent model prediction performance as follows: (1) gene scores were calculated for each individual in the Intermediate dataset following the Sum procedure; (2) the best IonTorrent LOOCV model architecture (i.e., a top-10-gene model) was used to construct the prediction model using the LRG/HRG IonTorrent dataset and features (i.e., genes) that were present in > 50% of the LOOCV model iterations; (3) individuals in the Intermediate set were evaluated to generate the prediction score. Note that our final model is not expected to be generalizable to predict individual phenotypes of other cohorts as it is overfit to our LRG/HRG cohort/dataset. As such, the classification accuracy for the LRG/HRG data is informative of the maximum achievable performance using these genes to classify this data set, while the Intermediate data set predictions elucidate model performance on individuals from the same cohort that were not used in training. Gene expression profiles of every candidate gene were collected for early preimplantation embryonic stages, including zygote, 4 cell, 8 cell, compacted morula, early blastocyst (ICM), and late blastocyst (epiblast and primitive endoderm) ( Stirparo et al. 2018 ). For gene function annotation, missense Z score (mis_Z) and putative Loss-Of-Function (pLoF) Z scores (lof_Z) were extracted from gnomAD v2.1.1 for each gene ( Karczewski et al. 2020 ). Experimentally verified meiosis genes in model organisms and predicted meiosis genes in mice were extracted from MeiosisOnline ( Jiang et al. 2021 ) ( https://mcg.ustc.edu.cn/bsc/meiosis/index.html ). Gene functions and related diseases were manually curated. Three databases were used to investigate protein–protein interaction (PPI) networks among candidate genes, ConsensusPathDB (CPDB) ( Herwig et al. 2016 ), STRING ( Szklarczyk et al. 2017 ), and GIANT_v2 ( Greene et al. 2015 ; Wong et al. 2018 ). For CPDB, human gene entrez ids from HGNC ( https://www.genenames.org/download/statistics-and-files/ ) were used in the “induced network modules” analysis using only high-confidence interactions and no intermediate nodes. The resulting PPI network was downloaded for the following analyses. For STRING, the v11 full human PPI network was downloaded from the STRING website ( https://stringdb-static.org/download/protein.links.full.v11.5/ ). Self-interactions of genes were removed, and a cutoff of 0.7 was set for the “combined score” to obtain high-confidence interactions, as defined by STRING database. For GIANT_v2, the full PPI network with global evidence was downloaded from the website ( http://giant-v2.princeton.edu/static//networks/global.dab ) and converted to text format with a Python script ( https://github.com/FunctionLab/flib/blob/master/dat.py ). A cutoff for interaction score of 0.15 was set to select ~25% high-confidence interactions from GIANT_v2 database. The final numbers of interactions are 240,970, 252,013, and 223,714 for CPDB, STRING and GIANT_v2, respectively. All types of interactions from the three databases were combined for the network construction, as previously described ( Cao et al. 2021 ). In addition to candidate genes, genes that connect to at least two candidate genes were included as intermediate nodes. Enrichment analyses were performed for 53 genes in the PPI network (including 16 candidate genes and 37 intermediate genes), using the overrepresentation analysis provided by CPDB ( Herwig et al. 2016 ). Enriched terms (e.g., Gene Ontology, pathway, or protein complex) containing at least two input genes were selected for further analysis. Enrichment p values were determined by CPDB using a hypergeometric test. Q values represent the adjusted p-values using the false discovery rate method.

Figure 1 describes the overall design of the project. After individual QC and filtering (see Methods for details), 139 individuals were selected from the IonTorrent dataset and assigned to the LRG/HRG groups ( Fig. 2A , 82 LRG and 57 HRG). The ages were similar between the LRG individuals [median age 37 years, interquartile range (IQR) 5 years] and HRG individuals [median age 36 years, IQR 4.5 years]. In the Illumina dataset, 142 individuals were selected, including 68 in LRG (median age 36 years, IQR 3 years) and 74 in HRG (median age 34 years, IQR 6 years) ( Fig. 2B ). Fifty patients were selected for testing in the Intermediate dataset, with median age of 34 years and IQR of 6 years. Following selection of patients, we filtered variants based on several QC matrices as specified in the Methods section in detail. After filtering, 9,082 genes in the IonTorrent dataset and 10,029 genes in the Illumina dataset were selected for further analysis, with 8,053 genes present in both datasets. Using the AVA,Dx pipeline, we calculated a gene score for each gene in each sample. Subsequent application of the K–S test on gene scores between LRG and HRG samples provided a ranked list of genes that were selected for the model evaluation in the leave-one-out cross-validation (LOOCV) analysis. We tested two gene score calculation methods (Sum or Product) and two classification algorithms (RF or SVM) in the AVA,Dx pipeline for each cohort (IonTorrent and Illumina; see Methods for details). Because the two gene scoring methods and two machine-learning algorithms had similar classification performance ( Supplemental Table 2 ), we elected to present an in-depth analysis using only the Sum gene score and the RF classifier in the following sections. For the IonTorrent dataset, using a top-10-gene model architecture (see Methods for detail) we achieved the highest ROC-AUC of 0.77 ( Fig. 3A ). A total of 16 unique candidate genes were included in at least one iteration of the LOOCV analysis ( Table 1 ). The median prediction score of HRG samples (i.e., the probability of a patient belonging to HRG) was 0.531, which was ~0.3 higher than the median prediction score of LRG samples (0.237) ( Fig. 3B ). By applying different prediction score cutoffs for HRG vs. LRG, we could adjust the sensitivity and specificity of the model ( Fig. 3B ). For example, if a high sensitivity is desired as a first-level screening to identify as many high aneuploidy risk patients as possible for secondary diagnosis, a lower prediction score cutoff can be applied. To this end, assuming that HRG samples are defined as the individuals with a prediction score larger than 0.3, we obtain a sensitivity of 0.75 and a specificity of 0.66 in this dataset. On the other hand, if high specificity is preferred to select patients with high confidence, a high prediction score cutoff for HRG can be applied. For example, when the prediction score for HRG patients is set to larger than 0.7, we achieve a specificity of 0.94 but at the expense of reducing sensitivity to 0.33 in this dataset. Next, we applied the same procedure to the Illumina dataset, for which the sequencing platform and variant calling procedure are different from the IonTorrent dataset. The top-5-gene architecture achieved a best predictive performance of 0.68 ROC-AUC ( Fig. 3C ) and included 7 unique candidate genes in all LOOCV iterations ( Table 1 ). The median prediction score of HRG samples was 0.684, which was ~0.3 higher than the median prediction score of LRG samples (0.373) ( Fig. 3D ). There was no overlap between the 16 candidate genes selected from the IonTorrent dataset and the 7 candidate genes selected from the Illumina dataset. To evaluate the model performance in samples not included in model training, we applied the RF classification model trained on the IonTorrent dataset to the Intermediate dataset representing samples from the same cohort not used in the model training. These individuals (Intermediate dataset) have intermediate aneuploidy rates ranging from 0.3 to 0.5 ( Supplementary Table S1 ). By applying the top-10-gene model trained using the IonTorrent dataset, the median prediction score of the Intermediate dataset was 0.42, which falls in between the median predicted HRG and LRG scores in the IonTorrent dataset ( Fig. 4A ). There was no correlation between the prediction scores and the aneuploidy rates among the samples in the Intermediate dataset ( Fig. 4B ). To understand the biological function of candidate genes selected by the classification models, we annotated the candidate genes, including their mutation burden in LRG vs HRG, expression pattern in early embryonic stages, known disease association, predicted constraint metrics, and their molecular functions ( Table 1 , Table S3 ). Among the 23 genes, 3 are implicated in meiosis-related functions ( MCM5 , FGGY , DDX60L ) and 15 are predicted to be meiosis genes ( Table 1 ). In addition, several of the genes are involved in embryonic developmental processes and are responsible for tumor proliferation, migration, and metastasis ( e.g., FOXA1 and RALGDS , Table 1 ). Below, we briefly describe the three candidate genes that were previously implicated in meiosis-related functions in model organisms ( Table 1 ). MCM5 is a member of the MCM family of chromatin-binding proteins and is a component of the MCM2-7 complex (MCM complex). The MCM complex is a putative replicative helicase that is essential for DNA replication initiation and elongation in eukaryotic cells ( Bochman and Schwacha 2008 , 2009 ). Our analysis of the IonTorrent data revealed that MCM5 has a higher average gene score in HRG than LRG samples. Of note, the score suggests the extent of protein function change rather than indicating deleterious or advantageous changes on the gene level. Moreover, MCM5 shows high expression during early embryonic development ( Table S3 ). In Drosophila, Mcm5 is involved in the maturation of DNA double-strand breaks into crossovers in the meiotic recombination pathway, which is separable from its role in mitotic replication ( Lake et al. 2007 ). FGGY is a member of the FGGY family of carbohydrate kinases ( Singh et al. 2017 ). In contrast to MCM5 , FGGY shows a lower average gene score in HRG than LRG in the IonTorrent dataset. In human, knockdown of FGGY promotes cell proliferation and invasion, facilitates tumorigenesis, and causes dysregulated cell energy metabolism and cytokine/chemotaxin transcription ( Zhang et al. 2019 ). mei-1 , the C. elegans homolog of FGGY , has distinct functions in meiosis: the loss-of-function mutant blocks meiotic spindle formation in embryos whereas a gain-of-function allele results in spindle defects during early mitotic cleavages in embryos ( Clark-Maguire and Mains 1994 ). MEI-1/MEI-2 microtubule-serving complex, katanin, is required for the organization or post nucleation processing of microtubules ( Johnson et al. 2009 ). Katanin is involved during meiotic spindle assembly to increase polymer amount from a relatively inefficient chromatin-based microtubule nucleation pathway ( Srayko et al. 2006 ). DDX60L is a member of the DExD/H-box helicase family and has conserved domains that mediate functions in RNA metabolism and acts as cytosolic sensors of viral nucleic acid ( Kato et al. 2006 ; Linder and Jankowsky 2011 ). DDX60L shows a lower average gene score in HRG than LRG in the IonTorrent dataset. mus301 , the Drosophila ortholog of human DDX60L , is required for the proper specification of oocytes and for progression through meiosis ( McCaffrey et al. 2006 ). Specifically, mus301 is involved in chromosome segregation in meiosis and in the repair of DNA double-strand breaks in both meiotic and mitotic cells ( McCaffrey et al. 2006 ). Following the gene function review, we constructed a PPI network using the 23 candidate genes from the two datasets. When including genes that are connected to at least two candidate genes in the PPI database, 16 out of 23 candidate genes were connected into a single network of 53 genes ( Fig. 5 , see Supplemental Table S4 for a full list of interactions). We then performed enrichment analysis on the 53 connected genes using the overrepresentation analysis provided by CPDB. Several enriched GO terms and pathways were found related to processes associated with aneuploidy in embryos, such as “microtubule organizing center organization (MTOC)” ( q = 0.016) and “DNA recombination” ( q = 0.009). Other GO terms that are related to cell cycle and cell division are also enriched, such as “mitotic cell cycle phase transition” ( q = 3.3 × 10 −4 ) and “regulation of cell division” ( q = 7.9 × 10 −4 ). A full list of enriched pathways and GO terms is shown in Supplemental Table S5 .

Female infertility affects an estimated 12% of women of reproductive age in the USA. Factors including genetic, anatomical, endocrine, and hemostatic alterations are reported as established risk factors of recurrent spontaneous abortions or multiple implantation failures ( Toth et al. 2018 ; RPL et al. 2018 ; Practice Committee of the American Society for Reproductive 2012 ; Carrington et al. 2005 ). Embryonic chromosomal abnormalities account for half of the infertility cases, with aneuploidy accounting for a significant proportion of IVF failures and early miscarriages ( Hassold and Chiu 1985 ). Recent studies show that genetic variants in several genes decrease the fidelity of chromosome segregation and predispose women to a higher incidence of egg aneuploidy (Reviewed in ( Biswas et al. 2021 )). This association suggests the potential to develop diagnostic tests where a list of genetic variants could be used to identify subfertile patients and provide individualized treatment options. For example, patients with fewer risk variants might have better possibilities of successful IVF treatment whereas patients with more risk variants would be counseled to have multiple cycles of IVF treatments at a younger age. In this study, we used two sets of whole-exome sequencing data obtained from IVF patients to identify candidate aneuploidy genes and to develop a model that could be used to predict aneuploidy risks for new patients. Both datasets showed promising results for classifying LRG vs HRG samples, with best models achieving ROC AUC of 0.77 and 0.68 for IonTorrent and Illumina datasets, respectively. This prediction performance is in agreement with the previously reported results of using our analysis pipeline in diseases; e.g., 58% of the over 3,000 patients have been correctly identified as having Crohn’s disease with 82% precision with an ROC AUC of 0.75 ( Wang et al. 2019 ). For our two datasets, candidate genes selected in the best prediction model do not overlap, despite the 28 shared patients ( Supplemental Table S1 ). Two factors could contribute to this difference: different gene scores ( e.g. , due to differences in exome enrichment kits and sequencing platforms, downstream variant calling and filtering methods, etc.) and heterogeneity among non-shared patients in the two cohorts. To elucidate this difference between the two models, we first compared the variants in the shared samples. These 28 samples have 86,176 and 262,485 variants in the Ion Torrent and Illumina datasets, respectively. Among the total 271,773 unique variants, only 28% are shared, a pattern that persists in the candidate genes selected for each dataset. For example, in the 7 genes of the Illumina model there are 117 variants among the 28 shared samples. However, only 45 of these (38%) are shared between Illumina and IonTorrent, while 65 (56%) are unique to the Illumina dataset. The large difference in the platform-specific variants across genomes deemed to be identical (same person) suggests exome enrichment kit, sequencing, and variant calling are major contributors to the gene score differences. In addition, cohort differences could also contribute different candidate gene selection. The K–S test, which we use for gene selection, targets maximum gene score differences between the HRG and LRG in the training samples; these differences, in turn, depend heavily on the balance of cohort samples. In our datasets, a total of 301 and 284 genes passed the K–S test to be used for the RF model construction in the IonTorrent and Illumina dataset, respectively. Among them, only eight genes were in common ( ALK , CCDC105 , CYP4A22 , KRTAP25-1 , PLB1 , PLCG2 , USP45 , WWC1 ) and only one was included in the final model ( PLCG2 in Illumina). The small overlap in the top genes that were used in the model construction explains the lack of shared genes in the final models for the two datasets. Despite the lack of shared genes, our network analysis showed extensive connection among the two sets of candidate genes, suggesting these genes are involved in shared underlying biological functions and pathways. Note that our models are specific to the cohort on which they are built. As such, the classification accuracy we report for the final model (see Methods for detail), trained on the entire training set for each platform, reflects the maximum achievable performance for this particular dataset. Indeed, the prediction performance dropped substantially when applying these models to different datasets, despite 28 samples present in both datasets. The AUC-ROC dropped to 0.59 when applying the IonTorrent model to the Illumina dataset, compared to 0.88 achieved in training; the AUC-ROC was 0.67 when applying the Illumina model to the IonTorrent dataset, compared to 0.75 achieved in training. Overall, these results suggest it is not yet possible to apply the current prediction models to patient samples whose data are not generated in the same format and outside of the current cohorts. Nevertheless, with the continuous improvement of the sequencing technology and the reduction of sequencing cost, we expect high-quality whole-genome sequencing will become standard practice. It is reasonable to expect the transferability of the prediction models based on the whole-genome sequencing data to improve dramatically with such advancement. Our study has several limitations. One is the relatively small sample sizes of the two cohorts we tested. The small sample sizes could partly contribute to the observation that the classifiers are not transferable to another dataset. Another limitation is we used hard cutoffs for the phenotype definition. We chose the current LRG and HRG definition because the sequencing data were generated in previous studies and this definition maximizes the samples to be included in the analysis ( Tyc et al. 2020a , b ; Tyc et al. 2021 ). Because age is a known risk factor for aneuploidy, in the future it would be useful to take the maternal age into account while defining groups in the training dataset. Third, to avoid selecting false-positive ancestry markers, we limited our data set to patients who were genetically inferred to be of European ancestry. It is unclear whether models developed with a uniform ancestry can be applied to other populations. In the future, disparate ethnicity-centered models could be built with more samples and a more comprehensive set of genetic variants. We expect the better prediction models will also uncover additional genes that contribute to female aneuploidy with greater statistical power. Although our sample size and data format (whole-exome sequencing) limit the transferability of our prediction models, genes selected in the two models provide insights into the biology of aneuploidy. Some candidate genes selected in our models, such as MCM5 , FGGY and DDX60L , are biologically relevant to meiotic processes. Established phenotypes in model organisms associated with individual gene mutations further support the role of many of the candidate genes identified here in processes related to meiosis and linked to aneuploidy ( Table 1 ). Furthermore, several candidate genes were previously implicated in the developmental processes. Mutations in several genes are also implicated in tumor proliferation and migration, such as FOXA1 and RALGDS , which are catalogued in Cancer Gene Census (CGC) ( Sondka et al. 2018 ). This result supports the close relationship between tumor development and aneuploidy ( Ben-David and Amon 2020 ). Aneuploidy, as a large-scale genetic variation, can allow cells to adapt in changing environments such as nutrient fluctuations and hypoxia, and contribute to cancer evolution ( Giam and Rancati 2015 ). In addition to their functional relevance, most of the candidate genes in the two datasets can be connected with a limited number of intermediate genes to form a single PPI network, suggesting these genes act in concert to regulate and/or contribute to different aspects of embryonic development ( Fig. 5 ). Notably, we identified an enrichment of candidate genes in MTOC and DNA recombination ( Supplemental Table S5 ). An interesting interaction detected is one between MCM5 and Aurora kinase A ( AURKA ). Although MCM5 is a DNA replication factor, it also functions to prevent centrosome reduplication, to ensure that spindles do not become multipolar ( Ferguson et al. 2010 ). In mouse oocytes, AURKA is essential for meiosis I because it establishes meiotic spindle pole number and structure ( Solc et al. 2012 ; Bury et al. 2017 ; Wang et al. 2020 ; Blengini et al. 2021 ). To better understand the biology of aneuploidy and embryonic development, the identified 23 candidate genes should be further pursued in the lab, ideally using genetically engineered mouse models. In summary, using IonTorrent and Illumina exome-sequencing datasets, we demonstrated the feasibility of developing a classifier for screening female IVF patients for risk of embryonic aneuploidy. The classification results can potentially provide better prognosis of IVF treatment success and guidance on individualized treatment options, such as recommendations for early or multiple IVF cycles, or for alternative family planning. In addition, our analyses identified candidate genes that can facilitate research to understand the biology of human aneuploidy and meiosis.

Embryonic aneuploidy, where an embryo has an abnormal number of chromosomes, is the predominant genetic abnormality leading to pregnancy failures in women of reproductive age ( Hassold and Chiu 1985 ). When failing to achieve pregnancy, many women undergo in vitro fertilization (IVF) treatment to select embryos with the correct number of chromosomes. Advanced maternal age, prior miscarriage, and failed IVF implantation cycles are currently used to indicate the risk of embryo aneuploidy during IVF treatments. Although aneuploidy risk increases with maternal age ( Hassold and Chiu 1985 ; Carp et al. 2001 ; Hassold and Hunt 2001 ; Hogge et al. 2003 ; Kuliev et al. 2011 ; Franasiak et al. 2014 ; Kubicek et al. 2019 ), the rate of producing aneuploid eggs varies among IVF patients for a given age ( Franasiak et al. 2014 ; McCoy et al. 2015a , b ; Tyc et al. 2020a , b ). Thus, current predictors ( i.e. , young age and good reproductive history) are not sufficient to ensure giving birth to healthy offspring. Recent studies have shown that maternal genetic variants are associated with the risk of embryonic aneuploidy (reviewed in ( Biswas et al. 2021 )), such as variants in Aurora kinase B and C ( AURKB , AURKC ) ( Nguyen et al. 2017 ), transducin-like enhancer (TLE) family member 6 ( TLE6 ) ( Alazami et al. 2015 ), and centrosomal protein 120 ( CEP120 ) ( Tyc et al. 2020a , b ). The association between maternal genetic variants and embryonic aneuploidy risk suggests the potential of using genomic data to predict embryonic aneuploidy risk in female IVF patients. However, because of the intricate and complex interactions among genes and genetic variants, a single variant or multiple variants in a single gene cannot be used to accurately diagnose infertility or predict aneuploidy risk in most patients. A model integrating multiple informative variants in a patient’s genome could potentially provide a more accurate prediction. Discovery of relevant biological patterns from high-throughput genomic data have advanced over the last several years via application of various machine learning tools ( Zhang et al. 2016 ; Telenti et al. 2018 ; Cai et al. 2020 ). By extracting complex feature patterns and constructing classification models based on case–control cohorts, machine-learning methods have revealed genetic causes and shown promising results in predicting individual predispositions for several disorders, such as Crohn’s disease ( Wang et al. 2019 ) and endometriosis ( Akter et al. 2019 ). Recently, high-throughput genotyping array and next-generation sequencing technologies have been used to identify genetic variants that are crucial for producing healthy eggs in female IVF patients ( McCoy et al. 2015a , b ; Nguyen et al. 2017 ; Tyc et al. 2020a , b ). These genetic variants can potentially be used to guide IVF patients in terms of the timing and scheme of treatment. In this study, we used machine-learning approaches to assess the power of classifying female IVF patients in term of their embryo aneuploidy rate based on their genomic variation. Specifically, we first measured the blastocyst ploidy status in embryos of women undergoing IVF treatment to calculate the aneuploidy rates and determined their risk of aneuploidy. Following the identification of genetic variants from patients’ whole-exome sequencing data, we used a derivative of the Analysis of Variation for Association with Disease (AVA,Dx) tool to construct machine learning models and evaluated the model for predicting aneuploidy status of unseen individuals.