Epistasis regulates genetic control of cardiac hypertrophy

Deep learning of UK Biobank cardiac imaging quantifies left ventricular hypertrophy

We accessed all cardiac magnetic resonance images from the UKBB substudy (44,503 people at the time of this analysis)³⁴. We focused on the largest ancestry subset of 29,661 unrelated individuals (summary characteristics in Supplementary Table 1) and analyzed the most recent image per individual. We leveraged a recent deep learning model²⁸ to quantify LV hypertrophy from these 29,661 multislice cine magnetic resonance images (Fig. 2a). A fully convolutional network had been previously trained for image segmentation and was evaluated on manual pixelwise-annotations of images from 4,875 UKBB participants²⁸. This fully convolutional network learns features across five different resolutions through sequential convolutional layers interspersed with non-linearities, and has displayed accurate performance compared to cardiac segmentation by human experts²⁸. Using this segmentation model, we extracted areas of the LV chamber wall in each slice of the short axis image at the end of diastole. Areas extracted from each image slice in the same image stack were then integrated to calculate the heart muscle volume, which we converted to the LV mass using a standard density of 1.05 g/mL³⁵. This was normalized by body surface area, estimated using the Du Bois formula³⁶, to obtain the LV mass index (LVMi, Extended Data Fig. 1). Details regarding this analysis can be found in Methods.

• Download figure
• Open in new tab

Fig. 2: Low-signal signed iterative random forest (lo-siRF) prioritizes risk loci and epistatic interactions for left ventricular hypertrophy

a-e, Workflow of low-signal signed iterative random forest (lo-siRF). a, Lo-siRF took in as input single-nucleotide variant (SNV) data and cardiac MRI-derived left ventricular mass indexed by body surface area (LVMi) from 29,661 UK Biobank participants. b, Dimension reduction was performed via a genome-wide association study (GWAS) to concentrate the analysis on a smaller set of SNVs. c, LVMi was binarized into high and low LVMi categories according to three different binarization thresholds (represented by the stacked boxes). d, For each of the three binarization thresholds, a signed iterative random forest (siRF) was fitted using the GWAS-filtered SNVs to predict the binarized LVMi phenotype. Other popular prediction methods including polygenic risk scores, machine learning (ML), and deep learning (DL) models were also trained and evaluated as baseline comparison methods. The validation prediction accuracy of siRF was shown to be on par or better than these baseline comparisons, prior to interpreting the model fit. e, SNVs used in the fitted signed iterative random forest were aggregated into genetic loci based on annotations using ANNOVAR⁹⁰. Genetic loci and pairwise interactions between loci were finally ranked according to their importance across the three signed iterative random forest fits, as measured by our proposed stability-driven importance score. f, Lo-siRF-prioritized risk loci and epistatic interactions. (1) Loci stably prioritized by lo-siRF as epistasis participants are highlighted in green. (2) nIndSigSNVs, the number of independent significant SNVs that are stably prioritized by lo-siRF across the three different LVMi binarization thresholds (panel c). (3) nSNVs, the number of candidate SNVs extracted by FUMA⁴⁴ (v1.5.4) in strong LD (r² > 0.6) with any of the lo-siRF-prioritized independent significant SNVs. (4) Lo-siRF p-value, the mean p value from lo-siRF, averaged across the three LVMi binarization thresholds. (5) Lo-siRF p-value (excl. hypertension), the mean p value from lo-siRF when excluding hypertensive individuals from the analysis, averaged across the three LVMi binarization thresholds. (6) Max CADD, the maximum CADD⁴⁷ score of SNVs within or in LD with the specific locus. A high CADD score indicates a strong deleterious effect of the variant. A threshold of 12.37 has been suggested by Kircher et al.⁴⁷. (7) Min RDB, the minimum RegulomeDB⁴⁶ score of SNVs within or in LD with the specific locus. RDB is a categorical score to guide interpretation of regulatory variants (from 1a to 7, with 1a being the most biological evidence for an SNV to be a regulatory element)^44,46. (8) The top-ranked SNV or SNV-SNV pair showing the highest occurrence frequency (Extended Data Fig. 4) averaged across lo-siRF fits from the three LVMi binarization thresholds. A full list of lo-siRF-prioritized SNVs and SNV-SNV pairs can be found in Extended Data 3. (8) Genomic location (hg38) and GWAS statistics information (using PLINK³⁷) of the top SNV for each lo-siRF-prioritized locus. Abbreviations: MAF, minor allele frequency; NEA/EA, non-effect-allele/effect-allele; SE, standard error. (10) nPartnerSNVs, number of partner SNVs that interact with the given SNV in lo-siRF. These SNV-SNV pairs interacted in at least one lo-siRF decision path across every LVMi binarization threshold (details in Methods).

Low-signal signed iterative random forests prioritize epistatic genetic loci

We developed low-signal signed iterative random forests (lo-siRF, Fig. 2a-2e) to prioritize statistical epistatic interactions from the extracted LV mass and single-nucleotide variants (SNVs) from UKBB. Given the inherent low signal-to-noise ratio and aforementioned challenges, lo-siRF aims to recommend reliable candidate interactions for experimental validation rather than directly assessing claims of statistical significance from data. This prioritization pipeline is guided by the PCS framework²⁶ and builds upon signed iterative random forests^22,23, a computationally-tractable algorithm to extract predictive and stable nonlinear higher-order interactions that frequently co-occur along decision paths in a random forest. More specifically, lo-siRF proceeds through four steps:

1. Dimension reduction (Fig. 2b): we combined the results of two initial genome-wide association studies, implemented via PLINK³⁷ and BOLT-LMM³⁸ (Extended Data Fig. 2, Extended Data 1) to reduce the interaction search space from 15 million imputed variants down to 1405 variants (Extended Data 2). Details can be found in the Methods section Lo-siRF step 1: Dimension reduction of variants via genome-wide association studies.

2. Binarization (Fig. 2c): we partitioned the LV mass measurements into high, middle, and low categories using three different partitioning schemes (Supplementary Table 2). The partitioning enabled us to transform the original low-signal regression problem for a continuous trait into a relatively easier binary classification task for predicting individuals with high versus low LV mass measurements (omitting the middle category). This transformation is necessary to obtain a sufficient prediction signal, ensuring that the model indeed captures pertinent information about reality (Supplementary Table 3). Further justification and details on the partitioning can be found in the Methods section Lo-siRF step 2: Binarization of the left ventricular mass phenotype.

3. Prediction (Fig. 2d): we trained a signed iterative random forest using the 1405 GWAS-filtered SNVs to predict the binarized LV mass measurements. The learnt model yields on average the highest (balanced) classification accuracy (55%), area under the receiver operator characteristic (0.58), and area under the precision-recall curve (0.57) compared to other common machine learning prediction algorithms (Supplementary Table 4). Details about the model and prediction check can be found in the Methods section Lo-siRF step 3: Prediction.

4. Prioritization (Fig. 2e): we developed a stability-driven feature importance score (Extended Data Fig. 3), which leveraged the fitted signed iterative random forest and a permutation test, to aggregate SNVs into genetic loci and prioritize interactions between genetic loci. This importance score provides the necessary new interpretable machine learning ingredient to complete the lo-siRF discovery pipeline. Details can be found in the Methods section Lo-siRF step 4: Prioritization.

Additional discussion of the philosophy and modeling decisions driving lo-siRF can be found in Supplementary Note 1, an interactive HTML webpage hosted at https://yu-group.github.io/epistasis-cardiac-hypertrophy/. The webpage also provides a comparison of lo-siRF to alternative epistasis detection methods, including an exhaustive regression-based pairwise interaction search^39,40 and MAPIT⁴¹, demonstrating the challenges and limitations of existing methods for analyzing low-signal, complex phenotypes.

Lo-siRF identified six genetic risk loci that exhibited stable and reliable associations with LV mass (Fig. 2f). Because these loci are either located within a gene body or in between two genes (Fig. 3a), for convenience we denote these loci by their nearest genes. Notably, out of the six loci, three (TTN, CCDC141, and IGF1R) were prioritized by lo-siRF as epistatic loci. These loci not only interact with other loci, but also marginally affect LV mass. The other three lo-siRF-prioritized loci are LOC157273;TNKS, MIR588;RSPO3, and LSP1. The LOC157273;TNKS locus is located within the intergenic region between genes LOC157273 and TNKS (semicolon indicates intergenic region). This locus was prioritized by lo-siRF to be hypostatic (i.e., effects are deemed stable by lo-siRF only when interacting with the CCDC141 locus). Interestingly, all three identified epistatic interactions involved the CCDC141 locus (Fig. 3a, green links in circle 1). Furthermore, while the MIR588;RSPO3 and LSP1 loci lacked evidence for epistasis by lo-siRF, they were each identified to be marginally associated with LV mass. The specific prioritization order of these loci can be found in Supplementary Table 5, and details regarding the direction or sign of the interactions can be found in Supplementary Note 1. In total, lo-siRF identified 283 SNVs located within the six loci (Extended Data 3, Extended Data Fig. 4). Ninety percent of the 283 SNVs have previously been shown to harbor multiple distinct cardiac function associations⁴² in phenome-wide analyses (e.g., pulse rate, Extended Data 3), suggesting a strong likelihood that these lo-siRF-prioritized loci contribute to determining cardiac structure and function.

• Download figure
• Open in new tab

Fig. 3: Lo-siRF finds epistatic interactions between genetic risk loci for left ventricular hypertrophy

a, Circos plot showing the genetic risk loci identified by lo-siRF (green, circle 2) and regions after clumping FUMA-extracted SNVs in LD (r² > 0.6) with lo-siRF-prioritized SNVs (black, circle 3). Circle 1 shows the top 300 epistatic SNV-SNV pairs with the highest frequency of occurrence in lo-siRF (green), SNV-gene linkages (FDR < 0.5) based on GTEx⁴⁸ V8 cis-eQTL information from heart and skeletal muscle tissues (purple), and 3D chromatin interactions⁴⁴ based on Hi-C data of left ventricular tissue obtained from GSE87112. Circle 5 and 6 show bar plots of the occurrence frequency and number of partner SNVs in epistasis (normalized by the maximum value of the corresponding locus) identified by lo-siRF, respectively. Circle 7 shows the ChromHMM⁴⁵ core-15 chromatin state for left ventricle (LV), right ventricle (RV), right atrium (RA), and fetal heart (Fetal). Circle 8 shows the GWAS Manhattan plot from PLINK³⁷ (circles) where only SNVs with p < 0.05 are displayed. The 283 lo-siRF-prioritized SNPs and their LD-linked (r² > 0.6) SNVs are color-coded as a function of their maximum r² value. A portion of these LD-linked SNPs (the outer heatmap layer in circle 8) are extracted from the selected FUMA reference panel (thereby with no GWAS p-values). SNVs that are not in LD (r² ≤ 0.6) with any of the 283 lo-siRF-prioritized SNVs are gray. Dashed line indicates GWAS p = 5E-8. Circle 9 shows the 21 protein-coding genes mapped by FUMA. b, Pie charts showing ANNOVAR enrichment performed for each of the 6 lo-siRF loci (circle 2 in Fig. 3a and Fig. 2f). The arc length of each slice indicates the proportion of SNVs with a specific functional annotation. The radius of each slice indicates log₂(E + 1), where E is the enrichment score computed as (proportion of SNVs with an annotation for a given locus)/(proportion of SNVs with an annotation relative to all available SNVs in the FUMA reference panel). The dashed circle indicates E = 1 (no enrichment). Asterisks indicate two-sided p-values of Fisher’s exact tests for the enrichment of each annotation. Details can be found in Extended Data 4 and 5.

Considering the correlations between LV hypertrophy and hypertension⁴³, we evaluated whether these identified variants affect LV mass through regulating blood pressure. Specifically, we repeated the lo-siRF analysis using only the subset of UKBB individuals without hypertension (details in Methods). All previously highlighted loci and interactions maintained priority in this non-hypertensive subset, except for the MIR588;RSPO3 locus (Fig. 2f) which was not stably prioritized across all three binarization thresholding schemes. Additionally, none of the lo-siRF-prioritized variants showed a strong marginal association with hypertension, failing to meet the genome-wide (p < 5E-8) and even the suggestive (p < 1E-5) significance level. However, the MIR588;RSPO3 locus with lead SNV rs2022479 gave the smallest p-value of 5E-5, which may suggest a possible pleiotropic effect of MIR588;RSPO3 on both LV hypertrophy and blood pressure. In brief, while we cannot completely rule out pleiotropy, the highly stable prioritization of all three epistatic interactions in both analyses with and without hypertensive individuals suggest that the identified epistases on LV mass is not solely driven by blood pressure (additional discussion in Supplementary Note 1).

Loci associated with left ventricular mass exhibit regulatory enrichment

We performed functional mapping and annotation (FUMA)⁴⁴ for the 283 lo-siRF-prioritized SNVs (Fig. 1, purple and Fig. 3). For linkage disequilibrium (LD), we used a default threshold of r² = 0.6 and chose the UKBB release 2b reference panel created for British and European subjects to match the population group used for lo-siRF prioritization. FUMA identified 572 additional candidate SNVs (Extended Data 4) in strong LD (r² > 0.6) with any of the 283 lo-siRF-prioritized SNVs, including 492 SNVs from the input GWAS associations (points in Fig. 3a, circle 8) and 80 non-GWAS-tagged SNVs extracted from the selected reference panel (heatmap tracks in Fig. 3a, circle 8). We then assigned these 572 FUMA-extracted candidate SNVs to a lo-siRF-prioritized locus (Fig. 2f) based on the corresponding lo-siRF-prioritized SNV (out of the 283 SNVs), which has the maximum r² value with the candidate SNV.

The two loci contributing to the top-ranked epistatic interaction by lo-siRF, the CCDC141 and IGF1R loci (Fig. 2f), both showed a significant enrichment of intronic variants relative to the background reference panel (Fig. 3b, Extended Data 5). Over 88% of the SNVs in or in LD with these two loci were mapped to actively transcribed chromatin states (TxWk) or enhancer states (Enh) in left ventricles based on the ChromHMM Core 15-state model⁴⁵ (Fig. 3a, circle 7). More than 47% and 76% of the identified SNVs in or in LD with the CCDC141 and IGF1R loci, respectively, showed the highest RegulomeDB^44,46 categorical score (ranked within category 1 from the 7 main categories). The Combined Annotation-Dependent Depletion (CADD) score⁴⁷ was used to judge the deleteriousness of prioritized variants (Extended Data 4). As expected, GTEx⁴⁸ data revealed that 82% of SNVs in or in LD with the IGF1R locus are expression quantitative trait loci (eQTLs) for the gene IGF1R. In contrast, of the SNVs in or in LD with the CCDC141 locus, only 14% are eQTLs for gene CCDC141 and 22% are splicing quantitative trait loci (sQTLs) for gene FKBP7. Furthermore, Hi-C data indicated that all SNVs identified in or in LD with the IGF1R locus are in 3D chromatin interaction with gene SYNM while more than 54% SNVs identified in or in LD with the CCDC141 locus are in 3D chromatin interaction with gene TTN. These known 3D chromatin interactions could suggest a possibility of higher-order interactions between more than two genes.

The CCDC141 and TTN loci exhibit genomic proximity (Fig. 3a). Their interaction, however, does not appear to stem from this proximity. Indeed, the CCDC141 and TTN genes have been individually associated with LV mass^49,50. Due to this proximity, previous studies^51,52 have assumed CCDC141 as a secondary gene that affects LV mass through the TTN gene expression. However, we found low LD (r² < 0.6) between any two of the 283 lo-siRF-prioritized SNVs, suggesting that the identified CCDC141-TTN interaction is unlikely driven by non-random LD associations between SNVs in these two loci. In addition, we compared all the epistasis-contributing SNVs that were aggregated to the TTN locus, including both lo-siRF-prioritized SNVs and their LD-linked variants, with the complementary set of TTN-annotated SNVs in lo-siRF. We found that the TTN locus showed a significant depletion of SNVs located close to (<10 kb) the gene CCDC141 (p = 2.38E-9, two-sided Fisher exact test). Similarly, the CCDC141 locus showed a substantially decreased enrichment of SNVs that are close to gene TTN (p = 0.02, two-sided Fisher exact test). These results suggest that although the CCDC141 and TTN loci are located close to each other in the genome, the prioritized epistatic SNVs are located farther apart relative to randomly selected SNVs from the two loci.

In contrast to the CCDC141 and IGF1R loci, the TTN locus showed a significant enrichment of exonic variants and intronic variants that are transcribed into non-coding RNA (ncRNA_intronic, Fig. 3b). Of those exonic variants, 62% are nonsynonymous. This differential enrichment of exonic variants for the TTN locus may suggest a potential epistatic contribution to the structural alterations in the titin protein. Over 90% of SNVs in or in LD with the TTN locus were mapped to actively transcribed states (Tx, TxWk) in left ventricles (Fig. 3a, circle 7). Interestingly, these SNVs were associated with a quiescent chromatin state (Quies) in the right atrium, indicating that the epistatic effects of the TTN locus may be specific to ventricular tissues. Nearly half of SNVs in or in LD with the TTN locus are eQTLs for the gene FKBP7. In addition, 83% of these SNVs are sQTLs for gene FKBP7 or TTN, suggesting a regulatory effect of the TTN locus on the expression and splicing of gene FKBP7. Moreover, the TTN locus was suggested to impact genes PDE11A, RBM45, PRKRA, and DFNB59 through 3D chromatin interactions.

The hypostatic locus LOC157273;TNKS showed a significant enrichment of variants within non-coding RNA regions of exons and introns (Fig. 3b). Over 95% of identified SNVs in or in LD with this locus were mapped to inactive chromatin states (ReprPCWk, Quies) in left ventricles (Fig. 3a, circle 7). This suggests that in the absence of an epistatic partner, the LOC157273;TNKS locus is epigenetically quiescent or repressed by polycomb group proteins. In addition, of all the SNVs in or in LD with this locus, 66% are eQTLs for MFHAS1 or CLDN23 and 22% are in 3D chromatin interaction with gene TNKS.

Functional annotations for the other two lo-siRF-prioritized loci that were marginally associated with LV mass can be found in Extended Data 4 and 5.

Epistatic loci functionally map to twenty-one protein-coding genes

Three strategies, positional, eQTL, and chromatin interaction, mapped the six LV hypertrophy risk loci to 21 protein-coding genes (Fig. 4a). Genes prioritized by eQTL and chromatin interaction mapping are not necessarily located in the corresponding risk locus, but they are linked to SNVs within or in LD with the locus (Fig. 3a). Among the 21 genes, CCDC141 and IGF1R were prioritized by all the three mapping strategies (Fig. 4a), suggesting that these two genes are very likely involved in determining LV mass. Interestingly, none of the SNVs mapped to IGF1R were statistically significant in our GWAS studies using BOLT-LMM and PLINK (Extended Data Fig. 2 and Extended Data 1). Set-based association tests using SKAT-O⁵³ and MAGMA⁵⁴ also did not identify the IGF1R locus (details in Methods and Supplementary Note 1). This reveals the potential of lo-siRF to identify risk loci that may be overlooked by GWAS. Based on the expression data from GTEx V8, TTN, TNNT3, and SYNM are up-regulated while CLDN23 and MFHAS1 are down-regulated in both heart and muscle tissues (Fig. 4b). In contrast, CCDC141 is up-regulated specifically in heart tissues whereas RSPO3 is down-regulated in heart but up-regulated in muscle tissues (Fig. 4b).

• Download figure
• Open in new tab

Fig. 4: Genes mapped from epistatic loci show strong correlations in multiple functional co-association networks

a, UpSet plot showing the number of lo-siRF-prioritized SNVs (dark blue) and their LD-linked (r² > 0.6) SNVs (light blue, circle 8 in Fig. 3a) that are functionally mapped to each of the 21 protein-coding genes by positional, eQTL, and/or chromatin interaction (CI) mapping using FUMA⁴⁴. CCDC141 and IGF1R (highlighted in red) are prioritized by all the three types of SNV-to-gene mapping. Details can be found in Extended Data 4. b, Heatmap of averaged expression (from GTEx) per tissue type per gene (50% winsorization, log₂(TPM + 1)) for these functionally mapped genes. c, Co-association network built from an enrichment analysis integrating multiple annotated gene set libraries for gene ontology (GO) and pathway terms from Enrichr⁵⁵. The co-association network connects top enriched GO and pathway terms with genes (green nodes in the network) functionally linked from lo-siRF-prioritized epistatic and hypostatic loci (Fig. 2f). Strengths (indicated by the edge width in the network) of the co-association between enriched terms and genes were measured and ranked by the empirical p-value from an exhaustive permutation of the co-association score for all possible gene-gene combinations in the network (Details in Methods and Extended Data 6). d, A comparison between the top 10 transcription factors (TFs) enriched from genes prioritized (top) and deprioritized (bottom) by lo-siRF. The lo-siRF-prioritized genes are genes functionally linked from lo-siRF-prioritized SNVs (panel a). The lo-siRF-deprioritized genes are genes functionally linked from SNVs that failed to pass the lo-siRF prioritization thresholds. For each of the two gene groups, enrichment results against nine TF-annotated gene set libraries from ChEA3⁵⁶ and Enrichr⁵⁵ were integrated and ranked by the number of significantly (FET p < 0.05) overlapped libraries (numbers in the nLibraries column) and the mean scaled rank across all libraries containing that TF (colored boxes in the nLibraries column). The balloon plot shows the lowest FET p-values for each TF (horizontal axis) and the proportion of overlapped genes (balloon size) between the input gene set and the corresponding TF-annotated gene set. e, Heatmap showing the TF co-association strength of gene-gene combinations among lo-siRF-prioritized genes relative to randomly selected gene pairs in the co-association network. More details are available in Extended Data Fig. 5 and Extended Data 7.

Ten of twenty-one genes mapped from epistatic loci show strong correlations in network analysis

We performed gene ontology (GO) and pathway enrichment analysis on the 21 genes mapped from lo-siRF loci. We adopted previously established approaches^55–57 and integrated enrichment results across libraries from multiple sources to establish a GO and pathway co-association network (Fig. 4c). To evaluate the correlation strength between any two genes in the network, we calculated a co-association score for every possible gene-gene combination (n = 72,771) from both genes prioritized and deprioritized by lo-siRF. Lo-siRF-prioritized genes are the 21 genes functionally mapped from the 283 lo-siRF-prioritized SNVs and their LD-linked SNVs (Fig. 4a). Lo-siRF-deprioritized genes are those functionally mapped from the SNVs that failed to pass the lo-siRF prioritization threshold. Compared to random gene pairs in the network, 10 genes that were functionally mapped from the lo-siRF-prioritized epistatic and hypostatic loci showed significant co-associations with multiple GO/pathways (Fig. 4c, Extended Data 6). Consistent with our hypothesized epistasis (Fig. 2f), gene CCDC141 showed a significant co-association to SYNM (functionally linked to the IGF1R locus) and PDE11A and PLEKHA3 (both functionally linked to the TTN locus) through the GO term of hyperactivity (excessive movement), which has been linked to increased risk of cardiac disease⁵⁸. Beyond that, TTN, IGF1R, and SYNM are co-associated with kinase activity and cardiac structure related GO terms, indicating that these genes may jointly affect cardiac structure by regulating the process of kinase activity.

Genes mapped from epistatic loci are co-associated with myogenic regulatory factors

We next performed an integrative enrichment analysis to assess transcriptional regulation of genes prioritized and deprioritized by lo-siRF. Due to assay-specific limitations and biases, we integrated the enrichment results across nine distinct gene set libraries^55,56 (Fig. 4d, Extended Data 7). We found that the lo-siRF-prioritized epistatic genes shared important myogenic regulatory factors, such as MYOD1, MYF6, and MYOG (Fig. 4d, top). These myogenic regulatory factors coordinate to regulate muscle development and differentiation. In contrast, Transcription factors enriched from lo-siRF-deprioritized genes display a less coordinated regulatory pattern (Fig. 4d, bottom). These analyses enriched transcription factors based on their associations to given sets of individual genes rather than co-association to gene pairs^55,56. To further evaluate the correlation strength between any two genes that share transcription factors, we calculated a transcription factor co-association score for all the 72,771 possible gene-gene combinations (see Methods). Compared with random gene pairs,16 gene-gene combinations from the lo-siRF-prioritized genes displayed a significant co-association (empirical p < 0.05, Fig. 4e). These co-associations were found in gene-gene combinations from both intra- and inter-lo-siRF-prioritized loci (Fig. 4e). In particular, pairwise combinations among TTN, TNNT3, CCDC141, and SYNM share a common splicing regulator, RBM20 (Extended Data Fig. 5). RBM20 has been reported to regulate the alternative splicing of genes important for cardiac sarcomere organization⁵⁹. This suggests that the splicing patterns of these four genes are likely to be co-regulated by RBM20, which is consistent with the exhibited enrichment of sQTLs by the CCDC141, TTN and LSP1 lo-siRF loci (Extended Data 4).

Genes mapped from epistatic loci exhibit strong co-expression and connectivity change in human heart failure transcriptomics

We proceeded to the fourth stage for experimental confirmation (Fig. 1, blue) and evaluated how the identified epistases contribute to the progression of heart failure (Fig. 5). We employed a series of weighted gene co-expression networks derived from human cardiac transcriptomic data from 177 failing hearts isolated at the time of heart transplant and 136 non-failing hearts harvested from cardiac transplant donors whose organs were not able to be placed⁶⁰ (Fig. 5a). We compared the molecular connectivity of genes identified as statistical epistatic interactors. We defined connectivity as the edge weights between two genes normalized to the distribution of all network edge weights, and compared this to the connectivity of all other available gene-gene combinations in the network. This revealed strong co-expression correlations between CCDC141 and genes functionally linked to the IGF1R locus (SYNM and LYSMD4) and TTN locus (TTN and FKBP7) in the healthy control network (Fig. 5b). In contrast, most of these gene pairs (except for CCDC141-TTN) no longer exhibit a strong connectivity in the heart failure network (Fig. 5c). All of these connectivities showed a significant decrease (indicated by the negative connectivity difference score and p < 0.05 in Fig. 5d) in the differential network, suggesting a declined co-expression correlation between these gene pairs relative to random gene pairs during the progression of failing hearts. This difference is potentially related to the rewired gene modular assignments between the control and heart failure networks⁶⁰ (Fig. 5e and Extended Data Fig. 6). For instance, CCDC141, SYNM, TTN, and TNNT3 are co-associated with the electron transport chain/metabolism module in the control network. In the failing hearts, SYNM and TTN rewire to the muscle contraction/cardiac remodeling module, whereas CCDC141 and TNNT3 remain associated with the metabolism module (Fig. 5e). In addition, other genes functionally linked to IGF1R and TTN lo-siRF loci are co-associated with the membrane transport or unfolded protein response module in healthy hearts and rewire to the muscle contraction/cardiac remodeling or cell surface/immune/metabolism module in failing hearts.

• Download figure
• Open in new tab

Fig. 5: Network analysis using transcriptomic data from 313 human hearts indicates strong correlations between statistical epistasis contributors

a, Control (blue) and heart failure (red) gene co-expression networks were established from a weighted gene co-expression network analysis (WGCNA) on transcriptomic data obtained from 313 non-failing and failing human heart tissues⁶⁰. b-c, The connectivity between lo-siRF-prioritized genes in this study was compared against the full connectivity distributions for all possible gene-gene combinations in the control (b) and heart failure (c) networks. CCDC141 showed a significant connectivity to SYNM and LYSMD4 (IGF1R lo-siRF locus) and TTN and FKBP7 (TTN lo-siRF locus) in the control network. d, Comparing the difference between the control and heart failure networks indicate a significant decrease in the in the connectivities of these gene pairs during the progression of failing hearts. e, A sanky plot demonstrating the rewired gene modular assignments for the lo-siRF loci-associated genes (middle column) in the control vs. heart failure networks. Names of the control (left column) and heart failure (right column) network modules were derived from KEGG and Reactome associations of genes within each module.

Perturbation confirms epistatic relationships in cardiomyocyte hypertrophy

We interrogated epistatic associations in a genetic model of cardiac hypertrophy (Fig. 1, blue): induced pluripotent stem cell cardiomyocytes derived from patients with and without hypertrophic cardiomyopathy caused by the cardiac myosin heavy chain (MYH7) p.R403Q variant³³ (Fig. 6a). Cardiac myosin heavy chain 7 is a key component of the cardiac sarcomere, and the most common cause of hypertrophic cardiomyopathy³³. The patient presented with typical symptoms, and echocardiography revealed severe LV hypertrophy and a small LV cavity³³. At the cellular level, cardiomyocytes exhibit an elevated mean cell size and non-Gaussian size distribution with a long tail relative to the unaffected control (Fig. 6d).

• Download figure
• Open in new tab

Fig. 6: CCDC141 non-additively interacts with TTN and IGF1R to modify cardiomyocyte morphology

a, Human induced pluripotent stem cell (iPSC)-derived cardiomyocytes with and without hypertrophic cardiomyopathy (carrying an MYH7-R403Q mutation) were transfected with scramble siRNA or siRNAs specifically targeting single (CCDC141, IGF1R, and TTN) or combined (CCDC141-IGF1R and CCDC141-TTN) genetic loci prioritized by lo-siRF. b, Gene-silenced cardiomyocytes were bifurcated into two focused streams of large and small cells using a spiral microfluidic device (cell focusing mechanism illustrated in Extended Data Fig. 7) to allow high-resolution single-cell imaging. c, Workflow of the image analysis process. Time-lapse image sequences of single cells passing through the top and bottom microchannel outlets (panel b) were fed into a customized MATLAB-based program that extracts cell size/shape features via a sequential process of bright field background correction, cell boundary detection, cell tracking and stuck cell removal, cell feature extraction, data quality control and postprocessing, and morphological feature analysis. Extracted single cell features for each gene-silencing condition were compared with their scrambled control values to validate the potential role of epistasis in the genetic regulation of cardiomyocyte hypertrophy. d, Violin plots of cell diameters of unaffected (blue) and MYH7-R403Q variant (red) cardiomyocytes. Solid and dashed lines in box plots represent median and mean values, respectively. Asterisks indicate significant difference (***p < 1E-36, Wilcoxon signed rank test). e, Gene-silencing efficiency in unaffected (blue, n = 5 to 9) and MYH7-R403Q variant (red, n = 3) cells based on RT-qPCR analysis (details in Methods). Error bars indicate standard deviations. f, Percent change in median cell diameter (relative size difference) of gene-silenced cardiomyocytes relative to scrambled control values due to monogenic and digenic gene knockdown effects. Relative size differences were averaged across data from two to four independent batches of cells. Error bars indicate standard deviations computed on 1000 bootstrap samples of these batches with the following sample size: n = 13147 (TTN), 19460 (IGF1R), 45304 (CCDC141), 19979 (CCDC141-TTN), and 26135 (CCDC141-IGF1R) for unaffected cells and n = 22134 (TTN), 33801 (IGF1R), 21158 (CCDC141), 39515 (CCDC141-TTN), and 52049 (CCDC141-IGF1R) for MYH7-R403Q variant cells. Asterisks indicate significant difference between gene-silencing and scrambled control conditions based on the maximum p-values of Wilcoxon signed rank test across all batches of cells (*p < 0.05, **p < 0.001, and ***p < 1E-4). g, Violin plots highlighting the magnitudes and directions of non-additive interaction effects () for unaffected (blue) and MYH7-R403Q variant (red) cells and each gene pair compared to marginal effects () and ), gray), estimated via a quantile regression model across 10,000 bootstrap samples (details in Methods). h, CCDC141 non-additively interacts with IGF1R (left) and TTN (right) to modify boundary and texture features of unaffected (blue) and MYH7-R403Q variant (red) cells. Cell boundary waveness and texture irregularity were measured by the roundness error (i, top) and normalized peak number (i, bottom), respectively. j, Representative single-cell images overlapped with detected cell boundaries (red lines) show that a higher roundness error indicates increased irregularity of the cell boundary. k, Representative single-cell images with detected peaks (blue plus signs) of the brightfield intensity distribution enclosed within the cell boundaries (red lines) indicate a varying level of cell textural irregularity. Scale bars: 10 μm. Detailed statistical information of cell morphology measurements and non-additivity analysis for the studied gene pairs can be found in Extended Data 8.

To determine if CCDC141 can act both independently and in epistatic interactions with other genes to attenuate the pathologic cellular hypertrophy caused by MYH7-R403Q, we silenced genes CCDC141, IGF1R, TTN, and gene pairs CCDC141-IGF1R and CCDC141-TTN using siRNAs in both diseased and healthy cardiomyocytes and compared them with cells transfected with scramble siRNAs (control) (Fig. 6a and 6e). Phenotypic consequences of these perturbations on cellular morphology were then evaluated in high-throughput using a spiral inertial microfluidic device (Fig. 6b) in combination with automated single-cell image analysis (Fig. 6c). The microfluidic device adopted the Dean flow focusing principle³¹ (details in Extended Data Fig. 7 and Methods) to mitigate the non-uniform cell focusing⁶¹, thereby enhancing the imaging resolution⁶² affected by the large variations in cardiomyocyte diameter (Fig. 6d).

We first assessed the knockdown effects of the CCDC141-IGF1R interaction on cardiomyocyte size (Fig. 6f). Bootstrapped hypothesis tests were performed, for which the p-values are capped below by p < 1E-4 (Extended Data 8). Silencing IGF1R alone reduces the median cell size by 5.3% ± 0.4% (p < 1E-4) in diseased cells compared to scrambled control and 6.6% ± 0.5% (p < 1E-4) in healthy cells. Silencing CCDC141 alone also decreases median cell size by 3.2% ± 0.5% (p < 1E-4) in diseased cells, but had no impact on healthy cells. Digenic silencing of CCDC141 and IGF1R reveals a synergistic effect on attenuating pathologic cell hypertrophy in diseased cells, resulting in an 8.5% ± 0.3% (p < 1E-4) decrease in the median cell size. This is consistent in healthy cells, where silencing CCDC141 alone fails to affect cell size, but digenic silencing of CCDC141 and IGF1R decreases the median cell size by 9.3% ± 0.5% (p < 1E-4). Moreover, according to our estimated quantile regression analysis (details in Methods), this interaction effect appears to be non-additive for both healthy and diseased cells ( < 0, Fig. 6g; p < 1E-4 for non-additivity, Extended Data 8), consistent with an epistatic mechanism. These findings serve to confirm the strongest epistatic association identified by lo-siRF (Fig. 2f).

We found a comparable non-additive effect for the CCDC141-TTN interaction. Digenic silencing of CCDC141-TTN leads to a pronounced reduction in median cell size (by 5.8% ± 0.6% for healthy cells and 3.3% ± 0.4% for diseased cells, p < 1E-4) relative to monogenic silencing (Fig. 6f). This interaction appears to be non-additive for both healthy and diseased cells (p values in Extended Data 8) yet demonstrating opposite epistatic directions in these two cell states (Fig. 6g). Additionally, CCDC141 and TTN show distinctive independent roles in repressing cardiomyocyte hypertrophy. In healthy cells, monogenic silencing of TTN leads to a larger cell size reduction compared to the case of silencing CCDC141. In contrast, diseased cells display a larger size reduction in response to monogenic silencing of CCDC141.

Furthermore, both CCDC141-IGF1R and CCDC141-TTN interactions show a stronger effect on rescuing larger cardiomyocytes over smaller ones in both cell lines (Extended Data Fig. 8 and 9). In contrast, monogenic silencing does not exhibit such a non-uniform effect on reshaping the cell size distribution, which reinforces the hypothesized non-additivity of these two epistatic interactions (details in Extended Data 8 and Supplementary Note 2).

Recent studies have shown that cellular morphological features, such as cell boundary and textural irregularities, are informative readouts of cytoskeletal structure, which is highly associated with disease state in hypertrophic cardiomyopathy^32,63. We analyzed relative changes in cell shape and texture (Fig. 6h) by measuring the counts of peak intensities normalized to the total number of pixels enclosed by the cell boundary (Fig. 6i). Cells with a high normalized peak number display a ruffled texture, which manifests in unevenly distributed 2D intensities (Fig. 6k). Our analysis shows that silencing both CCDC141 and IGF1R (circles in Fig. 6h, left) yields a larger increase in intensity peak number than silencing IGF1R alone (triangles in Fig. 6h, left) for both cell lines, exhibiting a synergistic epistasis between CCDC141 and IGF1R (p < 1E-4 for non-additivity). We also analyzed cell roundness error, a measure of how far radii measured on the cell outline deviate from a perfect circle (Fig. 6i). This parameter increases with an increasing cell boundary waviness or elongation (Fig. 6j). We show that the silencing of CCDC141 and IGF1R synergistically interact to increase roundness error of diseased cardiomyocytes (p < 1E-4 for non-additivity, Fig. 6h, left). In addition, CCDC141 and TTN display antagonistic epistasis and synergistic epistasis in their impact on roundness error for healthy and diseased cells (p < 1E-4 for non-additivity, Fig. 6h, right), respectively.

Study participants

The use of human subjects (IRB −4237) and human-derived induced pluripotent stem cells (SCRO −568) in this study has been approved by the Stanford Research Compliance Office. The UK Biobank received ethical approval from the North West -Haydock Research Ethics Committee (21/NW/0157).

The UK Biobank (UKBB) is a biomedical database with detailed phenotypic and genetic data from over half a million UK individuals between ages 40 and 69 years at recruitment⁷⁷. In this study, we restricted our analysis to the largest ancestry subset (i.e., the White British population) of 29,661 unrelated individuals who have both genetic and cardiac magnetic resonance imaging (MRI) data from the UKBB (Supplementary Table 1). More specifically, we considered only those individuals from the UKBB cohort who self-reported as White British and have similar genotypic backgrounds based on principal components analysis as described in prior work⁷⁷. We also identified related individuals (i.e., third-degree relatives or closer) via genotyping and omitted all but one individual from each related group in the analysis. Details regarding this cohort refinement have been described and implemented previously^77,78. This refinement resulted in a cohort of 337,535 unrelated White British individuals from the UKBB, of which 29,661 have both genetic and cardiac MRI data. We randomly split this data into training, validation, and test sets of size 15,000, 5,000, and 9,661 individuals, respectively.

Genotyping and quality control

For the study cohort of 29,661 individuals described above, we leveraged genotype data from approximately 15 million imputed autosomal SNVs. These have been imputed from 805,426 directly assayed SNVs (obtained by the UKBB from one of two similar Affymetrix arrays) using the Haplotype Reference Consortium and UK10K reference panels⁷⁷. Imputed variants were subject to several quality-control filters, including outlier-based filtration on effects due to batch, plate, sex, array, and discordance across control replicates. Further, we excluded variants due to extreme heterozygosity, missingness, minor allele frequency (< 10^-4), Hardy-Weinburg equilibrium (< 10^-10), and poor imputation quality (< 0.9). Further details can be found in previous studies^77,78.

Quantification of left ventricular hypertrophy

We retrieved cardiac MRI images from 44,503 UKBB participants, taken during their most recent imaging visit, and closely followed the method previously described by Bai et al.²⁸. A fully convolutional network²⁸ was previously trained using a dataset of 4,875 subjects with 93,500 pixelwise segmentations of UKBB short-axis cardiac MRI multi-slice images generated manually with quality control checks for inter-operator consistency⁷⁹. The cardiac MRI image resolution was 1.8 x 1.8 mm², with a slice thickness of 8.0 mm and gap of 2.0 mm, typically consisting of 10 slices. Each slice was converted to an image and cropped to a 192 x 192 square, and measurements were 0-1 normalized. The network architecture employed multiple convolutional layers to learn image features across five resolution scales. Each scale involved two or three convolutions with kernel size 3 x 3 and stride 1 or 2 (2 appearing every 2 or 3 layers), followed by batch normalization and ReLU transformation. Feature maps from the five scales were upsampled back to the original resolution, combined into a multi-scale feature map, and processed through three additional convolutional layers with kernel size 1 x 1, followed by a softmax function to predict the segmentation label for each pixel. For an exact description of the model architecture, we refer to the original publication on the model²⁸. Notably, each of the pixelwise annotations used for training and evaluation was hand-segmented by a human expert and validated for quality. Furthermore, the model was validated in the UKBB and demonstrated strong concordance with the human-generated gold standard²⁸, ensuring that model predictions in the same dataset are of high quality. To our knowledge, this is the only published model trained in the UKBB on gold standard labels. We thus applied this trained deep learning model to our entire dataset of 44,503 cardiac MRIs. This resulted in segmentations of the LV cavity and myocardium from each short axis frame, which allowed for both an area calculation of each segment as well as the application of quality control checks²⁸ based on consistency within and between slices and time steps. There were 44,219 segmentations that passed the quality control. Using the calculated areas, we computed the volume of the LV myocardium through simple integration over slices. This volume was then converted to a left ventricular mass (LVM) using a standard density estimate of 1.05 g/mL³⁵. LVMi was computed by dividing LVM by an estimate of body surface area based on height and body weight calculated using the Du Bois formula³⁶. From the 44,219 segmentations, we restricted the analysis to LVMi measurements for 29,661 unrelated White British individuals using the measurements from their most recent imaging visit if multiple imaging visits were recorded.

Lo-siRF step 1: Dimension reduction of variants via genome-wide association studies

As the first step in the lo-siRF pipeline, we performed a genome-wide association study (GWAS) on the training data for the rank-based inverse normal-transformed LVMi using two algorithms, PLINK³⁷ and BOLT-LMM³⁸, in order to filter the number of features from over 15 million SNVs to a more computationally-feasible size (Fig. 2b). This step is akin to typical screening phases in fine-mapping⁸⁰ and other tree-based epistasis detection methods^72,81. Since BOLT-LMM and PLINK rely on different statistical models, we chose to employ both implementations to mitigate the dependence of downstream conclusions on this arbitrary choice. Specifically, for the first GWAS run, we fitted a linear regression model, implemented via ‘glm’ in PLINK⁸². For the second GWAS run, we used BOLT-LMM³⁸, a fast Bayesian-based linear mixed model method. Each GWAS was adjusted for the first five principal components of ancestry, sex, age, height, and body weight. We then ranked the SNVs by significance (i.e., the GWAS p-value) for each GWAS run separately and took the union of the top 1000 SNVs (without clumping) from each of the two GWAS runs. This resulted in a set of 1405 GWAS-filtered SNVs that were used in the remainder of the lo-siRF pipeline. Here, we chose to use the top 1000 SNVs per GWAS method (without clumping) as it yielded the highest validation prediction accuracy compared to choosing other possible thresholds (500 and 2000 SNVs per GWAS with and without clumping). Though the GWAS is not the focus of this work, we provide a summary of the PLINK and BOLT-LMM GWAS results for completeness and for comparison in Extended Data 1. We also provide a list of the 1405 GWAS-filtered SNVs in Extended Data 2. We note that these 1405 GWAS-filtered SNVs strictly contain the SNVs that passed the genome-wide significance threshold (p = 5E-8).

Lo-siRF step 2: Binarization of the left ventricular mass phenotype

Next, we partitioned the raw (continuous) LVMi phenotype into a low, middle, and high LVMi group before fitting signed iterative random forest to classify individuals with low versus high LVMi (Fig. 2c). That is, for a given threshold x, we binned individuals within the top and bottom x% of LVMi values into two classes with the high and low LVMi values, respectively, while omitting the individuals in the middle quantile range. Due to the sex-specific biological variation of LVMi (Supplementary Note 1), we performed this partitioning for males and females separately. For males, low and high LVMi was considered under 43.8-46.0 g/m² and above 55.4-58.5 g/m², respectively, depending on our choice of binarization threshold (Supplementary Table 2). For females, low and high LVMi was defined as under 35.1-36.8 g/m² and above 43.8-46.1 g/m², respectively, depending on our choice of binarization threshold. We performed this binarization step in order to simplify the original low-signal regression problem into a relatively easier binary classification task: to distinguish between individuals with very high LVMi from those with very low LVMi. This binarization approach was motivated by the observation that the validation R² values from the original regression problem of predicting each individual’s raw (continuous) LVMi were smaller than 0 (Supplementary Table 3 and Supplementary Note 1), raising the question of whether the regression models were capturing anything relevant to reality. At a minimum however, the PCS framework for veridical data science²⁶ advocates the importance of ensuring that the model fits the data well, as measured by prediction accuracy, before trusting any extracted interpretations from that model. We will see in the next section that the binarization procedure not only strengthened the prediction signal but also helped us more readily interpret and assess the performance of prediction methods with respect to the prediction screening step of the PCS framework²⁶. We importantly note that the practical use of this approach depends heavily on whether understanding the differences between the high and low categories is relevant to the scientific goals. Here, we believe that the connection between cardiac hypertrophy and those with high LVMi helps to justify the binarization approach and that studying how individuals with high LVMi differ from those on the other end of the spectrum may yield relevant scientific insights. Since the specific threshold choice is arbitrary, we ran the remainder of the lo-siRF pipeline using three different binarization thresholds (15%, 20%, 25%) to balance the improvement in prediction signal and amount of data lost. In the end, we aggregated the results that were stable across all three binarization thresholds, described in the Method section Lo-siRF step 4.4: Ranking genetic loci and interactions between loci.

Lo-siRF step 3: Prediction

Lo-siRF step 4: Prioritization

To lastly interpret the siRF fit for LVMi, we developed a novel stability-driven importance score to prioritize genetic loci and more interestingly, interactions between loci for follow-up experimental validation (Fig. 2e). The assessment of importance at the level of genetic loci, instead of individual variants, is necessary since variant-level importances here are incredibly unstable (detailed in Supplementary Note 1). This is due to the high correlation between SNVs in LD and the weak phenotypic signal. Consequently, our new importance score aims to aggregate weak, unstable variant-level importances into stronger, more stable locus-level importances via three steps: (1) assigning each variant to a genetic locus, (2) evaluating the local (or per-individual) importance of each genetic loci or interaction between loci in the siRF fit via a stability-driven measure, and (3) conducting a permutation test to summarize the importance of the genetic locus or interaction between loci across all individuals. We provide details for each step next.

Lo-siRF: PCS documentation and additional stability analyses

We acknowledge that many human judgment calls were inevitably made throughout our veridical machine learning pipeline and that alternative choices could have been made (e.g., different dimension reduction techniques, binarization procedures, and prediction models). In an effort to facilitate transparency of these human judgment calls, we provide extensive documentation, discussion, and justification in Supplementary Note 1. In particular, Supplementary Note 1 includes a discussion of our reasoning and motivation behind the choice of phenotypic data as well as choices in the dimension reduction, binarization, prediction, and prioritization steps. We also performed additional stability analyses in accordance with the PCS framework²⁶, to ensure that our findings are stable and robust to these human judgment calls (e.g., the choice of GWAS method and binarization threshold) and to bolster the reproducibility of our findings. Supplementary Note 1 is an HTML document, which can be downloaded and displayed in a browser or found at https://yu-group.github.io/epistasis-cardiac-hypertrophy/.

Non-hypertensive cohort analysis

We defined hypertensive individuals as anyone with self-reported hypertension, high blood pressure as diagnosed by a doctor, or any ICD10 billing code diagnosis in I10-I16. Out of the 29,661 UKBB participants in the original lo-siRF analysis, 7,371 individuals had hypertension, leaving 22,290 individuals for the non-hypertensive analysis. Specifically, using the same set of 1405 GWAS-filtered SNVs as in the original lo-siRF analysis, we performed steps 2-4 of the lo-siRF analysis using only the non-hypertension cohort. We also assessed the marginal effect of each of the 1405 GWAS-filtered SNV on hypertension. Here, we fit logistic regression models, regressing hypertension (i.e., a binary indicator of whether or not one has hypertension) on each SNV marginally, while adjusting for the first five principal components of ancestry, sex, age, height, and body weight. A more detailed discussion of the non-hypertension analysis results can be found in Supplementary Note 1.

Implementation of existing epistasis detection methods

Implementation of existing set-based genome-wide association tests

To investigate the importance of the IGF1R locus using existing set-based association methods, we performed SKAT-O⁵³ using the subset of 1405 GWAS-filtered SNVs with minor allele frequency > 0.05 as input and the rank-based inverse normal-transformed LVMi as the response. We also adjusted for sex, age, height, body weight, and the first five principal components of ancestry in the SKAT-O null model. This analysis was carried out using the SKAT (v2.2.5) R package. In addition to SKAT-O, we also ran the gene-based test as computed by MAGMA⁵⁴ (v1.6) using the LVMi PLINK GWAS results as input. This MAGMA analysis was carried out using FUMA with the default settings. Results of these two analyses are detailed in Supplementary Note 1.

Functional interpretation of lo-siRF-prioritized variants

Functional interpretation step 4: Functional gene mapping

In SNP2GENE, we performed three functional gene mapping strategies – positional, eQTL, and 3D chromatin interaction mapping – using the lo-siRF-prioritized SNV set and the reference SNV set described in the Methods section Functional interpretation step 1: Extraction of candidate SNVs and LD structures. For positional mapping^44,46, a default value of 10 kb was used as the maximum distance between SNVs and genes. For eQTL mapping, cis-eQTL information of heart left ventricle, heart atrial appendage, and muscle skeletal tissue types from GTEx v8⁴⁸ was used. Only significant SNV-gene pairs (FDR < 0.05 and p < 1E-3) were used for eQTL mapping. For 3D chromatin interaction mapping, Hi-C data of left ventricle tissue from GSE87112 was chosen with a default threshold of FDR < 1E-6. A default promoter region window was defined as 250 bp upstream and 500 bp downstream of TSS^44,46. Using these three gene mapping strategies, we mapped the lo-siRF-prioritized SNV set to 21 protein-coding genes (Fig. 4a), of which 20 are HGNC-recognizable. Each of the 21 genes was also functionally linked to a specific lo-siRF-prioritized LV hypertrophy risk locus (Fig. 2f), to which the highest proportion of SNVs mapped to the given gene were assigned. A Circos plot (Fig. 3a) showing comprehensive information of the lo-siRF-prioritized epistatic interactions, FUMA-prioritized eQTL SNV-to-gene connections and 3D chromatin interactions, as well as LD structures and prioritized genes was created by TBtools⁹¹. We then submitted these 21 genes to the GENE2FUNC process in FUMA and obtained GTEx gene expression data for 19 (out of the 21) genes across multiple tissue types (Fig. 4b). In addition, we used the same approach and mapped the reference SNV set (mentioned in the Methods section Functional interpretation step 1: Extraction of candidate SNVs and LD structures) to a separate gene set that contains 382 HGNC-approved genes. The lo-siRF-prioritized gene set and the reference gene set were used for gene set enrichment analysis that are described in the following sections.

Gene ontology and pathway enrichment analysis

Genes that co-associate to shared gene ontology (GO) and pathway terms are likely to be functionally related. To assess the differential GO and pathway co-association among the lo-siRF-prioritized genes relative to their counterparts that were deprioritized by lo-siRF, we performed an integrative GO and pathway enrichment analysis followed by an exhaustive permutation of co-association scores between any possible gene-gene combinations found in the aforementioned 382 HGNC-approved genes (see the Methods section Functional interpretation step 4: Functional gene mapping).

In order to improve GO and pathway prioritization, we adopted the concept from Enrichr-KG⁵⁷ and ChEA3⁵⁶ to assess enrichment analysis results across libraries and domains of knowledge as an integrated network of genes and their annotations. We first queried the 382 HGNC-approved genes from the reference gene set against various prior-knowledge gene set libraries in Enrichr⁵⁵ (https://maayanlab.cloud/Enrichr/). We selected five representative libraries from the GO and pathway Enrichr categories as follows: GO biological process^92,93, GO molecular function^92,93, MGI Mammalian Phenotypes⁹⁴, Reactome pathways⁹⁵, and KEGG pathways⁹⁶. Other FUMA-extracted genes that were not approved by HGNC using synonyms or aliases were discarded. This enrichment analysis allowed us to search for a union of enriched GO or pathway terms and their correspondingly annotated gene sets, from which we built a co-association network. According to the method by Enrichr-KG⁵⁷, nodes in the co-association network are either the enriched GO and pathway terms or genes.

To measure the degree of co-association to specific GO and pathway terms for two given interactor genes, we computed a co-association score for each of the 72,771 possible gene-gene combinations (from the 382 queried genes). The co-association score was calculated by R = N_(A⋂B)⧸N_(A⋃B). Here, N_(A⋂B) denotes the number of GO or pathway terms that were significantly enriched for both gene A and gene B in the proposed gene-gene combination, and N_(A⋃B) is the number of GO or pathway terms that were enriched for either gene A or gene B. For cases where N_(A⋃B) = 0, we defined R = 0 to indicate that no GO or pathway terms were found to be co-associated with the respective gene-gene combinations. Of the 382 HGNC-approved genes, 20 genes were mapped to lo-siRF-prioritized loci by FUMA functional gene mapping (one of the 21 lo-siRF-prioritized genes is HGNC-unrecognizable and is discarded). We compared the co-association scores R for gene-gene combinations in the lo-siRF-prioritized gene set relative to the full distribution of R provided by an exhaustive permutation of all possible gene-gene combinations in the set of 382 HGNC-approved genes. The ranking of gene-gene combinations was determined by the two-sided empirical p-values. Fig. 4c displays significant co-associations (empirical p < 0.05) between enriched GO or pathway terms and genes functionally mapped to lo-siRF-prioritized epistatic and hypostatic loci (Fig. 2f). Further details can be found in Extended Data 6.

Transcription factor enrichment analysis

Owing to the limitations and biases of various specific assays, we performed an integrative transcription factor (TF) enrichment analysis against multiple annotated gene set libraries in ChEA3⁵⁶ and Enrichr⁵⁵. To preserve the variety of library types, we assembled 9 gene set libraries (Fig. 4d) from distinct sources as follows:

1. Putative TF target gene sets determined by ChIP-seq experiments from ENCODE⁹⁷;

2. Putative TF target gene sets determined by ChIP-seq experiments from ReMap⁹⁸;

3. Putative TF target gene sets determined by ChIP-seq experiments from individual publications⁵⁶;

4. TF co-expression with other genes based on RNA-seq data from GTEx⁴⁸;

5. TF co-expression with other genes based on RNA-seq data from ARCHS4⁹⁹;

6. Single TF perturbations followed by gene differential expression⁵⁵;

7. Putative target gene sets determined by scanning PWMs from JASPAR¹⁰⁰ and TRANSFAC¹⁰¹ at promoter regions of all human genes;

8. Gene sets predicted by transcriptional regulatory relationships unraveled by sentence-based text-mining (TRRUST)¹⁰²;

9. Top co-occurring genes with TFs in a large number of Enrichr queries⁵⁶.

Of the mentioned 9 gene set libraries, libraries 1, 3, and 5 were assembled by combining gene set libraries downloaded from both ChEA3⁵⁶ and Enrichr⁵⁵. Libraries 2, 4, and 9 were downloaded from ChEA3⁵⁶. Libraries 6, 7, and 8 were downloaded from Enrichr⁵⁵. According to the integration method by ChEA3⁵⁶, for libraries in which multiple gene sets were annotated to the same TF, the unique gene set with the lowest FET p-value was used. As mentioned in previous sections, we used separate FUMA SNP2GENE processes and functionally mapped lo-siRF-prioritized SNVs and all GWAS-filtered SNVs to a lo-siRF-prioritized gene set (20 HGNC-recognizable genes) and a reference gene set (382 HGNC-recognizable genes), respectively. Because the lo-siRF-prioritized gene set is a subset of the reference gene set, we considered the 362 genes complementary to the lo-siRF-prioritized gene set as the lo-siRF-deprioritized gene set. Taking the 20 lo-siRF-prioritized genes and 263 lo-siRF-deprioritized genes as two separate input gene sets, we performed enrichment analysis against the 9 gene set libraries. For each of the 9 libraries, we ranked the significance of overlap between the input gene set and the TF-annotated gene sets in that library by FET p-values. Those TFs with identical FET p-values were ranked by the same integer number. A scaled rank was then assigned to each TF by dividing the corresponding integer rank by the maximum integer rank in its respective library. We then integrated the 9 sets of TF rankings and re-ordered the TFs by two sequential criteria: (1) the number of libraries that display a significant overlap with the input gene set (FET p < 0.05) and (2) the mean scaled rank across all libraries containing that TF. Using this method, we prioritized two distinct sets of TFs for the lo-siRF-prioritized genes (Fig. 4d, top) and lo-siRF-deprioritized genes (Fig. 4d, bottom).

The above analyses aim to enrich TFs based on their associations to given sets of individual genes rather than co-associations to gene pairs^55,56. To evaluate the differential TF co-association between the lo-siRF-prioritized genes relative to the lo-siRF-deprioritized genes, we used the same approach described in the Methods section Gene ontology and pathway enrichment analysis. A TF co-association score was computed for each of the 72,771 possible gene-gene combinations from the 382 genes. Still, the TF co-association score was computed by R = N_(A⋂B)⧸N_(A⋃B), except that N_(A⋂B) and N_(A⋃B) denote the number of enriched TF terms instead of GO or pathway terms. Pairwise interactions between lo-siRF-prioritized genes were extracted and ranked by the empirical p-values (Fig. 4e) from an exhaustive permutation of TF co-association scores for the 72,771 possible gene-gene combinations. Details regarding TF enrichment from both lo-siRF-prioritized and lo-siRF-deprioritized genes and TF co-association strengths can be found in Extended Data 7.

Disease-state-specific gene co-expression network analysis

In order to evaluate the connectivity between genes and their potential roles in the transition from healthy to failing myocardium, we compared gene-gene connectivity and changes in the topological structure between gene co-expression networks for healthy and failing human heart tissues (Fig. 5). To construct gene co-expression networks, cardiac tissue samples from 177 failing hearts and 136 donor, non-failing (control) hearts were collected from operating rooms and remote locations for RNA expression measurements. We performed weighted gene co-expression network analysis (WGCNA) on the covariate-corrected RNA microarray data for the control and heart failure networks separately (Fig. 5a). Detailed steps for generating these co-expression networks, which included calculating the correlation matrix, TOM transformation, and Dynamic Tree Cut module finding, are described in our previous study⁶⁰, and data for these networks is available at https://doi.org/10.5281/zenodo.2600420. To evaluate the degree of connectivity between correlating genes in each of the networks, we compared the edge weights between the lo-siRF-prioritized genes demonstrated in this study relative to the distribution of all possible pairwise combinations of genes (Fig. 5b and 5c). We also evaluated the difference of edge weights (Z-score normalized) between the control and heart failure networks to understand how these gene-gene connectivities change between non-failing and failing hearts (Fig. 5d). The two-tailed empirical p-value represents the proportion of the absolute difference in edge weights of all gene pairs that exceed the absolute difference score for gene pairs of interest. We then compared the structure of modules derived from dendrograms on the WGCNA control and heart failure networks (Extended Data Fig. 6). Modules were labeled according to Reactome enrichment analysis of genes within each module. The full gene module descriptions and Benjamini-Hochberg-adjusted enrichment p-values can be found in the Supplementary Data 5 and 6 in the study by Cordero et al.⁶⁰.

Induced pluripotent stem cell cardiomyocytes differentiation

The studied patient-specific human induced pluripotent stem cells (hiPSCs) were derived from a 45-year-old female proband with a heterozygous MYH7-R403Q mutation. Derivation and maintenance of hiPSC lines were performed following Dainis et al.³³. Briefly, hiPSCs were maintained in MTeSR (StemCell Technologies) and split at a low density (1:12) onto fresh 1:200 matrigel-coated 12 well plates. Following the split, cells were left in MTeSR media supplemented with 1 μM Thiazovivin. The hiPSCs were maintained in MTeSR until cells reached 90% confluency, which began Day 0 of the cardiomyocyte differentiation protocol. Cardiomyocytes were differentiated from hiPSCs using small molecule inhibitors. For Days 0-5, cells were given RPMI 1640 medium + L-glutamine and B27 -insulin. On Days 0 and 1, the media was supplemented with 6 μM of the GSK3β inhibitor, CHIR99021. On Days 2 and 3, the media was supplemented with 5 μM of the Wnt inhibitor, IWR-1. Media was switched to RPMI 1640 medium + L-glutamine and B27 + insulin on Days 6-8. On Days 9-12, cells were maintained in RPMI 1640 medium + L-glutamine -glucose, B27 + insulin, and sodium lactate. On Day 13, cells were detached using Accutase for 7-10 minutes at 37 ℃ and resuspended in neutralizing RPMI 1640 medium + L-glutamine and B27 + insulin. This mixture was centrifuged for 5 minutes at 1000 rpm (103 rcf). The cell pellet was resuspended in 1 μM thiazovivin supplemented RPMI 1640 medium + L-glutamine and B27 + insulin. For the rest of the protocol (Days 14-40), cells were exposed to RPMI 1640 medium + L-glutamine -glucose, B27 + insulin, and sodium lactate. Media changes occurred every other day on Days 14-19 and every three days for Days 20-40. On Day 40, cardiomyocytes reached maturity.

RNA silencing in induced pluripotent stem cell-derived cardiomyocytes

Mature hiPSC-derived cardiomyocytes were transfected with Silencer Select siRNAs (Thermofisher) using TransIT-TKO Transfection reagent (Mirus Bio). Cells were incubated for 48 hours with 75 nM siRNA treatments. Four wells of cells were transfected with each of the six siRNAs: scramble, CCDC141 (ID s49797), IGF1R (ID s223918), TTN (ID s14484), CCDC141 and IGF1R, and CCDC141 and TTN. After 2 days, hiPSC-CMs were collected for RNA extraction.

RT-qPCR analysis for siRNA gene silencing efficiency

Following cell morphology measurement, all cells for each condition were centrifuged for 5 minutes at 1000 rpm (103 rcf). Cell pellets were frozen at −80 ℃ prior to RNA extraction. RNA was extracted using Trizol reagent for RT-qPCR to confirm gene knockdown occurred. Reverse Transcription of RNA was done using High-Capacity cDNA Reverse Transcription Kit (Thermofisher). qPCR of the single stranded cDNA was performed using TaqMan Fast Advanced MM (Thermofisher) with the following annealing temperatures: 95°C 20” and 40 cycles of 95°C 1” and 60°C 20”. qPCR of the silenced genes was performed using TaqMan® Gene Expression Assays, including CCDC141 (Hs00892642_m1), IGF1R (Hs00609566_m1), and TTN (Hs00399225_m1). For gene-silencing efficiency analysis, gene RPLP0 (Hs00420895_gH) was used as a reference gene. Data were analyzed using the delta-delta Ct method.

Cell sample preparation for cell morphology measurement

Following siRNA treatments, cells were detached for microfluidic single-cell imaging using a mixture of 5 parts Accutase and 1 part TrypLE, treated for 6 minutes at 37 ℃. Cells were then added to the neutralizing RPMI 1640 medium + L-glutamine and B27 + insulin. These mixtures were centrifuged for 5 minutes at 1000 rpm (103 rcf). For each gene-silencing condition, the four wells of cells were resuspended in 4 mL of the MEM medium, which is composed of MEM (HBSS balanced) medium, 10% FBS, and 1% Pen Strep (Gibco). Cells were filtered with 100 μm strainers (Corning) before adding into the microfluidic devices.

Microfluidic inertial focusing device

We developed a new spiral inertial microfluidics system on the basis of the study by Guan et al.³¹ to focus randomly suspended cells into separate single streams based on cell size for high-resolution and high-throughput single-cell imaging. The microfluidic device (Extended Data Fig. 7) contains 5 loops of spiral microchannel with a radius increasing from 3.3 mm to 7.05 mm. The microchannel has a cross-section with a slanted ceiling, resulting in 80 μm and 150 μm depths at the inner and outer side of the channel, respectively. The channel width is fixed to 600 μm. The 495 μm wide slanted region of the channel ceiling is composed of ten 7 μm deep stairs. This particular geometry induces strong Dean vortices in the outer half of the channel cross-section, leading to high sensitivity of size separation and cell focusing. The device has two inlets at the spiral center to introduce cell suspensions and sheath flow of fresh medium. At the outlet region, the channel is expanded in width and split into two outlet channels with a width of 845 μm for the top outlet and 690 μm for the bottom outlet. Depths of the two outlet channels are designed to create equal hydraulic resistance. The top and bottom outlet channels are connected to 80 μm and 50 μm deep straight observation channels for high-throughput cell imaging.

Microdevice fabrication

The spiral microchannel was fabricated by CNC micromachining a piece of laser-cut poly (methyl methacrylate) (PMMA) sheet, which was bonded with a PMMA chip machined only with the inlet channels and another blank PMMA chip using a solvent-assisted thermal binding process to form the enclosed channel²⁹. Before bonding, PMMA chips were cleaned with acetone, methanol, isopropanol and deionized water in sequence. Droplets of a solvent mixture (47.5% DMSO, 47.5% water, 5% methanol) were evenly spread over the cleaned chips. The PMMA chips were assembled appropriately and clamped using a customized aluminum fixture, and then heated in a ThermoScientific Lindberg Blue M oven at 96 ℃ for 2 hrs. After bonding, fluid reservoirs (McMaster) were then attached to the chips using a two-part epoxy (McMaster). Microchannels were flushed with 70% ethanol followed by DI water for sterilization.

High-throughput single-cell imaging

Before each experiment, microchannels were flushed with 3 mL of the MEM medium. Prepared cell samples and fresh MEM medium were loaded into 3 mL syringes, which were connected to the corresponding microchannel inlets using Tygon PVC tubing (McMaster). Both cells and the fresh MEM medium were infused into the microchannel using a Pico Plus Elite syringe pump (Harvard Apparatus) at 1.2 mL/min. Microscope image sequences of cells focused to the top and bottom observation channels were captured using a VEO 710S high-speed camera (Phantom) with a sampling rate of 700 fps and a 5 μsec light exposure.

Image analysis for cell feature extraction

For each gene-silencing condition of each biological repeat, 21,000 images were processed to extract cell morphology features. To analyze cell size and shape changes induced by gene silencing, we developed a MATLAB-based image analysis pipeline, which includes three major steps: image preparation, feature extraction, and image post-processing (Fig. 6c). In step one, image sequences were fed into the MATLAB program and subtracted from the corresponding background image to correct any inhomogeneous illumination. The program automatically generates background images, in which each pixel value is computed as the mode pixel intensity value among the same pixel of the entire corresponding image sequence. After illumination correction, step two detects cell edges by looking for the local maxima of the bright field intensity gradient, following which the program closes edge gaps, removes cells connected to the image borders, cleans small features (noise), and then fills holes to generate binary images and centroid positions for each single cell. Cell locations were then traced and stuck cells were removed by a double-counting filter if present. The double-counting filter excludes measurements collected around the same location with similar cell sizes using a Gaussian kernel density method (bandwidth = 0.09) when the estimated density for a certain location and size exceeds a particular threshold. The maximal density value for experimental runs where no repeated measurements were observed was used as the threshold. This procedure was manually validated using visual inspection of the removed cells. Binary images passing the double-counting filter were used to create coordinates (X, Y) of cell outlines, which leads to a range of cell size and shape parameters, including cell diameter and area, solidity, roundness error, circularity, and intensity spatial relationship enclosed within the cell boundary. Cell area was computed as the 2D integration of the cell outline, and the cell diameter was computed as . Solidity is the ratio of cell area to the area of the smallest convex polygon that contains the cell region. Roundness error was computed as the ratio between the standard deviation and mean of radii on the cell outline measured from the centroid. Circularity was calculated as 4Areaπ/Perimeter². The 2D intensity distributions within cell outlines were used to derive peak locations and count peak numbers, which is a measure of intensity spatial relationship and a gating parameter to remove clumped cells. In the post-processing step, data were cleaned using three filters with the following gating threshold. To remove large clumps, the peak-solidity filter removes data outside of the polygonal region defined by {(0.9, 0), (0.9, 3.2), (0.934, 8.26), (1, 28), (1, 0)} in the (solidity, peak No.) space. Then, the roundness filter removes cells with weird shapes by excluding data with a roundness error higher than 0.3 or a circularity lower than 0.6. Finally, the small size filter removes cell debris whose major diameter is lower than 15 μm (12 μm) or minor diameter is lower than 12 μm (10 μm) for images photographed at the top (bottom) outlet microchannels.

Statistical assessment of gene-silencing effects in high-throughput single-cell experiments

To analyze the experimental results, we began by examining how and where cell size distributions differ between each of the gene/gene-pair silenced groups and their respective scrambled control groups. We thus conducted various statistical analyses to investigate these size distribution disparities. First, we compared the difference in median cell size (i.e., diameter measured in μm) between the gene/gene-pair silenced cells and their scrambled controls. We performed two different tests – a Wilcoxon signed rank test (Fig. 6f) and a bootstrap quantile test at the 0.5 quantile level. In accordance with the PCS framework, we used two different tests to ensure that our findings are robust to this arbitrary modeling choice and that the underlying assumptions do not drive our results. We note that the difference in median cell size was of greater interest than the difference in mean cell size due to the heavy right-skewness of the cell size distribution. Still, as an additional stability check, we performed a bootstrap-t test for the difference in trimmed means with varying levels of trimming (ranging from 0-0.3). These differences in trimmed mean tests (data not shown) yielded similar results to the tests for difference in medians, providing further evidence for the robustness of our conclusions. Secondly, in addition to comparing differences in central behavior, we compared differences in upper quantiles of size distributions for the gene/gene-pair silenced cells versus the scrambled controls. Identifying cell size differences at these upper quantiles, which focus on the larger, hypertrophic cells, is particularly relevant for the pathologic phenotype of cardiac hypertrophy and its clinical implications. To assess these differences in cell sizes at upper quantiles, we performed a bootstrap quantile test at the 0.6, 0.7, 0.8, and 0.9 quantile levels (Extended Data Fig. 8). All tests are performed on each experimental batch separately. To be as conservative as possible when claiming a significant effect, the maximum p-value across batches is reported in the main text. Similar analyses were conducted for assessing differences in morphological features (i.e., cell roundness error and normalized peak number).

Statistical assessment of non-additivity in high-throughput single-cell experiments

We also assessed whether the silenced gene pairs (i.e., CCDC141-TTN and CCDC141-IGF1R) are interacting in an additive or non-additive (i.e., epistatic) way to affect cell size in the high-throughput single-cell experiments (Fig. 6g). More formally, to assess this (non-)additivity for a given gene pair, say gene 1 and gene 2, we fit the following quantile regression:

where y_i is the diameter (μm) of cell i, x_i1 is an indicator whether gene 1 was silenced in cell i, x_i2 is an indicator whether gene 2 was silenced in cell i, g_i encodes the batch identifier from which cell i came (so that µ_gi is a batch effect term), and ε_i is the random error or noise term for cell i. This regression was fitted using the scrambled control cells (x_i1 = x_i2 = 0), gene 1-silenced cells (x_i1 = 1; x_i2 = 0), gene 2-silenced cells (x_i1 = 0; x_i2 = 1), and gene 1 – gene 2 jointly silenced cells (x_i1 = x_i2 = 1). Under this regression model, we tested the null hypothesis of β₁₂ = 0 versus the alternative hypothesis of β₁₂ ≠ 0 via a percentile bootstrap t-test and a traditional t-test (Extended Data 8) for varying quantile levels (0.5, 0.6, 0.7, 0.8, 0.9). A small p-value suggest that the gene pair is a non-additive epistatic interaction under the above model. We note again that two different tests were performed to check the robustness of our conclusions against modeling assumptions associated with each statistical test. To further bolster the robustness of our conclusions, we repeated this assessment of epistasis using the rank-based inverse normal-transformed cell diameter as the response y and under an ordinary linear regression model, finding that both the p-values and the direction of the non-additive interaction effects are similar to the reported quantile regression results. We thus omit these results for brevity. Since these regression models require comparisons between gene-silencing conditions (e.g., silencing CCDC141 and TTN vs only silencing CCDC141 vs only silencing TTN), and gene-silencing efficiency varied across silencing conditions, we only include cells from experimental batches with high gene-silencing efficiencies (>60%) for each regression. This helps to mitigate the possibility that the gene-silencing efficiency differences are driving spurious epistatic signals. We also conducted a simulation study to better understand how differences in gene-silencing efficiencies across batches might impact conclusions in Supplementary Note 3 (available on the website: https://yu-group.github.io/epistasis-cardiac-hypertrophy/simulations_efficiency). In general, we found that differences in our observed gene-silencing efficiencies do not typically lead to high false positive rates under the low signal-to-noise regimes that we expect in reality. We defer further discussion to Supplementary Note 3. Similar analyses were conducted for assessing differences in morphological features (i.e., cell roundness error and normalized peak number).