The ABC model defines the score of a potential E/G in a cell type as the product of the activity of the potential enhancer E in this cell type, and the contact between E and gene G, divided by the sum of the same products but across all potential enhancers in a 5 Mb region from G. The ABC model starts by defining candidate regulatory regions E, as regions of open chromatin (defined by either DNAse-seq or ATAC-seq) in a cell type. It then quantifies the enhancer activity (A) of these regions E by computing the geometric mean of the read counts of chromatin accessibility (usually assessed using DNAse-seq or ATAC-seq) and H3K27ac ChIP-seq at E. The contact (C) between E and G is then computed either as the Knight-Ruiz (KR) matrix-balancing normalized Hi-C contact frequency between E and the promoter of gene G, if cell type specific Hi-C data are available, or simply as the inverse of the distance (fractal globule model) between E and G otherwise (Fulco et al. (2019)). In order to predict E/G only for expressed genes, the ABC model can either take cell type specific gene expression data in, or consider as a proxy of gene expression, the activity of its promoter as defined above using chromatin accessibility and H3K27ac ChIP-seq data ² .

The Average-Rank method defines the score of a potential E/G as the inverse of the average of the ranks provided by the Sheffield and the distance methods (Moore et al. (2020a)). The Sheffield method was introduced in 2013 and defines the score of a potential E/G as the Pearson correlation between the logarithm of the chromatin accessibility at E (assessed by DNAse-seq) and the logarithm of the expression of G across many cell types (Sheffield et al. (2013)). The distance method scores a potential E/G as the inverse of the distance between E and G. Here potential enhancers are all distal enhancer elements (distal enhancer like signature elements or dELS) of the ENCODE registry of candidate cis-regulatory elements (cCREs) (Moore et al. (2020b)).

TargetFinder defines true (ground positive) E/G based on 3D data (HiC) and learns features associated to those using gradient boosting. The learnt features are as diverse as open chromatin, methylation, histone marks or transcription factors, and can both be taken from enhancer and promoter regions and from the window between them (Whalen et al. (2016)). Indeed its authors showed that features located in enhancer-promoter windows (EPW) are also predictive of true E/G relationships and should be incorporated in the model.

The two first link methods mentioned above also propose their own reference/evaluation datasets.

Fulco et al. (2019)’s reference set is based on previous CRISPR-based experiments performed in K562 cells and on the output of a new genetic screening technique developed by the authors, called CRISPRi-FlowFISH. This technique was specifically designed to predict E/G in a cell type for a given small number of genes. As stated by the authors, it perturbs “hundreds of non-coding elements in parallel and quantifies their effects on the expression of an RNA of interest, combining CRISPR interference, RNA fluorescence in situ hybridization (FISH) and flow cytometry”. In this approach, they “deliver KRAB-dCas9 to many candidate regulatory elements in a population of cells by using a library of guide RNAs”. The results of this technique are then subjected to a statistical framework to determine the sets of E/G that are active and inactive in the cell type. The technique was then applied to thirty genes in five genomic regions (spanning 1.1–4.0 Mb) for which they tested all DNase I hypersensitive (DHS) elements (representing open chromatin regions) in K562 cells within 450 kb of the gene of interest. Together with previous CRISPR experiments, this approach yielded 109 ground positives (i.e., “positive in the evaluation set”) and 3,754 ground negatives (i.e., “negative in the evaluation set”) E/G, which are considered as a reference set for the evaluation of numerous methods of the field including the ABC model and the distance method (Fulco et al. (2019)). We will use this reference set here and call it CRiFF (Table 1). Note that the 30 selected genes had an RPKM expression level above 20 in K562, and that some of them were erythroid-specific while others were ubiquitous. No filtering on chromatin accessibility level was applied to open chromatin regions, however the sequences of the probes that were designed to target open chromatin regions through gRNAs had to be specific enough. Fulco et al. also garantee a 5% FDR to detect E/G and more than 80% power to detect a 25% effect on gene expression with CRiFF. Finally it has to be noted that the 109 ground positives and the 3,754 ground negatives of CRiFF exclude repressive elements and promoter-promoter interactions (interactions where the targeted element is located less than 500bp away from a TSS).

To complement CRiFF which is rather small, we decided to use another recent genetic screening set that differs from CRiFF in being enhancer-centric instead of gene-centric. It was also available on K562 and could be retrieved from Moore et al. (2020a)’s reference set. It is made of 651 ground positive and 24,576 ground negative E/G relationships, and we will call it CRISPRi (Table 1).

Moore et al.‘s reference set is entitled BENGI (Benchmark of candidate Enhancer-Gene Interactions) and is made of sets of E/G active and inactive in different cell lines according to different types of data (3D, genetic). We focus our evaluation on the GM12878 cell line, which has the largest amount of annotation data, with 6 sets of active and inactive E/G available. The active E/G sets result from the processing of four types of 3D data, Hi-C (Rao et al. (2014)) and promoter capture Hi-C (Mifsud et al. (2015)) data and ChiA-PET of polymerase II and CTCF (Tang et al. (2015)) data, and of eQTL data from two different studies, GEUVADIS (Lappalainen et al. (2013)) and GTEx (Consortium et al. (2015)). The sets of ground negatives are built by taking, for each enhancer of a positive set, all the genes not connected to it in the positive set and lying within the 95 percentile of the positive set distances from it. The number of ground positive and negative E/G obtained are indicated in Table 1. Since 3D and eQTL data are not specifically generated to identify E/G relationships, the BENGI reference sets are expected to be overall less reliable than the CRiFF and the CRISPRi reference sets. However, the fact that BENGI provides genome-wide information is an advantage over CRiFF.

Given all these data we proceeded to the evaluation of the ABC model, the Average-Rank method, the distance method, and TargetFinder, on all three reference datasets: BENGI, CRiFF and CRISPRi.

For the ABC model, the Average-Rank method and the distance method, we used the code provided by the authors (Fulco et al. (2019); Moore et al. (2020a)), with some adjustments, while for the last one we downloaded the predictions provided by the authors (Whalen et al. (2016)) ³ . The obtained results are presented on Figure 3 for BENGI, and on Figure 4 for CRiFF and CRISPRi (see Material and methods). Note that while Fulco et al. provided a code associated to their proposed method ⁴ , Moore et al. only provided a code for the evaluation of their proposed method on a given evaluation set, which is less generic ⁵ . In addition, the fact that we used TargetFinder’s already thresholded predictions only allowed us to compute a single pair of (precision, recall) values for each reference set, and explains the absence of AUPR curves for this tool. In fact, two different pairs of (precision, recall) values, a pessimistic one and an optimistic one, were computed and plotted for TargetFinder, leading to two different dots for TargetFinder on the plots, TargetFinder_pes and TargetFinder_opt (see Material and methods for details).

We replicated the results of the two evaluation papers (Fulco et al. (2019); Moore et al. (2020a)). The curves and AUPRs (Area Under the Precision-Recall curve) of the Average-Rank and the distance methods of Figure 3 are in agreement with Figure S2 of Moore et al. (2020a) derived for GM12878 cell line (all pairs, natural ratio). Similarly, the curves and AUPRs of the ABC model and the distance method of Figure 4A agree with Figure 3A of Fulco et al. (2019). Like Moore et al. (2020a), we also found that TargetFinder performs better than Average-Rank, except on eQTL reference sets. Altogether, these positive controls confirmed the validity of the pipeline we implemented.

Figure 3 further shows low AUPR values for the first three methods on the six BENGI datasets, and that TargetFinder performs best overall, followed by Average-Rank, distance and ABC model. Note that TargetFinder is much more precise than sensitive, and performs much better on HiC and CTCF sets. This last result can be explained by the fact that TargetFinder learns true E/G based on HiC data. Nevertheless, Figure 4A shows larger AUPRs for the three last methods on the CRiFF set, and that the ABC model performs best (AUPR = 0.63), before the distance method and finally the Average-Rank method. Contrary to its result on BENGI, TargetFinder does not perform well on CRiFF.

When comparing the performances of the methods on CRISPRi with respect to CRiFF, if the methods that perform best (ABC model) and worst (TargetFinder) are the same, it has to be noted that Average-Rank performs better than distance and that the AUPRs are globally lower (Figure 4B). It is also important to note that on half of the evaluation sets used here, the number of TargetFinder’s validated E/G relationships is lower or equal to 9 (Table 2), showing that this tool is not easily applicable and might only deliver predictions for small subsets of references.

Therefore, state-of-the-art E/G identification methods do not perform very well overall, and using 3D or genetic screening data as reference provides completely opposite answers to the question of the best performing E/G identification method.

In order to illustrate unsupervised/heuristic enhancer/gene identification methods, we chose the simplest one, the pairwise chromatin accessibility correlation across cell types, and represented it on Figure 2. The actions that led to this figure are the following: we first downloaded the ENCODE uniformly processed read alignments (bam files) of DNAse-seq data (single end) from 10 cell types: stomach, HepG2, K562, thymus, adrenal gland, small intestine, GM12878, IMR-90, heart and H1-hESC, with accession numbers provided in Table 3. We then called the chromatin accessibility peaks from the mapped reads in each cell type using macs2 ⁷ (Figure 2; Table 3).

We obtained from about 60,000 (GM12878) to about 200,000 (IMR-90) peaks per cell type. By concatenating, sorting and merging on the genome the peaks called in each cell type using bedtools merge, we then obtained 473,766 consensus peaks across all cell types. We then quantified the chromatin accessibility of the 473,766 consensus peaks in each cell type by simply counting the number of mapped reads of each cell type overlapping each consensus peak using bedtools intersect, and normalized the number of reads of each peak in each cell type by the total number of mapped reads in peaks for this cell type. Finally we computed the consensus peak pairwise Pearson correlation between the log10 of the normalized chromatin accessibility across the 10 cell types of these peaks for all pairs of peaks less distant than 500 kb using a script that we wrote: compute_correlations.py ⁸ . We then only considered as E/G relationships, the pairs of peaks with a correlation above 0.7 and for which one of the two peaks overlapped the most 5’ bp (TSS) of a Gencode v19 ⁹ gene (vertical green rectangles on Figure 2).

In addition to evaluating the ABC model on BENGI and the Average-Rank and the TargetFinder methods on CRiFF, and since Moore et al. provided the code of the Average-Rank method, and Fulco et al. the code to run the ABC model, we decided to try and reproduce the evaluation of the Average-Rank method on BENGI and of the ABC model on CRiFF. We also used the code of the distance method provided by Moore et al. (with some modifications), to evaluate the baseline distance method on BENGI, CRiFF and CRISPRi (Tables 1, 2; Figures 3, 4).

In total we evaluated four methods, the ABC model, the Average-Rank, the distance and the TargetFinder methods on three references sets, BENGI, CRiFF and CRISPRi (Figures 3, 4). It has to be noted that contrary to the other methods, TargetFinder’s predictions were downloaded directly from its authors’s website ¹⁰ , therefore only allowing us to compute a single pair of (precision, recall) values, and not AUPR curves. In fact, we used two different ways to compute TargetFinder’s recall, an optimistic and a pessimistic one, which led to two different dots for this tool in the evaluation plots (see sections below about TargetFinder). In addition, since the code to generate the Precision-Recall curves and the AUPRs was not provided in the papers, we generated our own R code to make these plots using existing R packages. The code used to perform all these analyses was stored in Jupyter notebooks that we provide below, together with additional details about these analyses.