Using a K-mer Based Approach with Machine Learning Classifiers for Enhancer Identification and Classification

doi:10.21203/rs.3.rs-4836518/v1

Download PDF

Research Article

Using a K-mer Based Approach with Machine Learning Classifiers for Enhancer Identification and Classification

https://doi.org/10.21203/rs.3.rs-4836518/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

Enhancers are genetic elements that modulate the level of gene expression in cells; because they are essential for cellular function, enhancer dysfunction is associated with many complex diseases including many types of cancers. Furthermore, current research highlights the challenge of enhancer identification due to variability in enhancer definitions and identification approaches. Developing a machine learning pipeline to distinguish enhancers from other DNA elements would greatly aid the ability to study enhancers and their role in disease. In this project, we developed a random forest machine learning model to distinguish between human liver enhancer sequences with low versus high levels of support across eight different enhancer identification methods.

Results

Enhancers were classified as “shared” if they had support from multiple methods and “unique” otherwise; the threshold for the number of methods that made an enhancer region “shared” or “unique” was adjusted as an independent variable. We also extracted a set of random sequences from the human genome to evaluate potential differences between shared/unique enhancer regions and genomic background. Our machine learning model used the frequency of k length sequence (k-mer) patterns as features to classify the shared versus unique enhancer regions. We hypothesized that the distribution of k-mer frequencies would correlate with transcription factor binding sites on the enhancer and distinguish them from other types of DNA sequence.

Conclusions

The final shared vs random models perform with 93.8–99.9% accuracy, the final unique vs random models perform with 93.4–99.8% accuracy, but the final shared model performs with only 56.8–61.4% accuracy. These results indicate that the k-mer frequency distributions for shared and unique enhancers are not significantly different although shared/unique sequences differ from the genomic background.

gene regulation

genomics

machine learning

enhancers

The central dogma of molecular biology states that DNA is transcribed into RNA, which is translated into protein. Regulatory elements known as enhancers can increase the rate at which DNA is converted into RNA, and ultimately alter the amount of protein that is produced [1]. Regular cellular fluctuations in protein requirements make enhancers essential for day-to-day cellular function, and enhancer dysfunction is involved in a variety of diseases such as obesity, Alzheimer’s disease, diabetes, and cancer [2]. Given their importance, enhancers have been a central focus of many biological studies.

Enhancers interact with transcription factors and cofactors to regulate gene expression in a distance- and orientation-independent manner. They are thought to be brought near their target gene(s) through three-dimensional chromatin looping [3]. Due to the variability for the timing of this interaction and the distance between enhancers and their target genes, identification of enhancers within the human genome remains a significant challenge. Several experimental assays and computational approaches have been used to identify enhancers within the genome, with some success, but no single “gold-standard” method exists [1], [4], [5].

To increase confidence in enhancer sets, researchers have combined results from various enhancer identification methods. Previous work hypothesized that enhancers identified by multiple methods – shared enhancers – are more likely to be true enhancer regions than unique enhancers [6]. In this study, we evaluate this hypothesis by comparing shared and unique enhancers from eight common enhancer identification methods using a machine learning pipeline. Each of these methods generates a set of enhancers using computational models, experimental techniques, or a combination of both. We also compare the shared and unique enhancers to a set of length-matched, randomly generated sequences to identify differences between putative regulatory sequences and the genomic background. Because we expect active enhancers to contain transcription factor binding sites, we trained our model using DNA sequence patterns as features. We hypothesized that there would be differences in sequence patterns between both enhancer and non-enhancer regions as well as between shared and uniquely predicted sequences.

We downloaded reference enhancer data for our analysis from the seven methods detailed in Benton et al. [1], focusing on liver enhancer predictions. We added an additional set of regions comprised of ENCODE Project candidate cis-regulatory elements (cCREs) identified in the liver [7]. Each of the eight total methods for enhancer identification represents a different collection of putative enhancers in the GRCh37/hg19 build of the human genome. All methods use data derived from either chromatin immunoprecipitation assays followed by sequencing (ChIP-seq) and DNase I hypersensitive site sequencing (DNAse-seq) from the ENCODE Project and Roadmap Epigenomics consortia, or cap analysis of gene expression (CAGE) from the FANTOM consortium. The data were integrated using computational pipelines to identify genome-wide sets of enhancer annotations, as previously described [1]. We performed data preprocessing with AWK (v4.2.1) and BEDtools (v2.30.0) on a RHEL8 Linux server [8]. Our downstream analysis and machine learning models were implemented in Python 3.6 using publicly available libraries, including sci-kit learn (v1.2.2) [9].

Identifying shared and unique enhancers

To reduce bias during feature extraction for our model, we trimmed or lengthened all putative enhancer regions to a standard length of 500 bp around the midpoint [10], [11], [12]. Although there is not a single length for enhancer sequences, previous work suggests that 500 bp is a reasonable approximation [13]. We then used the BEDtools mutiIntersect command to annotate the enhancers from all eight methods with the amount of overlap between methods. The resulting database organized genomic locations by the number of overlapping methods for each sequence (from 1, no overlap, to 8, overlap between all methods). We used a variable threshold to separate the shared enhancer regions (highly overlapping with other methods) from the unique enhancer regions (enhancers with low cross-method support). One potential confounder is that the BEDtools multiIntersect algorithm can produce enhancer fragments of varying sizes. To combat the discrepancy, we discarded overlapping regions that were less than 50 bp and extended/trimmed all other regions to 500 bp. This step also ensured that our frequency data would not be skewed by very short overlaps.

Generating Random Genomic Regions

We generated a set of random enhancer sequences to evaluate the difference between the shared and unique enhancer regions and genomic background. The additional random dataset allowed for comparison to genomic regions that were excluded from the 8 enhancer sets referenced above. We generated the datasetusing the BEDtools random algorithm. Each random region was specified to be 500 bp in length, and we generated a total number of random regions equivalent to the number of ‘shared’ enhancers.

Creating the Frequency Distribution

We hypothesized that the shared enhancers would have a characteristic frequency of k-length base pair sequences (k-mers) that could distinguish them from both the unique enhancers and randomly generated sequences. The k-mers produced from the sequences could serve as a proxy for transcription factor (TF) binding sites, where TFs serve as enhancer activators [3]. A schematic for this intuition is presented in Fig. 1. We chose two possible values of k: k = 3 and k = 4. Mathematically, there are 64 (4³^3) possible way to arrange the nucleotide bases ‘A’, ‘G’, ‘T’, and ‘C’ in groups of 3. Similarly, there are 256 (4) ways to arrange the nucleotide bases in groups of 4. A k of 2 or lower would yield an insufficient 16 (4²) features while a k of 5 or greater would yield a computationally expensive 1024 (4⁵) features. To create our frequency distribution, we iterated through each enhancer region k base pairs at a time and incremented the count for the matching k-mer. We used a sliding window approach for grouping k-mers to avoid missing a reading frame. For example, for the sequence “ACGT,” the 3-mer frequency counts for ACG and CGT would both be incremented. The frequency counts were then divided by the total number of k-mers in each enhancer, which is given by n-k-1, where n is the length of the enhancer region in base pairs and k is the length of the k-mer in base pairs. For a 3-mer frequency distribution of 500 bp length enhancers, this results in 498 k-mers within each region. The overall process of generating k-mer frequencies is detailed in Fig. 2. We use the two resulting frequency distributions (for k = 3 and k = 4) in the subsequent models.

Learning the Frequency Distributions with Random Forests

Using our newly generated frequency distributions, we applied machine learning to classify our sequences based on patterns in the data. Every enhancer – shared, unique, and random – was considered an observation in our dataset with their 64 or 256 k-mer frequencies as features. To compare model pattern recognition for unique and random enhancers, we designed two sets of models: a first set that classified shared versus unique enhancers and a second one that classified shared versus random enhancers. We selected a random forest (RF) model to carry out the pattern discovery due to its success in similar contexts [8]-[10].

The RF classifier is an ensemble learning model that creates several decision trees for different random sub-samples of the dataset [4]. Decision trees are flowcharts that divide the sub-sample into progressively smaller groups while minimizing variability for the predictor variable [14]. Randomly selecting the sub-samples protects the model from overfitting. Each decision tree has a vote on what category any new point should be assigned to; the categorical variable with the highest number of votes represents the overall prediction of the model. The RF classifier is a popular ML model utilized in the field of bioinformatics[4], [14], [15] due to its inherent resistance to over-fitting and non-parametric structure [15].

Training the random forest model

An important consideration for our feature set was that the number of unique enhancer regions represented an overwhelming majority of the data. This was expected given previous work that showed that enhancer annotations produced by identification strategies have relatively little overlap [1]. We mitigated this by down sampling the unique dataset to improve class balance. To avoid bias from a single down-sample of the unique set, we generated 5 different random down-samplings of the unique set for each shared set. The pairs of shared and down-sampled unique sets were used to train five versions of the shared vs unique and unique vs random model. This method minimized the error that would be present from one poorly representative down sampled set. We evaluated our model performance using the five resulting models.

Random forest models are trained to distinguish between shared enhancers and genomic background and shared versus unique enhancers

Once the random forest models had been trained, we assessed their accuracy on the held-out test set. For all conditions, we set a minimum threshold for the overlap required to be considered “shared” and maximum threshold for the overlap allowed to be considered “unique”. For example, a maximum unique threshold of 2 would require that no more than 2 methods overlap that a region before it was no longer “unique” to a method; a minimum shared threshold of 4 means that overlap between at least 4 methods is required to consider that region “shared” across methods. Model accuracies and cross validation scores organized by threshold, k-mer, and classification are presented in Table 1.

Average accuracy suggests strong ability for models to distinguish between enhancers and genomic background but not between enhancer groups

Across all five-model iterations, the shared vs random and unique vs random results produced significantly higher average accuracy and cross validation scores than the shared vs unique results. In Table 1, the highly desirable accuracy/cross validation scores of 0.90 or greater are depicted in green. The shared vs random and unique vs random accuracy and cross validation scores are depicted entirely in this category, indicating strong ability for the random forest algorithm to discern between groups. This is further bolstered by the fact that results were consistent for several iterations of manipulated threshold and k values (minimum average accuracy across all shared/random and unique/random sets is 0.938). However, the shared vs unique average accuracy and cross validation scores are entirely orange or red, indicating scores of less than 0.70 and greater than 0.25 and less than 0.25 respectively. The poorer scores indicate the random forest algorithm’s relative inability to distinguish between shared/unique enhancers.

Model performance is significantly different between the Shared/Unique and both the Shared/Random and Unique/Random models

We compared the cross-validation scores of shared/unique with each of the shared/random and unique/random models using the non-parametric Mann-Whitney-U test. The p-value between the shared/unique vs shared/random is \(\:3.11*{10}^{-7}\:\)and the p-value between shared/unique and unique/random is \(\:3.18*{10}^{-7}\). Both indicate strongly significant differences in model performance when classifying between enhancers and genomic background and enhancers themselves.

Enhancer identification is an important problem in the field of genomics. Previous work has been instrumental in cataloging and exploring the many complementary approaches for enhancer prediction and validation; however, it remains challenging to prioritize functional sequences [1], [6]. In this project, we sought to evaluate whether predictions made by multiple enhancer identification strategies, or shared enhancers, had systematically different sequences than those that were identified by a small number of approaches. We hypothesized that sequences backed by a large amount of data would correspond to higher confidence, functional sequences, while those that were identified by fewer methods were more likely to be false positives.

We trained several random forest classification models to distinguish between shared versus unique enhancer sequences, and to distinguish between predicted enhancer sequences and sequences from the genomic background. We found that the shared vs unique models performed with significantly worse accuracy than either the shared vs random or the unique vs random models. We suspect that this is because both the shared and unique enhancer regions are identified by at least identification strategy, suggesting that they have at least some characteristics of functional sequences. The random enhancer regions were not identified by any other approach, thus were less likely to overlap with regulatory sequence patterns. However, the poor performance of the shared versus unique classifier suggests that the commonly used method of overlapping regions to identify true enhancers may not have a significant benefit.

A critical step of our analysis was the standardization of all enhancer sequences to 500 bp, a chosen standard that lies between the usual 200–1000 bp length of enhancers [16]. Prior to this addition, our preliminary RF models utilized initially performed with a high accuracy (within the 90–95% range). However, an important confounder had been inflating the model performance: the sequence length. Extending shorter sequences and trimming longer sequences to the 500 bp standard allowed for us to negate the confounding length factor. However, limitations exist in the interpretation of our results as no other enhancer lengths were tested. Trimming and extending sequences was done about the midpoint and with reference to the genome to capture the core “signal” present that could aid the random forest model in identifying enhancer sequences. However, it is possible that the core signal was not captured upon trimming. Furthermore, the signal could be diminished by extending enhancers to the 500 bp mark. Although our model results are thorough for a 500 bp enhancer standard, further testing is required to determine whether our results hold true for enhancer base pairs between the usual 200–1000 bp range.

Though the model can well distinguish between random and actual enhancer regions, further analysis must be conducted to see what factors contribute to different classifications. Our aim is to utilize a set of software tools known as the MEME-suite. Specifically, SEA is a software that utilizes a database of known TFBSs and finds enrichment for each of them in a set of sequences [17]. Our stock set of sequences is the overall set of enhancer regions, both shared and unique. By running these regions through SEA, we can find enriched TFBSs within our enhancer regions. The overall k-mer frequencies for the enriched TFBSs can be calculated with our above algorithm to determine what distribution of 3-mers or 4-mers is highest within enhancer regions. Furthermore, having this analysis done will be instrumental in determining the TFBSs enriched within enhancers specific to liver cells. Malfunctioning liver enhancers can then be pinpointed using the TFBS frequency distributions and studied as a target for genetic therapy.

Within the genomics community, it is assumed that regions that are identified by multiple enhancer identification methods is more likely to be an actual enhancer.

However, our results show that for the 500 bp standard and the utilized random forest pipeline it is likely that there is not a significant difference in sequence composition between shared enhancers that are identified by numerous methods and unique enhancers identified by a few. Therefore, given that enhancers are notorious and often difficult to place within the genome, it is likely beneficial to classify enhancers determined by any of the above accredited eight methods as a true enhancer region. Practicing a more inclusive selection process for true enhancers could be beneficial to avoid omission of potential regulatory regions responsible for disease from available enhancer databases.

Competing interests:

The authors declare that they have no competing interests.

Funding:

Not applicable.

Author Contribution

M.L.B. supervised the study. C.P. and M.L.B. analyzed the data and performed all statistical analyses. C.P. generated all figures and tables. C.P. and M.L.B. wrote the manuscript. All authors read and approved the final manuscript.

Acknowledgement

The authors would like to thank the members of the Benton Lab for constructive conversations and feedback during this project. Figure 1, Figure 2a, and Figure 2b were created with BioRender.com

Data Availability

All data supporting the findings of this study are publicly available. Enhancer annotations were downloaded from Benton et al. 2019 and are available in the associated creDB (https://www.kostkalab.net/software.html). Additional liver enhancer annotations were derived from the SCREEN registry of cCREs by the ENCODE project (version 3).

Benton ML, Talipineni SC, Kostka D, Capra JA. Genome-wide enhancer annotations differ significantly in genomic distribution, evolution, and function. BMC Genomics. 2019;20(1):511. 10.1186/s12864-019-5779-x.
Claringbould A, Zaugg JB. Enhancers in disease: molecular basis and emerging treatment strategies, Trends Mol Med, vol. 27, no. 11, pp. 1060–1073, Nov. 2021, 10.1016/J.MOLMED.2021.07.012
Uyehara CM, Apostolou E. 3D enhancer-promoter interactions and multi-connected hubs: Organizational principles and functional roles, 10.1016/j.celrep.2023.112068
Butt AH, Alkhalifah T, Alturise F, Khan YD. A machine learning technique for identifying DNA enhancer regions utilizing CIS-regulatory element patterns. Sci Rep. 2022;12(1):15183. 10.1038/s41598-022-19099-3.
Blinka S, Reimer MH, Pulakanti K, Pinello L, Yuan G-C, Rao S. Identification of Transcribed Enhancers by Genome-wide Chromatin Immunoprecipitation Sequencing, 10.1007/978-1-4939-4035-6_8
Gasperini M, Tome JM, Shendure J. Towards a comprehensive catalogue of validated and target-linked human enhancers, Nature Reviews Genetics 2020 21:5, vol. 21, no. 5, pp. 292–310, Jan. 2020, 10.1038/s41576-019-0209-0
Dunham I, et al. An integrated encyclopedia of DNA elements in the human genome. Nature. 2012;489(7414):57–74. 10.1038/nature11247.
Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. Mar. 2010;26(6):841. 10.1093/BIOINFORMATICS/BTQ033.
Pedregosa F. FABIANPEDREGOSA Scikit-learn: Machine Learning in Python Gaël Varoquaux Bertrand Thirion Vincent Dubourg Alexandre Passos PEDREGOSA, VAROQUAUX, GRAMFORT ET AL. Matthieu Perrot, Journal of Machine Learning Research, vol. 12, pp. 2825–2830, 2011, Accessed: Apr. 18, 2024. [Online]. Available: http://scikit-learn.sourceforge.net
Colbran LL, Chen L, Capra JA. Short DNA sequence patterns accurately identify broadly active human enhancers. BMC Genomics. 2017;18(1):536. 10.1186/s12864-017-3934-9.
Colbran LL, Chen L, Capra JA. Sequence Characteristics Distinguish Transcribed Enhancers from Promoters and Predict Their Breadth of Activity, Genetics, vol. 211, no. 4, pp. 1205–1217, Apr. 2019, 10.1534/genetics.118.301895
Fish A, Chen L, Capra JA. Gene Regulatory Enhancers with Evolutionarily Conserved Activity Are More Pleiotropic than Those with Species-Specific Activity, Genome Biol Evol, vol. 9, no. 10, pp. 2615–2625, Oct. 2017, 10.1093/gbe/evx194
Rao S, Ahmad K, Ramachandran S. Cooperative binding between distant transcription factors is a hallmark of active enhancers, Mol Cell, vol. 81, no. 8, pp. 1651–1665.e4, Apr. 2021, 10.1016/J.MOLCEL.2021.02.014
Tyryshkina A, Coraor N, Nekrutenko A. Predicting runtimes of bioinformatics tools based on historical data: five years of Galaxy usage, Bioinformatics, vol. 35, no. 18, pp. 3453–3460, Sep. 2019, 10.1093/bioinformatics/btz054
Rajagopal N, RFECS. Mar., : A Random-Forest Based Algorithm for Enhancer Identification from Chromatin State, PLoS Comput Biol, vol. 9, no. 3, pp. e1002968-, 2013, [Online]. Available: https://doi.org/10.1371/journal.pcbi.1002968
Swarr DT, Wert SE, Whitsett JA. Molecular Determinants of Lung Morphogenesis, Kendig’s Disorders of the Respiratory Tract in Children, pp. 26–39.e4, 2019, 10.1016/B978-0-323-44887-1.00002-X
Bailey TL, Grant CE. SEA: Simple Enrichment Analysis of motifs, bioRxiv, p. 2021.08.23.457422, Jan. 2021, 10.1101/2021.08.23.457422
Moore JE, et al. Expanded encyclopaedias of DNA elements in the human and mouse genomes. Brenton R Graveley. 2020;583:25. 10.1038/s41586-020-2493-4.

Table 1 is available in the Supplementary Files section.

No competing interests reported.

Table1.jpg

Download PDF

Version 1

posted

You are reading this latest preprint version

Using a K-mer Based Approach with Machine Learning Classifiers for Enhancer Identification and Classification

Status:

Version 1

Abstract

Figures

Background

Methods

Identifying shared and unique enhancers

Generating Random Genomic Regions

Creating the Frequency Distribution

Learning the Frequency Distributions with Random Forests

Training the random forest model

Results

Model performance is significantly different between the Shared/Unique and both the Shared/Random and Unique/Random models

Discussion

Declarations

Competing interests:

Funding:

Author Contribution

Acknowledgement

Data Availability

References

Table 1

Additional Declarations

Supplementary Files

Status:

Version 1