Modeling Gene Expression Along Branches

Why Model Gene Expression Along Branches?

A lineage model describes cells moving along a shared trajectory before a branchpoint and into distinct branches afterward. A gene may follow the same expression pattern before that point but change differently along the emerging lineages. The statistical question is not simply whether expression differs between two groups of cells: it is whether the expression trend varies by branch after accounting for the shared developmental path.

Here, \(t_i\) denotes the pseudotime of cell \(i\), \(\tau\) the branchpoint, \(z_i\in\{1,2\}\) its branch, and \(y_{ig}\) the measured expression for gene \(g\). For clarity, the worked calculations use log-normalized expression and treat the branchpoint and pseudotimes as fixed. In practice, pseudotime is an estimated coordinate, and that estimation can affect downstream uncertainty.

A Continuous Model for Two Branches

A simple model uses straight-line trends after the branchpoint. Define \(x_i=t_i-\tau\) and \(x_{i+}=\max(0,x_i)\). For one gene, consider

The indicator \(I(z_i=2)\) equals 1 for branch 2 and 0 for branch 1. Before the branchpoint, \(x_{i+}=0\), so both branches have the same fitted mean, \(\beta_0+\beta_1x_i\). At \(x_i=0\), the post-branch terms are also zero, so the fitted curves meet continuously at \(\beta_0\). After the branchpoint, the branch 1 slope is \(\beta_1+\beta_2\), and the branch 2 slope is \(\beta_1+\beta_2+\beta_3\). Thus \(\beta_3\) is the difference in post-branch slopes, in expression units per pseudotime unit.

This model assumes that the conditional mean is adequately described by the specified piecewise-linear curves. For the usual exact small-sample \(t\) and \(F\) tests, it also assumes independent errors with mean zero, a common variance, and normally distributed errors conditional on the predictors. These assumptions concern residuals around the fitted means, not the marginal distribution of all observed expression values. Log-normalized expression can sometimes be modeled this way as an approximation; it is not a claim that raw RNA counts are Gaussian.

For curved trends, replace the single post-branch hinge with a set of basis functions. Let \(B_k(s)\) be post-branch basis functions that equal zero at \(s=0\). A model can include shared terms and branch-by-basis interactions:

The shared part includes the common pre-branch curve and any shared post-branch trend. The functions \(B_k\) allow branch 2 to deviate from branch 1 after the branchpoint, while their zero value at the branchpoint preserves continuity. For two branches, testing \(H_0:\gamma_1=\cdots=\gamma_K=0\) asks whether the branch-specific deviations improve fit. A linear trend comparison is the special case \(K=1\).

Testing the Branch-Specific Trend

For the linear model, the null hypothesis of no branch-specific slope difference is \(H_0:\beta_3=0\), against \(H_a:\beta_3\ne0\). A coefficient test uses \(t=\hat{\beta}_3/SE(\hat{\beta}_3)\). A confidence interval for \(\beta_3\) expresses uncertainty in the difference in slopes, with units of expression per pseudotime unit. For several post-branch basis functions, compare the full model with a reduced model that removes all branch-by-basis terms.

The nested-model \(F\) statistic follows from comparing residual sums of squares. If the full model has \(p_{\mathrm{full}}\) parameters and the reduced model removes \(q\) parameters, then, under the linear-model assumptions and the null hypothesis,

The numerator measures the decrease in residual variation per restriction released; the denominator estimates residual variance in the full model. A large value indicates that the branch-specific terms reduce residual variation substantially relative to the noise estimate. The test does not establish that a branch caused the expression change, nor does a small p-value measure the size or biological importance of the change.

For a curved model with three branch-deviation coefficients, a joint test is preferable to testing each coefficient separately: the null is that all three equal zero. Testing them jointly asks whether the branch-specific curve differs anywhere in the modeled basis space, while separate tests can miss a combined pattern and inflate the chance of false-positive findings if interpreted as a collection of unadjusted tests.

Diagnostics, Alternatives, and Limitations

A branch-specific test is informative only if the data contain information about the curves being compared. Inspect expression against pseudotime separately by branch, along with fitted curves and residuals. Check that both branches have observations across enough of the post-branch range to distinguish a trend difference from a few isolated cells. Look for residual patterns, changing residual spread, influential cells, and systematic curvature left by a straight-line fit. For spline models, check whether the fitted shape changes sharply with the number or placement of basis functions.

The branchpoint also matters. If \(\tau\) is estimated rather than fixed, treating its estimate as known can understate uncertainty. A reported coefficient test is then conditional on the chosen branchpoint unless the fitting and inference procedure accounts for its estimation. Repeating the fit over plausible branchpoint values is a useful sensitivity check. Likewise, if branch labels are uncertain, assigning each cell to its most likely branch and then treating that assignment as certain can understate uncertainty; use a procedure that carries the branch uncertainty into the fit.

The alternatives answer different questions. A single pooled curve ignores branch-specific trends and may conceal divergence. Fitting unrelated curves on each branch allows divergence, but does not enforce a shared pre-branch segment or continuity at the branchpoint. The constrained shared-plus-deviation model makes that biological structure explicit and permits a direct nested comparison. A more flexible spline can represent curvature, but uses more parameters and may fit noise when data are sparse. In each case, compare models that use the same response scale and observations, and interpret a test in light of the chosen basis and assumptions.

Finally, these calculations are gene-specific. If many genes are tested, the collection of p-values requires an appropriate multiple-testing procedure before making claims across genes. A significant trend difference is statistical evidence of branch-associated expression patterns under the model, not proof of a causal regulatory mechanism. Poor coverage, misspecified mean curves, nonconstant variance, dependence among cells, and uncertainty in pseudotime can all make nominal standard errors and p-values unreliable. In single-cell data, also inspect whether conclusions are robust to the expression scale and to reasonable modeling choices.

Common Mistake

Wrong reasoning: “The branch 1 test is significant and the branch 2 test is not, so the branches have different trends.” Significance in one subgroup and nonsignificance in another is not itself a test of the difference between subgroups. Correction: fit a model containing the branch-by-pseudotime term and test that term directly, or compare the nested shared-trend and branch-specific models. The interaction estimates the difference of interest and its uncertainty.