Why Simulate Integration Data?
A benchmark using observed datasets can show how methods behave on realistic data, but the biological and technical components of those data are not usually known separately. A controlled simulation addresses a different question: if the underlying cell states and batch effects were known, would an integration method remove the intended technical differences while retaining the intended biology? The answer depends on the simulation design. A simulation that omits sampling imbalance, count variation, or states found in only one batch may make a method look reliable in circumstances where it would fail.
The goal is not to create one “realistic” synthetic dataset. It is to specify a generative model, record its truth, and vary important design features in a controlled way. This tutorial constructs such a model for single-cell counts. Its central distinction is between a batch shift conditional on cell state and a difference in the mix of states sampled in each batch.
A Generative Model With Explicit Truth
Let \(Y_{gi}\) be the count for gene \(g\) in cell \(i\), let \(b_i\) be that cell’s batch, and let \(z_i\) be its cell state. Let \(s_i\) be a positive cell-specific size factor. A useful starting point is a negative-binomial count model:
Here \(\alpha_g\) is a gene-specific baseline, \(\theta_{gz}\) is the biological expression difference associated with state \(z\), and \(\delta_{gb}\) is a gene-specific batch shift. The optional interaction \(\gamma_{gzb}\) represents a batch shift that varies across states. Under this parameterization, the negative-binomial variance is \(\mu_{gi}+\phi_g\mu_{gi}^2\). Counts are random even when all the model parameters are fixed.
For each batch \(b\), specify state probabilities \(\pi_{bz}=P(z_i=z\mid b_i=b)\), with probabilities summing to one within that batch. Drawing states from different probability distributions creates sampling imbalance. Shared states have positive probability in more than one batch. A batch-specific state has positive probability in one batch and zero in the others. In practice, the simulation should retain the state label, batch label, size factor, and generating parameters for every cell.
The simple additive model is not the only useful design. Gene-specific shifts allow some genes to be unaffected and others strongly shifted. State-by-batch interactions test whether a single correction per gene and batch is too restrictive. These are deliberate design choices: the benchmark should state which terms are present, rather than treating every difference between batches as the same kind of effect.
Separate Batch Shifts From Sampling Imbalance
Even if the conditional expression of every shared state is identical across batches, batch-level averages can differ when the batches contain different proportions of those states. The decomposition is exact for expected log-expression values in the additive model. Write \(m_{gbz}=\alpha_g+\theta_{gz}+\delta_{gb}\) for a gene’s mean log-expression in batch \(b\) and state \(z\). Then
Thus, the difference between two batches’ mixture averages equals the difference in their conditional batch shifts plus the difference in their state-composition averages. This identity explains why a marginal batch comparison alone does not reveal how much of the difference is technical.
Worked Example: Unequal State Frequencies Mimic a Batch Shift
Suppose one gene has biological log-expression \(1\) in state S1 and \(3\) in state S2. Batch A contains 80% S1 and 20% S2; batch B contains 20% S1 and 80% S2. Assume the conditional batch shift is zero in A and \(0.5\) in B. Calculate the expected mixture averages and distinguish the composition difference from the batch shift.
For batch A, the expected mixture average is \(0.8(1)+0.2(3)=0.8+0.6=1.4\). For batch B, before adding its batch shift, it is \(0.2(1)+0.8(3)=0.2+2.4=2.6\). Adding B’s shift gives \(2.6+0.5=3.1\). The observed expected difference is \(3.1-1.4=1.7\).
Of that \(1.7\), \(1.2\) comes from the state proportions: \(2.6-1.4=1.2\). The remaining \(0.5\) is the conditional batch shift. If the two states are equally weighted in both batches, the biological average is \((1+3)/2=2\); the corresponding shifted average in B is \(2.5\). A benchmark should therefore retain the state labels and proportions rather than interpreting the whole mixture difference as a technical effect.
Worked Construction: From Parameters to Counts
Worked Example: A Known Shift With Count-Sampling Variation
Consider one cell, one gene, and a shared state whose target log-expression is \(2.0\). Set the known batch contribution to \(0.9\), and let the cell size factor be \(s_i=1.5\). The log mean count per unit size factor is \(2.0+0.9=2.9\), so the mean count is \(1.5\exp(2.9)\approx27.26\). Suppose a hypothetical draw from the count model is \(Y=31\). Calculate the size-factor-adjusted log count and subtract the known shift.
The observed size-factor-adjusted log count is \(\log(31/1.5)=\log(20.6667)\approx3.029\). Subtracting the known batch contribution gives \(3.029-0.9=2.129\). This is about \(0.129\) above the biological target of \(2.0\), a possible consequence of count-sampling variation. It is not a failure of the simulation’s truth: the known target describes the generating mean, not a guarantee that each random count equals its expectation.
This calculation is illustrative rather than a claim that log-normalized counts are an unbiased estimate of log-expression. In particular, applying a logarithm to a random count does not generally yield an unbiased estimate of the logarithm of its mean. The simulation records the generating parameters as truth; it does not define truth by the noisy observed count.
Represent Shared and Batch-Specific States
A sound benchmark includes both overlap and non-overlap when the intended use includes those conditions. Shared states make conditional comparisons possible: a state observed in two batches provides a setting in which the method can be asked to align corresponding biology while removing a known shift. A batch-specific state tests a different requirement. It has no observed counterpart to anchor it in another batch, so the simulation must not label its distinct biology as a shift that should be erased.
Worked Example: A Batch-Specific State Is Not a Batch Effect
For a gene, let the biological log-expression be \(1\) in shared state S1, \(3\) in shared state S2, and \(5\) in state U. State U occurs only in batch B. Batch A has proportions 80% S1 and 20% S2. Batch B has proportions 50% S1, 30% S2, and 20% U. Let B have a conditional shift of \(0.4\), with no shift in A. Calculate the mixture averages, then identify the valid shared-state comparison.
Batch A’s expected mixture average is \(0.8(1)+0.2(3)=1.4\). Batch B’s biological mixture average before its shift is \(0.5(1)+0.3(3)+0.2(5)=0.5+0.9+1.0=2.4\). With its shift, the expected mixture average is \(2.4+0.4=2.8\), so the marginal difference is \(2.8-1.4=1.4\).
For the shared states, equal weighting gives an A average of \((1+3)/2=2\) and a shifted B average of \(((1+0.4)+(3+0.4))/2=2.4\). The shared-state difference is the \(0.4\) batch shift. State U’s biological value of \(5\) is a distinct state, not evidence that B has a larger technical shift. There is no batch A observation of U that could identify a cross-batch counterpart from the data alone.
This example also exposes an identifiability limit. If a state occurs in only one batch, state and batch are confounded for that state in the observed data. A simulation can record their separate generating values because the designer set them, but an integration method cannot recover that separation from the observations without additional assumptions or information. A benchmark should test whether a method preserves unique biology; it should not imply that the method can infer an unobserved match with certainty.
Designing a Useful Simulation Study
Choose a baseline design that represents the question, then vary one or a few features at a time. This makes it possible to tell which conditions drive a change in performance. A practical sequence is:
Specify genes, shared states, state-specific expression values, and any states restricted to a batch. Record the cell-state label for every generated cell.
Specify the batch shifts and state probabilities separately. Include balanced and imbalanced state proportions if composition is part of the target conditions.
Draw cell-specific size factors and counts from the stated model. Save generating parameters, random seeds, labels, and the exact simulation settings.
Use independent simulation replicates for each setting. Variation across random draws is part of the evidence; one especially favorable draw is not a reliable benchmark.
The count model’s assumptions should be checked against the purpose of the simulation. The negative-binomial distribution is a model for nonnegative integer counts with variance exceeding the mean under \(\phi_g>0\); the model above assumes its stated mean and dispersion conditional on state, batch, and size factor. The cell states are drawn according to the specified batch-specific probabilities. If the simulation assumes independent cell draws, that assumption should be stated. If cells are clustered by donor or sample, generate that dependence explicitly rather than treating all cells as independent replicates.
The assumed size factors, dispersion, gene-expression distributions, state frequencies, and batch-shift magnitudes should be plausible for the benchmark’s intended setting, but no single parameter choice is universally realistic. Check generated data for empty or nearly empty states, extreme library sizes, implausible count ranges, and accidental parameter combinations that erase the intended overlap. Such checks assess whether the simulation implements its design; they do not establish that the design captures every feature of a real experiment.
Alternatives, Diagnostics, and Limitations
| Approach | What it provides | Main limitation |
|---|---|---|
| Fully parametric simulation | Direct control of known biological and technical components, including count noise and sampling imbalance. | Results depend on the chosen generative model; omitted features cannot challenge a method. |
| Semi-synthetic perturbation | Begins with observed data and adds a specified shift or resampling scheme, retaining some observed structure. | The added perturbation has known truth, but all biological and technical structure in the original data does not thereby become known. |
| Real-data benchmarking | Tests methods on observed complexity and practical data-processing choices. | The correct adjustment for each cell is generally not known directly. |
These approaches answer complementary questions. A parametric simulation makes the data-generating truth explicit, while semi-synthetic and real-data studies can provide forms of realism that a simple model lacks. Agreement across designs is more informative than relying on one simulation family. Earlier material on benchmark design provides the broader rationale for using multiple evidence sources; here the emphasis is on making the truth in a synthetic design explicit and auditable.
A simulation is going wrong when its implementation contradicts its design. For example, state probabilities may fail to sum to one, a nominally shared state may be absent from a batch after sampling, or a purportedly constant shift may accidentally vary by state. Check simulated conditional means against the generating equation, tabulate realized state counts by batch, and confirm that every truth label matches its intended role. Record both intended probabilities and realized counts: finite samples need not reproduce the probabilities exactly.
A second limitation is that a method can be favored by the data-generating assumptions. If the simulation has only additive shifts and the evaluated method assumes additive shifts, performance in that setting does not establish robustness to state-dependent shifts. Conversely, a method may perform poorly under a deliberately misspecified model without being poor in all settings. Vary effect sizes, imbalance, noise, and overlap, and report which conditions each conclusion concerns. The later work on benchmark metrics can then summarize correction and preservation against these known targets; the simulation itself must first define those targets unambiguously.
Common Mistakes
- Calling every between-batch difference a batch effect. The correction is to distinguish conditional shifts from differences in state proportions, as shown by the mixture decomposition.
- Expecting a simulated observation to equal its generating mean. Counts are random. Compare repeated draws or summaries with the known target, rather than demanding exact equality cell by cell.
- Removing a state because it appears in only one batch. A batch-specific state can be genuine biology. Its lack of a counterpart creates an identification limit, not proof that it should be corrected away.
- Reporting only one favorable random seed. Repeat the simulation with independent random draws and preserve the settings and seeds so results can be reproduced.
Check Your Understanding
Use the model and examples above to answer each question. Short solutions follow.
- Two batches have identical state-specific expression, but different state proportions. Can their marginal average expressions differ? Why?
- A state occurs only in batch B and has higher expression than all shared states. Does this alone show that batch B has a positive technical shift?
- With biological log-expression \(1.5\), batch shift \(-0.2\), and size factor \(2\), what is the expected count under the model?
- Why should a simulation study use independent replicates rather than report a result from one random seed?
Solutions. (1) Yes. The mixture average weights state-specific expression by each batch’s state proportions, so unequal weights can produce different averages despite identical conditional expression. (2) No. The state’s biological expression and batch are confounded in observed data because there is no counterpart in another batch; the known simulation parameters can distinguish them, but the observed difference alone cannot. (3) The log mean count is \(\log(2)+1.5-0.2=\log(2)+1.3\), so the expected count is \(2\exp(1.3)\approx7.339\). (4) A single draw may be unusually favorable or unfavorable because of count and sampling variation. Independent replicates show how results vary across draws under the same specified design.