PETScan: score-based genome-wide association analysis of RNA-Seq and ATAC-Seq data

TL;DR

PETScan is a new method for analyzing RNA-Seq and ATAC-Seq data to study gene regulation more efficiently and accurately.

Contribution

PETScan introduces a score-based approach using negative binomial models for faster and more accurate genome-wide association analysis.

Findings

PETScan is three orders of magnitude faster than Wald tests in real-world datasets.

The method identifies significant gene-peak pairs while accounting for RNA-Seq data's count-based nature.

PETScan improves computational efficiency using score tests and matrix calculations.

Abstract

High-dimensional sequencing data, such as RNA-Seq for gene expression and ATAC-Seq for chromatin accessibility, are widely used in studying systems biology. Accessible chromatin allows transcription factors and regulatory elements to bind to DNA, thereby regulating transcription through the activation or repression of target genes. The association analysis of RNA-Seq and ATAC-Seq data provides insights into gene regulatory mechanisms. Most existing analytic tools exclusively focus on cis-associations, despite regulatory elements being able to physically interact with distant target genes. Furthermore, conventional approaches often utilize Pearson or Spearman correlations, which ignore the count-based nature of RNA-Seq data.

To address these limitations, we introduce PETScan, a computationally efficient genome-wide PEak-Transcript Score-based association analysis, utilizing negative binomial models to better accommodate RNA-Seq data. We leverage score tests and matrix calculations for improved computational efficiency, and combine an empirical permutation method with genomic control to ensure valid p-value calculations in studies with limited sample sizes. In real-world datasets, PETScan achieved three orders of magnitude faster than Wald tests, while identifying similar significant gene-peak pairs.

The PETScan R package is available on GitHub at https://github.com/yajing-hao/PETScan.

Figures

Left: Empirical type I error rates, and Right: Computational time (minutes) for PETScan (without permutation, L=0, and with L=5,10,20 permutations), Wald and LRT across different sample sizes.

Significant gene-peak pairs across the genome from ENCODE human tissue-specific data are shown, with peak locations on the X-axis and gene locations on the Y-axis. Different levels of significance are indicated using different colors. The middle panel illustrates the number of genes associated with each peak, while the bottom panel displays the number of peaks linked to each gene.

Scatter plots of chromatin accessibility at chr12:124,871,497–124,872,709 and the expression levels of its associated genes.

Heatmap of chromatin accessibility of associated peaks for the RPL32 pseudogene (left) and ST8SIA5 (right).

Scatter plots of significant gene-peak pairs from ENCODE human tissue-specific data validated by GTEx cis-eQTLs.

Significant gene-peak pairs across the genome from 10x Genomics embryonic mouse brain data are shown, with peak locations on the X-axis and gene locations on the Y-axis. Different levels of significance are indicated using different colors. The middle panel illustrates the number of genes associated with each peak, while the bottom panel displays the number of peaks linked to each gene.

Scatter plots of chromatin accessibility at chr7:54,633,595–54,635,192 and the expression levels of its negatively associated genes.

Scatter plots of chromatin accessibility at chr7:54,633,595–54,635,192 and the expression levels of its positively associated genes.

Scatter plots of significant gene-peak pairs from 10x Genomics embryonic mouse brain data validated by PLAC-Seq data.

Equations

Full text

Introduction

Transcriptomics data analysis investigates RNA transcripts produced by the genome under specific biological conditions or in particular cell types, which plays a significant role in understanding complex traits and diseases (Casamassimi et al. 2017). Although every cell in the human body shares an identical set of genes, distinct gene expression patterns define the specialized functions of different cell types both physiologically and pathologically (Ralston and Shaw 2008, Zeng 2022). These patterns can be regulated by epigenetic mechanisms, such as DNA methylation (Dhar et al. 2021) and histone modifications (Gagnidze and Pfaff 2022), which regulate gene expression through modifications of DNA molecules or chromatin without altering the underlying DNA sequence. Additionally, changes to chromatin accessibility represent another important epigenetic mechanism that regulates gene expression (Tsompana and Buck 2014). Structurally, a chromatin region can be densely packed (heterochromatin), leading to gene silencing, while loosely packed chromatin (euchromatin) tends to be transcriptionally active (Huisinga et al. 2006). The dynamic nature of chromatin accessibility plays a crucial role in gene regulation by allowing transcription factors (TFs) and other regulatory elements to bind to DNA, thereby regulating transcription through the activation or repression of target genes (Klemm et al. 2019).

High-throughput sequencing technologies, such as RNA sequencing (RNA-Seq) (Nagalakshmi et al. 2008) and Assay for Transposase-Accessible Chromatin sequencing (ATAC-Seq) (Buenrostro et al. 2013), enable the study of gene expression and chromatin accessibility, respectively. When RNA-Seq and ATAC-Seq data are jointly collected from the same samples or cells, systematic investigations of the relationship between gene expression and chromatin accessibility can be conducted to identify chromatin regions that regulate gene expression and reveal specifically impacted genes. Such analysis sheds light on the regulatory mechanisms underlying chromatin accessibility (Starks et al. 2019). Understanding these gene regulatory mechanisms is crucial for interpreting transcriptomic differences across various biological conditions, unraveling the progression of diseases (Kundaje et al. 2015), and developing effective therapeutic strategies (Dai et al. 2024). For example, research has identified disease-specific genes and uncovered their regulatory mechanisms underlying chromatin accessibility in diabetes (Ackermann et al. 2016), hepatocellular carcinoma (Yang et al. 2021, Bai et al. 2024), cellular senescence (Song et al. 2022, Ding et al. 2023), and asthma (Zhang et al. 2023). These studies provide insights into how genetic and epigenetic factors contribute to disease development and progression.

To identify associations between gene expression and chromatin accessibility, most existing analyses rely on Pearson correlation (Zhu et al. 2019, Chen et al. 2019, Stuart et al. 2021) or more robust Spearman correlation (Ma et al. 2020, Kartha et al. 2022, Park and Rhee 2024). For instances, cisDynet calculates Pearson correlation between gene expression and peak accessibility (Zhu et al. 2023), while TRIPOD detects peak-TF-gene trio regulatory relationships using Spearman correlation (Jiang et al. 2022). More advanced statistical methods include partial least-squares regression (Staitieh et al. 2023), LASSO (least absolute shrinkage and selection operator) (Cao et al. 2018), adaptive elastic-net regression (Ledru et al. 2024), and a tree-based non-linear regression method called gradient-boosting machine (Bravo González-Blas et al. 2023). All these methods ignore the count-based nature of RNA-Seq data, which can result in suboptimal performance and biased results, necessitating specialized statistical methods that appropriately handle RNA-Seq data.

Furthermore, current association analysis of RNA-Seq and ATAC-Seq data primarily focuses on identifying cis-associations in cis-regulatory regions (Frankel 2012, Klemm et al. 2019). For instance, cisDynet evaluates the association between a gene’s transcriptomic expression and peak accessibility within 500 kb upstream and downstream of its transcription start site (TSS) (Zhu et al. 2023). However, research has shown that gene regulatory elements bound to accessible chromatin can physically interact with distant target genes through chromatin looping (Kadauke and Blobel 2009), thereby modulating genes located far along the same chromosome or even on different chromosomes (Miele and Dekker 2008, Dean 2011, Dekker and Misteli 2015, Van Den Heuvel et al. 2015, Maass et al. 2019). A recent study has generated a comprehensive meta-Hi-C chromatin contact map, providing both functional cis- and trans-chromosomal interactions (Lohia et al. 2022), where trans-interactions refer to pairs of chromatin regions that are either located on different chromosomes or positioned far apart on the same chromosome.

To fully understand how chromatin accessibility regulates gene expression, a holistic approach that considers both cis- and trans-associations, or genome-wide association analysis, is necessary. However, RNA-Seq data typically includes expression profiles for over 20,000 genes, while ATAC-Seq data contains peak measurements from approximately 100,000 chromatin regions. This results in roughly 20,000 × 100,000, or 2 billion, paired association analyses, requiring a computationally efficient method that is currently unavailable.

To address these issues, we introduce PETScan, a computationally efficient genome-wide PEak-Transcript Score-based association analysis. The algorithm PETScan employs negative binomial models to respect the count nature of RNA-Seq data with potential overdispersion (Robinson et al. 2010, Di et al. 2011, Love et al. 2014). To mitigate the computational challenge, PETScan adopts two strategies to ensure its practical usage. First, PETScan utilizes score tests that require only parameter estimates under the null hypothesis. This feature becomes particularly useful since, for a given gene, its restricted estimates under the null hypothesis remain the same across all tested ATAC-Seq peaks and need to be evaluated only once.

Second, to further enhance computational efficiency, PETScan mimics Matrix eQTL (Shabalin 2012), an ultra-fast expression quantitative trait loci (eQTL) analysis tool. eQTL analysis also involves billions of paired association tests between millions of single-nucleotide polymorphisms (SNPs) and tens of thousands of genes across the entire genome. While Matrix eQTL converts traditional linear models into correlation coefficients for efficient large matrix operations, we combine score tests with ATAC-Seq data transformation to facilitate efficient matrix calculations.

Furthermore, in studies with small sample sizes, score tests may fail to converge to their theoretical asymptotic distributions. To better control the family-wise error rate or false discovery rate, PETScan incorporates an empirical permutation method inspired by genomic control (Devlin and Roeder 1999), which has been successfully used for robust statistical inference in differential gene expression analysis with RNA-Seq data (Zou et al. 2014). We demonstrate the effectiveness of PETScan using ENCODE human tissue-specific data and 10x Genomics embryonic mouse brain data. Computationally, PETScan is three orders of magnitude faster than the equivalent Wald tests.

Methods

Assume there are N paired samples of RNA-Seq and ATAC-Seq data with G genes and P peaks, respectively. Let $[eqn]$ represent a set of covariates in the ith sample ( $[eqn]$ ), such as the top principal components derived from the RNA-Seq data, which are commonly used to account for batch effects. For gene $[eqn]$ and peak $[eqn]$ , let $[eqn]$ be the expression level of gene g, and $[eqn]$ denote the log-transformed chromatin accessibility for peak p in the ith sample. For simplicity of notation, we drop g and p out of the subscripts when describing the model for testing the association between the gene g and peak p pair in the following discussion, whenever there is no confusion.

To adjust for potential over-dispersion, we assume that $[eqn]$ follows a negative binomial distribution with the density function $[eqn]$ , where $[eqn]$ represents the mean and $[eqn]$ denotes the dispersion parameter with E( $[eqn]$ ) = $[eqn]$ and Var( $[eqn]$ ) = $[eqn]$ . Specifically,

$$[eqn]$$

with $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ .

The log likelihood function of the data is therefore

$$[eqn]$$

with the corresponding score function $[eqn]$ and the Fisher information matrix $[eqn]$ . For the gene and peak pair, the primary interest is to test $[eqn]$ versus $[eqn]$ , which assesses whether the chromatin accessibility at the peak affects the expression of the gene. Let $[eqn]$ be the restricted maximum likelihood estimate of $[eqn]$ under the null hypothesis. Note that the parameter estimate $[eqn]$ does not depend on the $[eqn]$ . Therefore, for gene g, it remains consistent across all P peaks and only needs to be evaluated once.

For testing $[eqn]$ , PETScan uses the score test $[eqn]$ , where $[eqn]$ denotes the $[eqn]$ th element of $[eqn]$ , and $[eqn]$ represents the $[eqn]$ th entry of the inverse of $[eqn]$ , both corresponding to the parameter $[eqn]$ . Specifically, $[eqn]$ and $[eqn]$ , where $[eqn]$ , $[eqn]$ , and $[eqn]$ = diag( $[eqn]$ ). Leveraging blockwise matrix inversion,

$$[eqn]$$

where $[eqn]$ . After constructing $[eqn]$ by applying a transformation to z, we can extract specific elements without computing the full inverse of a matrix.

Again, for gene g, $[eqn]$ and $[eqn]$ depend only on $[eqn]$ and remain consistent across all P peaks. Consequently, the score test for peak p can be extended easily to all peaks utilizing matrix calculations as described below: let $[eqn]$ , $[eqn]$ , then $[eqn]$ , where $[eqn]$ represents the Hadamard (element-wise) product and $[eqn]$ denotes the element-wise reciprocal of a vector or matrix A. Unlike the Wald and likelihood ratio tests, for a given gene, score tests require fitting the negative binomial model under the null hypothesis only once for all peaks, instead of P full models separately for each peak. This approach facilitates the proposed matrix operations to simultaneously calculate the score tests across all peaks, rather than repeated calculations for one peak at a time. This streamlines computations and significantly enhances computational efficiency.

For studies with finite samples, score tests may be artificially inflated, leading to inaccurate p-values derived from the asymptotic null distribution. To ensure robust inference, we combine an empirical permutation method with genomic control when there is concern about asymptotic results. For gene g, let the original score test statistics across the P peaks be denoted as $[eqn]$ . We then randomly sample the ATAC-Seq data without replacement for L times, and similarly denote the permuted score tests as $[eqn]$ , where $[eqn]$ . Following the genomic control approach, we empirically estimate the inflation factor $[eqn]$ by $[eqn]$ = max(1, $[eqn]$ ) and rescale the original score test statistics $[eqn]$ to $[eqn]$ for $[eqn]$ , which are then compared to $[eqn]$ for p-value calculations.

Results

Simulations

To evaluate PETScan in controlling type I error, we conducted simulations with synthetic data as follows. For each gene $[eqn]$ , we let its expression levels $[eqn]$ follow $[eqn]$ with $[eqn]$ , where the covariate $[eqn]$ . For each peak $[eqn]$ , its log-transformed chromatin accessibility values were simulated as $[eqn]$ . We set the dispersion parameter $[eqn]$ and varied the sample size $[eqn]$ . We also varied the number of permutations $[eqn]$ to investigate the impact of L on type I error control.

We compared the empirical type I error rates and computational efficiency of PETScan against the Wald test (Wald) and likelihood ratio test (LRT), two asymptotically equivalent alternatives to the score test, or PETScan. Because Wald and LRT require full model parameter estimation, they cannot be efficiently implemented using matrix operations across all P chromatin accessibility peaks, making them computationally prohibitive for large scale genome-wide association analysis. Figure 1 presents the empirical type I error rates of PETScan (without permutation and with varying numbers of permutations), Wald and LRT across different sample sizes. As expected, all three tests in general controlled type I error when the sample size was large ( $[eqn]$ ). For smaller sample sizes, the simulation results indicated that genomic control with $[eqn]$ permutations was sufficient to control type I error, although slightly conservative. For all real data analyses, we used $[eqn]$ to ensure robust results.

Left: Empirical type I error rates, and Right: Computational time (minutes) for PETScan (without permutation, L=0, and with L=5,10,20 permutations), Wald and LRT across different sample sizes.

In terms of computational efficiency, the runtime of PETScan scaled approximately linearly with the number of permutations, consistent with the increasing dimensionality of matrix computations. Overall, PETScan was over 1,000 fold faster than both Wald and LRT while achieving comparable performance.

ENCODE human Tissue-Specific data

We applied PETScan to paired RNA-Seq and ATAC-Seq data from 46 healthy individuals across 15 human tissues, obtained from the Encyclopedia of DNA Elements (ENCODE) project (ENCODE Project Consortium 2012), which were accessible at https://www.encodeproject.org/search/?type=Experiment. The RNA-Seq dataset contained 59,526 genes, and we applied filtering based on gene expression counts and selected variable genes, resulting in a final set of 11,637 genes. Similarly, the ATAC-Seq dataset included 805,169 peaks, which were filtered down to 132,670 peaks. We incorporated gender and the top two principal components derived from log-transformed gene expression data as covariates. P-values culminated in assessing associations across 1,548,656,910 gene-peak pairs. Leveraging parallel computation across 64 CPU cores (AMD EPYC 7313), the analysis was completed in approximately 40 minutes with 20 permutations, significantly outperforming Wald tests, which required over 1,000 times longer. We identified 53,113 significant gene-peak pairs with Bonferroni-adjusted p-values below 0.05, suggesting a widespread impact of chromatin accessibility on gene expression. For further analysis, we focused on the top 7,834 most significant gene-peak pairs with Bonferroni-adjusted p-values below $[eqn]$ , as visualized in Figure 2, which was adapted from (Sun 2009) and is available at http://www.bios.unc.edu/∼weisun/software/eMap.pdf.

Significant gene-peak pairs across the genome from ENCODE human tissue-specific data are shown, with peak locations on the X-axis and gene locations on the Y-axis. Different levels of significance are indicated using different colors. The middle panel illustrates the number of genes associated with each peak, while the bottom panel displays the number of peaks linked to each gene.

One notable result involved an ATAC-Seq peak at chr12:124,871,497–124,872,709, which was associated with 47 genes across the genome, represented as a vertical line in Figure 2. Among these genes, HNF1B, EGF, and SLC4A4 were linked to pathways for maturity onset diabetes of the young, pancreatic cancer, and pancreatic secretion, respectively. As shown in Figure 3, chromatin accessibility at this peak was negatively related to the expression of these associated genes. Notably, the pancreas exhibited lower chromatin accessibility at this peak yet higher gene expression, distinguishing it from other tissues. These findings highlighted the role of chromatin accessibility in tissue-specific gene regulation and shed light on the regulatory landscape of genes involved in pancreatic function and disease. For other associated genes, BARX2 encoded the transcription factor Barx2, which regulated GRHL2; in turn, the transcription factor Grhl2 may potentially target multiple downstream genes, including BICDL2, CA12, CLDN4, INAVA, PKP3, RHOV, and TRPV6, as per the Cistrome Data Browser, a comprehensive database for transcriptional and epigenetic regulation studies (Mei et al. 2017, Lambert et al. 2018, Zheng et al. 2019). This regulatory cascade illustrated how chromatin accessibility at a specific peak could influence the coordinated expression of multiple genes across the genome. Such co-association suggested a shared regulatory mechanism, emphasizing the dynamic interplay between chromatin structure and gene regulation.

Scatter plots of chromatin accessibility at chr12:124,871,497–124,872,709 and the expression levels of its associated genes.

In addition, PETScan identified a few genes whose expression levels were associated with a large number of peaks, represented by horizontal lines in Figure 2. The most prominent was the RPL32 pseudogene on chromosome 7, which originated from the retrotransposition of RPL32 mRNA. This association primarily reflected tissue-specific gene expression and chromatin accessibility. Specifically, this gene exhibited higher expression in the pancreas while being lowly expressed in other tissues, while the chromatin accessibility of its associated peaks in the pancreas was either significantly higher or lower compared to other tissues (Figure 4). This finding was consistent with results from the GENCODE project, indicating that the transcription and regulation of pseudogenes were tissue-specific (Pei et al. 2012). A similar pattern was observed for the gene ST8SIA5 on chromosome 18, which showed significantly higher expression in the adrenal gland, along with chromatin accessibility that was either much higher or lower in this tissue compared to others (Figure 4), suggesting a strong tissue-dependent regulatory mechanism. Existing research has shown that the expression of this protein-coding gene is linked to open chromatin regions, particularly in brain tissues (Ehrlich et al. 2024). Our study did not include samples from brain tissues, but our analysis suggested that this highly expressed gene in the adrenal gland was also associated with open chromatin regions in the adrenal gland.

Heatmap of chromatin accessibility of associated peaks for the RPL32 pseudogene (left) and ST8SIA5 (right).

To evaluate our findings, we assessed whether significant gene-peak pairs were validated by previously reported cis-eQTLs from the Genotype-Tissue Expression (GTEx) analysis V8 release (Lonsdale et al. 2013), as cis-eQTLs data are more comprehensive than trans-eQTLs data, even though we investigated both cis- and trans-associations. Two gene-peak pairs overlapped with GTEx cis-eQTLs, where variants from significant variant-gene associations were located within peak regions, both exhibiting higher chromatin accessibility and higher gene expression in the pancreas, setting it apart from other tissues (Figure 5). For instance, variant rs4444903, situated within a peak spanning chr4:109,911,610–109,916,356, explained a fraction of the genetic variance in EGF expression. Overexpression of EGF and its receptor EGFR is frequently observed in pancreatic cancer, leading to uncontrolled cell proliferation and tumor progression, and is associated with poor prognosis (Heby et al. 2020). This finding highlighted the tissue-specific accessibility of this region in modulating EGF expression, thereby contributing to pancreatic functionality and possibly to the pathological mechanisms underlying pancreatic cancer. Similarly, variant rs66817580 was linked to FAM3B expression, which was located within a peak at chr21:41,315,678–41,317,753. FAM3B, also known as PANDER (pancreatic-derived factor), is primarily expressed in the pancreas, particularly within beta cells. This gene is instrumental in the regulation of glucose metabolism and insulin secretion, and influences beta cell functionality and survival, making it relevant to diabetes and pancreatic disorders (Dong et al. 2023). Further investigation into this regulatory mechanism could provide insights into beta cell dysfunction in diabetes and uncover potential therapeutic targets. PETScan demonstrated its effectiveness in uncovering gene regulatory mechanisms, particularly pancreatic-specific associations between chromatin accessibility and targeted gene expression.

Scatter plots of significant gene-peak pairs from ENCODE human tissue-specific data validated by GTEx cis-eQTLs.

10x genomics embryonic mouse brain data

We next conducted PETScan on publicly available 10x Genomics embryonic mouse brain data, which included 4,881 cells from the fresh cortex, hippocampus, and ventricular zone of the embryonic mouse brain at day 18 (10x Genomics 2021). The scRNA-Seq and scATAC-Seq data were processed following the TRIPOD vignette (Jiang et al. 2022), with cell-type labels transferred from an independent scRNA-Seq reference (La Manno et al. 2021) using SAVERCAT (Huang et al. 2020). We retained 3,851 cells with consistent transferred labels from four major cell types: neuroblast, forebrain GABAergic neuron, cortical or hippocampal glutamatergic neuron, and glioblast. To address the sparsity issue in single-cell data, we partitioned cells into 84 metacells by clustering based on their scRNA-Seq profiles and assigned cell types to metacells according to the majority cell types within each cluster.

After filtering, we retained 5,935 genes and 42,048 peaks. Using parallel computation across 64 CPU cores, we analyzed 249,554,880 gene-peak pairs in just 7 minutes with 20 permutations, adjusting for the first two principal components of the log-transformed gene expression data. PETScan demonstrated computational efficiency that was three orders of magnitude faster than Wald tests. We identified 26,991 significant gene-peak pairs with Bonferroni-adjusted p-values below 0.05, providing insights into the regulatory landscape of the embryonic mouse brain. For further analysis, we focused on the top 6,858 significant gene-peak pairs with Bonferroni-adjusted p-values below $[eqn]$ , as illustrated in Figure 6.

Significant gene-peak pairs across the genome from 10x Genomics embryonic mouse brain data are shown, with peak locations on the X-axis and gene locations on the Y-axis. Different levels of significance are indicated using different colors. The middle panel illustrates the number of genes associated with each peak, while the bottom panel displays the number of peaks linked to each gene.

We performed Gene Ontology (GO) and KEGG enrichment analyses on genes from these pairs to identify over-represented GO terms and KEGG pathways using clusterProfiler. The top enriched GO terms included forebrain development, telencephalon development, regulation of neurogenesis, gliogenesis, central nervous system neuron differentiation, and neuron projection guidance. These terms highlighted key processes involved in the development and specialization of the nervous system, particularly regarding the formation of the brain and its functional components. These findings aligned with the fact that embryonic mouse brain cells were collected during the transition phase between neurogenesis and gliogenesis, further supporting the validity of our results. The top enriched KEGG pathways grouped into two main categories: (1) nervous system: including glutamatergic synapse, GABAergic synapse, and synaptic vesicle cycle, all critical for neurotransmission and communication between neurons; (2) signal transduction: comprising calcium signaling pathway and neuroactive ligand-receptor interaction, which played essential roles in regulating cellular responses to external stimuli and modulating various physiological processes. These results were consistent with the characteristics of glutamatergic and GABAergic neurons included in embryonic mouse brain cells.

Furthermore, an ATAC-Seq peak at chr7:54,633,595–54,635,192 was associated with 27 genes across the genome. Notably, open chromatin at this peak was linked to lower expression of Nrxn3, Gad1, and Gad2 (Figure 7), all of which were involved in the GABAergic synapse. Conversely, it was correlated with higher expression of Neurod6, Zbtb18, Nfix, Lhx2, and Id2 (Figure 8), all of which were essential for telencephalon development and forebrain development. The upregulated or downregulated genes exhibited functional coherence, indicating that changes in chromatin accessibility may co-regulate genes involved in related biological processes.

Scatter plots of chromatin accessibility at chr7:54,633,595–54,635,192 and the expression levels of its negatively associated genes.

Scatter plots of chromatin accessibility at chr7:54,633,595–54,635,192 and the expression levels of its positively associated genes.

To validate our findings, we compared significant gene-peak pairs with previously reported chromatin interactions, using PLAC-Seq data of the mouse fetal forebrain (Zhu et al. 2019). Our analysis identified nine significant gene-peak pairs in PLAC-Seq data, all exhibiting positive correlations between chromatin accessibility and gene expression, predominantly in GABAergic and glutamatergic neurons (Figure 9). We also performed PETScan to evaluate cis-associations for all genes and their corresponding peaks within 1,000 kb of the transcription start site, identifying 4,211 significant gene-peak pairs with Benjamini & Hochberg (BH) adjusted p-values below 0.05, of which 756 (18.0%) were validated by PLAC-Seq data. These results showed the biological relevance of our approach in capturing chromatin interactions that regulate gene expression in mouse brain development.

Scatter plots of significant gene-peak pairs from 10x Genomics embryonic mouse brain data validated by PLAC-Seq data.

Discussion

PETScan is a computationally efficient tool for providing a comprehensive perspective on tissue- or cell-type-specific gene regulation and its implications for diverse biological processes and diseases. Our findings align with previously reported cis-eQTLs or PLAC-Seq data, validating the robustness and biological relevance of the identified gene-peak pairs. This concordance underscores *PETScan’*s potential to effectively bridge chromatin accessibility and gene expression data.

Beyond cis-associations, PETScan identifies trans-associations, capturing interactions between regulatory elements and distal target genes to uncover long-range regulatory mechanisms that may not be detectable using traditional approaches. Genome-wide associations identified by PETScan offer important insights into regulatory mechanisms, providing experimental validation guidelines for confirming the functional roles of specific open chromatin regions in regulating the transcriptomic expressions of their targeted genes.

Although PETScan is developed for performing billions of association tests between high-dimensional RNA-Seq and ATAC-Seq data, it is also versatile for other paired omics data, provided that negative binomial models offer an appropriate framework. For example, PETScan can be a viable approach for eQTL analysis when gene expression profiles are measured using RNA-Seq data, where the normality assumption of Matrix eQTL is no longer appropriate. Additionally, PETScan can facilitate genome-wide association analysis between RNA-Seq and ChIP-Seq data, evaluating how transcription factor binding or histone modifications impact gene expression. Similarly, analysis of paired RNA-Seq and Bisulfite-Seq data with PETScan can provide insights into how DNA methylation regulates gene expression.

Less obviously, PETScan can be adopted for differential gene expression analysis of RNA-Seq with a given factor, such as a binary outcome. When the sample size is limited, necessitating a large number of permutations, say P, to ensure valid type I error control, the P sets of permuted outcomes can be treated analogously to chromatin accessibility measurements across P peaks, on which PETScan can be applied to derive an empirical null distribution. Genomic control can similarly be implemented through permutations to address potential limitations associated with small sample sizes, thereby enhancing the reliability and robustness of the association findings. Similar research direction can be found in (Barry et al. 2025).

A notable feature of PETScan is its reliance on score tests, which focus solely on assessing the presence of associations without estimating effect sizes and directionality (positive or negative) as traditional regression-based methods do. While this approach offers computational efficiency and stability, it limits the ability to determine whether a peak positively or negatively influences gene expression, a crucial factor in elucidating underlying biological mechanisms. Nevertheless, PETScan remains a powerful tool for large-scale screening in genome-wide association analysis, particularly when computational efficiency is paramount. By identifying candidate interactions for further validation, PETScan provides a foundation for subsequent in-depth investigations.

The use of negative binomial models provides the advantage of flexible covariate adjustment compared to simple correlations. Currently, PETScan applies the same set of covariates to all genes. However, some covariates may influence the relationship between chromatin accessibility and gene expression in a gene-specific manner, such as the cis-effect of DNA methylation and transcript factor expression levels. PETScan can be easily modified to include gene-specific covariates when relevant data is available, enabling more accurate inference.

When multiple paired RNA-Seq and ATAC-Seq samples are collected from the same subject, such as from different biological time points, the correlations between these samples must be accounted for. Our future plan is to extend PETScan from negative binomial models to negative binomial mixed models (Tsonaka and Spitali 2021, Sun et al. 2016, Rivellese et al. 2022), presenting an even more computationally daunting task. In addition, negative binomial models in PETScan may not adequately account for counts with excess zeros commonly observed in other omics data, such as single-cell RNA-Seq data and microbiome data, suggesting the need for zero-inflated Poisson or zero-inflated negative binomial models, which could be another useful extension direction for PETScan.

In conclusion, PETScan offers a computationally efficient and biologically meaningful framework for assessing genome-wide associations between chromatin accessibility and gene expression. It holds promise for advancing the study of gene regulation and guiding experimental validation efforts, ultimately contributing to a deeper understanding of the regulatory mechanisms underlying gene expression.