Decomposing co-expression in twin sample space

TwINFER uses recently divided sister cells, or twins. Right after division, two sisters share the molecules and the cell state they inherited from their mother cell; with time, the expression of each sister drifts away from that of the other. This page derives how, as a result, the co-expression of two genes X and Y, their correlation across cells, splits into a part that twins share and a part that differs between twins. The two parts tell apart two reasons why two genes are correlated:

• regulation: X acts on Y inside every cell, so the correlation also shows up in the differences between twins;

• heterogeneity: cells sit in different states (for example, two cell types) that set both genes; twins sit in the same state, so the correlation is shared by twins and is absent from their differences.

Section 4 shows how TwINFER turns the two parts into the z-scores and the inference scenarios it reports.

1. Definitions

Consider a sample \(S\) with cells expressing the two genes X and Y, with corresponding transcript counts \(x_i(t)\) and \(y_i(t)\) in the \(i\)-th cell at time \(t\). The co-expression of genes X and Y is defined as

\[ \rho := \rho(x,y) = \operatorname*{spearman}_{i\in S}\left[x_i(t),\, y_i(t)\right]. \tag{E1} \]

Now imagine that each cell in \(S\) has divided to yield a pair of twin cells at \(t = 0\) and let us define the ordered twin sample set \(\hat{S}\), where each element \((i,j)\in\hat{S}\) is a twin pair \((i\neq j)\). Note that \((i,j)\) and \((j,i)\) are both in the set. After division the \(N\) twin pairs represent \(2N\) individual cells, for which we compute the ranks \(r_{x_i}(t)\) and \(r_{y_i}(t)\) of the \(i\)-th cell as before. Since twins inherit the same initial condition and dynamics, the co-expression in the offspring population is

\[ \operatorname*{spearman}_{(i,j)\in\hat{S}}\left[x_i(t),\, y_i(t)\right] = \operatorname*{spearman}_{i\in S}\left[x_i(t),\, y_i(t)\right] = \rho. \tag{E2} \]

The fact that we can think of \(\rho\) in the twin sample space without modifications will later allow us to express the conventional co-expression in terms of components in twin sample space.

In the twin sample space \(\hat{S}\), we define the twin shared correlation component

\[ \hat{\rho}_s := \hat{\rho}_s(x,y) = \operatorname*{spearman}_{(i,j)\in\hat{S}}\left[x_i(t),\, y_j(t)\right], \tag{E3} \]

namely, the correlation across twin pairs of the expression of gene X in one cell with the expression of gene Y in its twin. For convenience we introduce the notation

\[ \hat{h}_x := \hat{\rho}_s(x,x), \qquad \hat{h}_y := \hat{\rho}_s(y,y). \tag{E4} \]

Next we consider the difference \(\Delta x = x_i - x_j\). We define the twin difference correlation

\[ \hat{\rho}_\Delta := \hat{\rho}_\Delta(\Delta x,\Delta y) = \operatorname*{pearson}_{(i,j)\in\hat{S}}\left[\Delta r_x(t),\, \Delta r_y(t)\right], \tag{E5} \]

where \(r_{x_i}(t)\) and \(r_{y_i}(t)\) are the ranks across the sample, and \(\Delta r_x = r_{x_i} - r_{x_j}\). A subtle point to note is that this is a correlation between a difference of ranks and thus not a Spearman correlation.

In TwINFER. For a gene pair measured in \(N\) twin pairs:

• \(\rho\) is the conventional co-expression over all \(2N\) cells, the correlation that a method ignoring the twin structure would use;

• \(\hat{\rho}_s\) compares gene X in one sister with gene Y in the other sister. The two sisters are connected only through what they inherited from the mother cell, so \(\hat{\rho}_s\) picks up the co-expression carried by this shared state;

• \(\hat{h}_x\) and \(\hat{h}_y\) measure how alike the two sisters are for each gene. They are close to 1 right after division and decrease as the sisters drift apart;

• \(\hat{\rho}_\Delta\) correlates the differences between sisters. What the sisters inherited cancels in \(\Delta r_x = r_{x_i} - r_{x_j}\), so \(\hat{\rho}_\Delta\) picks up how X and Y co-vary inside each cell after division.

2. Decomposition

We start by computing the correlation in Eq. (E5). The covariance is given by

\[ \begin{aligned} \operatorname{cov}(\Delta x,\Delta y) &= \operatorname{cov}(x_i,y_i) + \operatorname{cov}(x_j,y_j) - \operatorname{cov}(x_i,y_j) - \operatorname{cov}(x_j,y_i)\\ &= 2\left[\operatorname{cov}(x_i,y_i) - \operatorname{cov}(x_i,y_j)\right], \end{aligned} \tag{E6} \]

where for clarity we omitted the time dependence and wrote the equation for the counts (the ranks will be substituted at the end). Similarly,

\[ \begin{aligned} \operatorname{var}(\Delta x) &= \operatorname{var}(x_i) + \operatorname{var}(x_j) - 2\operatorname{cov}(x_i,x_j)\\ &= 2\left[\operatorname{var}(x_i) - \operatorname{cov}(x_i,x_j)\right]. \end{aligned} \tag{E7} \]

Plugging in the definition of the Pearson correlation we obtain

\[ \begin{aligned} \operatorname{pearson}(\Delta x,\Delta y) &= \frac{\operatorname{cov}(\Delta x,\Delta y)}{\sqrt{\operatorname{var}(\Delta x)\operatorname{var}(\Delta y)}}\\[4pt] &= \frac{\operatorname{cov}(x_i,y_i) - \operatorname{cov}(x_i,y_j)}{\sqrt{\operatorname{var}(x_i) - \operatorname{cov}(x_i,x_j)}\,\sqrt{\operatorname{var}(y_i) - \operatorname{cov}(y_i,y_j)}}\\[4pt] &= \frac{\operatorname{pearson}(x_i,y_i) - \operatorname{pearson}(x_i,y_j)}{\sqrt{1 - \operatorname{pearson}(x_i,x_j)}\,\sqrt{1 - \operatorname{pearson}(y_i,y_j)}}, \end{aligned} \tag{E8} \]

where in the last step we divided the numerator and the denominator by \(\sqrt{\operatorname{var}(x_i)\operatorname{var}(y_i)}\) and used \(\operatorname{var}(x_j) = \operatorname{var}(x_i)\) and \(\operatorname{var}(y_j) = \operatorname{var}(y_i)\) for \((i,j)\in\hat{S}\).

Substituting the ranks, and according to the definitions in the previous section, we obtain after rearrangement

\[ \rho = \hat{\rho}_s + \hat{\rho}_\Delta\sqrt{\left(1-\hat{h}_x\right)\left(1-\hat{h}_y\right)}. \tag{E9} \]

Eq. (E9) is an exact finite-sample algebraic identity: every step above uses only sample covariances over \(\hat{S}\), where the symmetry of the set makes \(\operatorname{var}(x_j) = \operatorname{var}(x_i)\) and \(\operatorname{cov}(x_j,y_i) = \operatorname{cov}(x_i,y_j)\) hold exactly. The stochastic interpretation of its two terms in Eqs. (E10)–(E13) is, in contrast, justified at the population level.

Thus, the co-expression can be decomposed into a shared component and an unshared (difference) component. To see this, we write \(x_i(t) = x_i^{s}(t) + x_i^{u}(t)\), where \(x_i^{s}(t)\) is the contribution of the initial condition inherited from the mother cell, identical in both twins, \(x_i^{s}(t) = x_j^{s}(t)\) for \((i,j)\in\hat{S}\), and \(x_i^{u}(t)\) is the private contribution that evolved in each cell during \(t\) independent of its twin. Similarly, we write \(y_i(t) = y_i^{s}(t) + y_i^{u}(t)\). For the ranks, the shared and unshared components are defined directly in rank space using conditional expectations,

\[ r_x^{s} = \mathbb{E}\left[r_x \mid Z\right], \qquad r_x^{u} = r_x - r_x^{s}, \]

and similarly for \(y\), where \(Z\) represents the shared initial state of the twins. We then obtain

\[ \rho = \underbrace{\frac{\operatorname{cov}\left(r^{s}_{x_i}, r^{s}_{y_i}\right)}{\sqrt{\operatorname{var}\left(r_{x_i}\right)\operatorname{var}\left(r_{y_i}\right)}}}_{=\,\hat{\rho}_s} + \underbrace{\frac{\operatorname{cov}\left(r^{u}_{x_i}, r^{u}_{y_i}\right)}{\sqrt{\operatorname{var}\left(r_{x_i}\right)\operatorname{var}\left(r_{y_i}\right)}}}_{=\,\hat{\rho}_\Delta\sqrt{\left(1-\hat{h}_x\right)\left(1-\hat{h}_y\right)}}, \tag{E10} \]

where we used \(\operatorname{cov}(r^{s}_{x_i}, r^{u}_{y_i}) = \operatorname{cov}(r^{s}_{x_i}, r^{u}_{y_j}) = 0\) and \(\operatorname{cov}(r^{u}_{x_i}, r^{u}_{y_j}) = 0\) (the same holds when exchanging \(x\) and \(y\)), see the discussion in the next section. Namely, the shared and unshared components enter the co-expression separately, through \(\hat{\rho}_s\) and \(\hat{\rho}_\Delta\), respectively.

In TwINFER. (E10) is what allows TwINFER to tell regulation from heterogeneity with twins:

• co-expression carried by a state that the sisters share (heterogeneity) enters \(\hat{\rho}_s\) and cancels in the twin differences, so it does not enter \(\hat{\rho}_\Delta\);

• co-expression built up by regulation inside each cell after division enters \(\hat{\rho}_\Delta\); it reaches \(\hat{\rho}_s\) only through what the sisters still share, which fades as they drift apart.

The weight \(\sqrt{(1-\hat{h}_x)(1-\hat{h}_y)}\) sets how much of \(\rho\) the twin differences can carry. Right after division \(\hat{h}_x\) and \(\hat{h}_y\) are close to 1, the weight is close to 0, and nearly all of \(\rho\) is shared; as the sisters drift apart, the weight grows and the unshared component becomes visible.

3. Interpretation via the law of total covariance

In the previous section we asserted that \(\operatorname{cov}(r^{s}_{x_i}, r^{u}_{y_i}) = \operatorname{cov}(r^{s}_{x_i}, r^{u}_{y_j}) = 0\) and \(\operatorname{cov}(r^{u}_{x_i}, r^{u}_{y_j}) = 0\). This follows directly from the definition of the shared and unshared parts: the shared part \(r^{s}_{x_i}(t)\) is determined by the shared initial state \(Z\), that is, the initial condition \(r_{x_i}(0), r_{y_i}(0)\) and any initial hidden variable that the twins may share, and the unshared part is independent between twins given \(Z\). Moreover, by the definition of \(r_x^{s}\), the unshared part has zero conditional mean,

\[ \mathbb{E}\left[r_x^{u} \mid Z\right] = 0. \]

Together with the conditional independence of the twins given \(Z\), this justifies the zero-covariance identities in Eqs. (E11)–(E12). By the law of total expectation, averaging first over the dynamics at fixed initial condition, and then over the initial conditions,

\[ \begin{aligned} \mathbb{E}\left[r^{s}_{x_i}(t)\, r^{u}_{y_j}(t)\right] &= \mathbb{E}\left[r^{s}_{x_i}(t)\, \mathbb{E}\left[r^{u}_{y_j}(t) \mid Z\right]\right]\\ &= \mathbb{E}\left[r^{s}_{x_i}(t)\right] \mathbb{E}\left[r^{u}_{y_j}(t)\right], \end{aligned} \tag{E11} \]

hence

\[ \operatorname{cov}\left(r^{s}_{x_i}, r^{u}_{y_j}\right) = \mathbb{E}\left[r^{s}_{x_i}(t)\, r^{u}_{y_j}(t)\right] - \mathbb{E}\left[r^{s}_{x_i}(t)\right] \mathbb{E}\left[r^{u}_{y_j}(t)\right] = 0. \tag{E12} \]

In a similar manner one can show the other two identities.

Using the law of total covariance,

\[ \begin{aligned} \operatorname{cov}\left(r_{x_i}, r_{y_i}\right) &= \operatorname{cov}\left(\mathbb{E}\left[r_{x_i} \mid Z\right],\, \mathbb{E}\left[r_{y_i} \mid Z\right]\right)\\ &\quad {}+ \mathbb{E}\left[\operatorname{cov}\left(r_{x_i}, r_{y_i} \mid Z\right)\right]. \end{aligned} \tag{E13} \]

Plugging in \(r_{x_i}(t) = r^{s}_{x_i}(t) + r^{u}_{x_i}(t)\) and \(r_{y_i}(t) = r^{s}_{y_i}(t) + r^{u}_{y_i}(t)\) and applying all the above identities, and after division by \(\sqrt{\operatorname{var}(r_{x_i})\operatorname{var}(r_{y_i})}\), we obtain Eq. (E9).

In TwINFER. The two properties used above are what sister cells provide: both sisters inherit the same state \(Z\) from the mother cell, and after division the expression of each sister evolves independently of the other given \(Z\). As with any correlation computed from \(N\) twin pairs, \(\hat{\rho}_s\) and \(\hat{\rho}_\Delta\) carry sampling error, so TwINFER tests them against zero with z-scores (section 4).

4. How TwINFER uses the decomposition

For a gene pair measured in \(N\) twin pairs, TwINFER turns the three correlations into analytic z-scores, reported as z_rho, z_rho_s and z_rho_delta:

\[ z_\rho = \frac{\rho}{\sigma_\rho}, \qquad z_s = \operatorname{sign}(\rho)\,\frac{\hat{\rho}_s}{\sqrt{(1 + \rho^2 + \hat{h}_x\hat{h}_y)/(2N)}}, \qquad z_\Delta = \sqrt{N}\,\hat{\rho}_\Delta, \]

where \(\sigma_\rho\), with \(\sigma_\rho^2 = (1 + \hat{h}_x\hat{h}_y)/(2N)\), is the standard deviation of \(\rho\) when X and Y are not correlated (negative values of \(\hat{h}_x\) and \(\hat{h}_y\) are set to 0 in these standard deviations). The twin pairs, not the \(2N\) cells, are the independent units: the more alike the sisters, the larger \(\sigma_\rho\), from \(1/\sqrt{2N}\) for \(\hat{h}_x\hat{h}_y = 0\) to \(1/\sqrt{N}\) for \(\hat{h}_x\hat{h}_y = 1\).

Dividing (E9) by \(\sigma_\rho\) splits \(z_\rho\) into a shared and an unshared contribution, reported as A_s and A_u:

\[ z_\rho = A_s + A_u, \qquad A_s = \frac{\hat{\rho}_s}{\sigma_\rho}, \qquad A_u = \frac{\hat{\rho}_\Delta\sqrt{\left(1-\hat{h}_x\right)\left(1-\hat{h}_y\right)}}{\sigma_\rho}. \]

In terms of (E10), a significant \(z_s\) means that the co-expression has a shared component, and a significant \(z_\Delta\) that it has an unshared component. The three z-scores are tested in this order, \(z_\rho\), then \(z_s\), then \(z_\Delta\), and assign each gene pair to one of four inference scenarios:

scenario	\(z_\rho\)	\(z_s\)	\(z_\Delta\)	reading
1 nothing	not significant			no correlation
2 regulation	significant	not significant		correlation not shared by twins
3 regulation + heterogeneity	significant	significant	significant	a shared part and a twin-difference part
4 heterogeneity	significant	significant	not significant	correlation only through a shared state

Significant means \(|z| > 1.96\) for \(z_\rho\) and \(z_\Delta\), and \(|z| > 2.326\) for \(z_s\); when many gene pairs are scored at once, as in real data, \(z_\rho\) and \(z_\Delta\) are instead tested with a Benjamini–Hochberg correction.

For lineage-tracing data, where a clone can contain more than two cells, the same quantities are computed from all pairs of cells within each clone, each clone given the same weight. With twins measured at two time points, the correlation between sisters measured at different times also gives the direction of regulation. Both are covered in the tutorials inference_real_data_lineage_tracing.ipynb and inference_simulation_two_time_points.ipynb.