Tutorials › Statistical Methods for Single-Cell RNA Sequencing › Canonical Correlation for Paired RNA and Protein

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

Canonical Correlation for Paired RNA and Protein

Learn how CCA chooses RNA and protein directions whose cell-level scores have maximal correlation, and when regularization or caution is needed.

Advanced 12 min read

What You'll Learn

  • Define paired RNA and protein matrices and their canonical variates.
  • Derive the correlation-maximization objective under unit-variance constraints.
  • Connect the canonical correlations to singular values of a whitened cross-covariance matrix.
  • Interpret canonical directions without treating them as causal effects or unique biological programs.
  • Identify singular covariance matrices and overfitting risks in high-dimensional single-cell data.
  • Compare CCA with PCA and partial least squares.

Why Relate RNA and Protein Through Directions?

In CITE-seq, each cell can have an RNA profile and a set of antibody-derived protein measurements. A question of interest is whether a weighted combination of RNA features varies across cells in step with a weighted combination of protein features. Looking at one gene and one protein at a time can miss a relationship that is distributed across several features. Canonical correlation analysis (CCA) addresses this by choosing one linear score from each modality to maximize their correlation across paired cells.

CCA is a method for finding linear associations in the observed paired measurements. It does not establish that RNA changes cause protein changes, and a high correlation does not by itself identify a biological mechanism. The earlier material on scaling modalities is relevant when selecting and preparing features, especially because regularized versions of CCA can depend on feature scale.

Data, Canonical Variates, and the Objective

Suppose \(n\) cells have measurements on \(p\) RNA features and \(q\) protein features. Let \(X\) be the \(n \times p\) matrix of centered RNA measurements and \(Y\) the \(n \times q\) matrix of centered protein measurements. Rows must refer to the same cells in the same order. Centering means subtracting each feature’s sample mean from its measurements. Define the sample covariance matrices

$$ C_{XX}=\frac{X^\top X}{n-1}, \qquad C_{YY}=\frac{Y^\top Y}{n-1}, \qquad C_{XY}=\frac{X^\top Y}{n-1}. $$

For coefficient vectors \(a \in \mathbb{R}^p\) and \(b \in \mathbb{R}^q\), the canonical variates are the cell-level scores \(u=Xa\) and \(v=Yb\). Their sample covariance is \(a^\top C_{XY}b\); their sample variances are \(a^\top C_{XX}a\) and \(b^\top C_{YY}b\). Thus the correlation being optimized is

$$ \operatorname{corr}(Xa,Yb) = \frac{a^\top C_{XY}b} {\sqrt{(a^\top C_{XX}a)(b^\top C_{YY}b)}}. $$
Definition: The first pair of canonical variates consists of linear combinations \(Xa\) and \(Yb\) that maximize their sample correlation over nonzero \(a\) and \(b\). The first canonical correlation is the resulting maximum correlation, conventionally reported as a nonnegative value. Later pairs maximize correlation subject to being uncorrelated, within each modality, with the earlier canonical variates.

Because multiplying either coefficient vector by a nonzero constant does not change the absolute correlation, the maximization is usually written with unit-variance constraints:

$$ \max_{a,b}\ a^\top C_{XY}b \quad \text{subject to} \quad a^\top C_{XX}a=1,\qquad b^\top C_{YY}b=1. $$

The sign of either variate is arbitrary: changing \(a\) to \(-a\) reverses \(Xa\) without changing the magnitude of its correlation with \(Yb\). A reported negative correlation can therefore be made positive by reversing one direction. The conventional canonical correlation is the nonnegative magnitude.

Worked Example: One RNA Feature and One Protein Feature

Consider four hypothetical paired cells. Their centered RNA measurements are \(x=(-1.5,-0.5,0.5,1.5)\), and their centered protein measurements are \(y=(-1,-1,1,1)\). With one feature per modality, there is only one possible direction in each feature space, so CCA reduces to the sample correlation.

The sample covariance is \(s_{xy}=\sum_i x_i y_i/(n-1)\). The cross-products sum to \(1.5+0.5+0.5+1.5=4\), so \(s_{xy}=4/3\). The sample variances are \(s_x^2=5/3\) and \(s_y^2=4/3\). Therefore,

$$ r=\frac{4/3}{\sqrt{(5/3)(4/3)}} =\frac{4}{\sqrt{20}} \approx 0.8944. $$

The first canonical correlation is approximately \(0.8944\). In this hypothetical set of four paired measurements, the centered RNA and protein scores move together strongly. This is a descriptive association among these cells, not evidence that one measurement causes the other.

Why the Solution Is a Singular-Value Problem

The constraints in the objective are expressed using each modality’s covariance. When both covariance matrices are positive definite, they can be used to put the scores on unit-variance scales. Set \(z=C_{XX}^{1/2}a\) and \(w=C_{YY}^{1/2}b\). The constraints become \(z^\top z=w^\top w=1\). Substituting \(a=C_{XX}^{-1/2}z\) and \(b=C_{YY}^{-1/2}w\) into the objective gives

$$ a^\top C_{XY}b = z^\top \left(C_{XX}^{-1/2}C_{XY}C_{YY}^{-1/2}\right)w. $$

Write the matrix in parentheses as \(R\). For any fixed unit vector \(z\), the largest value of \(z^\top Rw\) over unit \(w\) is \(\|R^\top z\|\), by the Cauchy–Schwarz inequality. Maximizing that quantity over unit \(z\) gives the largest singular value of \(R\). If \(R=U\Sigma V^\top\) is a singular value decomposition, its leading left and right singular vectors provide the maximizing directions.

$$ R=C_{XX}^{-1/2}C_{XY}C_{YY}^{-1/2} =U\Sigma V^\top, \qquad a=C_{XX}^{-1/2}u_1,\quad b=C_{YY}^{-1/2}v_1. $$

The first canonical correlation is the largest singular value \(\sigma_1\) of \(R\). Subsequent singular values give further canonical correlations, with corresponding directions that produce scores uncorrelated with earlier scores within each modality. This result explains the procedure: first account for each modality’s covariance, then find the strongest remaining cross-modal association.

The derivation requires positive-definite \(C_{XX}\) and \(C_{YY}\) for the ordinary inverse square roots. It also assumes paired, consistently ordered rows and finite sample covariances. The optimization itself is an algebraic description; treating cells as independent observations for uncertainty assessment would require a separate justification. In studies with multiple cells per donor, donor-level dependence can make cell-level inferential claims overconfident.

Worked Example: Two Directions From a Covariance Summary

Suppose hypothetical paired-cell measurements yield the following covariance matrices, with the rows and columns ordered as two RNA features and two protein features:

$$ C_{XX}= \begin{pmatrix}1&0\\0&1\end{pmatrix},\qquad C_{YY}= \begin{pmatrix}1&0\\0&1\end{pmatrix},\qquad C_{XY}= \begin{pmatrix}0.8&0\\0&0.3\end{pmatrix}. $$

Both within-modality covariance matrices are positive definite, so the ordinary whitening step is defined. Their inverse square roots are identity matrices, making \(R=C_{XY}\). Its singular values are \(0.8\) and \(0.3\). For the first pair, \(a=(1,0)^\top\) and \(b=(1,0)^\top\), giving RNA score \(X_1\), protein score \(Y_1\), and correlation \(0.8\). For the second pair, the directions are the second features and the correlation is \(0.3\).

Thus, among these two hypothetical feature combinations, the first paired directions have the stronger linear association across cells. This example uses covariance matrices as the given numerical summary of the hypothetical paired data; it does not imply that the original feature values are independent or normally distributed.

Interpreting Directions and Diagnosing Failure

A coefficient in \(a\) or \(b\) describes how a measured feature contributes to its canonical variate, conditional on the other features in that linear combination. It is not automatically a feature’s correlation with the variate, nor is it a causal effect. When features are strongly correlated, different coefficient vectors can yield similar scores, so individual coefficients may be unstable even when the resulting canonical correlation is high. Inspecting the score-feature relationships and the stability of directions across resamples or donors can help reveal that problem.

Classical CCA can fail computationally or become misleading when the number of features is large relative to the number of paired cells. If \(p \geq n\), the centered matrix \(X\) has rank at most \(n-1\), so \(C_{XX}\) cannot be full rank; the same issue applies to \(C_{YY}\) when \(q \geq n\). Nearly redundant features can also make a covariance matrix ill-conditioned. Inverting such a matrix can give disproportionate influence to low-variance directions and produce unstable coefficients or very high in-sample canonical correlations.

Worked Example: A Singular RNA Covariance Matrix

Suppose three hypothetical paired cells have two centered RNA feature columns that are identical. Their covariance summary is

$$ C_{XX}= \begin{pmatrix}1&1\\1&1\end{pmatrix}. $$

The eigenvalues are \(2\) and \(0\), since the vectors \((1,1)^\top\) and \((1,-1)^\top\) are eigenvectors with those eigenvalues. The zero eigenvalue means the inverse square root required for classical CCA does not exist. This is not fixed by reporting more decimal places: the two columns contain no independent variation.

One possible regularization replaces \(C_{XX}\) with \(C_{XX}+\lambda I\), where \(\lambda>0\). For \(\lambda=0.1\), the eigenvalues become \(2.1\) and \(0.1\), and the matrix is invertible; its condition number is \(2.1/0.1=21\), rather than infinite. This stabilizes the calculation, but the choice of \(\lambda\) affects the fitted directions. Regularized CCA is therefore a modified objective, not a way to recover information absent from the measurements. Its tuning and validation should be reported, and a regularized singular value should not be casually interpreted as the ordinary empirical correlation.

Useful diagnostics include checking the rank and conditioning of each covariance matrix, confirming that rows truly match the same cells, and examining whether the estimated directions are dominated by a small number of noisy features. Compare in-sample association with performance on held-out cells or, where dependence matters, held-out donors. A very large training correlation that collapses on held-out data is evidence of overfitting, not strong generalizable cross-modal structure. CCA describes linear association; nonlinear relationships may not be captured, and a high canonical correlation does not establish agreement between modalities or a shared cause.

How CCA Differs From Nearby Linear Methods

MethodWhat its directions prioritizeHow that differs from CCA
Principal component analysis (PCA)Variance within one modalityPCA can represent large RNA variation without finding the RNA direction most correlated with protein.
Partial least squares (PLS)Cross-modal covarianceCCA maximizes correlation after scaling each score to unit variance; PLS can favor directions with larger score variances.
Canonical correlation analysis (CCA)Correlation between one linear score from each modalityIt balances cross-modal association against the within-modality variance constraints.

These objectives suit different questions. PCA is useful for summarizing variation within a modality; PLS is a direct choice when covariance is the target; CCA is appropriate when the target is linear association between paired scores. Classical CCA can be sensitive to small or singular covariance matrices, while regularization can improve numerical stability at the cost of adding choices that affect the result. None of these methods, by itself, resolves confounding or establishes a causal relationship.

Key takeaway: CCA selects paired linear combinations of RNA and protein measurements to maximize their across-cell correlation. Its standard singular-value solution relies on invertible within-modality covariance matrices; in high dimensions, regularization and out-of-sample checks are essential, and the resulting association remains descriptive.

Check Your Understanding

Try each question, then compare your reasoning with the brief solutions.

  1. Why must the rows of the RNA and protein matrices refer to the same cells for paired CCA?
  2. If \(C_{XX}\) and \(C_{YY}\) are identity matrices and \(C_{XY}\) has singular values \(0.65\) and \(0.20\), what is the first canonical correlation?
  3. Why can classical CCA be undefined when there are more RNA features than paired cells?
  4. How does PCA’s objective differ from CCA’s objective?
  5. Does a canonical correlation of \(0.9\) show that RNA changes cause protein changes? Why or why not?

Solutions. (1) The cross-covariance pairs measurements cell by cell, so mismatched rows would calculate a different, invalid pairing. (2) The first canonical correlation is the largest singular value, \(0.65\). (3) A centered matrix with \(n\) rows has rank at most \(n-1\); with more features than cells its covariance matrix is singular, so the inverse square root in classical CCA is unavailable. (4) PCA maximizes within-modality variance, whereas CCA maximizes correlation between paired modality scores. (5) No. The value describes strong linear association among the paired scores; it does not establish causation or rule out shared causes.

Clinical Biostats