Genome Sequencing is Critical for Forecasting Outcomes following Congenital Cardiac Surgery

Human subjects

All patients were diagnosed, phenotyped, and recruited by PCGC centers and participating regional hospitals into the PCGC Congenital Heart Disease Network Study (CHD GENES: ClinicalTrials.gov identifier NCT01196182; [https://clinicaltrials.gov/]). Informed consent was obtained from all participants or the participants’ guardians. Approval for human subjects research was obtained by the institutional review boards of participating centers, including Boston’s Children’s Hospital, Brigham and Women’s Hospital, Great Ormond Street Hospital, Children’s Hospital of Los Angeles, Children’s Hospital of Philadelphia, Columbia University Medical Center, Icahn School of Medicine at Mount Sinai, Rochester School of Medicine and Dentistry, Steven and Alexandra Cohen Children’s Medical Center of New York, Lucile Packard Children’s Hospital Stanford, University of California-San Francisco, University of Utah, and Yale School of Medicine. Automated CHD phenotype classification was performed on 14,765 PCGC participants. A subset of these participants who had both exome sequencing and perioperative data (2,253) was used for network analyses.

Clinical phenotypes

Cardiac diagnoses were obtained from review of echocardiogram, cardiac MRI, catheterization, and operative reports at the time of enrollment into the PCGC^4,5. Detailed cardiac diagnoses for each patient were coded using the Fyler system¹⁵. Extra-cardiac anomalies (ECAs) were identified at the time of PCGC enrollment^4,5 (Supplemental Table 1). Any structural anomaly that was not acquired was classified as an extra cardiac anomaly (ECA).

Post-operative variables

For patients undergoing open heart surgery, surgical and hospitalization data were obtained from participating centers using the local data collected for submission to the Society of Thoracic Surgeons Congenital Heart Surgery Database (STS-CHSD)¹⁷. A total of 59 surgical complication variables were extracted for analysis. The size of the final data set was constrained to 2,253 patients, such that all patients had WES and surgical variables had no more than 10% missing data. Most patients had multiple cardiac surgeries. A patient was scored as having an adverse event or surgical complication (e.g., prolonged mechanical ventilation) if that event occurred for any surgery at any age (Supplemental Table 2).

Surgical complexity

Surgical complexity is a well-known driver of mortality and morbidity. In response, the STS-European Association for Cardio-Thoracic Surgery (STAT) has created risk assessment categories in which procedures are grouped based on similar mortality rates¹⁸. STAT categories range from 1 to 5, with STAT1 representing the procedures with the lowest mortality rates and STAT5 representing the procedures with the highest mortality rates.

CHD classification

The PCGC has classified cardiac diagnoses for over 14,000 CHD probands using the Fyler coding system, which describes the congenitally malformed heart using a vocabulary of more than 3,000 possible phenotypic descriptors¹⁵. While this system allows for highly granular descriptions of heart defects, we hypothesized that condensing these terms into a few clinically relevant phenotypic categories might render them more tractable for outcomes analyses. Thus, we sought to automate cardiac phenotype classification across the entire PCGC cohort, assigning each patient to a single category. To do so, we used five major cardiac categories derived from a previous PCGC study⁵: left ventricular outflow tract obstructions (LVO), laterality and heterotaxy defects (HTX), atrioventricular canal defects (AVC), conotruncal defects (CTD), and other defects (OTH), which includes simple atrial septal defects and more complex heart defects not assigned to the other four categories^5,15. Each participant was assigned uniquely to one of the five phenotypic categories.

A gradient-boosted decision tree model was built to automatically classify PCGC probands (14,765) into one of these five CHD categories (see Supplemental Methods, Supplemental Figures 1, 2). Model learning and classification was performed using an ensemble-based method using the XGBoost library (v1.5)¹⁹. We created an analysis toolkit to streamline XGBoost training, grid-based parameter optimization, and performance evaluation. The truth set for training the classifier included 3,000 CHD patients, 2,752 PCGC patients previously assigned into the five CHD categories⁵ and 248 randomly selected PCGC patients that were manually reviewed and assigned to a CHD phenotype category. A gradient-boosted probabilistic patient classifier was built with XGBoost¹⁹ using the 3,000 patients and 698 Fyler features (Supplemental Table 3). Model training was performed with five-fold cross-validation and replacement subsampling (Supplemental Table 4). Model accuracy was assessed by comparing the patient’s known phenotype category label derived from the literature⁵ to its predicted phenotype category label (Supplemental Figure 3, Supplemental Table 5). Single-class prediction accuracy for the training data was higher for HTX (98.9%), AVC (97.8%), and LVO (97.7%) than for CTD (95.3%) and OTH (91.9%), where a low level of ambiguity occurred. Overall classification accuracy was 97.7% with a specificity of 99.3%, and sensitivity of 97.7% (Supplemental Tables 6). We then applied the trained classifier to 14,765 PCGC CHD patients with Fyler descriptors (Supplemental Tables 7, 8). Final classification of the 3,000 training patients was identical to the original training predictions. The five phenotype categories remained generally proportional between the training data and the full data set, with a maximum observed difference of 4.2% for CTD patients.

AI-based scoring of predicted damaging genetic variants

AI-based identification of candidate disease-causing genotypes was performed using Fabric GEM¹⁶ (Fabric Genomics, Oakland, CA). GEM incorporates Human Phenotype Ontology (HPO) terms, sex, genotype frequency (gnomAD), evolutionary conservation, Online Mendelian Inheritance in Man (OMIM), GnomAD, and ClinVar information in a probabilistic AI framework to identify the most likely genetic variant or genotype that explains the patient’s disease phenotype. Because WES are difficult substrates for CNV calling, we restricted our analyses to SNVs and short indels. HPO terms utilized in the GEM analyses were based on each patient’s Fyler phenotypes, which were mapped to HPO terms using the Clinithink software package (Clinithink, London).

GEM’s gene scores are log₁₀ transformed Bayes factors²⁰ that summarize the relative support for the hypothesis that the prioritized genotype damages the gene in which it resides and explains the patient’s phenotype versus the hypothesis that the variant neither damages the gene nor explains the patient’s phenotype. We used a stringent GEM score of ≥ 1.0 to represent a likely pathogenic genotype. A recent genomic analysis of critically ill newborns¹⁶ showed that a GEM score of ≥ 1.0 identified 90% of all true positive damaging variants, with a median of two candidate variants per patient¹⁶. Gene penetrance for GEM calculations was set to 0.95 to enforce strict consideration of known dominant and recessive disorders. For downstream analyses, damaging genetic variants were classified as de novo, dominant, or recessive/biallelic variants based on their inheritance pattern in trios. Dominant and de novo damaging variants were required to have a frequency of < 1/10000 in gnomAD databases (v2.1, v3.1) and most variants were not observed. Overall, we identified damaging de novo or recessive genotypes in 10.56% of the study cohort (Supplemental Tables 9, 10), in line with previous studies that utilized different methods of defining pathogenicity^4,5,8,9. Damaging genetic variants were assigned to several functional gene pathways. Gene lists for gene pathways were obtained using the reactome pathway browser. Gene lists are shown in Supplemental Table 11 and have been previously described^5,8,9 There is overlap between gene lists, with some genes represented in more than one gene pathway/category (Supplemental Figure 4).

Probabilistic graphical models

Probabilistic graphical models (PGMs) provide a robust explainable AI methodology capable of discovering and quantifying additive and synergistic effects amongst broad classes of variables. For the work presented here, we used a form of PGMs known as a Bayesian networks²¹. Bayesian networks are fully transparent, and their graphical representation offers an intuitive and visual mechanism for understanding the relationships between variables and the impacts of multiple variables on outcomes of interest^22–29. Moreover, Bayesian networks offer practical advantages over regression approaches, by capturing the entire joint probability distribution of the data, encompassing all interrelationships among the variables incorporated in the non-linear model. Thus, a single network can be used to explore any combination of variables as a target outcome in one query and then as a risk factor for a different target outcome, all within the same model. For more on these points see ^21,30,31.

Feature selection

Single variables, such as damaging genetic variants in chromatin-modifying genes, were tested for conditional dependency with phenotype variables using exact Bayesian Networks. Each conditional probability [e.g., the probability of LVO given a damaging de novo chromatin variant: P (LVO | chromatin dGV) was estimated as the median conditional probability from 1000 independent networks. Conditional probability estimates were divided by the baseline probability for each respective phenotype to obtain absolute risk ratios. Associations with absolute risk ratios ≥1.0 were selected for further analyses. Surgery-related variables associated with gene categories and CHD phenotypes were identified in a similar way. Each surgical feature was tested individually as a conditional variable with genetic variables and each of the five specific CHD phenotypes. All conditional variables, individually or as composite variables (e.g., mortality and ECA), were required to have at least six events.

Network construction

Bayesian networks were created for each CHD phenotype category. Each network included genetic and surgical variables identified in the feature selection stage. All-cause mortality was also included in each network. All input conditional variables were coded as presence/absence. A small number of missing surgical values (< 10% for any single variable) was imputed using a K-nearest neighbors approach (k = 10). The structure of each network was learned with the Silander-Myllymaki exact algorithm with Bayesian information criterion (BIC) scoring³² or by a greedy hill-climbing method for the large networks with more than 15 nodes. Posterior probabilities were network propagated using exact inference. Network structure learning and belief propagation were performed with the bnstruct and gRain R packages^33,34. To improve convenience and functionality, we created the BayesNetExplorer.jl package which implements network structure learning and belief propagation methods and provides tools for feature selection, risk estimation, graphics, and other network analysis tasks.

Risk Calculations

We used two risk ratios (RR) to summarize our risk estimates, absolute and relative RR. For example, the absolute RR of a LVO phenotype given a damaging genotype in a chromatin-modifying gene is as follows: . The relative RR estimates the relative change in mortality risk for LVO patients with damaging mutations in chromatin-modifying genes, compared to similar patients without a damaging chromatin genotype: . Final risk ratios and their confidence intervals are reported as the median and 95% confidence interval from an empirical distribution of risk ratio estimates. The empirical distributions are created by randomly resampling the data set with replacement and recreating 1000 independent networks and risk estimates. A Laplace correction (k = 1) or network smoothing value was used to prevent zero-state probability estimates during bootstrapping.

Refining the genetic architecture of CHD

For genetic and outcomes analysis, the study population consisted of 2,253 PCGC probands (1992 trios, 12 duos, 245 singletons) with both exome sequencing and surgical outcomes data, classified into five CHD phenotype categories. The AI-based genome analysis tool GEM¹⁶ identified predicted damaging de novo genotypes in 238 participants (10.6% of the cohort). A total of 131 damaging de novo / dominant and 198 damaging recessive/biallelic genetic variants were discovered (Supplemental Tables 9, 10). There were 17 genes with damaging de novo variants in two or more patients. The most commonly recurrent de novo variants were in known CHD-related genes such as KMT2D (11), CHD7 (6), RAF1 (3), JAG1 (3), and TAB2 (3). Biallelic damaging genotypes were observed in multiple patients for several genes including DYNC2H1 (3), DNAH5 (3), LAMA2 (3), GDF1 (2), and IFT140 (2).

We discovered that CHD phenotype categories were enriched for damaging genetic variants in specific gene pathways/categories (Table 1). For example, the LVO class was enriched 1.61-fold (CI 1.41-1.81) for damaging de novo genotypes in chromatin-modifying genes, with this signal driven primarily by patients with hypoplastic left heart syndrome (HLHS), a subset of LVO where these genotypes were enriched 1.92-fold (CI 1.32-2.52). While previous studies implicated damaging chromatin-modifying gene variants in CHD cohorts at large^4,5,8,9, our analyses here help to define the specific CHD subtypes most influenced by damaging variants in chromatin-modification genes. The LVO phenotype class was also enriched for de novo genotypes in WNT genes (2.13-fold, CI 1.86-2.40), signal transduction genes (1.51-fold, CI 1.01-2.01), and a curated list of genes known to cause CHD (1.22-fold, CI 1.11-1.33; see Supplemental Tables 9 and 10). Notably, damaging genotypes in these pathways were not enriched in HLHS patients, further underscoring the complex genetic landscape underlying CHD.

The HTX phenotype class was enriched for damaging recessive/biallelic variants in cilia-related genes (2.63-fold, CI 2.06-3.20) and showed proportionally higher enrichment in the subset of motile cilia genes modulated by FOXJ1 (6.89-fold, 3.30-10.36), findings consistent with previous reports^5,8,35. The OTHER phenotype class was enriched for damaging de novo variants in chromatin-modifying genes (1.85-fold, CI 1.44-2.26) and the curated CHD genes (1.59-fold, CI 1.38-1.80) lists. We were underpowered to detect enrichment in damaging genotypes in the AVC and CTD phenotype classes.

Damaging genotypes impact surgical outcomes

To further explore the relationships between genetic and clinical variables, we utilized Bayesian networks, a powerful statistical framework that can model complex dependencies, including non-linear relationships and indirect associations, in a probabilistic manner. In Bayesian networks, variables are depicted as nodes in a graph and conditional dependencies between variables are represented by the edges connecting those nodes. Once a network is constructed, the impact of any combination variables on any selected outcome can be quantified, while controlling for the effects of other variables incorporated into the network. The networks describing the conditional dependencies between damaging genotypes, CHD phenotypes and post-operative variables are shown in Figure 1a, b.

• Download figure
• Open in new tab

Figure 1. Damaging chromatin and cilia genotypes predict adverse post-operative outcomes in the context of CHD phenotypes and surgical complexity.

a, An exact Bayesian network depicting the relationship among damaging de novo genetic variants in chromatin-modifying genes (green), phenotypes: LVO, HLHS, and ECAs (blue), surgical STAT4 or STAT5 category (red), and surgical outcomes (orange). b, An exact Bayesian network depicting the relationship among damaging recessive genetic variants in cilia genes (green), phenotypes: laterality defects (HTX) and extra cardiac anomalies (ECAs) (blue), surgical STAT4 category (yellow), and surgical outcomes (orange). Directed acyclic graphs were moralized and displayed as non-directional networks. c, Relative risk ratios for adverse post-operative outcomes and CHD phenotypes or surgical complexity, comparing probands with and without damaging genotypes. Empirical ninety-five percent confidence intervals (CI 5, 95) are based on 1000 resampled network-based probabilities (see Methods). Abbreviations: Chromatin dGV: de novo damaging genotypes in chromatin-modifying genes, LVO: left ventricular outflow tract obstruction, HLHS: hypoplastic left heart syndrome Cilia dGV: biallelic damaging genotypes in cilia-related genes, ECA: extra cardiac anomaly, HTX: heterotaxy/laterality defects, MORT: mortality, STAT4: surgical STAT4 category, STAT4-5: surgical STAT 4 or STAT5 category, VENT: post-operative ventilation time >7 days.

Damaging genotypes in chromatin-modifying and cilia-related genes (defined by a GEM¹⁶ score > 1.0) increased the probability of severe adverse clinical outcomes following congenital cardiac surgery, including mortality, cardiac arrest, and prolonged mechanical ventilation (> 7 days post-surgery). For example, damaging de novo chromatin genetic variants increased the probability (relative risk) of mortality 1.81-fold (CI 1.50-3.21), cardiac arrest 1.74-fold (CI 1.40-2.94) and prolonged ventilation 1.65-fold (CI 1.41-2.27). Likewise, damaging recessive/biallelic cilia genotypes increased the probability of mortality 1.40-fold (CI 1.09-2.10), cardiac arrest 1.50-fold (CI 1.13-2.34) and prolonged ventilation 1.43-fold (CI 1.09-1.96). Reciprocally, the absence of a damaging genotype was protective for these adverse post-operative outcomes. Thus, for a proband without a damaging de novo chromatin genotype, the relative risk ratio for mortality was 0.55 (CI 0.31-0.69), 0.58 (CI 0.34-0.72) for cardiac arrest and 0.61 (CI 0.44-0.72) for prolonged ventilation. For a proband without at damaging recessive/biallelic cilia genotype, the relative risk ratio for mortality was 0.72 (CI 0.48-0.91), 0.63 (CI 0.43-0.88) for cardiac arrest and 0.70 (CI 0.51-0.92) for prolonged ventilation.

Damaging genotypes impact surgical outcomes in the context of CHD phenotypes

More broadly, considering mortality in the context of CHD phenotypes, LVO patients who died were 2.26-fold (CI 1.22-2.94) more likely to harbor a damaging de novo chromatin genotype, while HTX patients who died were 2.80-fold (CI 2.45-3.15) more likely to harbor a damaging recessive/biallelic cilia genotype (Figure 1c). Similarly, damaging chromatin or cilia genotypes were overrepresented in probands with LVO, HLHS and HTX who experienced cardiac arrest or prolonged post-operative ventilation (Figure 1c). Specifically, HTX patients who arrested post-operatively were 3.34-fold (CI 1.44-5.49) more likely to harbor a damaging recessive/biallelic cilia genotype, compared to similar patients without a damaging cilia genotype. Collectively, these findings demonstrate that genome sequencing data are critical for predicting severe post-operative events in the context of specific CHD phenotypes and the highest risk congenital heart surgeries.

Damaging genotypes impact surgical outcomes in the context of extracardiac phenotypes

Given the recognized impact of ECAs on outcomes following congenital cardiac surgery^10,14,36, we also explored the relationship between ECAs and adverse post-operative outcomes in the context of genotypes and CHD phenotypes. ECAs increased the probability (relative risk) of mortality 2.85-fold (CI 1.47-2.86) and prolonged ventilation 1.72-fold (CI 1.63-1.72) following congenital cardiac surgery. Consistent with previous findings¹⁰, damaging de novo genetic variants, including de novo variants in chromatin-modifying genes, were enriched in probands with ECAs 1.47-fold (CI 1.44-1.48) and 2.09-fold (CI 2.08-2.11), respectively. By contrast, damaging recessive/biallelic cilia genotypes were not enriched in probands with ECAs (0.84-fold; CI 0.85-1.13).

We also examined reciprocal effects, i.e., impact of predicted damaging genotypes on ECAs and adverse outcomes. For example, probands with damaging de novo chromatin genotypes identified by GEM were 2.49-fold (CI 2.17-4.99) more likely to have an ECA and die, compared to probands without a damaging chromatin genotype, and 2.44-fold (CI 1.97-3.41) more likely to have an ECA and prolonged ventilation (Figure 2). Moreover, a damaging recessive/biallelic cilia genotype identified by GEM increased the probability of mortality in probands with an ECA 1.48-fold (CI 1.02-2.85), compared to similar probands without a damaging cilia genotype, and increased the probability of prolonged ventilation in the presence of an ECA 1.52-fold (CI 1.06-2.51). Additionally, a damaging cilia genotype increased the probability of prolonged ventilation in HTX patients with an ECA 4.01-fold (CI 1.67-10.6), compared to similar patients without a damaging cilia genotype (Figure 2). Taken together, these findings demonstrate that damaging genotypes in chromatin and cilia genes predict severe post-operative events in the setting of ECAs.

• Download figure
• Open in new tab

Figure 2. Damaging chromatin and cilia genotypes predict adverse post-operative outcomes in the context of extracardiac anomalies.

Relative risk ratios for adverse post-operative outcomes and extracardiac anomalies (ECAs), comparing probands with and without damaging genotypes in chromatin-modifying or cilia-related genes. The \\ symbol represents CI95 that exceeds the x-axis range.

The number of probands experiencing adverse outcomes and harboring damaging gene pathway variants in the AVC, CTD, and OTHER categories was too low to warrant generation of Bayesian networks for outcomes prediction in these CHD phenotypes. Consequently, larger cohorts are necessary to adequately predict the impact of genetics on outcomes for these CHD phenotypes. However, damaging genotypes in a number of gene pathways/categories, such as FOXJ1-controlled genes, high murine heart expression genes (HHE), WNT signaling genes, NOTCH signaling genes, and genes in a curated CHD gene list were predictive of mortality for the most complex surgical categories (Figure 3). Damaging genotypes in signal transduction and TGF-β pathways were not predictive of mortality. Taken together, these findings highlight the value of genomic data for predicting adverse outcomes following congenital cardiac surgery, especially in the context of CHD phenotypes, ECAs and surgical complexity.

• Download figure
• Open in new tab

Figure 3. Damaging genotypes in various gene categories/pathways are predictive of mortality for the most complex surgical procedures.

Relative risk ratios for adverse post-operative outcomes and surgical complexity, comparing probands with and without damaging genotypes in various gene pathways or categories. Gene lists are described in Supplemental Table 11 and have been previously published.^5,8,9 There is overlap between gene lists, with some genes represented in more than one gene pathway/category (Supplemental Figure 4).