Deciphering the digenic architecture of congenital heart disease using trio exome sequencing data

Abstract

Congenital heart disease (CHD) is the most common congenital anomaly and a leading cause of infant morbidity and mortality. Despite extensive exploration of the monogenic causes of CHD over the last decades, ∼55% of cases still lack a molecular diagnosis. Investigating digenic interactions, the simplest form of oligogenic interactions, using high-throughput sequencing data can elucidate additional genetic factors contributing to the disease. Here, we conducted a comprehensive analysis of digenic interactions in CHD by utilizing a large CHD trio exome sequencing cohort, comprising 3,910 CHD and 3,644 control trios. We extracted pairs of presumably deleterious rare variants observed in CHD-affected and unaffected children but not in a single parent. Burden testing of gene pairs derived from these variant pairs revealed 29 nominally significant gene pairs. These gene pairs showed a significant enrichment for known CHD genes (p < 1.0 × 10^-4) and exhibited a shorter average biological distance to known CHD genes than expected by chance (p = 3.0 × 10^-4). Utilizing three complementary biological relatedness approaches including network analyses, biological distance calculations, and candidate gene prioritization methods, we prioritized 10 final gene pairs that are likely to underlie CHD. Analysis of bulk RNA-sequencing data showed that these genes are highly expressed in the developing embryonic heart (p < 1 × 10^-4). In conclusion, our findings suggest the potential role of digenic interactions in CHD pathogenesis and provide insights into unresolved molecular diagnoses. We suggest that the application of the digenic approach to additional disease cohorts will significantly enhance genetic discovery rates.

Keywords: biological proximity; congenital heart disease; digenic interaction; enrichment analysis; trio exome sequencing.

Graphical abstract

Figure 1

Study workflow (A) Presumably deleterious rare variants were extracted from the CHD-control trio ES data. (B) Gene pairs for each CHD-affected individual and control were generated using variant pairs that were not observed in their parents. (C) Gene pairs with >5 observations were employed in the burden testing, resulting in 29 significant gene pairs enriched in CHD. (D) Ten final gene pairs containing known genes and/or genes prioritized through three complementary biological relatedness approaches.

Figure 2

Burden testing of the digenic interactions in the combined CHD and control cohort (A) PCA plot showing the distributions of genetically determined ancestries of CHD-affected probands and unaffected siblings from the SPARK dataset. (B) Distribution of presumably deleterious variants used in the analysis across different predicted variant impact categories and gnomAD MAF bins. The high-impact category includes frameshift, stop gain, stop loss, and essential splice site variants. (C) Histograms displaying the number of gene pairs per sample carried by CHD-affected probands (left) and controls (right). Gene pairs with ≥10 observations were binned together. (D) QQ plot of the burden testing results. The red dashed line shows the Bonferroni-corrected p-value threshold (p = 5.59 × 10⁻⁵), whereas the green dashed line indicates the nominal significance level (p = 0.05).

Figure 3

Assessment of known CHD gene enrichment, biological distance to known CHD genes, and network analysis of the 29 significant CHD gene pairs (A) Density plot showing the distribution of the number of gene pairs with at least one known CHD gene from 29 random gene pairs in 10,000 resampling iterations. The dashed red line represents the number of CHD gene pairs with a known CHD gene. (B) Distribution of the average distances of 29 random gene pairs (D_{random_pair}) to known CHD genes obtained from 10,000 resampling iterations. The dashed red line represents the average distance of the significant CHD gene pairs (D_{candidate_pair}, n = 29) to known CHD genes. (C) Subnetwork identified by IPA software. Known CHD genes and candidate CHD genes from significant CHD gene pairs are shown in purple and red, respectively, whereas the blue color indicates known CHD genes identified in the extended network but not observed in the CHD gene pairs.

Figure 4

Prioritization of the final candidate CHD gene pairs (A) Sankey plot displaying the prioritization steps of the final 10 CHD gene pairs. Forty genes from the significant CHD gene pairs were prioritized if they have been previously reported as known human CHD genes or prioritized by all three complementary methods (HGC, network analysis with IPA, and Enrichr-KG and ToppGene). (B) PCA plot showing the distribution of the samples with the final 10 CHD gene pairs. (C) Percentiles (pth) of developing heart expression levels of the 40 genes from the 29 CHD gene pairs. The dashed vertical line indicates the average percentile of expression levels of the 475 known human CHD genes.