Tutorials › Statistical Methods for Single-Cell RNA Sequencing › Linking Accessible Peaks to Gene Expression

CITE-seq, ATAC and Multiome Data · Tutorial 612 of 700

Linking Accessible Peaks to Gene Expression

You will learn how to calculate and test peak–gene correlations, incorporate genomic distance into candidate selection and null comparisons, and interpret links cautiously when cells share latent states.

Advanced 11 min read

What You'll Learn

  • Define peak–gene candidates using genomic coordinates and a prespecified distance rule.
  • Calculate Pearson correlations from paired accessibility and expression measurements.
  • Derive and apply the classical test of zero correlation, including its assumptions.
  • Distinguish distance-based eligibility and calibration from evidence of association.
  • Diagnose correlations induced by shared cell states and use conditional analyses cautiously.
  • Interpret statistical links as associations rather than proof of regulatory effects.

Why Link Peaks to Genes?

An accessible chromatin peak may lie near a gene whose RNA abundance varies across the same cells. A correlation between peak accessibility and gene expression can prioritize a possible regulatory link, but it does not by itself show that the peak regulates the gene. The correlation may reflect a shared cell state, technical effects, or other features of the cells rather than a direct peak-to-gene relationship.

The basic unit of analysis is a candidate pair: one genomic peak and one gene. For each cell in which both modalities were measured, the pair contributes an accessibility value and an expression value. The introductory material on constructing gene activity scores concerns aggregating accessibility near genes; here the question is narrower: whether the measured accessibility of a particular peak is associated with the measured expression of a particular gene across paired cells.

Definition: For a candidate peak \(j\) and gene \(g\), let \(A_{ij}\) be the accessibility measurement for peak \(j\) in cell \(i\), and let \(Y_{ig}\) be the expression measurement for gene \(g\) in that same cell. A peak–gene association summarizes how \(A_{ij}\) and \(Y_{ig}\) vary together across the analyzed cells. A positive association means that larger values tend to occur together; it does not establish a causal regulatory effect.

Define Candidate Links Before Testing

Distance is useful first as a rule for defining which pairs will be considered. A common distance measure is the number of base pairs between a peak and a gene’s transcription start site (TSS). For a peak interval, one can define the distance from the TSS to the nearest point of the interval, or use another clearly specified coordinate convention. The convention should be applied consistently, especially for peaks that overlap a TSS.

For example, an analysis might consider a peak eligible for a gene if its interval is within 100 kilobases of that gene’s TSS. The window is an analysis choice, not a biological boundary: a peak outside it might still be relevant, and a nearby peak is not necessarily relevant. A peak can be eligible for more than one gene, and a gene can have multiple eligible peaks. These facts should be retained when interpreting results rather than hidden by assigning each peak to only its nearest gene without justification.

Once candidate pairs are defined, measure the association across the same paired cells. If accessibility and RNA measurements are highly skewed or dominated by zeros, a transformation or normalization may be needed before using a linear correlation. For example, an analysis may use normalized, log-transformed measurements. The choice changes what the correlation describes; it should be made before examining candidate-specific results and applied consistently.

Worked Example: Calculating a Peak–Gene Correlation

Suppose five hypothetical paired cells have normalized accessibility values \(A\) and normalized expression values \(Y\) for one peak–gene pair. The peak is 18 kilobases from the gene’s TSS, and the prespecified candidate window is 100 kilobases.

CellAccessibility \(A\)Expression \(Y\)
101
212
322
434
546

The pair meets the distance eligibility rule because \(18\) kilobases is less than \(100\) kilobases. This says only that the pair is included for analysis. The means are \(\bar{A}=2\) and \(\bar{Y}=3\). The centered values are \((-2,-1,0,1,2)\) for accessibility and \((-2,-1,-1,1,3)\) for expression. Thus the sum of cross-products is \(12\), while the sums of squared deviations are \(10\) and \(16\).

$$ r=\frac{\sum_i(A_i-\bar{A})(Y_i-\bar{Y})} {\sqrt{\sum_i(A_i-\bar{A})^2\sum_i(Y_i-\bar{Y})^2}} =\frac{12}{\sqrt{10\cdot16}} =0.9487. $$

The observed values have a strong positive linear association in these five cells. The small, hypothetical data set does not establish that the association is reproducible, and distance eligibility does not turn the correlation into evidence of regulation.

Correlation Testing and Its Assumptions

Pearson’s sample correlation \(r\) describes the direction and strength of a linear association. Under the classical test of \(H_0:\rho=0\), where \(\rho\) is the population correlation, the test statistic is

$$ t=r\sqrt{\frac{n-2}{1-r^2}}, \qquad df=n-2. $$

For the usual exact \(t\)-test interpretation, the paired observations must be independent across cells, and the paired measurements should follow a bivariate normal distribution under the model. In practice, the test is also sensitive to strong outliers, and Pearson correlation describes linear rather than general nonlinear association. With enough independent observations it can be used as an approximation under less restrictive distributions, but this does not resolve dependence between cells.

These assumptions require particular care in single-cell data. Cells from the same donor are not automatically independent biological replicates, and cells with related states may have correlated measurements. A large cell count does not necessarily provide an equally large number of independent units of biological evidence. Inspect scatterplots, distributions, and patterns by donor or cell state; a single extreme cell or a mixture of states can dominate the apparent association.

Worked Example: Testing a Correlation Under the Classical Model

In a hypothetical analysis of \(n=20\) independent paired cells, a candidate pair has sample correlation \(r=0.30\). Consider a two-sided test of \(H_0:\rho=0\) against \(H_a:\rho\ne0\) at \(\alpha=0.05\). For this calculation, suppose the classical assumptions are reasonable.

Substitute the observed correlation and sample size into the statistic:

$$ t=0.30\sqrt{\frac{20-2}{1-0.30^2}} =0.30\sqrt{\frac{18}{0.91}} \approx1.334,\qquad df=18. $$

The two-sided \(p\)-value is approximately \(0.199\). Because \(0.199>0.05\), this test does not provide sufficient evidence against zero population correlation at the chosen significance level. It does not prove that the true correlation is zero. The test’s classical interpretation depends on its assumptions; if the 20 cells share substantial state or donor dependence, this \(p\)-value may not be appropriately calibrated.

What Distance-Aware Testing Means

Genomic distance can enter an association analysis in more than one way, and those uses answer different questions. A prespecified distance window defines eligible candidates. Distance can also be used when comparing a candidate’s correlation with a background of other peak–gene pairs. Neither use changes the arithmetic definition of \(r\), and neither automatically corrects a correlation test for dependent cells.

For a distance-aware background comparison, one can compare a candidate pair with non-target or otherwise suitable peak–gene pairs at similar genomic distances. The comparison can also account for features that affect how correlations are measured, such as expression level or detection rate. If a candidate’s correlation is larger than those of suitable distance-matched pairs, that is evidence that it stands out relative to that chosen background. It is not the same as testing \(H_0:\rho=0\), and its interpretation depends on how the background pairs were selected.

A different approach is to calculate an association test for each eligible pair and then account for the many tests. If \(m\) candidate pairs are tested, false discovery rate procedures such as Benjamini–Hochberg adjustment can help control the expected proportion of false discoveries among rejected hypotheses under their applicable conditions. A small adjusted \(p\)-value is not a measure of regulatory effect size, and distance eligibility alone does not account for multiple testing.

1
Specify the candidate set.
Choose the peak-to-gene distance definition and eligibility rule before inspecting the correlations.
2
Calculate pairwise associations.
Use paired accessibility and expression measurements, with transformations and cell inclusion criteria stated in advance.
3
Choose the inferential comparison.
Use a correlation test for a zero-correlation question, or a distance-matched background for a relative-to-background question. State which question is being answered.
4
Check dependence and multiplicity.
Assess cell-state and donor structure, then account for the number of tested pairs using an appropriate procedure.
5
Interpret as a candidate link.
Report the association, distance, uncertainty or adjusted evidence, and relevant limitations without claiming causality.

Shared Cell States Can Create Correlation

Suppose cells from one state tend to have both low accessibility and low expression, while cells from another state tend to have both high accessibility and high expression. Pooling the states can produce a positive correlation even if accessibility and expression are negatively associated within each state. In that case, the pooled correlation partly describes differences between states, not a consistent within-state peak–gene relationship.

One diagnostic is to plot or summarize the pair within known cell types, donors, or other relevant groups. Another is to estimate the association after adjusting for measured covariates or after calculating residuals from models that represent those covariates. The resulting conditional association answers a different question from the pooled correlation. It also depends on whether the covariates adequately capture the shared structure. An unmeasured continuous cell state can remain a source of dependence even after adjustment for broad cell-type labels.

Worked Example: A Pooled Link Reverses Within States

Consider eight hypothetical cells in two equally sized latent states. In state A, accessibility values are \(0,1,2,3\) and expression values are \(3,2,1,0\). In state B, accessibility values are \(8,9,10,11\) and expression values are \(11,10,9,8\). Within each state, expression decreases exactly as accessibility increases, so each within-state correlation is \(-1\).

Across all eight cells, both overall means are \(5.5\). The centered accessibility values are \(-5.5,-4.5,-3.5,-2.5,2.5,3.5,4.5,5.5\). The centered expression values are \(-2.5,-3.5,-4.5,-5.5,5.5,4.5,3.5,2.5\). Their cross-product sum is \(118\), and each centered sum of squares is \(138\). Therefore:

$$ r_{\text{pooled}}=\frac{118}{\sqrt{138\cdot138}} =\frac{118}{138}=0.8551. $$

The pooled association is positive even though the within-state associations are perfectly negative. The upward shift in both measurements between states dominates their negative relationship within states. A pooled positive peak–gene correlation in real data could likewise reflect a shared state pattern. State-stratified plots or a conditional analysis are needed before interpreting it as a within-state association.

Diagnostics, Alternatives, and Limits

A useful diagnostic sequence is to examine the paired scatterplot, label or facet points by donor and cell state, and check whether the association is driven by a few cells or by differences between groups. Compare the pooled correlation with correlations within sufficiently represented groups. Check whether the pattern is plausibly linear; a near-zero Pearson correlation can coexist with a strong nonlinear relationship. If measurements contain many zeros or are heavily skewed, inspect how the transformation and detection pattern affect the result rather than treating a correlation as automatically robust.

The classical Pearson test is simple and interpretable when its assumptions are reasonable. Rank correlations can be useful for monotone associations that are not linear, but they do not remove cell dependence or confounding. A permutation test can avoid a particular parametric reference distribution, but its validity still requires a defensible exchangeability scheme. For example, permuting RNA measurements across all cells may destroy cell-state structure; permuting only within suitable donor or state groups tests a more conditional null, provided observations are exchangeable within those groups. These alternatives answer different questions and are not universal fixes.

If there are many candidate links, also examine the number of tests, the multiplicity procedure, and whether the candidate set was defined independently of the observed association. Report the correlation and its direction, the genomic distance and coordinate convention, the cells or groups analyzed, and the testing or background comparison. Statistical association can prioritize candidates for further study, but correlation alone cannot distinguish direct regulation from shared state, indirect pathways, or technical effects.

Common Mistake

Common mistake: “This peak is within 100 kilobases of the gene, and its correlation is significant, so the peak regulates the gene.” The distance rule only makes the pair eligible, while a significant correlation supports an association under the chosen model and assumptions. The correction is to describe it as a candidate peak–gene link, investigate cell-state and donor structure, account for multiple testing, and avoid causal language unless supported by additional evidence.

Check Your Understanding

Answer each question, then compare with the brief solutions.

  1. A peak is 140 kilobases from a gene’s TSS. The prespecified candidate window is 100 kilobases. Is the pair eligible, and what does that decision imply about association?
  2. For \(n=20\) independent paired cells with \(r=0.30\), calculate the classical test statistic and degrees of freedom.
  3. Why can a pooled correlation be positive when within-state correlations are negative?
  4. What is the distinction between a zero-correlation test and comparison with a distance-matched background?
  5. Name two reasons a small correlation-test \(p\)-value may not justify a claim of direct regulation.

Solutions. (1) No; \(140>100\) kilobases, so the pair is outside the stated eligibility window. Eligibility is not itself evidence of association. (2) \(t=0.30\sqrt{18/0.91}\approx1.334\), with \(df=18\); this has its classical interpretation only when the test assumptions are appropriate. (3) Between-state shifts can dominate the within-state pattern, as in the worked example. (4) The zero-correlation test asks whether the population correlation is zero under its model; a distance-matched comparison asks whether the candidate is unusual relative to selected pairs at similar distances. (5) Shared cell-state or donor structure can invalidate the assumed independence, and many tested pairs require multiplicity control; even a well-calibrated association test does not establish causality.

Clinical Biostats