Worked Multiome Analysis From Peaks to Regulatory Links

From Matched Modalities to a Regulatory Hypothesis

A multiome experiment measures RNA and chromatin accessibility in the same cells. Matching modalities cell by cell is useful, but it does not by itself establish that an accessible peak regulates a particular gene. A defensible analysis has several distinct steps: summarize accessibility near genes, define candidate peak–gene pairs, estimate associations using an appropriate independent unit, account for testing many pairs, and weigh statistical evidence alongside biological context.

The earlier material on constructing gene activity scores and linking accessible peaks to gene expression establishes the basic quantities used below. This worked analysis emphasizes how to connect those quantities into a reproducible chain of evidence. All data and results here are hypothetical.

Build the Analysis Around the Experimental Design

Begin by checking that RNA and ATAC barcodes refer to the same cells, that genomic coordinates use the same reference assembly, and that the peak-to-gene assignment rule is specified before testing. Restricting analysis to a relevant cell type or state can make an association easier to interpret, but the restriction should be chosen for scientific reasons rather than because it gives a smaller p-value. Record donors, batches, cell types, and any other design variables that could affect both modalities.

For a simple activity score, let \(A_{ij}\) be the accessibility measurement for cell \(i\) at peak \(j\). Let \(\mathcal{P}_g\) be the selected set of peaks overlapping the promoter or gene body of gene \(g\). An unweighted score is the sum of those measurements. In practice, normalization for library depth and transformation may be needed before comparing cells. The exact genomic window, peak inclusion rule, and normalization are part of the score’s definition and should be reported.

For a candidate regulatory link, it is often clearer to test accessibility at one peak against expression of one gene than to test a broad activity score alone. Define the candidate set in advance using genomic proximity, prior biological knowledge, or another stated rule. A peak near a gene is a plausible candidate, not necessarily its target; distal peaks may act on more distant genes, and proximity can misassign targets.

Paired measurements mean the ATAC and RNA values belong to the same cell. They do not mean that thousands of cells from a few donors are thousands of independent biological replicates. If the target of inference is a population of donors, a practical strategy is to aggregate within donor and a prespecified cell type, then model the resulting donor-level values. This reduces cell-level pseudoreplication. When there are enough donors and the design warrants it, a model can instead account explicitly for donor-level clustering and relevant covariates.

Estimate a Candidate Association Across Donors

Consider a simple linear model for one candidate peak and one gene. Let \(X_d\) be normalized peak accessibility and \(Y_d\) be normalized expression for donor \(d\), each summarized within the same prespecified cell type. The model is \(Y_d=\beta_0+\beta_1X_d+\varepsilon_d\). The slope \(\beta_1\) is the expected change in the chosen expression scale for a one-unit increase in accessibility, under this model. It describes association, not a causal effect.

The usual small-sample t test for a slope assumes independent donor-level errors, a linear mean relationship, constant error variance, and normally distributed errors for exact inference. These assumptions apply to the model residuals, not to the raw single-cell counts. Check a scatterplot for curvature and influential donors, and inspect residuals for unequal spread or extreme deviations. With very few donors, such diagnostics have limited power; uncertainty should be reported rather than concealed.

Adjust for Many Candidate Links and Integrate Evidence

A workflow commonly tests many peak–gene pairs. Testing each at level \(\alpha=0.05\) without adjustment increases the chance of false discoveries across the set. The Benjamini–Hochberg procedure controls the false discovery rate under its standard conditions: sort the \(m\) p-values from smallest to largest, find the largest rank \(i\) satisfying \(p_{(i)}\le(i/m)\alpha\), and reject hypotheses at that rank and below. An equivalent reporting practice is to give adjusted p-values, often called q-values in this context, computed as the monotone adjusted values \(\min_{k\ge i}\{m p_{(k)}/k\}\), capped at 1. State the tested universe: changing which pairs enter the family changes the adjustment.

Statistical evidence is only one component of prioritization. For each retained candidate, consider the direction and size of the association, its uncertainty, consistency across donors, genomic distance, whether the peak contains a relevant motif, and whether expression of the proposed target is measurable in the analyzed cells. Motif presence can support a regulatory hypothesis, but it does not prove that the motif is occupied or that the peak targets the gene. Likewise, proximity and a small adjusted p-value are not independent causal validation.

Diagnostics, Alternatives, and Limitations

Potential warning signs include a relationship driven by one donor, a trend that disappears after accounting for batch, residual curvature, strong heteroskedasticity, and associations that change substantially when cell-type definitions or peak boundaries change. Report donor-level points whenever possible. Also report the number of donors, the number of cells contributing to each summary, the normalization and transformation, the candidate-generation rule, and the multiple-testing family. More cells per donor can improve measurement precision, but do not replace independent donors for population-level inference.

The ordinary linear model is interpretable and convenient when donor-level summaries are reasonably well behaved and its assumptions are plausible. A rank-based association can be less sensitive to outliers and nonlinearity in scale, though it does not solve donor dependence or confounding. A donor-aware mixed model can retain cell-level observations while accounting for within-donor dependence and covariates, but requires an appropriate model and enough independent donors to estimate its structure. These approaches answer related but not identical questions; choose and report the model based on the design rather than selecting whichever produces the smallest p-value.

A gene activity score depends on its chosen windows and normalization, and may miss regulatory elements outside those windows. Association can reflect shared cell state, donor composition, sequencing depth, or other common influences rather than direct regulation. Tests across many candidates require multiplicity control, and sparse measurements can make effects imprecise. Functional perturbation or other independent validation is needed to make a stronger causal claim.

A Compact Analysis Sequence