3DGenBench: a web-server to benchmark computational models for 3D Genomics

Abstract Modeling 3D genome organisation has been booming in the last years thanks to the availability of experimental datasets of genomic contacts. However, the field is currently missing the standardisation of methods and metrics to compare predictions and experiments. We present 3DGenBench, a web server available at https://inc-cost.eu/benchmarking/, that allows benchmarking computational models of 3D Genomics. The benchmark is performed using a manually curated dataset of 39 capture Hi-C profiles in wild type and genome-edited mouse cells, and five genome-wide Hi-C profiles in human, mouse, and Drosophila cells. 3DGenBench performs two kinds of analysis, each supplied with a specific scoring module that compares predictions of a computational method to experimental data using several metrics. With 3DGenBench, the user obtains model performance scores, allowing an unbiased comparison with other models. 3DGenBench aims to become a reference web server to test new 3D genomics models and is conceived as an evolving platform where new types of analysis will be implemented in the future.


INTRODUCTION
Understanding the mechanisms driving chromatin folding is among the most challenging questions in modern genomics. Robust experimental techniques allow nowadays to obtain various datasets describing chromosome organisation and dynamics quantitatively. The plethora of genomewide experimental data provides opportunities to infer relationships between 1D (epi)genomic and 3D genome folding properties, as well as the underlying molecular mechanisms, using statistical analysis (1,2) or biophysical modelling (3,4).
Statistical methods are being actively developed, and besides exploration of interdependencies between genomic datasets, dozens of recent models allow predictions of chromatin architecture in unseen cell types, species, and genetic conditions. Although each study contains certain benchmarks, a direct comparison of developed models is difficult due to differences in the model's input information, statistical engine, design of validation experiments and metrics characterising the model performance.
In contrast to the statistical approaches which are agnostic to the mechanisms underlying the observed dependencies, biophysical models can be used to (in)validate specific biological mechanisms or reconstruct 3D genome structures from data (3,5,6). In particular, biophysical models have predictive power and can be employed to understand the consequences of (epi)genetic perturbations (7,8) without the need to perform time-consuming and laborious experiments. However, as for statistical inference, standardised benchmarks and metrics are missing to evaluate and compare model performances and to quantify similarities between predictions and experiments.
To sum up, although multiple modeling approaches have been developed, the fundamental differences among the methodologies and validation datasets render the comparison challenging. To address this challenge, the International Nucleome Consortium (https://inc-cost.eu/) has recently started a benchmarking initiative, providing tools for the uniform assessment of computational models for 3D-Genomics. This initiative includes the collection and processing of a reference dataset based on previously published experimental data and the implementation of metrics modules scoring the similarity between experimental data and modeling results. The central component of this work is a web server, 3DGenBench, allowing users to submit their predictions and access performance metrics.

Datasets collection
3DGenBench performs two kinds of analysis. The first type aims at quantifying how good a method predicts one experimental dataset. We defined this analysis as horizontal benchmarking. The second analysis, defined as vertical benchmarking aims to estimate how well a method can predict changes in chromosome folding caused by structural perturbations of the genome (see results text and Figure 1A for details of each analysis type). For each of these analyses, we collected specific datasets.
The data for the vertical benchmark were collected from nine studies from 2016 to 2019 (7,(9)(10)(11)(12)(13)(14)(15) describing genomic mutations (deletions, inversions, duplications) of specific regions associated with chromatin architecture (Figure 1B,  Supplementary Table S1). We collected and processed capture Hi-C (cHi-C) reads for wild-type (WT) and mutated conditions. Since CTCF is particularly important for chromatin architecture in vertebrates and several examples in our dataset include mutations of CTCF binding sites, we collected CTCF ChIP-seq data for relevant cell types. For this aim, we either reprocessed CTCF ChIP-seq reads obtained from the same study as cHi-C reads or, if not available, we searched for CTCF ChIP-seq experiments on matching cell type from other studies in PubMed and EN-CODE datasets.
Data for the horizontal benchmarking include GM12878 and K562 human cell lines (16), mouse embryonic stem cell (mESC) line (17), and one Drosophila cell line Kc167 (18). We chose five regions for prediction for each cell type (Supplementary Table S2). Since epigenetics data are important as features for several modelling tools, we composed a table  (Supplementary Table S3) with download links to the most widely used formats for several epigenetics marks for each cell type.
In addition, we used CTCF ChIA-PET dataset for human GM12878 cell line available as GEO accession GSM1872886 (19) for horizontal ChIA-PET benchmarking.
ChIA-PET data processing. PET clusters file was downloaded from GSM1872886, single interaction for mitochondrial chromosome was removed, the remaining data were used without any further processing.
ChIP-seq data processing. All CTCF ChIP-seq data were processed using ENCODE ChIP-seq pipeline v1 (https: //github.com/kundajelab/chipseq pipeline) with parameters default for transcription factors analysis mode. Links to raw data and reference genomes are provided in Supplementary  Table S1.

Performance metrics
Horizontal benchmark metrics.
• Stratum adjusted correlation coefficient (SCC) metric. SCC is a special measure to quantify the similarity between Hi-C matrices. This algorithm stratifies all Hi-C contacts according to their genomic distance and computes a Pearson correlation coefficient for each stratum. We used SCC metric (23) implemented via the hicreppy package (https://github.com/cmdoret/hicreppy). SCC algorithm has a tuning parameter h controlling the smoothing level. Smoothing the Hi-C matrix minimizes the effect of noise and biases. We used h = 2 for all comparisons. • Insulation score (IS) metric. IS metric reflects the segmentation of the genome into domains. We used calculate insulation score function from coolltools.api.insulation package (version 0.5.0) with window size parameter equals to 5× bin size (resolution) to compute IS for all bins in the experimental and predicted data. We then computed Spearman's correlation coefficient between IS tracks obtained from experimental and predicted data. • Compartment strength (CS) metric. CS track for experimental and predicted data was calculated as in (24). CS for each bin (b) quantifies as the ratio between the average OoE (observed over expected) value of b with bins of the same compartment and the average OoE value of b with any other bin in the same chromosome. We then used Pearson's correlation coefficient between obtained tracks as the resulting metric.
W6 Nucleic Acids Research, 2022, Vol. 50, Web Server issue • Contacts scaling metric. We calculated the average number of contacts at separate genomic distances for experimental and predicted data. The obtained arrays of values were compared using Pearson's correlation coefficient.
• The percentage of captured ChIA-PET interactions (i.e. true-positive predictions); • The percentage of ChIA-PET interactions predicted by the model but not observed experimentally (i.e. falsepositive predictions); • The percentage of ChIA-PET observed experimentally but not predicted by the model (i.e. false-negative predictions); • The mean-squared error of PET counts: for N overlapping PET clusters, each characterized by observed count C Exp i and predicted count C Pr ed i , we computed MSE as Vertical benchmark metrics.
• Insulation score changes metric. We used calculate insulation score function from coolltools.api.insulation package (version 0.5.0) to compute IS for all bins in the WT and mutation (Mut) conditions for experimental and predicted data. Then the Mut IS was divided by WT IS both for experimental and predicted data. Finally, we used Pearson's correlation coefficient to compare Mut/Wt experimental and Mut/Wt predicted ratio tracks.
• Ectopic interactions area under the curve (AUC) metric. Ectopic interactions for experimental and predicted data were calculated as in (7). Normalization coefficients computation accounts for the effects of the copy number changes caused by structural variants. Specifically, I. We normalized Mut matrix by coverage excluding contacts inside of deleted or duplicated regions. Let M be a Hi-C matrix for Mut condition with elements m i j , and W be a matrix for WT condition with elements w i j . Than we define normalized matrix for Mut condition as:M In case of duplication, we computed an additional normalization coefficient as: wherem i j are elements ofM. Next, we calculated normalized mutated matrix as: where D coe f = 1 for deleted and inverted regions.
These normalization steps were performed both for experimental and predicted data. II. At the second stage, we subtracted the W matrix from theM matrix. III. Then we normalized the obtained differences matrix, dividing each entry by the average value at the corresponding diagonal. IV. Finally we calculated Z-scores at each diagonal. We filtered out elements above the 96th percentile of all diagonal values in the calculation of standard deviations.
Ectopic interactions were defined as interactions exceeding a Z-score of 2 for experimental data. These ectopic interactions for the experimental data were used as true labels for precision recall curve function from sckit-learn python library. Matrix with standard deviations of predicted data was a second argument of precision recall curve function. Different thresholds of standard deviation defined predicted ectopic interactions, precision, and recall. AUC was calculated using sckit-learn python library as well.

Generation of baselines
To provide baselines for comparison, we simulated fake predictions with varying quality either by adding different amounts of noise into experimentally observed data, or randomly by permuting observed contacts on each diagonal of a Hi-C matrix. Metrics were then computed between these baselines and the original experimental datasets.
Noise perturbation was implemented by adding a random value to each interaction of the experimental contact matrix. This random value was sampled from a normal distribution with the mean equal to the observed contact frequency and standard deviation ranging from 0.001 (the lowest noise level) to 50 (the highest noise level). We used nine control models with standard deviations equal to 0.001, 0.01, 0.01, 0.5, 1, 2, 10, 20, 50.
Diagonal shuffling was performed by grouping contacts of loci with the same genomic distance, and by randomly shuffling the observed contact values within each group. The Hi-C matrices obtained after this shuffling procedure were used as a baseline.

Generation of case study datasets
To demonstrate the metrics performance implemented in the horizontal benchmark challenge, we chose two genomic regions for GM12878 cell line at 10 Kb resolution. We used predictions generated by the DRAGON (25) and PRISMR (7) models as an example of computational predictions of chromatin architecture. DRAGON uses polymer simulations to predict Hi-C contact matrices and threedimensional (3D) chromosome structures at five kilobase resolution from 12 key epigenetic marks and the genomic locations and orientations of CTCF binding sites. It was parameterized with a maximum entropy optimization algorithm using Hi-C data from GM12878 cells. PRISMR modelling was performed at 10 kb resolution as described before (7,26). In brief, PRISMR infers, through a recursive procedure based on polymer physics, the minimal polymer model that, at equilibrium, best describes the input contact matrix.
Nucleic Acids Research, 2022, Vol. 50, Web Server issue W7 The inferred set of binding sites of the model correlates with key epigenetic factors and defines a code linking epigenetics to architecture, as detailed in (26). By using this code, the contact matrices of entire chromosomes can be predicted from only epigenetics (26).
To demonstrate the metrics performance for the vertical benchmark, we simulated predictions of WT and Mut data for three rearrangements (structural variants ids according to the Supplementary Table S1: deletion MusDel-B, inversion Inv1, and duplication dup-C) using 3DPredictor (27) and 3DPolyS-Fit (28). 3DPolyS-Fit infers pair-wise interaction strengths of a self-interacting polymer model from WT data and then use them to predict the MUT data by simulating the corresponding structural rearrangement as in (29).

Web site development
The web site was developed with the following tools: Visualization was performed with HiGlass JS API (30).

Benchmarking datasets and challenges
We implemented two types of analysis in the 3DGenBench benchmark platform, which we defined as horizontal and vertical benchmarking ( Figure 1A). The horizontal benchmarking challenge aims at quantifying the predictive power of computational models that predict Hi-C data from locusspecific 1D-epigenetic data. For this challenge, we chose four high-resolution Hi-C datasets describing well-studied human, mouse, and Drosophila cell lines (human K562 and GM12868; mouse embryonic stem cells; Drosophila Kc167 cells). In addition, we selected CTCF ChIA-PET data for human GM12868 cells to allow scoring ChIA-PET data prediction.
To standardise the 1D-epigenetic information, we provided a list of available epigenetic tracks for each of the cell lines included in the benchmarking, such as ChIP-seq profiles, transcriptome data, and other genomic features. We limited the challenge to six 20-Mb long regions from different chromosomes, which should allow most of the models to participate in the benchmarking. Model parameters can be learned from whole genome data from which the benchmarking validation datasets are excluded. The submission to the platform consists of one predicted Hi-C map for one of the six regions of the benchmarking dataset. ChIA-PET data submission consists of a single file of PET clusters.
We implemented several metrics measuring model performance in the horizontal benchmarking challenge (see Materials and Methods). For Hi-C data, these included Spearman's correlation and stratum adjusted correlation coefficient (SCC) (23) to measure overall similarity of experimen- tal and predicted Hi-C maps. SCC is advantageous over correlation because it is not biased by the strong distance dependence of chromatin interaction data. We measured the accuracy of the predicted distance dependence separately as the correlation between observed and predicted contact scaling functions. We used Pearson's correlation between predicted and observed compartmentalization profiles (24) to show how well the model captures preferences in interactions between chromatin compartments. Finally, we used Spearman's correlation of insulation scores to estimate how the model predicts boundaries of topologically associating domains (TADs). Since there are multiple models predicting TAD boundaries without inferring contact maps, we provided an option to submit insulation scores tracks only. To support visual assessment of the predicted data, we integrated a HiGlass view in the 3DGenBench output pages.
For ChIA-PET data prediction the overlap between predicted and reference lists of PET clusters are estimated (see Materials and Methods) and the mean-squared error of PET counts is computed to estimate the accuracy of cluster counts.
The second type of analysis implemented in the 3DGen-Bench platform is called vertical challenge. It aims at quantifying how accurate models are in predicting changes of chromatin architecture caused by genetic perturbations. For this challenge, we collected capture Hi-C-datasets available for CRISPR-edited mouse cells and tissues with matched wild-type controls. The resulting dataset includes 49 wildtype--mutation pairs from nine loci and nine cell types ( Figure 1B). The input data for this type of analysis consists of a pair of predicted Hi-C datasets matching experimental wild-type and mutation data. For each wildtype--mutation pair, 3DGenBench will extract the fraction of chromatin interactions that differ between both conditions. 3DGenBench then compares this predicted set of ectopic interactions with experimentally observed differences.
We implemented several previously-published metrics, specifically measuring how the model captures differences in chromatin architecture between two conditions (see Materials and Methods). These metrics score the overlap between observed and predicted ectopic chromatin interactions (7) and predictions of TAD boundaries changes.

3DGenBench workflow, input and output description
We expect the following usage of the 3DGenBench services. First, modellers can explore the 3DGenBench datasets available under the 'Datasets' tab to find relevant experimental data for benchmarking (Figure 2A). Each Hi-C dataset is available in the most commonly used formats: hic, cool (5, 10, 20 and 50 kb resolutions), and pairs. Also, there are supplementary tracks describing CTCF binding sites for most of the datasets. All the data can be downloaded via hyperlinks in the table. Alternatively, the complete metadata table can be obtained from the server in JSON or text format to fetch and process download URLs automatically.
After exploring and modeling experimental data, predictions for the corresponding datasets can be uploaded to the server via FTP or the submission page ( Figure 2B). Predicted Hi-C contacts should be formatted as 4-column files, each row containing the chromosome name, two genomic coordinates, and contact count (see Supplementary Data 1 for a representative example); ChIA-PET predictions as 6-column files, with the same format as the pairs format utilized by 4DNucleome portal (https://www.4dnucleome. org/). Example files are available at the 3DGenBench tutorial pages. Users can upload multiple files that will appear in the drop-down list at the next (submission) step.
Once the data is uploaded, users can follow the 'Compute metrics' tab from the home menu, choose the type of prediction (single for horizontal benchmarking or paired for vertical benchmarking) and fill all the fields related to the predicted sample ( Figure 2C).
The status of the submission is then available at the 'Submissions list' tab of the home menu ( Figure 2D). Blue status of submission indicates that the job is queued; yellow status of submission indicates the server is currently being processed the job; green status means that metrics computation was successfully completed, and red indicates an error with data processing.
When data processing is over, clicking on the ID link (e.g. bmXXXXXXXX) redirects the user to the page with the job results ( Figure 2E). This page contains metrics describing the prediction accuracy of the model. In addition, the result pages shows metrics obtained for baselines, simulating predictions of varying quality (see Materials and Methods for details). It is also possible to compare with previouslysubmitted jobs by choosing any available IDs with a green colour of the status icon.

Case study
To evaluate how 3DGenBench scores the similarity between Hi-C datasets, we generated fake predictions by injecting different amounts of noise into experimental Hi-C W10 Nucleic Acids Research, 2022, Vol. 50, Web Server issue data or by randomly shuffling contacts on each diagonal of Hi-C matrices (see Materials and Methods) and submitted all generated datasets to the 3DGenBench platform. We also generated predictions for contact maps of two loci in GM12878 cells based on 1D epigenomic data (horizontal challenge) using the DRAGON and PRISMR chromatin models (7,25,26), and predicted three pairs of WT and Mut contact maps (vertical challenge) using 3DPredictor, and 3DPolyS-Fit algorithms (25,27,28).
As shown in Figure 3A for horizontal benchmarking, the 3DGenBench metrics are concordant with the level of noise injection. Noisier models have lower SCC, compartmentalization and insulation score similarity metrics. Interestingly, metrics show differ sensitivity to noise injections. The correlation coefficient computed from all contacts decreases even when a small amount of noise was added to contacts, whereas compartments strength and insulation scores were relatively stable at these noise levels. Diagonal shuffling expectedly preserve contacts scaling, but result in strong changes of the insulator score.
According to our metrics, the DRAGON model perfectly captures contact scaling, shows good metrics of compartment prediction, but is less efficient regarding the insulation landscape ( Figure 3A and C). PRISMR model captures insulation landscape significantly better than DRAGON model, but is not as good in predicting contacts scaling. This toy example shows how 3DGenBench can be used to discover strengths and weaknesses of different predictive models.
Similarly, 3DGenBench correctly scores noise injection models in vertical benchmarking ( Figure 3B). For example, precision-recall AUC gradually decreases from 1 to ∼0.2 when increasing the amount of noise. Our comparison of 3DPredictor and 3DPolyS-Fit algorithms shows that in most of cases 3DPolyS-Fit reconstructs chromatin architecture and its changes better than 3DPredictor, although in few cases (i.e. for prediction of insulation score changes in the inversion Inv1) 3DPredictor outperforms 3DPolyS-Fit. Overall, this shows that the 3DGenBench metrics module provides a user-friendly interface to score computational predictions of genome architecture.

CONCLUSIONS AND FUTURE PROSPECTS
Here, we describe 3DGenBench, a web-platform for benchmarking of computational models of 3D chromatin architecture. There are two types of analysis implemented in 3DGenBench currently, but this platform can be easily extended to include additional species, cell types and benchmarks. It can be used by model developers to improve their methods based on robust metrics but also by model users to choose the best methods to tackle their questions. One promising direction of 3DGenBench development is the integration with 4DNucleome datasets (31), which is currently the largest available collection of 3D genomic data. In addition, different types of data could be integrated, e.g. from microscopy based experimental approaches (32), to test the models at the single molecule level (33) beyond the population averaged contact maps. Another future direction is the development of online services implementing the most popular algorithms for chromatin modelling.
Currently, the vast majority of computational models do not have user-friendly implementations, and only few of them have online-tools allowing simple access (27,34). A web-implementation of these modeling tools alongside with benchmarking results provided by 3DGenBench platform will allow the 4D Genomics community to easily integrate chromatin models in their research.