Contribution of each group of features to cancer breakpoint prediction

We built machine learning models to investigate the impact of every feature group on prediction of breakpoint hotspots (S7 Table) by modeling breakpoint hotspots based on each feature group. Each feature group was ranged according to the model performance (see Methods), and the results are presented in Fig 1A. The results show that the best feature group significantly (by 0.25) outperforms others for almost all cancer types, and this value even reaches 0.5 for 3 cancers (blood, pancreas and prostate). This outperforming group included TFs for 5 cancers (liver, skin, prostate, ovary, breast), non-B DNA structures for 2 cancers (brain and bone), genomic regions (pancreas), HMs (blood), and HDNase (blood).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Feature group ranging by model performance using two different metrics–lift of recall and ROC AUC.

A. Mean lift of recall for 0.03 probability percentile for each cancer type and feature group scaled and averaged for 99% and 99.5% labelling types. B. Mean ROC AUC for each cancer type and feature group scaled and averaged for 99% and 99.5% labelling types.

https://doi.org/10.1371/journal.pcbi.1008749.g001

Non-B DNA structures always enter the top 3 important groups for all types of cancers, TFs are present in the top3 list in 8 types of cancer (except for blood and brain), genomic regions–in 5, and HMs in 4 cancer types. A similar ranging is observed when using ROC AUC (Fig 1B) as the model evaluation metric: non-B DNA structures represent the first important group in 7 cancer types and rank among the top-3 most important groups for all cancers. For the remaining three cancer types the most important groups are HMs (blood), genomic regions (pancreas), and TFs (liver). TFs ranked in the top-3 list for 9 cancers, genomic regions for 6 cancers, and methylation for 3 cancers (brain, breast and skin).

When analyzing group feature importance for all cancers combined, the non-B DNA structures are in the first place followed by transcription factors and genomic regions (S2A and S2B Fig) with the mean lift of recall (3.18 and 3.02) and the mean ROC AUC (0.66 and 0.62) slightly higher for non-B DNA than for TF. The results indicate that large-scale chromatin organization does not contribute much to the breakpoint formation: HDNAse and TADs occupy the last two places and HDNase is the one before last in the list. The two remaining in-between groups–HMs and methylation have nearly equal influence.

We also performed an analysis to define the best individual features inside the best feature groups–non-B DNA and TFs. We established that among all non-B DNA features G-quadruplexes and direct repeats are the most influential factors for all cancer types, while short tandem repeats and Z-DNA also make a valuable contribution to model performance (S3 Fig). Moreover, the importance of the best feature is several times higher than that of the remaining features emphasizing the key role of the selected candidate. Breakpoint hotspots locations on density coverage plots for 20Mb region on the chromosome 13 (S4 Fig) clearly demonstrate a correlation between cancer breakpoint hotspots and different non-B DNA features. Concerning transcription factors, sets of the most important features are less similar to each other for different cancer types but nevertheless CTCF, GABPA, RXRA, SP1, MAX and NR2F2 are more frequently included in top features than other transcription factors (S5 Fig). Interestingly, for the majority of cancer types distant features are of higher importance than the local features, and this is a common observation for two feature groups (see Results below).

Additionally, we investigated the influence of individual features with the Boruta feature selection method (see Methods) considering individual contributions of every feature from every group (S6 Fig and S8 Table). Totally, the method selected 50 important features with 5–23 features per cancer type (Fig 2). Similarly to the feature group importance results presented earlier, the Boruta selection method also highlighted features mostly from non-B DNA, TFs and genomic regions groups (at least one feature of each group is included in a set of best features for 8 cancer types). Nevertheless, all groups except for TADs contribute to the model prediction power.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Boruta selected feature sets by cancer type and feature group.

The most important features for all cancer types selected by Boruta feature selection method for each cancer type separately.

https://doi.org/10.1371/journal.pcbi.1008749.g002

Performing importance analysis on Boruta feature selection procedure results, we observed that features of only four feature groups (non-B DNA, TFs, genomic regions and HDNase) are included in the top list according to the criteria of a feature ranking important in at least 300 of 3000 analyzed datasets for all cancer types (see Methods and S7 Fig). In the resulting list the top five features are represented by direct repeats and G-quadruplexes, both local and distant (see Results below), together with distant binding sites of SP1 transcription factor. Other individual features that were selected less frequently include Z-DNA, short tandem repeats, mirror repeats, transcription factors RXRA, NR2F2, GABPA, CTCF, genomic regions such as 5’ UTR, coding exons, 3’ UTR, promoters and downstream areas, and also HDNase, but their influence differs for cancer types (S8 Fig). Interestingly, the majority of the features were captured at 1Mb scale, not at 100kB scale, revealing the importance of taking into account a wider landscape for cancer breakpoint hotspots identification (see Results below).

Independently from the methods, either group or Boruta feature importance,two best performing feature groups have invariably been those of non-B DNA and TFs following by genomic regions, of which 5’UTR and coding exons were the most contributing. This result emphasizes the key role of non-B DNA structures and transcription factors for cancer breakpoint hotspots prediction, however the contribution of individual features from these two groups for each cancer type hotspots prediction is unequal. The importance of genomic regions such as 5’UTR and coding exons provides a hint about the major role of transcription processes in generating cancer breakpoints.

A comparison of model performance built on one group of features or on all features showed that an all-features model always outperforms the one-group-of-features model (either with the lift of recall or ROC AUC) with the exception of blood cancer (Figs 3 and S9 and S8 Table). We also report PR AUC as specifically informative for class-imbalance problems (Figs 3 and S9 and S8 Table). As we can see from Fig 3, according to PR AUC, the top three feature groups coincide with those determined according to ROC AUC or lift of recall for all types of cancer (except for the blood cancer for 99% labeling type), and the ranking of cancer types according to the model performances is generally preserved among three metrics. This result confirms the hypothesis that information from omics data improves the model prediction power [20]. However, it is of interest to note that the performance of the all-features model is only slightly higher than the performance of the model built on the best group of features. This result suggests the existence of structural variation signatures that would be very important to detect in the future, but is out of the scope of the present study.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Feature group contribution.

Contribution of different feature groups separately and all together. Distribution of the lift of recall (first panel) for 0.03 probability percentile, PR AUC (second panel) and ROC AUC (third panel) for each of one feature group-based models and models on Boruta selected features (for 99% labeling type).

https://doi.org/10.1371/journal.pcbi.1008749.g003

To evaluate goodness-of-fit of the all-features and single feature group models, we calculated McFadden’s R-squared (pseudo R-squared) for 99% and 99.5% labeling type (S10 Fig and S9 Table) by building logistic regression models. We found that for a half of cancers (breast, ovary, pancreas, skin, and blood) non-B DNA-based models produced the maximal pseudo R-squared while for the rest the all-features model performed the best. In absolute values, the maximal McFadden’s R-squared over all cancer types are observed for breast, prostate and ovary cancers with 0.087, 0.066 and 0.061 respectively.

We investigated whether the different feature group models identify the same or different breakpoint hotspots regions. Towards that end we calculated the mean Spearman correlation between predictions for each pair of feature groups (S11A Fig). The highest mean correlation over all cancer types is observed for histone modification-based models that are positively associated with TF- and methylation-based models (with correlation of 0.47 and 0.46 respectively). On the opposite, TAD-based models generate the most specific predictions with maximal correlation of 0.045 with TF-based predictions. Secondly, for each genome window in each test set we calculated the number of hotspots predicted with different number (1, 2, or 3) of feature groups and analyzed the proportion of true hotspots. A comparison of median quality of the best feature group model (corresponds to label "1") to overlapped predictions of 2- or 3-feature group models demonstrates that for the majority of cancer types a single feature group model is more precise than the overlapped predictions (S11B Fig). For several cancer types (uterus, breast, bone, pancreas) the median precision is only slightly higher with the increased standard deviation and the resulting number of identified hotspots (true positives) is lower. This fact, together with the observed uplift in combined models compared to the single feature group models, leads to a conclusion that in order to account for multiple factor contributions in cancer breakpoint hotspots formation we need to use a more complex model than the overlap between different several single feature group models.

The inclusion of distant features improves prediction power of the models

For cancer breakpoint hotspot modeling we used 100 kB coverage for all considered features as predictors. We decided to check whether model performance could be improved by the addition of different feature transformations to a current feature set (hereinafter referred to as local features). We created binary flags of feature presence, indicators of local/global coverage maxima and distant features (coverage aggregates on 1Mb window). We observed that a binary feature addition made no contribution to model quality while maximum indicators gave only slight performance increase (S12A Fig). On the contrary, distant features being added to local features improved ROC AUC by 0.03 on the average. However, the effect differs between cancer types being clearly positive for pancreatic, ovary, breast and liver cancer and slightly negative for blood. Precision and recall metrics are also improved with the inclusion of distant features (see S1 File). Based on these results we included distant features in the overall feature set (designated with the prefix “upper_”).

Feature importance analysis from Boruta selection results (S7 Fig) has shown that all top non-B DNA features (G-quadruplexes, direct repeats, Z-DNA, short tandem repeats, mirror repeats), HDNase, and 5`UTR are presented in both scales (local and distant). For the remaining features (all five TFs and four genomic regions) only distant aggregates ranked in the top list. These observations indicate that the information about broader landscape of the genomic regions is useful for breakpoint hotspots prediction. Moreover, it is especially meaningful when the distant are combined with the local features, because it enables inside-model comparison of densities and identification of anomalies. Feature importance distribution for non-B DNA-based, TF-based, and genomic region-based breast cancer hotspots prediction model (Fig 4) shows that the model considers distant features as the most powerful factors while taking into account local features leads to further gain.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Feature importance analysis inside non-B DNA, TF and genomic region groups.

A. Feature importance of models built on non-B DNA only for breast cancer (99% labeling type). B. Feature importance of models built on TFs only for breast cancer (99% labeling type). С. Feature importance of models built on genomic regions only for breast cancer (99% labeling type).

https://doi.org/10.1371/journal.pcbi.1008749.g004

Additionally, we checked whether binary flags alone or these combined with the indicators of local/global maxima could help to predict cancer breakpoint hotspots (S12B Fig). The results showed that for models using only binary features the quality dropped by 0.13 ROC AUC on average (compared to local features), while the addition of indicators of maxima reduced the difference to 0.03 ROC AUC. This means that in itself the presence o a feature in a genomic area around a hotspot or the presence of a feature together with indicators of its maxima is less informative than the exact feature coverage in given genomic region.

Non-predictable individual breakpoints and highly predictable breakpoint hotspots

Previous attempts to build models using full breakpoint distribution did not achieve good prediction performance (regression task, R2 of 10% for all cancers except for 18% for the breast cancer [9]). The modeling presented above was done for predicting breakpoint hotspots that correspond to genomic regions with recurrent mutations. Here we address the question whether individual breakpoints are just random noise or they also have general mutation determinants that can be captured by a machine learning approach. We performed a series of experiments for the comparison of models predicting breakpoint hotspots, individual breakpoints (S10 Table) and models that distinguish breakpoint hotspots from breakpoints (Fig 5 and S11 Table). The results show that for the majority of cancers the breakpoints are indeed unrecognizable from random locations (ROC AUC around 50% or slightly higher and the mean lift of recall ranging from 0 to 2.5), however for breast, ovary and prostate cancers even individual hotspots can be predicted with ROC AUC 65–75% (nevertheless the lift of recall is very low). In addition, models that can distinguish hotspots from breakpoints can be as good as models predicting hotspots, and reach 85% ROC AUC for breast cancer and more than 70% for brain, liver, pancreas, and prostate confirming that hotspots locations are completely different from individual breakpoints locations by considered features.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Comparison of models predicting individual breakpoints and hotspots.

Distributions of model prediction power are presented for all cancer types using ROC AUC (upper panel) and lift of recall (lower panel). Blue colour signifies models distinguishing between breakpoints and random genomic positions (“breakpoints vs all”), yellow colour signifies models distinguishing between breakpoint hotspots and and random genomic positions, grey colour signifies models distinguishing between individual breakpoints and breakpoint hotspots. The results are presented for 0.03, 0.05 and 0.1 probability percentiles.

https://doi.org/10.1371/journal.pcbi.1008749.g005

Besides it was established that for the majority of cancer types (all except bone, blood and brain) both ROC AUC and the median lift of recall increase with an increase in hotspots labelling threshold meaning that hotspots of higher breakpoints density differ more from other genomic regions than those of lower density (S13 Fig and S12 Table). As an example, decreasing the breakpoint density threshold for hotspot labeling by only 5% leads to 2 times lower mean lift of recall and 15% lower mean ROC AUC for all cancer types on average. Moreover, for considered low labeling thresholds ROC AUC exceeds 60% only for bone cancer while on average for all cancer types ROC AUC equals 54% with mean lift of recall of 1.6.

The results of these experiments are summarized in Table 1 (the median test ROC AUC) and Table 2 (the median test lift of recall for 0.03 probability percentile).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Median test ROC AUC for different models.

https://doi.org/10.1371/journal.pcbi.1008749.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Median test lift of recall for different models for 0.03 probability percentile.

https://doi.org/10.1371/journal.pcbi.1008749.t002

Randomness of breakpoint formation and breakpoint recurrence degree influence overall prediction model quality, and our results indicate that only breakpoints with high recurrence could be well explained by genomic and epigenomic features. Nevertheless, for some cancer types we achieved moderate quality of individual breakpoints prediction, which means that even on an individual level of breakpoints there exist structural mutation signatures underlying these processes, and it would be of interest to further research, distinguish and quantify the contribution of those structural mutation signatures in cancer mutagenesis.

Comparison with the state-of the art approach

The latest comprehensive study that explored the contribution of different factors in cancer mutagenesis with machine learning models was the study of Georgakopoulos-Soare et al. [9], in which the authors investigated contribution of non-B DNA structures and histone modifications, either individually or in combination. As of note, Georgakopoulos-Soare et al. focused on somatic point mutations, and reported the results of the modeling for point mutations, with breakpoint analysis reported in Supplementary Materials. We compared our approach that uses more groups of factors and considers long-distant features to the data analyzed by Georgakopoulos-Soare et al. The results are presented in Fig 6. First we processed the raw breakpoints data (genomic coordinates) provided by Georgakopoulos-Soare et al. with our preprocessing pipeline and get breakpoints density for 500 kb length genome windows. We added all features used in our models (see Methods) and built a Random Forest regression following the model in [9] (prediction of breakpoint densities based on densities of different features) and using only non-B DNA and histone modification groups. We further built models using all of the features, including all long-distant features (Fig 6A). As it was done in [9], we used log2(1+x) for feature densities while we applied log2(mean(x)+x) transformations for the target variable. The model quality was assessed through a 10-fold cross-validation procedure. As a result, even on two feature groups–non-B DNA and histone modifications–the model achieved a better prediction power, partially because the size of the data sets is different since the data are constantly growing in ICGC, and partially because we used all available data on histone modifications and, additionally, the annotations by short stem-loop structures for non-B DNA group (see Methods).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Comparison with the state-of the art approach.

The state of the art is the model in Georgakopoulos-Soare et al. [9]. To enable comparison the designation of the types of cancer is given as in Georgakopoulos-Soare et al. The horizontal line mark the median performance of the model of Georgakopoulos-Soare et al. We built models using different groups of predictors: restricted (non-B and histone groups), full (all groups used in this study), and full and higher-level (all groups considering locations up to 1Mb) A. Comparison of regression models (prediction of breakpoints density). B. Comparison of classification models (prediction of 99% labeling type hotspots).

https://doi.org/10.1371/journal.pcbi.1008749.g006

Next we tested the prediction power of classification tasks that we used in the presented study–by selecting breakpoint hotspots based on the distribution thresholds (Fig 6B). Qualitatively we obtained the same results when the best model was the model that used all feature groups with additionally included distant features. The higher performance of the models trained in our study is explained by the increase in the data accumulated in ICGC since the publication of Georgakopoulos-Soare et al. [9].

The data

To form a target variable we used the data for 628,126 cancer genome breakpoints from 10 cancer types (release 28 from ICGC [36]). Compared to previous analysis [14] performed on an earlier data release, the sample was enriched with more than 165,000 new breakpoints, which became available mostly for prostate cancer (more than 128,000 new breakpoints) (S1 Table). The breakpoint data were pre-processed prior to the analysis. First, the data were filtered so that breakpoints with a high position inaccuracy (a start/end location interval is wider than 10 bases) were removed from the sample (S2 Table). Then the genome was split into non-overlapping windows 100 kb wide, and further analysis was performed at the level of genome window.

To control for genome instability introduced by natural topological constrains some windows were removed. Specifically, we excluded the windows intersecting with centromeres and telomeres, the first and the last window of each chromosome. Additionally, we took into consideration the DAC blacklist (https://www.encodeproject.org/annotations/ENCSR636HFF/) comprising a set of genome regions recognized as low mappable or anomalous and recommended to be removed from genome studies.

Concerning the blacklisted regions, it is worthy of notice that the percent of breakpoints in these regions is very small. So it is minimal for liver cancer with 0.09% of breakpoints located in the regions, maximal for brain cancer (6,8%) and is in the range of 0.4% - 1.5% for other cancer types. Finally, Y and MT chromosomes were not included into the study as they contain relatively small number of breakpoints. Finally, we selected 29288 100-kb windows.

To test whether low mappability regions affect the set of top features determined with machine learning algorithms we filtered breakpoints falling into Duke 35-mer track with mappability less than 0.5 (https://genome.ucsc.edu/cgi-bin/hgFileUi?db=hg19&g=wgEncodeMapability).

For each window located in a specific chromosome, breakpoint density was calculated as a number of breakpoints overlapping with the window divided by the total number of breakpoints in the chromosome separately for each cancer type. Given breakpoint density we highlighted the hotspots–the regions of the genome significantly enriched with breakpoints. The window was designated as a hotspot if its breakpoint density is greater than the specified density distribution percentile (99%, 99.5%, 99.9%) for corresponding cancer type. These three thresholds are referred to as labeling types (labeling thresholds). The total number of breakpoints sorted by cancer type and labelling type is given in S3 Table. Additionally, binary labels were assigned to genome windows indicating the presence/absence of breakpoints.

Genomic features data were downloaded from several projects such as the Encode [18], DNA Punctuation (www.dnapunctuation.org)[17], Non-B DB projects [16], UCSC Genome Browser, and from materials of research papers [19]. The source data for secondary structure annotations for all types except for stem-loops were collected from Non-B DB, while the latter were provided by the DNA punctuation project. Then we used genomic regions markup from UCSC Genome Browser. Topologically associating domain (TAD) data were published in [19] (S3 Table “TAD boundary annotations”). The epigenetic features data were downloaded from the Encode project. We selected experiments with released status, a number of biological replicates greater than one, with no treatment applied. We used TF ChIP-seq for transcription factors,—DNase-seq for DNA accessibility,—Histone ChIP-seq for histone modifications, and RRBS for DNA methylation–.

Feature data were aggregated in chromosome windows by calculating the coverage as the summary length of all regions in a window occupied by a feature (avoiding double account of overlapping regions) divided by the window length, so that the coverage reflects how many bases in a window are covered by the feature. Feature densities were calculated as a percent of feature coverage of a genomic region of 100kB (local features) or 1 Mb (distant features). Based on these densities we obtained binary flags (1, if coverage of the genomic characteristic is greater than 0, and 0 otherwise) and local/global maximum indicators (a binary indicator of the feature value exceeding 90%/ 95%/ 99% percentiles of the feature distribution; a binary indicator of the feature value exceeding values of all 1/5/10 preceding and following adjacent genomic regions; the relative difference between the feature value and the maximum of 10 preceding and following neighbouring regions).

A list of groups and features available for each cancer type is presented in S5 and S6 Tables. The total analysis was performed for 30 datasets (10 cancer types with 3 different target labeling types) comprising hotspot binary labels and genomic feature coverage. The total number of features for each cancer type is given in S6 Table.

Modeling

We used Random Forest [37] model for the machine-learning approach, with hyperparameters (maximal number of terminal nodes, number of trees, minimal number of observations in a terminal node, number of sampled variables for a tree) selected to be optimal when averaging ROC AUC statistics across all cancer types using data for 99% labeling type.

To attain reliable model performance estimates we used a 30-times repeated and stratified train-test split with 30% of data in test sample, with the splits stratified by chromosome and the inside position The quality of models was measured on each of 30 test sets with ROC AUC and lift of recall/lift of precision metrics. These metrics show the number of times recall/precision of the model is greater than a random choice using specific probability percentile as a threshold where each probability percentile is given by the percentile of the model output probability distribution. In particular, the lift of recall for a specific probability threshold (given by a percentile of a probability distribution) is calculated as the recall divided by this percentile of the probability distribution. For example, since it is supposed that in case of a random choice recall approximates the probability percentile (labeling n% of examples gives n% recall), by dividing the model recall by the probability percentile we get the number of times a model is better compared to a random choice.

After training and evaluating hotspot prediction models based on all features (S15 Fig and S18 Table) we performed Boruta feature selection [38] to choose the most important features for hotspot prediction task for each type of cancer. Boruta feature selection is a method specifically designed for Random Forest which for various random initializations iteratively checks whether a considered feature is more important than all of the random (“shadow”) features and updates the feature set for the next iteration by removal of unimportant features. Based on this approach we selected top features for each type of cancer (based on results for 99% labeling type) and evaluated models on these feature sets (S16A Fig and S19 and S20 Tables). It was established that for three cancer types (pancreatic, prostate and breast cancer) the quality of models dropped in comparison with the quality of models on all features. To eliminate the difference we performed forward feature selection and found 1, 1 and 2 features respectively which addition to feature sets lead to comparable quality (S16B Fig and S21 Table). These final features sets for each type of cancer are referred to as “all” in Fig 3.

Feature group importance analysis was done as follows. We built models separately on features of each feature group for each cancer type. To obtain feature group performance rankings we fixed the maximum achievable mean lift of recall at 0.03 probability percentile for each labeling type and scaled the values of the mean lift of recall for all feature groups to that fixed value and then averaged for the cancer type for 99% and 99.5% labeling types (Fig 1).

For the analysis of the degree of randomness of cancer breakpoints we performed a set of experiments. First, we built models to predict hotspots with increasing hotspot labeling thresholds to 75%/90%/95% of breakpoint density distributions. Secondly, we estimated the quality of models discriminating hotspots from breakpoints and models discriminating breakpoints from the remaining genomic areas to compare it to performance of hotspot prediction models. In both cases for modeling we used all available features. At the end we compared all performance metrics and identified cancer types for which random nature of breakpoints formation with predictable recurrent locations is the best visible (S17 Fig).