Improved quality metrics for association and reproducibility in chromatin accessibility data using mutual information

Abstract

Background

Correlation metrics are widely utilized in genomics analysis and often implemented with little regard to assumptions of normality, homoscedasticity, and independence of values. This is especially true when comparing values between replicated sequencing experiments that probe chromatin accessibility, such as assays for transposase-accessible chromatin via sequencing (ATAC-seq). Such data can possess several regions across the human genome with little to no sequencing depth and are thus non-normal with a large portion of zero values. Despite distributed use in the epigenomics field, few studies have evaluated and benchmarked how correlation and association statistics behave across ATAC-seq experiments with known differences or the effects of removing specific outliers from the data. Here, we developed a computational simulation of ATAC-seq data to elucidate the behavior of correlation statistics and to compare their accuracy under set conditions of reproducibility.

Results

Using these simulations, we monitored the behavior of several correlation statistics, including the Pearson’s R and Spearman’s $\rho$ coefficients as well as Kendall’s $\tau$ and Top–Down correlation. We also test the behavior of association measures, including the coefficient of determination R$^2$, Kendall’s W, and normalized mutual information. Our experiments reveal an insensitivity of most statistics, including Spearman’s $\rho$, Kendall’s $\tau$, and Kendall’s W, to increasing differences between simulated ATAC-seq replicates. The removal of co-zeros (regions lacking mapped sequenced reads) between simulated experiments greatly improves the estimates of correlation and association. After removing co-zeros, the R$^2$ coefficient and normalized mutual information display the best performance, having a closer one-to-one relationship with the known portion of shared, enhanced loci between simulated replicates. When comparing values between experimental ATAC-seq data using a random forest model, mutual information best predicts ATAC-seq replicate relationships.

Conclusions

Collectively, this study demonstrates how measures of correlation and association can behave in epigenomics experiments. We provide improved strategies for quantifying relationships in these increasingly prevalent and important chromatin accessibility assays.

Background

Epigenetic modifications play an important role in regulating multiple cellular processes ranging from DNA replication to gene expression. These covalent additions to DNA and histone proteins do not alter the underlying DNA sequence, but rather, help modulate chromatin structure resulting in distinctive phenotypes. Genome-wide epigenetic modifications can be determined using several techniques: the gold-standard is chromatin immunoprecipitation followed by sequencing (ChIP-seq) [1,2,3]. Chromatin accessibility, or the analysis of the regions that are available for DNA:protein interactions potentially resulting in gene expression, is measured using an enzyme-driven assay called transposase-accessible chromatin via sequencing (ATAC-seq) [4]. These two methods have distinct advantages in probing the state of the epigenome, and both approaches generate paired-end sequencing libraries. These data are mapped to the genome to determine the loci that are occupied with a particular epigenetic modification or the loci that are localized within an open, accessible region. Epigenetic modifications and chromatin accessibility are visualized as peaks resulting from the aggregation of sequencing reads [5]. As such, many software platforms used for analysis of ChIP-seq and ATAC-seq data sets use ‘peak calling’ to determine locations of epigenetic modifications or accessible chromatin regions [6,7,8,9].

To ensure significance and consistency of identified peaks, best practices have been defined for quantifying reproducibility across experimental replicates [8, 10]. These include several quality control metrics and workflows that standardize analysis and enable comparison among different experiments [10]. These standards apply to the total number of sequenced reads, total number of identified significant peaks, and concentration of sequenced reads within said peaks. For example, pseudo-replication was developed for ChIP-seq analysis to assess the amount of variation between biological replicates [8]. In this protocol, synthetic replicates are created from true, experimentally derived data: to do this, aligned reads are merged from two true replicates and randomly reassigned into new alignments to create two synthetic replicates. This permutation practice homogenizes (and splits) signals present within the true, observed replicates, generating the null hypothesis of near perfect correlation between pseudo-replicates. Peak calling is then also conducted on pseudo-replicates, and the read counts of peaks conserved between the two pseudo-replicates are compared to the observed peaks in the true replicates. Landt et al. proposed that experiments, whose number of observed peak counts (among true replicates) divided by the total number of pseudo peaks (between pseudo-replicates), which nears a value of one, are broadly reproducible [8]. The ENCODE project has since extended this practice to ATAC-seq experiments [11, 12].

To better understand experimental reproducibility, many studies also conduct correlation analysis on binned signals between ATAC-seq replicates [13,14,15]. In such analyses, for each replicate, the genome is binned into smaller, contiguous regions, for example using windows of ten kilobase pairs [13]. The number of mapped sequenced fragments (defined by a pair of mapped reads) that overlap these bins are counted and standardized to fragments per kilobase pair per million reads (Fpkm) [16]. These Fpkm counts are then compared between replicates using correlation and association statistics such as Pearson’s R or the coefficient of determination (R$^2$), respectively. Values from these statistics trending toward a value of one generally indicate a reproducible experiment [17].

Correlation analysis is a useful tool, not singularly purposed for the analysis of reproducibility in ATAC-seq experiments. Such analysis can be found within studies of chromosome accessibility in cancer, ageing of human stem cells, cellular diversity, or new ATAC-seq protocols [18,19,20,21,22,23]. Furthermore, correlation analyses are ubiquitous, found in the fields of genetics, RNA-seq experiments, and in studies of 3D chromatin architecture [16, 24,25,26,27,28,29,30]. Given their popularity and use in genomic and epigenetic studies, software suites—for example deeptools and HiCExplorer—have developed methods and tools for calculating correlation metrics between replicates and experiments [13, 31,32,33,34].

The natural properties of data from genomic and epigenomic experiments make the application of commonly used correlation and association statistics, for example Pearson’s R and R$^2$, potentially problematic as none of these data (ATAC-, ChIP-, or Hi-C seq) are normally distributed [35]. Both ATAC- and ChIP-seq experiments are defined by numerous, loci-specific peaks of signal generated by the accumulation of sequencing reads [3, 4]. Mapped sequenced fragments may overlap contiguous genomic bins used in analysis, producing non-independent data points [24]. Conversely, regions lacking assayed modifications or with inaccessible chromatin will have little to zero signal for ChIP-seq or ATAC-seq data, respectively. Furthermore, during correlation analysis, several genomic bins may overlap an inaccessible chromatin region that is reproducible, appearing in both the ATAC-seq replicates (or experiments) being compared. As such, each of these bins will acquire zero Fpkm and within the bi-variate distribution formed between the replicates. These data points, which appear as zero Fpkm in both replicates, are referred to here as co-zeros. Some analysis programs, like deeptools, HiCExplorer, and HiCcompare, offer options to remove co-zeros prior to analysis [29, 31, 34]. However, there is no published guidance on this practice, and while the co-zero values are a feature common across genomic and epigenomic data sets [36], the effect of removing such features on correlation statistics has not been explored. Despite the known features of genomic and epigenomic data, and the underlying assumptions of statistical tests, there have been few studies that explore their expected behavior, accuracy, and use of alternative statistics determining reproducibility of such data [26, 27].

Here, we present a computational approach to generate synthetic ATAC-seq replicates to explore the behavior of various correlation and association metrics for epigenomics datasets. These synthetic ATAC-seq replicates are generated from eight true data sets to capture features uniquely present within ATAC-seq experiments. We have developed a random subsampling strategy to generate synthetic replicates with varying portions of shared peaks, as a proxy for reproducibility. Across our simulations, we apply the Pearson’s R [37,38,39] and Spearman’s $\rho$ [40] and monitor their behavior, including the effect of removing co-zeros. Additionally, we demonstrate the behavior of other statistics, including non-parametrics such as Kendall’s $\tau$ [41,42,43,44] and an information theoretic approach, normalized mutual information [45, 46], to determine their utility in assessing epigenomics data. Finally, we build a random forest model [47] using the normalized mutual information and R$^2$ coefficient between experiments to predict the biological relationships between replicates. Overall, our results demonstrate an improvement in the expected behavior of all statistics after removing co-zeros and normalized mutual information emerges as a promising statistic for measuring association between ATAC-seq samples.

Results

ATAC-seq data characteristics and subsamping strategy for synthetic replicate generation

Fig. 1

ATAC-seq profiles of chromosome 9 form A549 cells. A TN5 binds to open chromatin, cutting DNA and adding primers to generate a paired-end sequencing library. B A549, ATAC-seq replicates along chromosome 9. Samples were generated using fresh cells (green) and previously cryo-preserved cell cultures (orange and brown). Positively (black) and negatively oriented genes are annotated along the bottom. C Pair-wise, bi-variate scatter plots of whole fragments per kb per million values (x- and y-axis) using 10 kb genomic bins between A549, ATAC-seq replicates. Sample names are annotated along the diagonal. Pair-wise Pearson’s correlation statistic is annotated within subplots

Table 1 ATAC-seq experiments used, mapped reads, peak counts and FrIP scores

To study the behavior of correlation measurements between ATAC-seq replicates (Fig. 1A), we analyzed data from three experiments using the A549, human lung cell line and implemented a subsampling paradigm to generate synthetic replicates. Across these experiments, the total number of reads mapped to the human reference genome varied from 15 million to nearly 43 million (Table 1). The number of genome-wide peaks found in the ATAC-seq samples varied across experiments and between replicates, ranging from approximately 80 to 130 thousand (Table 1). The fraction of sequenced read-pairs mapped in peaks (i.e. the FrIP score as defined by the ENCODE project [8, 11]), was greater than 0.34 for all of the A549 ATAC-seq samples (Table 1). These samples displayed high spatial correlation of peaks across replicates (Fig. 1B). Counting all whole fragments per kilobase per million (WFpkm), every ten kilobases, we observed a high statistical correlation between replicates, with average Pearson’s R of 0.86, 0.87, and 0.94 (p-values $< 0.05$) between the technical replicates of the three biological replicate experiments (Fig. 1C).

Fig. 2

Synthetic replicate generation via peak down-sampling. A An example region along chromosome 17 of true, A549 ATAC-seq data. Real ATAC-seq signal (brown lines) is used to initialize two synthetic replicates. Red and black horizontal bodies depict negatively and positively oriented genes, respectively. B A portion of the genome-wide significant peaks (ranging from 0 to 1) are chosen randomly between the two synthetic replicates. Within one of the replicates, 85% of paired reads (blue and orange rectangles connected by grey dotted line) are removed to down-sample signal within that locus. C Example of two synthetic replicates with a known portion of peaks varying between them

Fig. 3

Synthetic replicate bivariate plots and statistical profiles. A Scatter plots displaying counts per genomic bin (10 kb in size) of whole fragments per kilobase per million (WFpkm) between two synthetic replicates (x- and y-axis) generated in process Fig. 2A–C. The percentage of shared peaks decreases between the two simulated replicates from top to bottom. B Correlation values (y-axis) as a function of percentage of shared peaks between synthetic replicates (x-axis). C Association scores (y-axis) as a function of the percent of shared peaks between synthetic replicates (x-axis). In B and C, red and grey curves depict the mean and 95% confidence interval (respectively) across simulations. A dashed line marks a one-to-one relationship between the x- and y-axis. Left and right columns display change in values as a function of removing co-zeros

For simulations, synthetic replicates were generated using the paired-end read alignment profiles from the eight ATAC-seq samples we generated. For each simulation, two synthetic replicates were initiated by duplicating a given true ATAC-seq experiment (Fig. 2A). Within the true ATAC-seq data set, reproducible, significant peaks were identified (see Methods). From these, a random portion of peaks was chosen to vary between the two synthetic replicates. This was accomplished by subsampling and removing 85% of the aligned sequenced fragments within each of the randomly chosen peaks between the two synthetic replicates (Fig. 2B, C). This process was repeated fifteen times for each of the eight real ATAC-seq samples, randomly varying the common peaks from 5 to 99% of peaks between the two synthetic replicates. Finally, across all 120 simulations, for each pair of synthetic replicates, the WFpkm values were calculated in ten kilobase windows and used in statistical comparisons (Fig. 3A).

Top–down correlation displays best behavior in correlation analysis across simulations

Across the 120 down sampling simulations, correlation and association statistics were calculated between each pair of synthetic replicates. The Wfpkm counts were used between synthetic replicates in statistical analysis (Fig. 3A). The values of correlation (Fig. 3B) and association (Fig. 3C) statistics were calculated within each simulation as a function of the number of shared peaks between synthetic replicates. For each examined statistic, the area under the curve (AUC), formed by the statistical values calculated across portions of shared peaks, for each simulation was used in comparisons (Additional file 1: Fig. S1). Of the correlation statistics, the Top–Down correlation statistic had the smallest average AUC of 0.6881 (95% CI 0.6860–0.6906) and was significantly smaller than the average AUC of the Pearson’s R, at 0.8284 (95% CI 0.8237–0.8335, p-value $= 0$, bootstrapped difference of mean AUC). Both the two non-parametric statistics, Spearman’s $\rho$ and Kendall $\tau$, had significantly larger average AUCs compared against the Pearson’s R (p-values $= 0$, bootstrapped difference of mean AUC). However, they demonstrated nearly identical AUC profiles compared to each other, with average AUC of 0.9140 (95% CI 0.9118–0.9162) and 0.9096 (95% CI 0.9074–0.9120) respectively (p-value $= 0.037$, bootstrapped difference of mean AUC).

Across the metrics of association, Kendall’s W, normalized mutual information, and the R$^2$ coefficient, between replicates, the R$^2$ coefficient exhibited the greatest sensitivity to the change in portion of shared peaks between synthetic replicates (Fig. 3C). Across simulations, the average AUC of the R$^2$ coefficient was 0.7026 (95% CI 0.6951–0.7102). This average AUC was significantly smaller than the average AUC of the Kendall’s W and normalized mutual information, with values of 0.957 (95% CI 0.9559–0.9581) and 0.8197 (95% CI 0.8153–0.8241), respectively (p-value $= 0$, bootstrapped difference of mean AUC).

Removal of co-zeros improves estimates of correlation and associations

Using this simulation paradigm, we evaluated the efficacy of removing co-zeros from the analysis to determine the impact on correlation and association statistics. Co-zero values were defined as value counts in ATAC-seq experiments that appeared to have zero aligned fragments in a genomic bin of ten kilobases between two replicates (Additional file 2: Fig. S2). On average, these values can make up nearly 5% of a given bi-variate distribution formed between real ATAC-seq replicates (Additional file 3: Fig. S3). Across all the correlation and association statistics examined here—except for Top–Down correlation—removing the co-zero values significantly reduced the average AUC (Table 2, Fig. 3B, C, Additional file 1: Fig. S1). The large reduction observed in the AUC after removing co-zeros from analysis was unexpected, as co-zeros are a modest portion of the bi-variate distribution formed between two replicates.

Table 2 Mean area under the curve across simulations

After removing co-zeros, all the correlation statistics, Top–Down correlation, Pearson’s R, Spearman’s $\rho$, and Kendall’s $\tau$, displayed nearly identical sensitivity to the change in shared peaks between replicates across simulations (Fig. 3B). However, the Pearson’s R had the largest average AUC of 0.6965 (95% CI 0.6946–0.6984) followed by the Top–Down statistic (AUC of 0.6872, 95% CI 0.685–0.6895, p-value $= 0$, bootstrapped difference of mean AUC). The Spearman’s $\rho$ (mean AUC: 0.6686, 95% CI 0.6665–0.6705) and Kendall’s $\tau$ (mean AUC: 0.6673, 95% CI 0.6654–0.6691) statistics had the smallest and identical average AUC after removing co-zeros (p-value $= 0.208$, bootstrapped difference of mean AUC). Furthermore, the AUC of the Top–Down correlation statistic was unaltered by the exclusion of co-zero values between synthetic replicates (Fig. 3B, Additional file 1: Fig. S1, Table 2, p-value $= 0.635$, bootstrapped difference of mean AUC). This observation was not surprising given how Top–Down correlation places emphasis on larger values, down-weighting smaller values, such as co-zeros [48].

Normalized mutual information best estimates difference between replicates

Removing co-zero values had a similar effect on association metrics, attenuating and improving the average AUC across the portion of shared peaks between synthetic replicates (Fig. 3C, Additional file 1: Fig. S1). Apart from Kendall’s W, the R$^2$ coefficient and normalized mutual information, on average, displayed a nearly one-to-one relationship with the portion of shared peaks between replicates (Fig. 3C). The average AUC of normalized mutual information was 0.5055 (95% CI 0.5045–0.5065) and was smaller than the average AUC of the R$^2$ coefficient, with a value of 0.5346 (95% CI 0.5324–0.5368, p-value $= 0$, bootstrapped difference of mean AUC). This difference in average AUC indicates that normalized mutual information better follows the designed proportion of shared peaks between synthetic replicates across our simulations, compared to the R$^2$ coefficient.

As introduced earlier, one parameter in this simulation is the removal of a percentage of aligned read-pairs from within randomly selected peaks (Fig. 2B). Initially set at 85%, this parameter was altered to simulate ATAC-seq replicates that are nearly reproducible (at 50%) at every selected peak or broadly unreproducible (at 95%) across all selected peaks. Comparing the results between the two simulation sets with 85 and 95% of reads removed, we observed no significant difference between the two simulations (Additional file 4: Fig. S4, Additional file 5: Fig. S5). This is somewhat expected when considering the small difference in magnitude between removing 85 and 95% of reads from within peaks. In simulations with only 50% of read pairs removed from selected peaks, after removing co-zeros, the two statistics that showed the largest response in our simulation were the R$^2$ coefficient and normalized mutual information.

Co-zeros inflate estimates of correlation and association in epigenomic assays

Fig. 4

Correlation and association statistics across epigenomic experiments. For samples from (A) ATAC-seq and ChIP-seq (assays for (B) H3K27ac and (C) H3K4me3 modifications) experiments, the Spearman’s $\rho$, Pearson’s R, R$^2$ coefficient, and normalized mutual information (y-axis of columns left to right, respectively) were calculated on WFpkm counts between replicates, with and without co-zeros (x-axis). Red and black dumbbells represent calculations between replicates or non-replicates, respectively, and connect calculations across the co-zero handling strategy

After successfully implementing normalized mutual information between simulated replicates, we next examined how this statistic behaves when used on replicates from real epigenomic experiments. We also monitored how dropping co-zeros affects estimates of correlation (and association) between samples. For these analyses, additional experiments were downloaded from the ENCODE project public repository [11]. These included additional ATAC-seq experiments (Table 1) as well as ChIP-seq experiments, specifically twenty and ten assays for the H3K27ac and H3K4me3 modifications, respectively (Additional file 10: Table S2). The Spearman’s $\rho$, Pearson’s R, R$^2$ coefficient, and normalized mutual information were calculated between these replicates. Correlation and association values were also calculated between non-replicates within each of the three assays. We then repeated these comparisons, eliminating co-zeros from calculations. With this design, we were able to gauge the effect of masking co-zeros within replicates, between non-replicates (within the same assay), and across different types of epigenomic data. Between real experiments, excluding co-zeros from analysis significantly decreased the computed correlation and association statistics (p-value $< 1^{-10}$, Wilcoxon signed-rank test). This reduction is seen in the distributions of the Spearman’s $\rho$, Pearson’s R, R$^2$ coefficient, and normalized mutual information (Additional file 6: Fig. S6) across all three assays, ATAC-seq (Fig. 4A), H3K27ac (Fig. 4B), and H3K4me3 (Fig. 4C). Further investigation revealed that omitting co-zeros primarily alters estimates of correlation and association between non-replicates (black lines and dots in Fig. 4), which were significantly decreased (p-value $< 1^{-19}$, Wilcoxon signed-rank tests). The correlation and association values between true replicates from H3K27ac and H3k4me3 assays were unaltered when ignoring co-zeros (p-value $> 0.001$, Wilcoxon signed-rank tests). However, a significant alteration in correlation (and association) estimates between replicate ATAC-seq experiments (red lines and dots in Fig. 4A) was detected (p-values $< 0.0009$, Wilcoxon signed-rank tests). Moreover, excising co-zeros expanded the difference in the average estimates of correlation and association between replicates and non-replicates, across all three assays (Fig. 4, Additional file 11: Table S3). Thus, removing co-zeros produced lower estimates of correlation and association for most samples, and overall improves the ability to differentiate pairs of replicates from non-replicates in real epigenomic data.

A random forest prioritizes mutual information for predicting replicate relationships

Fig. 5

Random forest prediction of experimental relationships. A Distributions of the coefficient of determination (R$^2$) and normalized mutual information scores calculated on binned counts of WFpkm between ATAC-seq experiments. Blue, orange, and green dots mark comparisons between independent experiments, independent experiments using the same cell line, and true experimental replicates, respectively. B Example confusion matrix from a random forest model using R$^2$ and normalized mutual information as features to predict experimental relationships (y-axis) presented in A (x-axis). The confusion matrix depicts results of model on a hold-out set (40% of data, accuracy = 95.12%). Light to dark colors depict the number of counts per class. C Bi-variate plot displaying the change of paired importance scores from ten-fold cross validation between the normalized mutual information (x-axis) and R$^2$ (y-axis) features. Dashed lines depict the uni-variate means of the normalized mutual information and R$^2$ scores. Blue and yellow colors depict the level of accuracy for each fold

After removing co-zeros, the R$^2$ and normalized mutual information metrics performed best in simulation. Furthermore, when tested on real ATAC-seq experiments, these two metrics produced the largest difference (on average) between replicates and non-replicates. Given the comparable behavior of normalized mutual information and the R$^2$ coefficient on true ATAC-seq data we set out to further assess their usefulness in predicting the relationships between experiments. To do this we combined our A549 ATAC-seq samples with ATAC-seq samples from the ENCODE project. These included additional biological replicates of the A549 cell line, as well as ATAC-seq experiments in the HepG2, RWPE2, GM12878, IMR-90, K562, and WTC11 cell lines (Table 1). With these combined data, comparisons between any two ATAC-seq experiments were classified into one of three discrete, replicate classes; (1) between independent ATAC-seq experiments in different cell lines, (2) between independent experiments using the same cell line, or (3) between true replicates. Plotting the normalized mutual information and R$^2$ coefficient calculated between ATAC-seq experiments with the above classifications revealed clustering of the classes between replicates (Fig. 5A). Between the two statistics, the R$^2$ coefficient displayed the largest variation across compared experiments (Table 2) and the most mixing of the three class labels (Fig. 5A). We also observed a strong co-linear relationship between the calculated R$^2$ coefficient and normalized mutual information scores (Additional file 7: Fig. S7, Pearson’s R $= 0.96$, p-value $< 1^{-10}$). This finding is not surprising given that both metrics attempt to measure the same relationship between samples.

To quantify which statistic (the R$^2$ coefficient or normalized mutual information) better estimates experimental relationships between replicates, we built a random forest model. This model uses the reported values of the R$^2$ coefficient and normalized mutual information between ATAC-seq experiments as features to predict the replicate class (as defined above). We utilized ten-fold cross validation, stratifying on the replicate class to build our random forest. An example confusion matrix from one of these folds demonstrates the model had difficulty distinguishing between independent experiments using the same cell line and true experimental replicates (Fig. 5B). This difficulty also manifested as lower f1-scores and recall for this class (Additional file 8: Fig. S8). The accuracy across these folds ranged from 88 to 98% (Fig. 5C). Across the folds, the feature importance score of the R$^2$ coefficient was inverted with that of normalized mutual information (Fig. 5C). Overall, we observed a greater feature importance score for normalized mutual information, with a significant average pair-wise difference between the R$^2$ coefficient and normalized mutual information of 6.78% (p-value $< 0.05$, Wilcoxon signed-rank test).

Discussion

To improve the assessment of reproducibility in epigenomic data sets, we sought to investigate the use of several correlation and association statistics on binned genomic signals. Our findings suggest that best practices should include analyzing association between compared replicates (or experiments) via normalized mutual information with binned, Fpkm counts rounded to the nearest whole integer, after the removal of co-zero values as input. In choosing a correlation statistic, after removing co-zero values, our results indicate little difference in the outputs from the Pearson’s R, Spearman’s $\rho$, Kendall’s $\tau$, or Top–Down correlation statistics. Notably, from simulations, we observed that the Top–Down correlation statistic was unaffected by the removal of co-zeros values. As such, this statistic should serve as an alternative for investigators if binned co-zeros values between replicates are retained.

As part of this study, we generated highly correlated, new ATAC-seq experimental replicates of the A549 cell line. Our data highly correlates with previously published ATAC-seq data of the A549 cell line generated by the ENCODE project. Using these data, we generated a novel simulation that utilizes down sampling to generate replicates with known varying signals. While similar simulation studies have been conducted on Hi-C sequencing data [30], to our knowledge, no prior study has examined the behavior of statistical metrics on ATAC-seq data. That said, there are several statistics and methodologies that may be used to analyze this data type, such as Poisson regression [49]. Improving on this simulation design could help generate a framework that allows researchers to develop new statistical tools for hypothesis testing.

In our simulations, we observed that most statistics overestimate the correlation of signal between replicates. One specific strategy we investigated to reduce this inflation was the removal of co-zeros, which is an option present in several bioinformatic software suites [29, 31, 34]. Our analysis demonstrated that removal of these values can provide a more accurate estimate of correlation between replicates as measured by the known number of peaks between replicates. Interestingly, we never observed a correlation value that perfectly trends with the designed number of peaks between synthetic replicates. We also did not observer negative correlation values between the replicate Fpkm counts. The first of these observations can be explained by background autocorrelation still present within our synthetic replicates. The second of these observations may point to a limitation in the design of our simulation, as negative correlation values have been observed in true ATAC-seq profiles [20, 31].

In epigenomics and chromatin accessibility data sets, biological interpretation of the data is dependent upon visualization of “peaks” where accumulation of sequenced reads denotes the presence of a modification or an accessible region. Regions with zero (or nearly zero) aligned sequenced reads are deemed unmodified or inaccessible and largely ignored when interpreting data. Correlation statistics should provide biologists with the confidence that replicates are truly comparable. As stated above, the inclusion of co-zeros seems to inflate values of most correlation and association statistics. Thus, removal of co-zeros formed by the genomic bins that overlap and account for inaccessible regions may be warranted.

Using our simulation, we also examined the behavior of three association statistics, which we distinguish from the set of correlation statistics as those metrics ranging in value from zero to one. These association statistics were the R$^2$ coefficient, normalized mutual information statistic, and Kendall’s W. Prior to the removal of co-zeros, the only association statistic that displayed any sensitivity to the change in shared peaks between replicates was the R$^2$ coefficient. Co-zeros inflate the value of this statistic by reducing the total summed error between data points during calculation. Similarly, co-zeros increase the information gained between replicates when calculating the normalized mutual information score. In other words, knowing a replicate has a value of zero at a given genomic bin provides information that there is a zero at the corresponding bin within the other replicate. After removing co-zeros, we saw a large improvement in the sensitivity of both these statistics.

Curiously, Kendall’s W displayed the least sensitivity to the designed peak counts between synthetic replicates. This statistic was of interest given Kendall’s W is capable of simultaneously examining the ranks of more than two input samples [41, 50]. This would have provided researchers with a statistical tool capable of examining correlation among a full set (triplicate) of replicates within a single test, rather than multiple pair-wise comparisons. Removing co-zeros did little to improve the sensitivity of this statistic. The other statistic from Kendall, Kendall’s $\tau$, displayed similar performance to the other non-parametric statistic, Spearman’s $\rho$. This finding is contrary to other studies of Kendall’s $\tau$ conducted in the fields of signal processing and psychology [43, 44]. For analysis of genomic data, the Spearman’s $\rho$ is standard in deeptools’ correlation functions [13].

We also examined the effect of dropping co-zeros when estimating correlation between real ATAC-seq and ChIP-seq samples. Much like our results from analysis on simulated ATAC-seq replicates, expunging co-zero values from correlation (and association) calculations reduced the value of the reported statistics between real samples. These effects were primarily seen in the correlation and association scores calculated between non-replicates within the examined assays. In particular, the correlation and association values between true H3K27ac and H3K4me3 replicates were unaffected by eliminating co-zeros. These correlation scores were high and remained high after excising co-zeros from calculations. This may be due to higher overall signal in these assays. Importantly, omitting co-zeros from analysis produced a larger difference in the average correlation between groups of replicates and non-replicates. Thus, excluding co-zeros from analysis is an important step for quality control procedures looking to identify errant samples.

Of the statistics examined here, the R$^2$ coefficient and normalized mutual information score were the most sensitive to the change in shared peaks between replicates (when co-zeros were removed). Comparison of these two statistics revealed that normalized mutual information was the better-behaved statistic. This behavior manifested as smaller AUC within simulations, less variation in values, and larger differences in values between groups of replicates and non-replicates. Similarly, the computational evidence provided by our random forest model suggests that normalized mutual information was better at estimating experimental relationships between true ATAC-seq replicates. Taken together, these results indicate that of the two metrics, normalized mutual information may be the stronger association metric for ATAC-seq data. Information theoretic approaches, such as normalized mutual information, have been utilized in several other biological fields, ranging from cancer genomics to fungal genetics [51,52,53,54,55,56,57]. Regarding ATAC-seq data, a handful of other studies have specifically used mutual information in data integration, analysis, and deep-learning of single-cell ATAC-seq profiles [58, 59]. For those investigator interested in using information theoretic approaches, several of these functions are made easily available within the python, scikit learn library [46].

To perform correlation and association analyses as seen here, we have generated python code and an executable for public use. These software, installation instructions, and a tutorial written as a jupyter notebook are hosted on the Github listed within the data availability section. We hope these tools will benefit investigators and students in their exploration of the mutual information statistic and the effect of excluding co-zero values in their epigenetic data.

Sparsity and zero mapped sequenced reads are not unique properties of ATAC-seq data. These extend to genomic, Hi-C, ChIP-seq, and RNA-seq data sets. Imputation along with modified zero-inflated models have been used with success for studying RNA sequencing data sets and detecting regions with differential expression [60]. Simulations and models of sampling zero-genomic count data have been developed to understand the effects of these values, particularly in the context of differential analysis [36]. Previous simulation studies of ATAC-seq have been focused on generating ATAC-seq data, for pipeline development, or single-cell ATAC-seq samples, to examined different approaches in their analysis [61, 62]. To our knowledge, this is the first example of using a simulation approach for studying reproducibility and association of ATAC-seq samples. Adapting strategies from these previous studies will help improve our simulation and expand it to other genomic and epigenomic sequencing data. The current results of our study strongly suggest that normalized mutual information is an appropriate metric for measuring reproducibility in chromatin accessibility assays.

Conclusions

For this study, we produced eight ATAC-seq experiments using the A549 Cancer cell line. Across replicates, these ATAC-seq samples are well correlated and reproducible. For investigations of chromatin accessibility (particularly in the A549 cell line), these experiments are an additional resource for developing analysis pipelines, peak detection algorithms, and machine learning approaches.

Leveraging the A549 ATAC-seq experiments, we designed a computational simulation to generate simulated replicates. Specifically, synthetic replicates were coded that share a known, fixed portion of significantly enriched loci. Using these replicates, correlation metrics—the Pearson’s R, Spearman’s $\rho$, Top–Down, and Kendall’s $\tau$—and association statistics (ranging from zero to one)—the R$^2$ coefficient, Kendall’s W, and normalized mutual information—were tested for accuracy. Overall, the reported value of these statistics was inflated and much larger than the fixed portion of shared, significant loci between replicates.

Removing specific outliers from ATAC-seq data, specifically the removal of co-zeros, improved estimates of correlation and association. We estimate that co-zero values, when comparing WFpkm counts between two real ATAC-seq experiments, occupy nearly five percent of a bi-variate distribution. While only a small portion of the total data, filtering these values from analysis greatly improves the measurements of most correlation and association statistics between samples, in simulation. Applied to real ATAC-seq and ChIP-seq data, removing co-zero values from comparison significantly reduced the reported correlation and association statistic, matching results from simulation.

One of the association statistics examined here is normalized mutual information, an information theoretic approach that is less well known across the (epi)genomics field. After removing co-zero values, normalized mutual information displayed the lowest inflation relative to the similarity between simulated replicates. The R$^2$ coefficient also performed well in simulations (after removal of co-zeros), displaying good sensitivity to differences between simulated replicates. Of these two association metrics, a random forest model selected normalized mutual information as the stronger feature when estimating experimental relationships between real ATAC-seq experiments. From these results we conclude that normalized mutual information is a powerful, non-parametric approach for estimating association between ATAC-seq experiments.

Methods

Construction of A549 ATAC-seq libraries

ATAC-seq experimental libraries were generated using A549 human lung carcinoma epithelial cells (ATCC, VA, catalog #CCL-185) [63,64,65]. Three biological replicate libraries were prepared from freshly harvested cells using an ATAC-seq kit (Active Motif, 53150) following the manufacturer’s protocol. The remaining five libraries were prepared using cryopreserved cells following methods outlined in Milani et al.  with modifications [18]. Briefly, A549 cells were cultured in T75 flasks and harvested by trypsinization using 0.25% (w/v) Trypsin-EDTA (0.5%) solution (Gibco, 15400054). Harvested cells were centrifuged and pellets resuspended in freezing media containing DMEM (Gibco, 11885-084), 10% FBS (Corning, 35-015-CV), and 10% DMSO (ATCC, 4-X). Pellets were frozen using an isopropyl alcohol chamber (Thermo Fisher Scientific, 5100-0001) at $-\,80\,^{\circ }$C. After 24 h, frozen cells were transferred to liquid nitrogen for long term storage. To perform experiments, cryopreserved cells were transferred to $-\,80\,^{\circ }$C for several days, and the tube was immersed in 37$^\circ$C water bath for approximately two minutes on the day libraries were prepared. Thawed cells were resuspended in 1X PBS with protease inhibitor cocktail (Thermo Fisher Scientific, 78430). Cell counts and viability were assessed and aliquots containing 80,000 cells per sample were processed into ATAC-seq libraries.

Sequencing, alignment and filtering

ATAC-seq libraries were sequenced at the sequencing facility at Los Alamos National Laboratory on an Illumina NextSeq2000 sequencer in paired end mode (PE151) using P3 chemistry. With Fastp, raw reads were trimmed and filtered to remove Nextra adaptors and reads with repetitive sequences [66]. Additionally reads were also filtered to remove bases with low quality scores (q $< 15$). These processed reads were aligned to the new, telomere-to-telomere human reference genome, version 2 [67] via bwa [68]. After alignment, duplicate sequenced pairs were marked via samblaster and removed from analysis [69]. Read pairs mapping to the mitochondria were also removed (see Additional file 9: Table S1).

Other data used

Raw ATAC-seq data, in the form of paired fastq.gz files, was downloaded from the ENCODE project for the A549, HepG2, RWPE2, GM12878, IMR-90, K562, and WTC11 cell lines [11, 70]. The ENCODE file experiment and replicate accession numbers are included in Table 1. For alignment, these data were passed through the same pipeline described above for ATAC-seq samples generated here, and aligned to the human, telomere-to-telomere, reference genome [67].

For ChIP-seq experiments, twenty and ten assays of the H3K27ac and H3K4me3 epigenetic modifications (respectively) were downloaded (also) from the ENCODE project as raw alignments in bam file format. The ENCODE accession numbers of these files are listed in Additional file 10: Table S2. Each of these raw alignments were made with the GRCh38 (v1.5.1) human reference genome. Filtered bam files were generated via samtools view command and the following flags: -F 4 -F 256 -F 512 -F 1024 -F 2048 -q 30. When comparing differences between correlations and association values within experiments, between calculations with and without co-zeros, and between groups of replicates and non-replicates, a Bonferroni correction was used to calculate the adjusted p-value $= 0.05/(12 \times 3) \sim 0.00139$, for establishing significance.

Peak calling, peak filtering and reproducibility

After filtering, sample alignments were analyzed to identify loci displaying significant enrichment of paired-end reads. This peak calling was conducted using MACS2 [6, 71]. Specifically, after removing duplicates and mitochondrial mapped reads, samples were further filtered using samtools with the following flags:

-F 4 -F 256 -F 512 -F 1024 -F 2048 -q 30 and then passed to MACS2 in BAMPE mode [72, 73].

Between true, biological replicates, reproducible peaks were identified via irreproducible discovery rate thresholding [74]. Using ChIP-R, replicate narrow peak files were filtered to retain only those peaks that were consistent across all replicates; in ChIP-R, where command line parameter, m $=$ number of biological replicates [75]. In addition to this setting the ’-fragment’ option was also invoked. These sets of final peak counts were retained for further analysis.

Genomic down-sampling and simulation design

For each of the eight ATAC-seq experiments of A549 cells generated in this study, synthetic replicates were generated by duplicating a given sample into two copies and then randomly, varying the total number of shared peaks between them. Specifically, for a given ATAC-seq experiment, a set portion of peaks was chosen at random, such that within one of the synthetic replicates, a given selected peak was depleted, randomly removing a portion of the alignments within the peak bounds (as defined by MACS2). These sets of peaks were randomly selected from the set of reproducible peaks for that sample and its associated biological replicates (see above). Three sets of simulations were conducted, removing 50, 85 and 95% of reads within selected peaks. This procedure results in two synthetic ATAC-seq replicates, generated from a single, true parent ATAC-seq data set. These synthetic ‘sister’ ATAC-seq data sets have identical genome-wide alignments except within a sub-set of loci that vary between them. From each true ATAC-seq data set, synthetic sister replicates were generated by varying the total percentage of shared peaks from 99 to 5%, with a delta of 5%. For each simulation, across the change in portion of shared peaks, a common random seed was used to preserve autocorrelation across this axis. This process was repeated fifteen times for each of the eight, A549 ATAC-seq samples, totaling a one hundred and twenty simulations.

Genomic binning, fragment counts, and standardization

On both synthetic samples from simulation studies or replicates from (true) ATAC-seq experiments, a genomic binning approach was used to estimate correlation and association statistics between samples. For each chromosome, contiguous bins were established 5’–3’, every ten kilobases. Within each of these bins, the number of sequenced fragments is counted and standardized to fragments per kilobase per million. These counts were rounded up to their nearest whole integer generating standardized counts of whole fragments per kilobase per million (WFpkm).

Calculating correlation and association metrics

In python scripts, using the scipy-stats module [76], the Pearson’s R, Spearman’s $\rho$, and Kendall’s $\tau$ were calculated on the WFpkm counts between pairs of ATAC-seq replicates. Functions for the Top–Down correlation metric [48] and Kendall’s W rank statistic [41, 50] were also developed using custom python code. The R$^2$ coefficient was calculated using the square of the Pearson’s R. The normalized mutual information statistic from pythons sklearn module [46] was used in association studies. Between any pair of WFpkm counts, the bi-variate distribution was examined to identify instances were both profiles contained a value of zero WFpkm. For studies of the effects of co-zero inflation, these co-zero values were removed, and the correlation (or association) statistics recalculated on these filtered distributions.

For correlation analysis on ATAC-seq experiments conducted here using A549 cells, the Pearson’s R correlation statistic was calculated on WFpkm values between replicates with co-zeros removed. Similarly, co-zeros were removed prior to calculating correlation and association statistics between replicates of ATAC-seq data downloaded from the ENCODE project public repository.

Statistical tests on area under the curve

Across simulations, values of correlation and associations statistics were calculated as a function of the designed portion of peaks between synthetic replicates. For each statistic tested, the 95% confidence interval of the average area under the curve was calculated via bootstrapping, with a thousand iterations. This was done for statistical profiles from simulations with and without co-zero values. For comparisons of the average area under the curve, either between statistics or within the same statistic after removing co-zeros, one thousand permutations were used to calculate the null distribution of the difference between the mean area under the curve [77]. The proportion of these differences greater than or equal to the true observed difference was used as the p-value. A significance level of 0.05 was used to reject the null hypothesis, H$_0$: no difference in mean area under the curve, in favor of our alternative hypothesis, H$_1$: difference of mean area under the curve.

Design of random forest model

A random forest model was built in python using the scikit learn module [46, 47]. Association statistics from the ATAC-seq data generated in this study on A549 cells and additional ATAC-seq data downloaded from the ENCODE project was used as input (see Table 1). As features in this random forest, the R$^2$ coefficient and normalized mutual information were calculated between every pair of ATAC-seq experiments using WFpkm counts, across ten kilobase pair, genomic bins and removing co-zero values. The comparison of each unique pair of experiments (totaling 276) were discretized as (1) between independent ATAC-seq experiments in different cell lines, (2) independent experiments using the same cell line, and (3) between true replicates. The total number of comparisons distributed among these three classes was 213, 45, and 18 (respectively). Given the over-representation of comparison between independent ATAC-seq experiments in different cell lines, 39 of the 213 comparisons were chosen randomly to represent the total, unique comparisons between experiments with unique cell lines. This down sampling resulted in 39, 45, and 18 comparisons between independent experiments in different cell lines, independent experiments using the same cell line, and true replicate experiments, respectively.

For the testing and training of the model, test and training sets of the classes defined above were selected using a stratified, 40:60 split of the data. Additionally, ten-fold, stratified cross validation was used to train and test the model [78]. A hundred estimators with the entropy selection criterion were used along with default settings in the python random forest classifier function within scikit learn [46].