Human cardiac cis-regulatory elements, their cognate transcription factors, and regulatory DNA sequence variants

Experimental detection of cardiac CREs is incomplete

A schematic overview of our approach is shown in Supplemental Figure S1. We first performed DNase-seq to identify open chromatin regions through DNase I hypersensitive sites (DHSs) as potential CREs in two adult human hearts (left ventricles) together with three publicly available heart DNase-seq data sets (Fig. 1A). Peak calling with MACS2 (Zhang et al. 2008) identified varying numbers of DHSs (50,000–110,000) across the five samples (Methods; Supplemental Fig. S2A). We defined DHSs as 600-bp regions centered at the identified summits. DHSs frequently cluster in small regions, and 22%–32% of extended DHSs overlap their neighbors, forming larger DHSs with multiple summits (Supplemental Fig. S2B). We next compared heart DHSs with other tissue DHSs using the top 50,000 DHSs from each sample and the Jaccard index (number of bases in the intersection over the union for each pair). These comparisons show higher similarity (∼50%) between the adult cardiac samples than with other tissues (∼30%), including fetal heart (Methods; Supplemental Fig. S2C). Thus, DHS maps do uncover tissue-relevant CREs, although many regions are open across many tissues. Nevertheless, technical and biological variation affects the adult cardiac DHS maps. To maximally detect CREs, we aggregated all DHSs across the five adult heart replicates and identified ∼160,000 distinct regions (“observed DHSs”) covering ∼4% of the genome. Using multiple samples is important because ∼40% of DHSs were detected only once across the replicates, while only ∼22% detected across all five (Fig. 1B).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 1.

A machine learning algorithm (gkm-SVM) accurately predicts cis-regulatory elements. (A) An example of heart DNase-seq signals (raw data) and peaks (MACS2) at the GATA4 locus across multiple human samples. (B) A genome-wide heat map of DNase-seq read densities in 1000-bp windows centered at heart DHSs. Randomly sampled 1000 regions were used. Regions were grouped based on the configuration of the DHS peaks across the five samples with at least one observed DHS. (C) ROC curves of gkm-SVM models for five replicates against reserved test sets. (D) Comparisons of the fraction of DHSs overlapping predicted regions (precision) and fraction of predicted regions overlapping DHSs (recall).

Sequence-based models can predict additional CREs

To elucidate the sequence code for these cardiac CREs, we used the gapped k-mer support vector machines (SVM) method, gkm-SVM (Ghandi et al. 2014), widely accepted as one of the state-of-the-art sequence-based methods (Kreimer et al. 2017). Also, extracting predictive sequence features from gkm-SVM is easier than from other methods. We used an improved software LS-GKM (Lee 2016) that enabled training on larger data sets (∼100,000 DHSs) (Supplemental Fig. S3A). The gkm-SVM attained considerable accuracy in classifying CREs from non-CREs on reserved Chromosome 9 data (Methods; Fig. 1C). Of note, the choice of a specific chromosome did not affect the method's performance as its accuracy was unchanged with fivefold cross-validation (Supplemental Fig. S3B). Importantly, our model trained on a single data set can predict DHSs observed in the other data sets and systematically assigns higher scores to DHSs experimentally missed than random genomic regions, demonstrating that sequence-based predictions do indeed identify true DHSs (Supplemental Methods; Supplemental Fig. S4A).

The number of samples used for model building determines a model's accuracy. Our tests show that using all five samples identified an additional 5%–10% true DHSs than any single sample, although the maximum accuracy was rapidly attained (Supplemental Methods; Supplemental Fig. S4B). Further, the randomness of negative sets for training is also a major source of variation. To minimize this effect, we averaged 10 different models, trained on independently generated negative sets, and used the combined model for DHS prediction at balanced precision and recall rates (∼43%) (Methods; Fig. 1D).

Scanning the entire genome yielded an additional 88,026 distinct “predicted” heart DHSs after excluding regions that overlapped observed ones (Methods; Fig. 2A). Thus, these predicted DHSs do not suffer from overfitting by construction. For validations, we used several functional annotations to compare them to observed regions. First, we assessed sequence conservation against randomly permuted regions (Supplemental Fig. S5A): Both observed DHSs and, to a lesser extent, predicted DHSs significantly overlap genomic conserved elements (binomial test, P < 2.2 × 10⁻¹⁶) (Davydov et al. 2010). Second, we asked how frequently predicted DHSs were in open chromatin in other tissues. We defined two different DHS sets: (1) a universal set by combining DHSs from all ENCODE and Roadmap data sets, except from adult heart; and (2) a heart-related set of DHSs from cardiac-related samples only (Methods). Both observed and predicted DHSs significantly overlapped both DHS classes (Supplemental Fig. S5B,C): 94.7% of observed and 57.7% of predicted DHSs are open in heart-related tissues, with an additional 30.7% in other cell types (Fig. 2A). Third, we compared H3K27ac histone modification marks in heart tissues to show the same pattern: 52.3% of observed but only 10.9% of predicted DHSs overlapped H3K27ac marked regions (Supplemental Fig. S5D). Yet, we identified ∼7500 additional regions, under the stringent criteria that predicted regions overlap H3K27ac marks and are heart-related DHSs, not detected by DNase-seq alone. Predicted DHSs are largely restricted to specific tissues and cells (Fig. 2B; Supplemental Fig. S5E) while nearly half of observed DHSs are open across many different noncardiac cell types. Moreover, predicted DHSs show systematically weaker DHS signals than observed DHSs (Supplemental Fig. S5F).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

Learned sequence features predict additional CREs. (A) Venn diagram of observed and predicted DHSs and their proportions in noncardiac cell types. (B) Heat map of DHS signal intensities in heart-related cell types and tissues at observed (left) and predicted (right) DHSs. Randomly sampled 5000 regions from observed and predicted DHSs were used. Observed DHSs are mostly located in ubiquitously open chromatin regions in cardiac-relevant cells and tissues, while predicted DHSs show greater cell-type specificity.

Varying degrees of overlap with functional annotations led us to further investigate potential sources of variation in CRE identification. We divided predicted DHSs into four distinct categories based on overlap: heart H3K27ac histone modifications and heart-related DHSs (tier 1); heart-related DHSs only (tier 2); universal DHSs only (tier 3); and the remainder (tier 4). For each of these categories, we considered the (1) proportion of ambiguous bases, (2) average SNP (single nucleotide polymorphism) frequency per base (common and rare SNPs, 1% MAF as a threshold), (3) proportional overlap with H3K9me3 histone modifications in heart tissues, a representative heterochromatin mark (Nakayama et al. 2001), and (4) proportional overlap with FANTOM5 enhancers (Methods; Lizio et al. 2015). As a positive control, we included observed DHSs for comparison (Supplemental Table S1). First, the predicted DHSs in the lower tiers (3 and 4) have poorer mappability: Predicted DHSs in tier 4 have 2–3.5× more ambiguous bases than observed DHSs, depending on read length. Since the functional annotations (H3K27ac and heart-related DHSs) were also identified by sequencing, the variability in CRE detection can be explained at least in part by differential mappability. SNP frequencies also follow the expected pattern; higher SNP frequency in the lower tiers (2 & 3) imply that these predicted DHSs are more variable, making read mapping difficult. Predicted DHSs in tier 4 show the lowest SNP frequency because their extremely poor mappability reduces SNP identification (Nielsen et al. 2011). DHSs predicted in the lower tiers are more enriched in heterochromatin, suggesting that chromatin organization, a feature difficult to predict using sequence-based models only, influences CRE detection as well. Lastly, we used FANTOM5 enhancers (Lizio et al. 2015) as an orthogonal enhancer validation set and compared them to predicted DHSs. Significant proportions of predicted DHSs in tiers 1, 2, and 3 (12.3%, 5.9%, and 2.7%) were FANTOM5 enhancers (binomial test, P < 2.2 × 10⁻¹⁶ for all cases), but tier 4 had only 0.5% FANTOM5 enhancers (P < 0.97).

Another potential cause of variation in CRE detection lies in peak calling, with some predicted DHSs being missed due to being below our detection threshold. Thus, we recalled DHS peaks with relaxed thresholds (false discovery rate <0.1, 0.15, 0.2, and 0.25) and identified larger numbers of DHSs (193, 205, 242, and 280 k, respectively). We calculated the proportion of the original predicted DHSs overlapping these newly identified DHSs (Supplemental Table S2) to show that predicted DHSs, especially in the higher tiers, overlap many of these less stringent DHSs. Thus, a significant fraction of predicted DHSs are likely true. Taken together, our results suggest that sequence-based predictions of an additional ∼55% cardiac CREs are highly complementary to their experimental detection and necessary for obtaining comprehensive CRE maps.

Sequence-based models can also predict cardiac TFs

We next tested whether the sequence features of gkm-SVM allow identification of the cognate TFs through their TF binding sequences (TFBS) in the corresponding tissue/cell types (Lee et al. 2011; Gorkin et al. 2012). First, we used the distribution of all 11-mer weights in the SVM model and compared it to those that match TFBSs for known heart TFs. Indeed, TFBSs for CTCF and MEF2A as well as other known cardiac factors are significantly enriched in the top fifth percentile of the SVM weight distribution, while an exemplar noncardiac factor such as POU5F1 is not (Fig. 3A). Thus, to search for cardiac TFs, we used the Cis-BP database (Weirauch et al. 2014) augmented by 93 new position weight matrices (PWMs) of C2H2 Zinc-Finger (ZF) TFs (Schmitges et al. 2016), to show that 11-mers matching 54% of PWMs (473/868) associated with 506 TFs were enriched in the top fifth percentile of the distribution (Supplemental Data S1). As expected, 11-mers matching these predicted cardiac TFs systematically have larger SVM weights than those that don't (Supplemental Fig. S6). Since not all TFs are expressed in any given tissue, we further restricted attention to the 334 TFs expressed in heart left ventricles using gene expression profiles from the GTEx project (Supplemental Methods); Fisher's exact test confirmed that predicted TFs were significantly associated with cardiac gene expression (Fig. 3B).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 3.

Learned sequence features identify cardiac TFs. (A) SVM weight distributions of PWM-matched 11-mers for the top two scoring TFs (CTCF and MEF2A), five well-known cardiac-specific TFs (GATA4, NKX2-5, HAND1, HAND2, SRF, and HEY2), and a cardiac-irrelevant TF (POU5F1). The fifth percentile is shown as a green line. The enrichment is defined as the fraction of TFBS-matching 11-mers in the top fifth percentile of all 11-mers compared to the expected fraction (5%). (B) 2 × 2 contingency tables comparing TFs enriched in heart DHSs with TFs expressed in heart left ventricles (top) and atrial appendages (bottom). (C,D) −log₁₀(P) of one-sided Fisher's exact test for every tissue tested using all TFs (C) and after removing commonly expressed TFs (D); the two tissues from adult heart (atrial appendage and left ventricle) are highlighted in red with the Bonferroni-corrected P-value threshold (0.05) shown as a dashed line.

These results are physiologically relevant since left ventricles and atrial appendages from adult hearts were the two most significant tissues among the 53 GTEx tissue gene expression profiles we examined (Fig. 3C). Many other tissues also showed significant association because many TFs are expressed across multiple tissues. Thus, enrichment tests after removing 353 commonly (>90% of tissues) expressed TFs dramatically reduced the significance for every tissue except for the two heart tissues (Fig. 3D). Thus, cardiac CRE functions are activated by a large number of both commonly expressed and cardiac-specific TFs. Our analyses revealed 78 selectively expressed and CRE-enriched TFs in heart tissues, of which 29 are C2H2 ZFs with diverse binding specificities, with 44 others belonging to seven other major classes; 16 basic-helix-loop-helix (bHLH), nine Nuclear Receptor (NR), five Sox, four Homeodomain (HD), four basic leucine zipper (bZIP), three GATA, and three ETS transcription factors (Fig. 4). Many expressed TFs in the same family have similar sequence specificities, explaining some part of the expected functional redundancy.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 4.

A list of 78 TFs enriched in heart DHSs as well as selectively expressed in cardiac tissues. The 78 TFs were grouped based on their DNA binding domain families. These PWMs were also clustered based on their sequence specificity. For each TF, cluster number (Clus), name, gene expression (FPKM) in left ventricle, fold enrichment (nfolds), and PWM are shown, respectively.

We validated these predicted cardiac TFs also using ENCODE ChIP-seq data (Supplemental Methods) based on two classes: potential cardiac TF ChIP-seq data (n = 1054 for 176 TFs) and the remainder (n = 313 for 126 TFs). For each, we calculated the overlap between ChIP-seq peaks and our observed and predicted DHSs to show that candidate cardiac TF-bound regions overlap heart DHSs significantly more often than noncardiac TFs (P < 2.2 × 10⁻¹⁶ for observed DHSs, P < 2.84 × 10⁻¹⁶ for predicted DHSs, with one-tailed two-sample Kolmogorov-Smirnov tests) (Supplemental Fig. S7).

Common cardiac DHSs are mostly promoters and CTCF sites

Many DHSs are accessible across diverse cell types (Xi et al. 2007; Song et al. 2011; Thurman et al. 2012). These common DHSs are typically enriched in transcription start sites (TSSs) or CTCF binding sites or both, with the remainder enriched in cell-type–specific distal enhancers. To discriminate the cardiac DHSs, we defined common DHSs as regions that are open in ≥30% of all ENCODE/Roadmap tissue samples (Supplemental Methods; Lee et al. 2015). Consistent with previous observations (Xi et al. 2007; Song et al. 2011; Thurman et al. 2012), ∼47% of observed heart DHSs are common DHSs, of which ∼63% overlap TSSs or CTCF bound regions (Supplemental Fig. S8). To understand the sequence features of these classes of DHSs we separately trained a gkm-SVM model after removing these common DHSs. In this heart-specific model, consistent with the CTCF ChIP-seq peak overlap analysis, most of the predictive 11-mers that showed a major decrease in SVM scores (Z-score differences > 6) were CTCF binding sites. In the top five, motifs showing a major decrease were additionally ZBTB4, PLAG1, SP4, and KLF10 binding sites, known to be promoter-binding TFs (Supplemental Fig. S8C,D). These analyses provide a principled, statistical way to distinguish TSSs, CTCF binding sites, and distal enhancers from DHS data, in a tissue- or cell-type–specific manner.

Predicted cardiac regulatory variants affect chromatin accessibility and gene expression

Our sequence-based model allows systematic predictions of whether a DNA sequence variant within a CRE is likely to affect its function, thus making detection of (cardiac) regulatory variants possible. To assess this function, we used the deltaSVM method (Lee et al. 2015; Beer 2017; Kreimer et al. 2017) to first assess a small set of high-confidence cardiac regulatory variants associated with allele-biased DHSs (Methods). We tested two different gkm-SVM models, one trained on all heart DHSs (“generic model”) and the other on cardiac-specific DHSs after removing common DHSs (“specific model”): deltaSVM scores from both models were significantly correlated with their allele-biased chromatin accessibility (Fig. 5A; Supplemental Fig. S9A). Both models achieved comparable precision—46% and 52% at 40% recall—using 4× larger sets of control variants with no allele bias (Fig. 5B). Although their performance is similar, they do not predict the same variants since ∼30% of variants are unique to each model (Supplemental Fig. S9B). This result is consistent with the k-mer weight distribution differences between the models (Supplemental Fig. S8C) and suggests that the specific model is better at detecting cardiac-specific TFBSs by ignoring CTCF and general promoter TFBSs.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 5.

Identification of common, cardiac-specific regulatory sequence variants. (A) deltaSVM scores from the cardiac-specific model as compared to allele-biased chromatin accessibility. (C) Pearson correlation coefficient, (n) number of variants, (P) t-distribution P-value. (B) Precision-recall curves of deltaSVM scores of allele-biased DHSs against 4× larger control SNP sets; the dashed green line indicates the recall rate 40%. Error bars are the standard errors calculated from 10 independently sampled control SNP sets. (C,D) Statistical significance (P) of deltaSVM SNPs from the generic (C) and specific (D) cardiac models as compared to eQTLs using the χ² test; the two heart tissues are highlighted. (AA) Heart atrial appendages, (LV) heart left ventricles. The Bonferroni-corrected P-value threshold (0.05) is shown as a dashed line.

These results prompted the question of whether deltaSVM predicted variants affected local gene expression. We identified all high scoring variants within heart DHSs (Supplemental Data S2) and compared them to GTEx expression quantitative trait loci (eQTLs) from 44 different tissues (Methods): deltaSVM variants predicted by the generic model are significantly associated with eQTLs in many tissues (Fig. 5C). This significance is strongly correlated with the number of eQTLs in CREs, which is also correlated with the sample size used for gene expression studies (The GTEx Consortium 2015). On the other hand, variants predicted by the specific model are mostly associated with eQTLs from heart tissues, although their overall statistical significance is lower (Fig. 5D). We concluded that, despite the pleiotropy of many regulatory variants, the specific model does identify cardiac-specific regulatory variants with increased specificity but at the cost of missing some regulatory variants common to many tissues. Further, inclusion of predicted DHSs enhanced detection by increasing the statistical significance of eQTL associations (Supplemental Fig. S9C,D).

Predicted cardiac regulatory variants explain cardiac phenotypes

We hypothesize that these predicted causal variants, within observed and predicted CREs explain a significant fraction of cardiac phenotype heritability, which we tested using QTi GWAS (Arking et al. 2014), an intermediate trait involved in long QT syndrome and sudden cardiac death (Tomaselli et al. 1994). We used our published QTi meta-analysis on 76,061 European ancestry subjects and ∼2.7 million SNPs and performed Q-Q analysis. First, variants within heart CREs (Fig. 6A, red dots) are significantly enriched in QTi GWAS SNPs in comparison to all common SNPs (black dots). A subset of these heart CRE variants, predicted to be causal by deltaSVM, especially by the specific model (green dots), are further enriched in QTi-associated variants (Fig. 6A).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 6.

Cardiac regulatory variants explain QT-interval heritability. (A) Q-Q plots of QTi GWAS results using different subsets of common genome-wide sequence variants (numbers of variants in parentheses); the dashed gray line indicates the genome-wide significance threshold (P < 5 × 10⁻⁸). (B) Comparison of enrichment P-values between 10 predefined cell-type group functional annotations. (C) Enrichment values, estimated as the fraction of the heritability (h_G²) explained by variants over the fraction of the SNPs, for various heart DHS-based annotations. (Cardiovascular) the predefined cell-type group functional annotation for cardiovascular tissues, (H3K27ac) H3K27ac marks in heart, (Dickel) heart CRE map from Dickel et al. (2016), (obsDHS 600 bp) observed DHSs, (obsDHS 1 kb/2 kb) observed DHSs with 1-kb/2-kb extension, (allDHS 2 kb) observed and predicted DHSs with 2-kb extension, (dSVM) deltaSVM predicted variants. A combination of multiple annotations indicates intersection of them. For example, “obsDHS 2 kbp & dSVM” means a set of variants predicted by deltaSVM and overlapping observed DHSs with 2-kb extension. The proportion of SNPs are in parentheses; error bars denote standard errors estimated by a block Jackknife method; no enrichment is shown as a dashed line.

Taken together, these variants contribute substantially to QTi heritability, as shown for some functional annotations of the genome relative to other complex phenotypes (Yang et al. 2011). To quantify this contribution, we used linkage disequilibrium (LD) score regression methods on GWAS summary statistics (Bulik-Sullivan et al. 2015; Finucane et al. 2015) to evaluate the heritability contribution from CRE causal variants (Supplemental Methods). All common autosomal variants (1KGP SNPs with MAF > 5% in European ancestry subjects) explained 11.2% of QTi variation. Using predefined cell-type group functional annotations (Finucane et al. 2015), the cardiovascular cell-type group contributed to the most significant enrichment as compared to other tissues (Fig. 6B); this cell-type group includes lung tissues, so that the specificity for QTi may be diluted. Consequently, we tested our heart DHS-based annotations by comparing the enrichment under various definitions of causality (Fig. 6C). We first evaluated different DHS lengths, discovering that a higher heritability (61%) could be explained by extending DHSs to 2-kb lengths, although higher enrichment (7.9×) was achieved at the original definition of 600 bp. We surmise that narrower DHSs truncate some CREs, missing true variants. As a comparison, we also evaluated a recently published cardiac CRE map based on EP300 and H3K27ac bound regions from multiple cardiac developmental stages (Dickel et al. 2016). This enrichment is comparable to that achieved by our 2-kb heart DHSs (5.4× vs. 5.3×; one-tailed Z-score P < 0.51). The average size of Dickel's cardiac CREs is bigger than our 2-kb heart DHSs (3211 bp versus 2467 bp), but the number of elements is smaller (82,119 versus 133,102). The majority of Dickel's CREs (49,506 out of 82,119, or 60%) overlap our DHSs, but many regions are still detected by just one approach, suggesting that each method detects somewhat different types of CREs. Next, adding H3K27ac marks to DHSs further increases (from 5.7× to 8.5×) the enrichment by effectively filtering out less informative DHSs. While the difference between these two enrichments is only marginally significant (one-tailed Z-score P < 0.06), H3K27ac marks alone capture more variants than H3K27ac overlapping DHS variants (7.7% vs. 5.8%) but explains a smaller heritability fraction (44.3% vs. 48.8%) (Fig. 6C; Supplemental Fig. S10). Thus, some H3K27ac peaks do not capture CRE variants, consistent with the fact that H3K27ac peaks are typically located at CRE boundaries. Adding deltaSVM predictions greatly increased enrichment but explained a smaller fraction of heritability, likely due to false negatives in deltaSVM predictions. Nonetheless, our heart DHS-based annotations, combined with deltaSVM predictions and H3K27ac marks, achieved a highly significant 33.1× enrichment. Tier 1 predicted DHSs consistently increase the explained heritability by 5%–10% (i.e., 1∼3% additional heritability) without compromising enrichment (Fig. 6C; Supplemental Fig. S10). These results confirm that regulatory sequence variation is the major source of QTi phenotypic variation and that we can identify a majority of such causal variants. Of note, including predicted DHSs with weaker functional support (Tiers 2, 3, and 4) decreased overall SNP heritability (Supplemental Table S3). This is probably due to the assumption in the LD score regression method that effect sizes are normally distributed. Under this assumption, adding a large number of SNPs with very small effect sizes can introduce a downward bias. It is also consistent with our observation that many predicted DHSs are weaker compared to observed ones, making them less contributory to the phenotype.

To demonstrate our method's broad applicability, we also analyzed three additional cardiac phenotypes (systolic blood pressure, diastolic blood pressure, and pulse rate) and one noncardiac phenotype (BMI) using the UK Biobank project public data (Methods; Sudlow et al. 2015). Similar to the QTi result, predicted cardiac regulatory variants significantly contributed to all three cardiac-relevant phenotypes but not BMI (Supplemental Fig. S11). The heritability of blood pressure explained by the cardiac variants is consistently less than that for QTi and pulse rate, suggesting that regulatory variants active in other tissues, such as kidney and blood vessels (Hoffmann et al. 2017), also contribute to blood pressure phenotypes.

The heart CRE map we generated and the deltaSVM predictions of causality for specific variants within these CREs can be used to probe each of the known QTi GWAS loci in much greater detail (Supplemental Data S3). The QTi GWAS meta-analysis (Arking et al. 2014) identified 35 independent loci with genome-wide significant (P < 5 × 10⁻⁸) common variants. By restricting attention to significant SNPs and their LD proxies (r² > 0.9), we identified 149 variants as potentially regulatory (Supplemental Table S4). Of these, only ∼50% (n = 72) are highly associated (r² > 0.6) with one of the 67 sentinel SNPs (cf. some of the 35 loci have multiple independent sentinel SNPs). On the other hand, 108 variants have alternative measures of high association (|D′| > 0.9) suggesting that many index SNPs may tag multiple regulatory variants within their haplotypes. Note that these regulatory variants are not uniformly distributed across the associated loci: 75% are located within eight loci (PLN, NOS1AP, ELP6, LITAF, CNOT1, SCN5A, LAPTM4B, and KCNH2). These variants are significantly enriched in GTEx eQTLs as well: 104 of 149 variants are eQTLs in at least one tissue (2.7× enrichment; binomial test, P < 3.6 × 10⁻³⁰), of which 57 are eQTLs in the two heart tissues (7.5× enrichment; P < 2.2 × 10⁻³⁴), implying that these specific variants affect the QTi phenotype by perturbing gene expression in the heart.

As an independent validation of deltaSVM predictions, we finally analyzed a small number of GWAS variants obtained from a recently published study (Wang et al. 2016). In this report, 18 putative functional SNPs were selected based on subthreshold GWAS significance (5 × 10⁻⁸ < P < 1 × 10⁻⁴) for the QTi and QRS duration, additionally supported by other functional annotations (histone modification marks, DHSs, or eQTLs). These were tested by luciferase assays in human iPSC-derived cardiomyocytes. Of these, 14 were found within our heart DHSs, of which five were predicted to be potentially regulatory by deltaSVM. All of these deltaSVM-predicted SNPs showed differential enhancer activities by the luciferase assay (100% specificity), but the study also identified five additional regulatory SNPs (50% false negative rate) (Supplemental Table S5). Although the small sample size explains the lack of statistical significance (Fisher's exact test P < 0.125), the result is consistent with our predictions.