ShapeGTB: the role of local DNA shape in prioritization of functional variants in human promoters with machine learning

Datasets

To obtain the positive dataset we used single-nucleotide variants (SNVs) annotated as regulatory mutations in The Human Gene Mutation Database (HGMD^®) professional version (release 2016.2) within five kilobases (kb) upstream from the annotated transcription start sites (TSS) and provided sequences (Stenson et al., 2014). The total number of experimentally validated disease-related variants in our dataset is equal to 1,772. The control dataset contains SNVs from the 1000 Genomes Project (The 1000 Genomes Project Consortium et al., 2015) with a global minor allele frequency ≥1%. The overlapping elements of both sets were removed. Only variants lying within five kb upstream of TSS were selected for further analysis (Rosenbloom et al., 2015). The sequences of neutral motifs (not associated with disease phenotype) were retrieved from Ensembl with BioMart (Kinsella et al., 2011). The total number of negative examples in our dataset is equal to 3,806. We ensured that positive and negative motif sets are matched in their basic properties (Kolmogorov–Smirnov two sample test results for GC-content distributions are as follows D-statistic = 0.02, p-value = 0.48, null hypothesis of identical distributions retained). Distributions of TSS distances in the two sets differed, but we made sure that it does not affect obtained results (see Supplementary Material S5).

Machine learning pipeline

We split the available data into training and test sets randomly keeping the ratio 8:2. Full training set contained 1,417 positives and 3,045 negatives, full test set contained 355 positives and 761 negatives. Training set was used to build feature ranking, train classifiers, and optimize their parameters, while test set was left for final validation and for comparison with other prediction methods. To validate our methods internally on the training set we used a cross-validation strategy in which in each fold SNPs from a single chromosome formed test set and SNPs from other chromosomes formed training set. This eliminated possibility of overfitting during parameter tuning and feature selection procedures, and additionally demonstrated whether our method generalizes across different chromosomes.

We applied the Monte Carlo feature selection (MCFS) algorithm (Draminski et al., 2008) to perform feature importance ranking. It is a universal feature selection strategy combining random subspace methods with decision trees. A random subset of the original features is drawn in each iteration of the algorithm and an equivalent of random forest is induced using the selected variables. Feature importance ranking is constructed based on all induced trees. Additionally, meaningful interdependencies between features are discovered by calculating how often two features are used together to predict the class value. MCFS aims at finding all features relevant for the classification task, and it guarantees that with sufficient number of iterations all features can be tested. Following general guidelines by the authors of the algorithm, we set the number of iterations to 1,000 and the subset of original features considered in each iteration to 0.25.

In the classification task gradient tree boosting was used (GTB)—a popular tree-based ensemble algorithm (Friedman & Meulman, 2003). It is known to perform very well in many domains, often outperforming methods such as random forest, support vector machines or neural networks (Sheridan et al., 2016; Ladds et al., 2016; Babajide Mustapha & Saeed, 2016). The key idea behind GTB is to build trees sequentially, training a tree at each step to explain the prediction error made by the combination of existing trees. Usually the trees are regularized to prevent overfitting. We used the state-of-the-art implementation provided by XGBoost library (Chen & Guestrin, 2016). Through cross-validation performed on the training set we selected optimal parameter values (number of trees—300, maximal tree depth—8, learning rate—0.1).

Comparison with existing approaches

Presently, the field of prediction and prioritization of human noncoding regulatory variants still lacks a large, independent and publicly available gold-standard dataset for training, testing, and validating existing in silico approaches. The comparison of our method to the current state-of-the-art methods is hampered even further by different aims and objectives. To our best knowledge all available tools were designed for genome-wide, regulatory variants prioritization and there are no computational methods focused on promoter regions. Nonetheless, we compared performance of our algorithm with other tools on our own hold-out test set and on independent high-quality data from ClinVar database (January 5, 2017 release) after excluding variants present in our training data (Landrum et al., 2016). Our hold-out test set contained 355 positives from HGMD and 761 negative examples from 1000 Genomes Project. External validation set contained 32 positive examples labeled as pathogenic in ClinVar database and 761 negative examples from 1000 Genomes Project (not present in our train set).

Features groups

We used the following feature groups to annotate each SNV in our pathogenic and control datasets (more detailed description can be found in Supplementary Material S1 and S4):

1. DNA sequence (52 variables): 9-nt sequence motifs centered on the mutated nucleotide. The sequence was encoded using 4-bits binary coding. Additional 12 binary (4-nt by three mutations) variables indicated what type of mutation occurred (e.g., A → C, G → T, etc.).

2. Local DNA shape features (88 variables): HelT, MGW, ProT, roll values in span of 9 nt. Differences (_diff) between reference and mutated scores were added as additional features (Chiu et al., 2016).

3. GC-content (8 variables): GC-content in span of 7- and 9-nt for reference and mutated sequences separately. Differences between the reference and mutated scores were added as additional features.

4. Histone modifications (38 variables): ChIP-seq data for histone 3 lysine 9 acetylation (H3K9ac) and histone 3 lysine 4 trimethylation (H3K4me3) across 16 cell lines from ENCODE (Ram et al., 2011). For H3K9ac, H3K4me3, or either modification mean values over all cell lines and binary variables indicating modifications in any cell line were added.

5. Transcription Factor Binding Sites (12 variables): TFBS ChIP-seq clusters (V3) from ENCODE data retrieving binding sites of top 10 TFs with the highest binding site coverage. Mean value over all TFs and 0–1 indicator of any TF occurrence were added in addition (ENCODE Project Consortium, 2012).

6. Transcription factor binding disruption (one variable): P-value of disrupting putative strongest transcription factor binding site due to mutation was calculated with Annotation of Regulatory Variants using Integrated Networks (ARVIN) algorithm (Gao et al., 2018) using Cis-BP database (Weirauch et al., 2014).

7. Maximum transcription factor binding log-odds ratio score (one variable): Maximum TF binding log-odds ratio score for reference and mutated sequences among scores calculated with ARVIN algorithm (Gao et al., 2018, Weirauch et al., 2014).

8. DNase I hypersensitivity (one variable): ENCODE DNase clusters (V3) from 125 cell line types (John et al., 2011; Thurman et al., 2012; Rosenbloom et al., 2013).

9. Evolutionary conservation (10 variables):

   1. GERP ++: Genomic Evolutionary Rate Profiling scores (Davydov et al., 2010).

   2. PhastCons: PhastCons conservation score by vtools (San Lucas et al., 2012).

   3. Z-score: recalculated Z-score values defined in our previous work (Wyrwicz et al., 2007) on whole genome human–mouse alignments (genome builds hg19 and mm9 (Chiaromonte, Yap & Miller, 2002; Kent et al., 2003; Schwartz et al., 2003) from UCSC Genome Browser (Kent et al., 2002) for the reference and mutated sequence and for window length 7 and 9. Differences of Z-scores for the reference and mutated sequence were added.

10. Dinucleotide content (16 variables): Observed vs. expected frequencies of 16 possible pairs of nucleotides appearing in the short sequence motif.

Feature importance

From MCFS we obtained the ranking of all 227 features according to their relative importance in the classification problem. Each feature group contained multiple individual features with different ranks in the overall ranking. In the context of machine learning task, usefulness of a particular group should be determined by the best performing features from this group.

Figure 1 presents detailed feature ranking including all features from each group. Generally, features that contribute to the correct classification mostly belong to GC-content group, shape group, and sequence group. Other feature groups were of lesser importance (the full ranking is included as Supplementary Material S2, feature names glossary as Supplementary Material S4). The most important feature was the difference in GC-content between the reference and the mutated sequence fragment (rank 1). Features describing raw nucleotide sequence and dinucleotide content appeared in the middle of the ranking. Among the shape features, those describing the closest neighborhood of the mutated nucleotide were the most important. This is not surprising because differences in shape are expected to have local effects on DNA properties. Among the shape features attributes concerning ProT were ranked as the most important, attributes concerning HelT and roll followed, and attributes concerning MGW occurred lower in the ranking. What is notable, most of the features appearing among the top 20 concerned differences in shape properties between SNP and wild type. Features derived from transcription factors were less important than sequence-based features. Histone modifications, conservation scores and DNase I hypersensitivity score were not identified as particularly informative features.

To investigate the role of individual features we calculated Welch’s t-score capturing the relationship between particular feature and class value. The decrease of GC-content between the reference and the mutated sequence correlated negatively with functionality (t-score −8.2088 for decrease for motif length 7, t-score −11.3710 for decrease for motif length 9), while the increase of ProT value correlated positively (t-score 9.7417 for increase immediately before the modified nucleotide, t-score 5.5047 for increase immediately after the modified nucleotide).

The role of each feature in a classification task lies not only in its correlation with class value, but also in how well it complements with other features. For example, Fig. 2 presents joint distributions of the two most important features in the two classes (difference of GC-content between the reference and the mutated sequence, difference of ProT at the 3rd position between mutated variant and wild type). For non-functional SNPs the features are uncorrelated, but there is a visible negative correlation for functional SNPs. MCFS allows studying that kind of dependencies through its interdependency discovery function. The full list of feature interdependencies and their relative strength is included as Supplementary Material S3. Figure 3 presents graph of the strongest interdependencies among the top selected features (GCSCORE, GC composition; SEQ, sequence feature; ROLL, roll; HELT, helix twist; PROT, propeller twist). Difference in GC-content acts as a central hub and interacts strongly with all groups of shape features except MGW. The simplified intuition is that functional SNPs should increase GC-content of the motif, and at the same time increase rotation of the DNA strand accordingly.

Classifier performance

Obtained feature ranking suggests that a large portion of information is contained in features derived from the DNA sequence, and features describing evolutionary conservation and functional properties play less significant role. To verify this hypothesis, we performed a cross-validation experiment (with folds determined by chromosomes) on the train set by training GTB classifier on different combinations of feature groups. Calculated values of multiple performance measures are presented in Table 1.

Classifier based on all available features performed better than the classifier using only 25 best ranked features. Among individual feature groups GC-content produced classifier with the largest AUC ROC (0.78). Combining GC-content with shape features and sequence features allowed achieving AUC ROC 0.98. No other combinations of features performed better. These results show that shape features are more meaningful when combined with another feature. In further experiments classifier trained on sequence, shape, and GC-content was used. We named this classifier ShapeGTB.

We compared final ShapeGTB classifier with more general SNP prioritization methods, which did not focus specifically on promoter regions: CADD, FATHMM-MKL, and DeepSEA. Figure 4 present precision-recall curves calculated on the hold-out test set constructed from our data (HGMD and 1000 Genomes Project) and for smaller experimental dataset (ClinVar and 1000 Genomes Project). Area under precision-recall curve can be interpreted as average precision (AP), and is an aggregated measure of classifier performance. It is preferred over AUC ROC when the problem is characterized by a large class imbalance. On the hold-out test set, ShapeGTB outperformed general-purpose methods by a large margin (AP 0.93 vs. 0.47–0.55). On the external validation set, the ShapeGTB aggregated performance was comparable with FATHMM-MKL (AP 0.47 vs. AP 0.42). However, shapes of precision-recall curves for those methods were very different: FATHMM-MKL displayed high precision only for small subset of examples, while ShapeGTB precision was relatively stable even for large values of recall. Differences between results obtained for the two datasets suggest that ClinVar-derived positives have different characteristics and pose a greater challenge. We speculated that the gap between ShapeGTB and reference tools on the hold-out test is due to inclusion of shape features and their interactions with GC-content. To verify this, we randomly permuted these features in our test set and evaluated performance of ShapeGTB again on permuted data sets. AP of ShapeGTB with GC-derived features permuted was 0.80, with shape-derived features permuted 0.44, and with both kinds of features permuted 0.35 (Fig. 5). This once more corroborates the hypothesis that shape features together with GC-content provide important information for distinguishing functional SNPs in our data set.