A systematic evaluation of Hi-C data enhancement methods for enhancing PLAC-seq and HiChIP data

Overview of the evaluation framework

In this study, we mainly applied three existing methods (HiCNN2, HiCPlus and DeepHiC) designed for Hi-C data to enhance the sequencing depth of HP data (Figure 1). We also explored HiCNN and VEHiCLE but chose not to include them for most assessments because HiCNN has a highly similar performance as HiCNN2 (Supplementary Figures S1 and S2, see Supplementary Data available online at https://academic.oup.com/bib); and VEHiCLE’s specific features tailored for Hi-C data render it suboptimal for HP data (Supplementary Figure S3, see Supplementary Data available online at https://academic.oup.com/bib). First, we generated low-depth HP (‘baseline’) datasets from mESC H3K4me3 PLAC-seq data and GM12878 Smc1a HiChIP data with different downsampling ratios (Data preprocessing in Methods section). We similarly generated low-depth mESC CTCF PLAC-seq data and H3K4me3 PLAC-seq data for four human brain cell types with downsampling ratio of 0.125. Next, we split each dataset into training and testing datasets with the training dataset consisting of chromosomes 1, 2, 3, 5, 7 and 9 (chromosome 2 was used as the validation data, as part of the training procedure, to select the best model); and the testing dataset containing all the other chromosomes (chromosomes 4, 6, 8, 10–19 for mESC or chromosomes 4, 6, 8, 10–22 for human cell types). Then on the training datasets, we applied each method to train models using the low-depth (i.e. baseline) input data and the high-depth (i.e. full data without downsampling) target data and subsequently applied the trained models to the low-depth testing data to obtain an enhanced high-depth data (Figure 1). Finally, we calculated the similarity between the enhanced HP data and the high-depth HP data (i.e. full data without downsampling for the testing chromosomes, which serves as the working truth). Specifically, we assessed similarity using four metrics: Pearson’s correlation coefficient, Spearman’s rank correlation coefficient, Brownian distance covariance (or distance correlation) [20] and HPrep, a new method to assess the reproducibility of HP data [21]. For presentation brevity, we have decided to only show the Pearson correlation coefficient results in the main text as the other statistics reach qualitatively the same conclusions. In addition, we performed 3D peak calling before and after enhancement to assess the impact of enhancement on the detection of chromatin interactions.

Next, we compared the performance of each method using the model trained on HP data to that trained on Hi-C data (detailed in later section Hi-C or HP data for training). Lastly, we evaluated the transferability of each method by enhancing HP data across different cell types detailed in the later section Model transferability).

Since HP data measure protein-mediated chromatin interactions, we evaluated enhancement results only for bin pairs where at least one bin contains the protein of interest. Specifically, we defined bin pairs where both bins contain the protein of interest as the ‘AND’ set, and bin pairs where only one bin contains the protein of interest as the ‘XOR’ set, following our previous work [18, 21]. We removed bin pairs where neither of the two bins contains the protein of interest from our downstream analysis, and only applied abovementioned similarity metrics (Pearson’s correlation, Spearman’s correlation, Brownian distance and HPrep) to bin pairs in the ‘AND’ and ‘XOR’ sets. Noticeably, in HP data, bin pairs in the ‘AND’ set usually show higher contact frequency than bin pairs in the ‘XOR’ set due to double ChIP enrichment. We thus evaluated similarity for the ‘AND’ and ‘XOR’ sets separately.

Performance comparison of different methods

We benchmarked the performance of the three methods (HiCPlus, HiCNN2 and DeepHiC) in terms of enhancing low-depth HP data at 10 Kb resolution. Note that we decided not to include HiCNN for most evaluations because HiCNN and HiCNN2 perform highly similarly (Supplementary Figures S1 and S2, see Supplementary Data available online at https://academic.oup.com/bib). For HiCNN2, which is an ensemble method with three models (HiCNN2-1, HiCNN2-2 and HiCNN2-3), we present only HiCNN2-1 for the rest of the manuscript because the three HiCNN2 models perform almost indistinguishably (Supplementary Figures S1 and S2, see Supplementary Data available online at https://academic.oup.com/bib). Evaluations of the three methods (namely, HiCPlus, HiCNN2-1 and DeepHiC) suggest that they perform reasonably, all significantly outperforming the low-depth HP data (Figure 2, Supplementary Figures S4 and S5, see Supplementary Data available online at https://academic.oup.com/bib). For example, when using HiCNN2-1 to enhance the GM12878 HiChIP data by 25x (i.e. from 0.04 depth to full), Pearson correlation coefficients for the ‘AND’ set are 0.70–0.81, which are 0.09–0.24 higher than the low-depth data, when the genomic distance is 20–250 Kb (Figure 2A). Similarly, when using HiCPlus to enhance the mESC PLAC-seq data by 25x, Pearson correlation coefficients for the ‘XOR’ set are 0.44–63, 0.16–0.21 higher than the low-depth data, when the genomic distance is 50–500 Kb (Supplementary Figure S5, see Supplementary Data available online at https://academic.oup.com/bib). VEHiCLE shows inferior performance when enhancing GM12878 HiChIP data (Supplementary Figure S3, see Supplementary Data available online at https://academic.oup.com/bib) possibly due to certain features tailored for Hi-C data that are no longer suitable for HP data. For instance, VEHiCLE performs Knight-Ruiz (KR) normalization for Hi-C data. However, because HP technologies enrich chromatin contacts at the region mediated by the protein of interest, the equal visibility assumption made in KR normalization is invalid for HP data. In addition, VEHiCLE requires all diagonal bin pairs to have nonzero counts. Therefore, we decided not to pursue further with VEHiCLE enhancement.

In addition, the three methods perform similarly for most of the seven HP datasets evaluated, with Pearson correlation differences largely within a difference of 0.1 (Figures 2 and 3, Supplementary Figures S4 and S5, see Supplementary Data available online at https://academic.oup.com/bib). For example, when downsampling ratio is 1/25 and the distance is 20 Kb–1.25 Mb for the GM12878 HiChIP data, HiCNN2-1 improves Pearson correlation by 0.024–0.048 and 0.012–0.059 on the ‘XOR’ set, compared with HiCPlus and DeepHiC, respectively (Figure 2D). For another example, when downsampling ratio is 1/16 and the distance is 250 Kb–1.5 Mb for the mESC PLAC-seq data, HiCNN2-1 and HiCPlus show highly similar performance and improve Pearson correlation by 0.01–0.09 on the ‘AND’ set, compared to DeepHiC (Supplementary Figure S5B, see Supplementary Data available online at https://academic.oup.com/bib). When enhancing some cell types, for instance, radial glia (RG) and intermediate progenitor cells (IPCs), DeepHiC is substantially worse than HiCNN2-1 and HiCPlus with a difference of 0.2 in Pearson correlation (Figure 3C, D and H, I). The inferior performance of DeepHiC may be due to mode collapse issues [22, 23] for GAN models (more details in the Discussion section).

3D peak calling

3D peak calling, or the detection of statistically significant long-range chromatin interactions, is one of the important downstream analyses for various types of chromatin conformation data, including HP data. To evaluate the impact of HP data enhancement, we further applied our model-based analysis of PLAC-seq and HiChIP (MAPS) pipeline [18] to identify significant chromatin interactions before and after enhancement, and compared them with chromatin interactions detected from the full data. We treated 3D peak calling results derived from the full data (without any downsampling) as the truth. Specifically, we defined true peaks as bin pairs with MAPS false discovery rate (FDR) <1%, contacts ≥12, and signal to noise ratio (SNR, i.e. the ratio of observed count over expected count) ≥2; and we defined true background bin pairs as those with MAPS FDR >10% and contacts ≥12. We found that even high-depth input HP data (e.g. 0.5 down-sampled GM12878 HiChIP data in Figure 4 or 0.5 downsampled mESC PLAC-seq data in Supplementary Figure S6, see Supplementary Data available online at https://academic.oup.com/bib, where the raw sequencing depth is ~322 million and ~568 million, respectively) benefit from enhancement in that enhanced datasets can improve the power of 3D peak calling. For example, for 0.5 downsampled GM12878 HiChIP data (Figure 4), the baseline (i.e. downsampled data before enhancement) has a low sensitivity of 0.32, while the enhanced data (using HiCPlus or HiCNN2) improve the sensitivity to 0.71–0.72, while maintaining the desired FDR 1%. DeepHiC increases sensitivity even more drastically but fails to maintain the desired 1% FDR. Similarly, we observed clear improvement with enhanced datasets for 0.5 downsampled mESC PLAC-seq data (Supplementary Figure S6, see Supplementary Data available online at https://academic.oup.com/bib). Observing that the FDRs from baseline is essentially 0, we relaxed the MAPS-FDR threshold to 0.2 for the baseline, which led to an actual FDR of 0.02, comparable to that after enhancement. With the relaxed FDR threshold, the power of the baseline increased substantially, from 0.32 to 0.6 but still clearly lower than 0.71 after enhancement (Figure 4). Similar patterns were observed for 0.5 downsampled mESC data (Supplementary Figure S6, see Supplementary Data available online at https://academic.oup.com/bib).

We additionally examined SNR (again defined as the ratio of observed count over expected count) before and after enhancement, both compared to SNRs from the full data without enhancement. We found that SNR estimates from baseline data without enhancement are significantly lower than those from full data (Supplementary Figure S7B, see Supplementary Data available online at https://academic.oup.com/bib). Treating the estimates from full data as the working truth, these results indicate that baseline data tend to underestimate the magnitude of 3D peaks. Data enhancement mitigates the underestimation issue, with enhanced data producing SNR estimates more closely approaching the working truth (Supplementary Figure S7C, see Supplementary Data available online at https://academic.oup.com/bib). Although we observed a significant difference in SNR estimates both at 3D peaks (Supplementary Figure S7E, see Supplementary Data available online at https://academic.oup.com/bib) and at background bin pairs (Supplementary Figure S7F, see Supplementary Data available online at https://academic.oup.com/bib), we noticed that the absolute difference is more pronounced among 3D peaks. Specifically, mean and median SNR at 3D peaks are 4.43 and 4.07 at baseline, 4.90 and 4.37 after HiCNN2-1 enhancement and 4.99 and 4.38 when using the full data (Supplementary Figure S7E, see Supplementary Data available online at https://academic.oup.com/bib). In contrast, the mean and median SNR at background bin pairs are 0.97 and 0.94, 0.98 and 0.97, and 1.01 and 0.99, respectively, with only ≤0.05 absolute difference. The statistical significance at background bin pairs is driven primarily by the huge number of background bin pairs (Supplementary Figure S7F, see Supplementary Data available online at https://academic.oup.com/bib).

Encouraged by the power improvement in 3D peak calling genome-wide, we proceeded to examine two specific loci in mESCs where previous studies [24, 25] have established enhancer–promoter interactions. These two loci are Med13l and Mtnr1a loci (Figure 5). From Figure 5, we observe that baseline without enhancement fails to identify many 3D peaks, including the most important EPIs. After HP data enhancement, we were able to rescue some of the EPIs. For example, for the two bin pairs corresponding to EPIs at the Med13l locus (illustrated with black arrows), full data identified both; baseline identified only one, while every enhanced data were able to rescue the missed one (Figure 5A–E). Similarly, for the bin pair corresponding to EPI at the Mtnr1a locus (illustrated with black arrows), full data identified it; baseline failed to detect it, while again every enhanced data managed to rescue the signal (Figure 5G–K). Interestingly, simply multiplying the baseline matrix with a constant of 8 can rescue the EPIs at the Med13l locus (Figure 5F), suggesting that simple amplification of the contact frequency matrix may help 3D peak calling when the input data are of low depth. However, this simple strategy still fails to detect the EPI at the Mtnr1a locus (Figure 5L), showcasing the advantage of data enhancement.

Finally, we assessed whether the identified chromatin interactions relate to gene expression. As shown in Supplementary Figure S8 (see Supplementary Data available online at https://academic.oup.com/bib), we found that genes with promoters involving 3D peaks show significantly higher expression levels than genes whose promoters do not involve any 3D peaks. The fact that genes with 3D peaks identified from baseline data without any enhancement are expected as they tend to be the lower hanging fruits with stronger magnitude of chromatin interactions that can be detected by low-depth data.

Hi-C or HP data for training?

We then evaluated the robustness and relative performance of HP depth enhancement for each method when training models with different assays, specifically Hi-C or HP. For enhancing the GM12878 HiChIP data, HiCPlus, HiCNN2 and DeepHiC, all showed comparable or improved enhancement by using models trained on HP data than trained on Hi-C data (Figure 6; Supplementary Figures S9–S11, see Supplementary Data available online at https://academic.oup.com/bib). For example, when enhancing the GM12878 HiChIP data by 8X (i.e. enhancing 1/8 downsampled data to full data), Pearson correlation coefficients using the HiCPlus model trained on HP data improved by up to 0.11 for the ‘AND’ set and 0.04 for the ‘XOR’ set, compared with models trained on Hi-C data (distance: 250 Kb–2 Mb, Figure 6A and D). Similarly, HiCNN2 and DeepHiC models trained on HP data demonstrated overall improved or comparable performance than those trained on Hi-C data, with more obvious improvement than HiCPlus. For example, within 250 Kb–2 Mb distance, HiCNN2 (Figure 6B and E) and DeepHiC (Figure 6C and F) improve Pearson correlation coefficient by 0.15 and 0.18 for the ‘AND’ set and 0.08 and 0.11 for the ‘XOR’ set, compared with the aforementioned 0.11 and 0.04 for HiCPlus (Figure 6A and D).

However, when enhancing the mESC PLAC-seq data, we observed mixed results using models trained on HP data versus those trained on Hi-C data, (Supplementary Figures S12–S14, see Supplementary Data available online at https://academic.oup.com/bib). Specifically, HiCPlus showed similar performance using HP data (light yellow) or Hi-C data (dark yellow) for training (Supplementary Figures S12–S14left panels, see Supplementary Data available online at https://academic.oup.com/bib); HiCNN2-1’s HP trained models (light red) outperformed its Hi-C-trained models (dark red) (Supplementary Figures S12–S14middle panels, see Supplementary Data available online at https://academic.oup.com/bib), while DeepHiC’s HP-trained models (light blue) were inferior to its Hi-C-trained models (dark blue) (Supplementary Figures S12–S14right panels, see Supplementary Data available online at https://academic.oup.com/bib).

One possible explanation for DeepHiC’s better performance of mESC Hi-C-trained models is the much higher sequencing depth of mESC Hi-C data relative to mESC PLAC-seq data. Specifically, mESC Hi-C data is 4.59x that of mESC PLAC-seq data, in terms of informative reads (Table 1). Such drastic depth difference could render the models trained on Hi-C data more advantageous than those trained on HP data for all three methods (Supplementary Figures S12–S14, see Supplementary Data available online at https://academic.oup.com/bib). Particularly in DeepHiC, we observed an obvious advantage of Hi-C-trained models over HP-trained models. To reduce the impact of the different sequencing depths, we downsampled mESC Hi-C data so that its informative reads are comparable to those in mESC PLAC-seq data (Generating HiC_Downsampled data in Methods section). After downsampling, we obtained 59.2 million informative reads in downsampled Hi-C data (HiC_downsampled), matching that (also 59.2 million) in PLAC-seq data (Table 1). We then retrained DeepHiC models using the downsampled Hi-C data. With a comparable number of informative reads, models trained on the downsampled Hi-C data showed worse or comparable performance than those trained on HP data (Supplementary Figure S15, see Supplementary Data available online at https://academic.oup.com/bib). Although worse performance is expected, the magnitude of performance impairment is drastic. For example, when enhancing by 8x (middle panel) within a distance 500 Kb–1 Mb, the Pearson correlation is 0.22–0.37 with DeepHiC models trained on the downsampled Hi-C data, compared with 0.50–0.67 with models trained on HP data and 0.55–0.76 with models trained on full Hi-C data. These results suggest that the DeepHiC method is more sensitive to the sequencing depth of Hi-C data than HiCNN2 and HiCPlus.

Model transferability

Although many HP datasets have been generated recently [16–19], deeply sequenced HP datasets are only available to limited cell types, making it infeasible to train models separately for each cell type. One potential solution is to use models pretrained on available datasets from other cell type(s).

For enhancing the GM12878 HiChIP data, HiCPlus performed similarly with either GM12878-trained or mESC-trained models (the left panel of Figure 7; the left panels of Supplementary Figures S16–S18, see Supplementary Data available online at https://academic.oup.com/bib). Comparatively, HiCNN2-1 and DeepHiC showed slightly higher or higher accuracy using GM12878-trained models than mESC-trained models (the middle and right panels of Figure 7), with the difference more obvious for DeepHiC (the right panel of Figure 7).

For enhancing the mESC PLAC-seq data, the performance of HiCPlus, HiCNN2 and DeepHiC using models trained on the GM12878 HiChIP data was comparable to or even slightly better than using models trained on the mESC PLAC-seq data (Supplementary Figures S19–S21, see Supplementary Data available online at https://academic.oup.com/bib). Specifically, for HiCPlus and HiCNN2-1, the two sets of models were nearly indistinguishable in terms of enhancing the mESC PLAC-seq data (left and middle panels of Supplementary Figures S19–S21, see Supplementary Data available online at https://academic.oup.com/bib). Interestingly, DeepHiC achieved even slightly better performance when using models trained on the GM12878 HiChIP data (right panel of Supplementary Figures S19–S21, see Supplementary Data available online at https://academic.oup.com/bib). One plausible reason is that the models for all three methods are originally developed and fine-tuned for the GM12878 Hi-C data.

Encouraged by the promising transferability results between mESC H3K4me3 PLAC-seq data and GM12878 Smc1a HiChIP data, we proceeded with transferability assessment across more cell types. Since HiCNN2-1 and HiCPlus achieved similarly best transferability performance, we presented only HiCNN2-1 results for brevity. Specifically, we enhanced 0.125 downsampled HP data for each of the six cell types [namely, GM12878 and mESC, and four human brain cell types (19]) including RG, IPCs, excitatory neurons (eN) and interneurons (IN)] using models trained from the corresponding cell type as well as using models trained from each of the other five cell types. For training and testing, we used the same chromosome splitting as illustrated in Figure 1. Results shown in Figure 8 and Supplementary Figures S22 and S23 (see Supplementary Data available online at https://academic.oup.com/bib) further support that the models learned are transferable across cell types. Specifically, Pearson correlation coefficients are almost indistinguishable when enhanced with models trained from the matching cell type or from other cell types (Figure 8 and Supplementary Figure S22, see Supplementary Data available online at https://academic.oup.com/bib), and all the models lead to similar performance in 3D peak calling (Supplementary Figure S23, see Supplementary Data available online at https://academic.oup.com/bib). For example, to detect chromatin interactions in IN, 0.125 downsampled PLAC-seq data (before any enhancement) had essentially no power at all (sensitivity to detect IN of IN-specific 3D peaks is 0.00, left most bars labeled ‘Baseline’ in Supplementary Figure S23C and G, see Supplementary Data available online at https://academic.oup.com/bib); in contrast, enhanced RG data using models trained with IN data resulted in a sensitivity of 0.49 (or 0.48) for IN (or IN-specific) 3D peaks (magenta bars in Supplementary Figure S23C and G, see Supplementary Data available online at https://academic.oup.com/bib); similarly and importantly, enhanced IN data using models trained with data in any of the other five cell types resulted in a comparable sensitivity of 0.47–0.59 (or 0.43–0.63) for IN (or IN-specific) 3D peaks (blue bars in Supplementary Figure S23C and G, see Supplementary Data available online at https://academic.oup.com/bib).

We further evaluated transferability in terms of capturing cell type-specific features, examining gene expression and open chromatin status in the corresponding cell types. Specifically, we compared the distribution of gene expressions for three groups of genes: (i) genes with 3D peak(s) looping to their promoters identified at ‘baseline’ (low-depth data without enhancement); (ii) genes without any 3D peaks at ‘baseline’ but with 3D peak(s) after enhancement, separately for enhanced data using models trained with each of the six cell types and (iii) genes without any 3D peaks even with the full data (‘background’). Not surprisingly, as shown in Supplementary Figure S24A–D (see Supplementary Data available online at https://academic.oup.com/bib), ‘baseline’ identified only few lower hanging fruits, and thus expression levels are the highest; genes with 3D peaks identified only after enhancement [whether using models trained with the matching cell type (yellow boxplots) or different cell types (nonbaseline and nonbackground cyan boxplots)], reassuringly, had only slightly lower expression levels, drastically higher than those ‘background’ genes. Similar patterns are observed when restricting only to cell type-specifically expressed genes (Supplementary Figure S24E–H, see Supplementary Data available online at https://academic.oup.com/bib). Following similar logic, we assessed 3D peaks in terms of their overlap with cell type-specific assay for transposase-accessible chromatin using sequencing (ATAC-seq) peaks, observing similar patterns (Supplementary Figure S24I–L, see Supplementary Data available online at https://academic.oup.com/bib). In particular, models trained with matching cell types (yellow bars) resulted in a similar proportion of overlap with cell type-specific ATAC-seq peaks as those trained with different cell types (nonbaseline and nonbackground cyan bars), suggesting that models trained from unmatching cell types can similarly retain cell type-specific features.

Finally, we explored model transferability across different proteins of interest by cross-applying models learned from H3K4me3 and CTCF PLAC-seq data in mESC. We used the same mESC CTCF PLAC-seq data as in Juric et al. [18]. We observed almost indistinguishable performance when using models trained with the same protein of interest or the other protein (Supplementary Figure S25, see Supplementary Data available online at https://academic.oup.com/bib), both visibly better than without enhancement (i.e. baseline). These results suggest that the enhancement models learned are likely transferable also across different proteins of interest, with the caveat that our assessments only involved three different proteins: Smc1a above for GM12878, CTCF here for mESC and H3K4me3. In the future, more high-depth HP data with different proteins of interest will allow us to perform a more comprehensive assessment across various proteins.

Model robustness

Throughout the manuscript so far, we have used chromosomes chr1, 3, 5, 7 and 9 as training; chromosome 2 as the validation (part of the training procedure to select the best model); and the remaining chromosomes as testing. In addition, when creating low-depth input data, we performed downsampling only once. We evaluated model robustness by swapping training and testing, by using leaving-one-chromosome-out and by performing downsampling five times. Results presented in Supplementary Figures S26 and S27 (see Supplementary Data available online at https://academic.oup.com/bib) show that the models trained are robust, resulting in highly similar Pearson correlation coefficient decay profiles.

Data preprocessing

All our assessed deep learning methods require training data, testing data and validation data as input. In our study, for the mESC PLAC-seq data, we assigned chromosomes 1, 3, 5, 7 and 9 as the training data, chromosome 2 as the validation data and chromosomes 4, 6, 8, 10–19 as the testing data. For the GM12878 HiChIP data, we assigned chromosomes 1, 3, 5, 7 and 9 as the training data, chromosome 2 as the validation data and chromosomes 4, 6, 8, 10–22 as the testing data. Here, validation data was used to select the best model (details in The principle of deep learning in Methods section). Both mESC PLAC-seq data and GM12878 HiChIP data consist of two parts: high-depth data and low-depth data, referring to the original/full HP data without downsampling and downsampled data, respectively. We used the low-depth data as the input for each deep learning method, and the high-depth data to calculate the loss function (details in The principle of deep learning in Methods section).

Specifically, we applied the following steps to generate low-depth data:

• (i) Converting read pairs into bin pairs. We followed our previous study [18, 21] to preprocess the HP data to retain only long-range read pairs (intrachromosomal contacts >1 Kb). We then randomly selected a subset of the read pairs with downsampling ratios 1/25, 1/16 or 1/8. Downsampling was implemented using command ‘SAMtools views ratio’ [28]. Next, we binned the original read pairs and downsampled read pairs into 10 Kb bin pairs, resulting in high-depth data and low-depth data, respectively.

• (ii) Converting high-depth and low-depth bin pairs into contact matrices. For each chromosome, based on whether the bins containing the protein of interest [H3K4me3 ChIP-seq peaks for mESC PLAC-seq data (18) and Smc1a ChIP-seq peaks for GM12878 HiChIP data (17)], we further grouped bin pairs into three categories: the ‘AND’ set (bin pairs where both bins contain the protein of interest), the ‘XOR’ set (bin pairs where only one bin contains the protein of interest) and the ‘NOT’ set (bin pairs where neither bins contains the protein of interest). Since HP technologies measure protein-mediated long-range chromatin interactions, we only focused on the ‘AND’ and ‘XOR’ sets for downstream analysis. In addition, we filtered out bin pairs with either end overlapping with the encyclopedia of DNA elements (ENCODE) blacklist regions [29] or with low mappability (mappability score < 0.9) [30]. After filtering, we created a 10 Kb bin resolution contact matrix for each chromosome. In this work, we only used 10 Kb bin pairs with 1D genomic distance less than 2Mb in our analysis.

• (iii) Converting contact matrices into training data, testing data and validation data. According to the required format of each deep learning method, we split the contact matrix for each chromosome into multiple submatrices. Different deep learning methods adopt different splitting strategies as their default configuration. DeepHiC splits the high-depth and low-depth contact matrices into nonoverlapping 40 × 40 submatrices with stride size 40 × 40. In contrast, HiCNN, HiCNN2 and HiCPlus partition the low-depth data with overlapping 40 × 40 submatrices with a stride size 34 × 34 (the overlapping region between two consecutive submatrices is 6 × 40). Next, HiCNN, HiCNN2 and HiCPlus partition the high-depth data into nonoverlapping 28 × 28 submatrices with stride size 28 × 28. The overlapping submatrices split by HiCNN, HiCNN2 and HiCPlus imply that all inferred regions (i.e. the 34 × 34 core regions) have flanking information. We applied each method with its default matrix splitting strategy. With those submatrices, we constructed three types of tensors: for training data, testing data and validation data, respectively. Here, training data is the tensor concatenating data from five chromosomes (1, 3, 5, 7 and 9), validation data is a tensor of chromosome 2, and testing data contains data from chromosomes 4, 6, 8, 10–19 or chromosomes 4, 6, 8, 10–22, for the mESC PLAC-seq data or the GM12878 HiChIP data, respectively.

The principle of deep learning

All three deep learning methods (HiCNN2, HiCPlus and DeepHiC) evaluated in this study are supervised learning algorithms, which can be formulated as the following:

(1)

where represents the low-depth data in the training dataset (see Data preprocessing), represents the enhanced data and represents the parameters of neural network , which approximates the mapping → by learning from the training dataset. Each parameter (with being the total number of parameters of the neural network) can be optimized by the gradient descent algorithm (e.g. stochastic gradient descent or Adam [31]) in Equation (2)

(2)

where is the learning rate, which controls the step size of gradient descending, is the predefined loss function and is the gradient of which is calculated by backpropagation algorithm [10]. In HiCNN2 and HiCPlus, is the mean squared error as specified in Equation (3)

(3)

where is the sample size and it is the product of batch size (hyperparameter), the width of the submatrix and the height of the submatrix ; is a batch of target data [high-depth submatrices, Data preprocessing (3) in Methods section], is the high-depth matrix; represents the low-depth matrix; is the enhanced matrix; is the neural network which represents the mapping of to ( → ).

In HiCNN2 or HiCPlus, the training loss (represented by Equation 3) is optimized by the gradient descent algorithm (equation 2) iteratively. Each method trains the network using multiple epochs, with the default being 500 and 40 000 for HiCNN2 and HiCPlus, respectively. One epoch involves passing all batches completely through the neural network. In each epoch, HiCNN2 or HiCPlus uses validation loss to evaluate whether to retain the current trained model or not. Specifically, the algorithm calculates the validation loss between full data (viewed as the target data) and the enhanced data using Equation 3 and updates to the current model only when the validation loss decreases.

Generating HiC_downsampled data

In addition, we conducted additional experiments for evaluation of transferability, where models were trained on the mESC Hi-C data with comparable sequencing depth as the mESC PLAC-seq data, which we referred to as HiC_downsampled data. We generated the HiC_downsampled data by downsampling read counts within 2Mb genomic distance of mESC Hi-C data to 59.2 million, matching the total number of reads (59.2 million) in the ‘AND’ and ‘XOR’ sets of the corresponding mESC PLAC-seq data (Table 1).

Yuchen Yang is an associate professor in the School of Ecology at Sun Yat-sen University.

Gang Li is a PhD student in the Department of Statistics and Operations Research at the University of North Carolina at Chapel Hill.

Minzhi Jiang is a PhD student in the Department of Applied Physical Sciences at the University of North Carolina at Chapel Hill.

Jia Wen is a postdoctoral researcher in the Department of Genetics at the University of North Carolina at Chapel Hill.

Armen Abnousi was a postdoctoral researcher in the Department of Quantitative Health Sciences at Lerner Research Institute, Cleveland Clinic Foundation, and now is a senior innovation software engineer at NovaSignal.

Jonathan D. Rosen is a PhD student in the Department of Biostatistics at the University of North Carolina at Chapel Hill.

Ming Hu is an assistant staff in the Department of Quantitative Health Sciences at Lerner Research Institute, Cleveland Clinic Foundation.

Yun Li is a professor in the Departments of Genetics, Biostatistics and Computer Science at the University of North Carolina at Chapel Hill.