# DeCovarT: A Probabilistic and Multidimensional Framework for Cellular Deconvolution in Heterogeneous Biological Samples

Bastien Chassagnol¹,²,*, Grégory Nuel², Etienne Becht¹

1 Institut De Recherches Internationales Servier (IRIS), FRANCE

2 LPSM (Laboratoire de Probabilités, Statistiques et Modélisation), Sorbonne Université,
4, place Jussieu, 75252 PARIS, FRANCE

* bastien_chassagnol@laposte.net

## Abstract

Although bulk transcriptomic analyses have greatly contributed to a better understanding of complex diseases, their sensibility is hampered by the highly heterogeneous cellular compositions of biological samples. To address this limitation, computational deconvolution methods have been designed to automatically estimate the frequencies of the cellular components that make up tissues, typically using reference samples of physically purified populations. However, they perform badly at differentiating closely related cell populations.

We hypothesised that the integration of the covariance matrices of the reference samples could improve the performance of deconvolution algorithms. We therefore developed a new tool, DeCovarT, that integrates the structure of individual cellular transcriptomic network to reconstruct the bulk profile. Specifically, we inferred the ratios of the mixture components by a standard maximum likelihood estimation (MLE) method, using the Levenberg-Marquardt algorithm to recover the maximum from the parametric convolutional distribution of our model. We then consider a reparametrisation of the log-likelihood to explicitly incorporate the simplex constraint on the ratios. Preliminary numerical simulations suggest that this new algorithm outperforms previously published methods, particularly when individual cellular transcriptomic profiles strongly overlap.
# 1 Introduction

The analysis of the bulk transcriptome provided new insights on the mechanisms underlying disease development. However, such methods ignore the intrinsic cellular heterogeneity of complex biological samples, by averaging measurements over several distinct cell populations. Failure to account for changes of the cell composition is likely to result in a loss of specificity (genes mistakenly identified as differentially expressed, while they only reflect an increase in the cell population naturally producing them) and sensibility (genes expressed by minor cell populations are amenable being masked by highly variable expression from major cell populations).

Accordingly, a range of computational methods have been developed to estimate cellular fractions, but they perform poorly in discriminating cell types displaying high phenotypic proximity. Indeed, most of them assume that purified cell expression profiles are fixed observations, omitting the variability and intrinsically interconnected structure of the transcriptome. For instance, the gold-standard deconvolution algorithm CIBERSORT [New15] applies nu-support vector regression (ν-SVR) to recover the minimal subset of the most informative genes in the purified signature matrix. However, this machine learning approach assumes that the transcriptomic expressions are independent.

In contrast to these approaches, we hypothesised that integrating the pairwise covariance of the genes into the reference transcriptome profiles could enhance the performance of transcriptomic deconvolution methods. The generative probabilistic model of our algorithm, DeCovarT (Deconvolution using the Transcriptomic Covariance), implements this integrated approach.

# 2 Model

First, we introduce the following notations:

* $(\mathbf{y} = (y_{gi}) \in \mathbb{R}_{+}^{G \times N}$ is the global bulk transcriptomic expression, measured in $N$ individuals.

* $\mathbf{X} = (x_{gj}) \in \mathcal{M}_{\mathbb{R}^{G \times J}}$ the signature matrix of the mean expression of $G$ genes in $J$ purified cell populations.

* $\mathbf{p} = (p_{ji}) \in ]0, 1]^{J \times N}$ the unknown relative proportions of cell populations in $N$ samples

As in most traditional deconvolution models, we assume that the total bulk expression can be reconstructed by summing the individual contributions of each cell population weighted by its frequency, as stated explicitly in the following linear matricial relationship (Equation (1)):

$$ \mathbf{y} = \mathbf{X} \times \mathbf{p} \qquad (1) $$

In addition, we consider unit simplex constraint on the cellular ratios, $\mathbf{p}$ (Equation (2)):

$$ \left\{ \begin{array}{l} \sum_{j=1}^{J} p_j = 1 \\ \forall j \in \tilde{J}, \quad p_j \ge 0 \end{array} \right. \qquad (2) $$

## 2.1 Standard linear deconvolution model

However, in real conditions with technical and environmental variability, strict linearity of the deconvolution does not usually hold. Thus, an additional error term is usually considered, and without further assumption on the distribution of this error term, the usual approach to retrieve the best of parameters is by minimising the squared error term between the mixture expressions predicted by the linear model and the actual observed response. This optimisation task is achieved through the ordinary least squares (OLS) approach (Equation (3)),

$$ \hat{\mathbf{p}}_{i}^{\text{OLS}} = \arg \min_{\mathbf{p}_{i}} \|\hat{\mathbf{y}}_{i} - \mathbf{y}_{i}\|^{2} = \arg \min_{\mathbf{p}_{i}} \|\mathbf{X}\mathbf{p}_{i} - \mathbf{y}_{i}\|^{2} = \sum_{g=1}^{G} \left( y_{gi} - \sum_{j=1}^{J} x_{gj} p_{ji} \right) \qquad (3) $$
If we additionally assume that the stochastic error term follows a *homoscedastic* zero-centred Gaussian distribution and that the value of the observed covariates (here, the purified expression profiles) is determined (see the corresponding graphical representation in Figure 1a and the set of equations describing it Equation (4)),

$$y_{gi} = \sum_{j=1}^{J} x_{gj}p_{ji} + \epsilon_{i}, \quad y_{gi} \sim \mathcal{N}\left(\sum_{j=1}^{J} x_{gj}p_{ji}, \sigma_{i}^{2}\right), \quad \epsilon_{i} \sim \mathcal{N}(0, \sigma_{i}^{2}) \quad (4)$$

then, the MLE is equal to the OLS, which, in this framework, is given explicitly by Equation (5):

$$\hat{p}_{i}^{\text{OLS}} = (\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{y}_{i} \quad (5)$$

and is known under the the Gauss-Markov theorem.

## 2.2 Motivation of using a probabilistic convolution framework

In contrast to standard linear regression models, we relax in the DeCovarT modelling framework the *exogeneity* assumption, by considering the set of covariates $\mathbf{X}$ as random variables rather than fixed measures, in a process close to the approach of DSection algorithm and DeMixt algorithms. However, to our knowledge, we are the first to weaken the independence assumption between observations by explicitly considering a multivariate distribution and integrating the intrinsic covariance structure of the transcriptome of each purified cell population.

To do so, we conjecture that the $G$-dimensional vector $\mathbf{x}_j$ characterising the transcriptomic expression of each cell population follows a multivariate Gaussian distribution, given by Equation (6):

$$\mathrm{Det}(2\pi\mathbf{\Sigma}_j)^{-\frac{1}{2}} \exp\left(-\frac{1}{2}(\mathbf{x}_j - \mathbf{\mu}_{.j})\mathbf{\Sigma}_j^{-1}(\mathbf{x}_j - \mathbf{\mu}_{.j})^{\top}\right) \quad (6)$$

and parametrised by:

* $\mu_{.j}$, the mean purified transcriptomic expression of cell population j

* $\Sigma_j$, the covariance matrix of each cell population, constrained to be *positive-definite* (see Appendix A.1). Precisely, we retrieve it from inferring its inverse, known as the precision matrix, through the gLasso [Maz11] algorithm. We define $\Theta_j \equiv \Sigma_j^{-1}$ the corresponding *precision matrix*, whose inputs, after normalisation, store the partial correlation between two genes, conditioned on all the others. Notably, pairwise gene interactions whose corresponding off-diagonal terms in the precision matrix are null are considered statistically spurious, and discarded.

To derive the log-likelihood of our model, first we *plugged-in* the mean and covariance parameters $\zeta_j = (\mu_{.j}, \Sigma_j)$ estimated for each cell population in the previous step. Then, setting $\zeta = (\mu, \Sigma)$, $\mu = (\mu_{.j})_{j \in \tilde{J}} \in \mathcal{M}_{G \times J}$, $\Sigma \in \mathcal{M}_{G \times G}$ the known parameters and $\mathbf{p}$ the unknown cellular ratios, we show that the conditional distribution of the observed bulk mixture, conditioned on the individual purified expression profiles and their ratios in the sample, $\mathbf{y}|(\zeta, \mathbf{p})$, is the convolution of pairwise independent multivariate Gaussian distributions. Using the *affine invariance* property of Gaussian distributions, we can show that this convolution is also a multivariate Gaussian distribution, given by Equation (7).

$$\mathbf{y}|(\zeta, \mathbf{p}) \sim \mathcal{N}_{G}(\mu\mathbf{p}, \Sigma) \text{ with } \mu = (\mu_{.j})_{j \in \tilde{J}}, \quad \mathbf{p} = (p_1, \dots, p_J) \text{ and } \Sigma = \sum_{j=1}^{J} p_j^2 \Sigma_j \quad (7)$$

The DAG associated to this modelling framework is shown in Figure 1b).

In the next section, we provide an explicit formula of the log-likelihood of our probabilistic framework, its gradient and hessian, which in turn can be used to retrieve the MLE of our distribution.
<figure><img src="image_1.png" /><figcaption>**Figure 1.** We use the standard graphical convention of graphical models, as depicted in RevBayes webpage. For identifiability reasons, we conjecture that all variability proceeds from the stochastic nature of the covariates.</figcaption></figure>

## 2.3 Derivation of the log-likelihood

From Equation (7), the conditional log-likelihood is readily computed and given by Equation (8):

$$ \ell_{y|\zeta}(\mathbf{p}) = C + \log \left( \text{Det} \left( \sum_{j=1}^{J} p_j^2 \Sigma_j \right)^{-1} \right) - \frac{1}{2} (\mathbf{y} - \mathbf{p}\mu)^{\top} \left( \sum_{j=1}^{J} p_j^2 \Sigma_j \right)^{-1} (\mathbf{y} - \mathbf{p}\mu) \quad (8) $$

## 2.4 First and second-order derivation of the unconstrained DeCovarT log-likelihood function

The stationary points of a function and notably maxima, are given by the roots (the values at which the function crosses the x-axis) of its gradient, in our context, the vector: $\nabla \ell : \mathbb{R}^J \to \mathbb{R}^J$ evaluated at point $\nabla \ell(\mathbf{p}) : [0, 1]^J \to \mathbb{R}^J$. Since the computation is the same for any cell ratio $p_j$, we give an explicit formula for only one of them (Equation (9)):

$$ \begin{aligned} \frac{\partial \ell_{y|\zeta}(\mathbf{p})}{\partial p_j} &= \frac{\partial \log(\text{Det}(\Theta))}{\partial p_j} - \frac{1}{2} \left[ \frac{\partial (y-\mu p)^\top \Theta (y-\mu p)}{\partial p_j} + (y-\mu p)^\top \frac{\partial \Theta}{\partial p_j} (y-\mu p) + (y-\mu p)^\top \Theta \frac{\partial (y-\mu p)}{\partial p_j} \right] \\ &= -\text{Tr}(\Theta \frac{\partial \Sigma}{\partial p_j}) - \frac{1}{2} \left[ -\mu_j^\top \Theta (y-\mu p) - (y-\mu p)^\top \Theta \frac{\partial \Sigma}{\partial p_j} \Theta (y-\mu p) - (y-\mu p)^\top \Theta \mu_{.j} \right] \\ &= -2p_j \text{Tr}(\Theta \Sigma_j) + (y-\mu p)^\top \Theta \mu_{.j} + p_j (y-\mu p)^\top \Theta \Sigma_j \Theta (y-\mu p) \end{aligned} \quad (9) $$

Since the solution to $\nabla (\ell_{y|\zeta}(\mathbf{p})) = 0$ is not closed, we had to approximate the MLE using iterated numerical optimisation methods. Some of them, such as the Levenberg-Marquardt algorithm, require a second-order approximation of the function, which needs the computation of the Hessian matrix.
Deriving once more Equation (9) yields the Hessian matrix, $\mathbf{H} \in \mathcal{M}_{J \times J}$ is given by:

$$ 
\begin{aligned}
\mathbf{H}_{i,i} &= \frac{\partial^2 \ell}{\partial^2 p_i} = -2 \operatorname{Tr}(\Theta \Sigma_i) + 4p_i^2 \operatorname{Tr}((\Theta \Sigma_i)^2) - 2p_i(y - \mu p)^\top \Theta \Sigma_i \Theta \mu_{.i} - \mu_{.i}^\top \Theta \mu_{.i} - \\
&\qquad 2p_i(y - \mu p)^\top \Theta \Sigma_i \Theta \mu_{.i} - (y - \mu p)^\top \Theta (4p_i^2 \Sigma_i \Theta \Sigma_i - \Sigma_i) \Theta (y - \mu p), \quad i \in \tilde{J} \\
\mathbf{H}_{i,j} &= \frac{\partial^2 \ell}{\partial p_i \partial p_j} = 4p_j p_i \operatorname{Tr}(\Theta \Sigma_j \Theta \Sigma_i) - 2p_i(y - \mu p)^\top \Theta \Sigma_i \Theta \mu_{.j} - \mu_{.i}^\top \Theta \mu_{.j} - \\
&\qquad 2p_j(y - \mu p)^\top \Theta \Sigma_j \Theta \mu_{.i} - 4p_i p_j (y - \mu p)^\top \Theta \Sigma_i \Theta \Sigma_j \Theta (y - \mu p), \quad (i,j) \in \tilde{J}^2, i \neq j
\end{aligned}
\quad (10) $$

in which the coloured sections pair one by one with the corresponding coloured sections of the gradient, given in Equation (9).

Matrix calculus can largely ease the derivation of complex algebraic expressions, thus we remind in Appendices A.1 and A.2 relevant matrix properties and derivations. The numerical consistency of these derivatives was asserted with the **numDeriv** package, using the more stable Richardson's extrapolation ([For81]).

However, the explicit formulas for the gradient and the hessian matrix of the log-likelihood function, given in Equation (9) and Equation (10) respectively, do not take into account the simplex constraint assigned to the ratios. While some optimisation methods use heuristic methods to solve this problem, we consider alternatively a reparametrised version of the problem, detailed comprehensively in Appendix A.3.

# 3 Simulations

## 3.1 Simulation of a convolution of multivariate Gaussian mixtures

To assert numerically the relevance of accounting the correlation between expressed transcripts, we designed a simple toy example with two genes and two cell proportions. Hence, using the simplex constraint (Equation (2)), we only have to estimate one free unconstrained parameter, $\theta_1$, and then uses the mapping function, defined in Appendix A.3 to recover the ratios in their original space.

We simulated the bulk mixture, $\mathbf{y} \in \mathcal{M}_{G \times N}$, for a set of artificial samples $N = 500$, with the following generative model:

* We have tested two levels of cellular ratios, one with equi-balanced proportions ($\mathbf{p} = (p_1, p_2 = 1 - p_1) = (\frac{1}{2}, \frac{1}{2})$) and one with highly unbalanced cell populations: $\mathbf{p} = (0.95, 0.05)$.

* Then, each purified transcriptomic profile is drawn from a multivariate Gaussian distribution. We compared two scenarios, playing on the mean distance of centroids, respectively $\mu_{.1} = (20, 22)$, $\mu_{.2} = (22, 20)$ and $\mu_{.2} = (20, 40)$, $\mu_{.2} = (40, 20)$) and building the covariance matrix, $\Sigma \in \mathcal{M}_{2 \times 2}$ by assuming equal individual variances for each gene (the diagonal terms of the covariance matrix, Diag($\Sigma_1$) = Diag($\Sigma_1$) = $\mathbf{I}_2$) but varying the pairwise correlation between gene 1 and gene 2, Cov[$x_{1,2}$], on the following set of values: $\{-0.8, -0.6, \dots, 0.8\}$ for each of the cell population.

* As stated in Equation (1), we assume that the bulk mixture, $\mathbf{y}_{.i}$ could be directly reconstructed by summing up the individual cellular contributions weighted by their abundance, without additional noise.

## 3.2 Iterated optimisation

The extremum, and by extension the MLE, is a root of the gradient of the log-likelihood. However, in our generative framework, the inverse function cancelling the gradient of Equation Equation (8) is non-closed. Instead, iterated numerical optimisation algorithms that consider first or second-order approximations of the function to optimise are used to approximate the roots.
<figure><img src="image_2.png" /><figcaption>**Figure 2.** We used the package **ComplexHeatmap** to display the mean square error (MSE) of the estimated cell ratios, comparing the NNLS output, as implemented in the DeconRNASeq algorithm ([Gon13]), in Panel A, with our newly implemented DeCovarT algorithm, in Panel B. The lower the MSE, the least noisy and biased the estimates. In addition, we added the two-dimensional density plot for the intermediate scenario, for which each population is parameterised by a diagonal covariance matrix, and the most extreme scenarios (those with the highest correlation between genes). The ellipsoids represent for each cell population the 95% confidence region and the red spherical icon and the green triangular icon represent respectively the centroids (average expression of gene 1 and gene 2) of cell population 1 and cell population 2.</figcaption></figure>

The *Levenberg-Marquardt* (*LM*) algorithm ([Lev44]) bridges the gap between between the steepest descent method (first-order) and the Newton-Raphson method (second-order) by inflating the diagonal terms of the Hessian matrix. Far from the endpoint, a second-order descent is favoured for its faster convergence pace, while the steepest approach is privileged close to the extremum since it allows careful refinement of the step size. Specially, we used the LM implementation of R package *marqLevAlg* to infer estimates of the cellular ratios from the bootstrap simulations ([Phi+21]). It notably includes an additional convergence criteria, the relative distance to the maximum (RDM), that sets apart extrema from spurious saddle points.

## 3.3 Results

We compared the performance of DeCovarT algorithm with the DeconRNASeq algorithm ([Gon13]).

Even with a limited toy example including two cell populations characterised only by two genes, we observe that the overlap was a good proxy of the quality of the estimation: the less the overlap between the two cell distributions, the better the quality of the estimation Figure 2.

The package used to generate the simulations and infer ratios from virtual or real biological mixtures with the DeCovarT algorithm is implemented on my personal Github account DeCovarT.

# 4 Perspectives

The new deconvolution algorithm that we implemented, DeCovarT, is the first one based on a multivariate generative model while enforcing explicitly the simplex constraint. Hence, it provides a strong basis to further derive statistical tests to assert whether the proportion of a given cell population differs significantly between two distinct biological conditions.
**Extend the Simulation Framework** To evaluate the biological and statistical interest of DeCovarT, we need to expand the simulation framework, by encompassing a larger number of cell types, genes, and testing the sensitivity of the model by voluntarily including noise in the benchmark evaluation.

The next phase of our evaluation involves real-world experiments, encompassing both blood and solid tumoral samples. To that end, we could start from the Kassandra benchmark, by [Zai22]. This large-scale project evaluates the performance of five established gold-standard and signature-based deconvolution algorithms, including EPIC [Rac17], CIBERSORT [New15], CIBERSORTx [New+19], quanTIsq [Fin19], and ABIS [Mon+19]. The evaluation involves deconvolving six publicly available datasets annotated with both flow cytometry and bulk RNA-seq expression.

**Enhanced Inference and Integration of Co-Expression Networks** All the popular differential gene expression analyses, such as *limma* + *voom* ([Rit+15]), *EdgeR* ([RMS10]) and *Deseq2* ([Var+16]), tend to overlook gene-gene interactions, comparing independently the expression between two conditions for each gene. The usual univariate approach of DGEAs additionally underlies the need of adjusting *p*-values, as numerous genes are examined simultaneously, and without accounting for interactions between them, the probability of observing false positives increases.

To account for correlations among observations, two consecutive papers, by [CL23] and [CCB22], present an innovative Bayesian framework which models proteomic expression across diverse biological conditions as multivariate Gaussian distributions. Insightful discussions with the main author, Marie Chion, suggest a straightforward extension of the method to transcriptomic expression, given the close relationship between two kinds of omics, both depicting counts.

While the methodology was originally designed to delineate differentially expressed genes between two conditions, the method can be readily extended to incorporate a *one-vs-all strategy*. This extension allows for the identification of markers specific to a particular cell population in comparison to all others. Furthermore, the generative model aligns closely with our deconvolution framework, leveraging the same distributions to describe cellular omic profiles. Alternatively, differential network approaches, such as INDEED, by [Zuo16], implement heuristic and dual-optimisation approaches, finding the sweet spot between maximising the mean differences between purified expression profiles and differentiating the neighbourhood network structure.

The gLasso algorithm used to derive the precision matrix associated to each purified cell profile is subjected to *parameter shrinkage*, like any penalty regularisation approach. Notably, in our setting, shrinkage tends to systematically underestimate the non-zero partial correlations of the precision matrix.

To mitigate this issue, one approach is to incorporate the *support* (indicating non-null inputs), derived from the gLasso output, into a conventional Maximum Likelihood Estimation (MLE) framework. The general concept is to utilising the true "zeros" to impose topological constraints on the final Gaussian Graphical Model (GGM). However, it's important to note that unless the undirected Markov network obtained from the gLasso output is a *chordal graph*, there is usually no straightforward mapping between the two topological spaces.

Finally, the inclusion of prior biological knowledge, such as the strength of relationships between transcription factors, retrieved from Protein-Protein Interaction (PPI) networks, can help reduce the exponential space of undirected graphs to explore.

**Enhanced Inference and Integration of Co-Expression Networks** All the methods outlined in Section 4 yield a subset of genes that distinguish a particular cell population from all others. However, when we combine these gene subsets, we often end up with a non-scalable signature matrix, presenting strong multicollinearity resulting from the redundancy between the gene markers identified.

To further refine the final set of genes able to delineating any cell population included in the signature matrix, *AutoGeneS*, by [Ali21], introduces a greedy genetic approach coupled with a dual optimisation approach <fnref n="1" />. Precisely, the loss function involves minimising inter-population correlation

<footnote n="1">¹In standard approaches that rely on linear regression, the condition number serves as the gold-standard metric for assessing the level of precision of the linear model achievable with the design matrix</footnote>
while simultaneously maximising the distance of the centroids.

We propose instead of this dual optimisation approach the minimisation of the *global overlap* between the concatenated distributions of the cellular profiles. Indeed, this metric not only captures in a single criterion the combined influence of mean inter-cluster distance and differential network structure in delineating purified cellular expression profiles, but supplies a straightforward score easy to interpret. The *overlap* metric precisely measures the shared probability mass and the degree of concurrence in probability densities. In simpler terms, it quantifies the global probability of incorrectly assigning an expression profile to the wrong cell subtype when utilising a maximum a posteriori approach, with the knowledge of each cellular profile's individual parameters.

**Joint Estimation of purified Expression Profiles and Cellular Ratios** The generative model underlying the DeCovarT framework (Figure 1b) assumes that both the ratios and the purified cellular expression profiles are unobserved and need to be inferred from our model. However, we derived explicit formulas for the Gradient (eq. (9)) and Hessian (eq. (10)) of the associated log-likelihood function as if the purified expression profiles had been observed, by heuristically replacing the unknown and sample-specific purified expression profiles $X_{.,i}$ with their averaged counterparts $\mu$. However, jointly optimising the cellular ratios and the purified expression profiles results in a non identifiable problem exhibiting an infinite number of solutions, without strong prior assumptions or regularisation of the unknown parameters to estimate. Finally, it's a highly intractable analytical task, and it is quite likely that no explicit form of the Gradient, nor the Hessian, could be derived.

We detail in Appendix B a Gibbs sampler to approximate the target distribution, here the joint value of the purified profiles and the cellular ratios. In addition, MCMC sampling allows for straightforward incorporation of prior knowledge, and streamlines the derivation of Maximum a Posteriori (MAP) estimates and *credible intervals*.

Precisely, by coupling Gibbs and Metropolis Hasting samplers, we ensure at each iteration that the estimated parameters adhered to the “balance condition”, an essential property guaranteeing the convergence of MCMC chains to a stationary distribution identifiable to the target distribution.

## References

[Lev44] Levenberg, Kenneth. "A method for the solution of certain non-linear problems in least squares". In: *Quarterly of Applied Mathematics* (1944). ISSN: 0033-569X, 1552-4485. DOI: 10.1090/qam/10666. URL: https://www.ams.org/qam/1944-02-02/S0033-569X-1944-10666-0/.

[For81] Bengt Fornberg. "Numerical Differentiation of Analytic Functions". In: *ACM Trans. Math. Softw*. (1981). DOI: 10.1145/355972.355979. URL: https://doi.org/10.1145/355972.355979.

[RMS10] Mark D. Robinson, Davis J. McCarthy, and Gordon K. Smyth. "edgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data". In: *Bioinformatics* (Jan. 1, 2010). ISSN: 1367-4803. DOI: 10.1093/bioinformatics/btp616. URL: https://doi.org/10.1093/bioinformatics/btp616.

[Maz11] Mazumder, Rahul and Hastie, Trevor. "The Graphical Lasso: New Insights and Alternatives". In: *Electronic Journal of Statistics* (2011). DOI: 10.1214/12-EJS740.

[Gon13] Gong, Ting and Szustakowski, Joseph D. "DeconRNASeq: A Statistical Framework for Deconvolution of Heterogeneous Tissue Samples Based on mRNA-Seq Data". In: *Bioinformatics* (Oxford, England) (2013). DOI: 10.1093/bioinformatics/btt090.

[New15] Newman, Aaron and Liu, Chih and others. "Robust Enumeration of Cell Subsets from Tissue Expression Profiles". In: *Nature methods* (2015). DOI: 10.1038/nmeth.3337.
[Rit+15] Matthew E. Ritchie et al. "Limma Powers Differential Expression Analyses for RNA-sequencing and Microarray Studies". In: *Nucleic Acids Research* (Apr. 20, 2015). ISSN: 0305-1048. DOI: 10.1093/nar/gkv007. URL: https://doi.org/10.1093/nar/gkv007.

[Var+16] Hugo Varet et al. "SARTools: A DESeq2- and EdgeR-Based R Pipeline for Comprehensive Differential Analysis of RNA-Seq Data". In: *PLOS ONE* (June 9, 2016). ISSN: 1932-6203. DOI: 10.1371/journal.pone.0157022. URL: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0157022.

[Zuo16] Zuo, Yiming and Cui, Yi and others. "INDEED: Integrated Differential Expression and Differential Network Analysis of Omic Data for Biomarker Discovery". In: *Methods (San Diego, Calif.)* (2016). DOI: 10.1016/j.ymeth.2016.08.015.

[Rac17] Racle, Julien and de Jonge, Kaat and others. "Simultaneous Enumeration of Cancer and Immune Cell Types from Bulk Tumor Gene Expression Data". In: *eLife* (2017). Ed. by Alfonso Valencia. DOI: 10.7554/eLife.26476.

[Fin19] Finotello, Francesca and Mayer, Clemens and others. "Molecular and Pharmacological Modulators of the Tumor Immune Contexture Revealed by Deconvolution of RNA-seq Data". In: *Genome Medicine* (2019). DOI: 10.1186/s13073-019-0638-6.

[Mon+19] Gianni Monaco et al. "RNA-Seq Signatures Normalized by mRNA Abundance Allow Absolute Deconvolution of Human Immune Cell Types". In: *Cell Reports* (2019). DOI: 10.1016/j.celrep.2019.01.041.

[New+19] Aaron M. Newman et al. "Determining cell type abundance and expression from bulk tissues with digital cytometry". In: *Nature Biotechnology* (2019). DOI: 10.1038/s41587-019-0114-2.

[Ali21] Aliee, Hananeh and Theis, Fabian J. "AutoGeneS: Automatic Gene Selection Using Multi-Objective Optimization for RNA-seq Deconvolution". In: *Cell Systems* (2021). DOI: 10.1016/j.cels.2021.05.006.

[Phi+21] Viviane Philipps et al. "Robust and Efficient Optimization Using a Marquardt-Levenberg Algorithm with R Package marqLevAlg". In: *The R Journal* (2021). ISSN: 2073-4859. DOI: 10.32614/RJ-2021-089. URL: http://arxiv.org/abs/2009.03840.

[CCB22] Marie Chion, Christine Carapito, and Frédéric Bertrand. "Accounting for Multiple Imputation-Induced Variability for Differential Analysis in Mass Spectrometry-Based Label-Free Quantitative Proteomics". In: *PLOS Computational Biology* (Aug. 29, 2022). ISSN: 1553-7358. DOI: 10.1371/journal.pcbi.1010420. URL: http://arxiv.org/abs/2108.07086.

[Zai22] Zaitsev, Aleksandr and Chelushkin, Maksim and others. "Precise Reconstruction of the TME Using Bulk RNA-seq and a Machine Learning Algorithm Trained on Artificial Transcriptomes". In: *Cancer Cell* (2022). DOI: 10.1016/j.ccell.2022.07.006.

[CL23] Marie Chion and Arthur Leroy. "A Bayesian Framework for Multivariate Differential Analysis Accounting for Missing Data". In: (July 18, 2023). DOI: 10.48550/arXiv.2307.08975. arXiv: 2307.08975 [stat]. URL: http://arxiv.org/abs/2307.08975 (visited on 02/01/2024). preprint.
# Appendix of Article 4

We recall for readability motivations the log-likelihood of DeCovarT's non-constrained generative model, conditioned on the purified and global bulk expression profiles, along with its gradient and its Hessian.

The conditional log-likelihood is readily computed and given by Equation (D.1):

$$ \ell_{y|\zeta}(\boldsymbol{p}) = C + \log \left( \text{Det} \left( \sum_{j=1}^{J} p_j^2 \boldsymbol{\Sigma}_j \right)^{-1} \right) - \frac{1}{2} (\boldsymbol{y} - \boldsymbol{p}\boldsymbol{\mu})^{\top} \left( \sum_{j=1}^{J} p_j^2 \boldsymbol{\Sigma}_j \right)^{-1} (\boldsymbol{y} - \boldsymbol{p}\boldsymbol{\mu}) \quad (\text{D.1}) $$

The Jacobian is given by Equation (D.2):

$$ \frac{\partial \ell_{y|\zeta}(\boldsymbol{p})}{\partial p_j} = -2p_j \operatorname{Tr}(\boldsymbol{\Theta}\boldsymbol{\Sigma}_j) + (\boldsymbol{y} - \boldsymbol{\mu}\boldsymbol{p})^{\top}\boldsymbol{\Theta}\boldsymbol{\mu}_{.j} + p_j(\boldsymbol{y} - \boldsymbol{\mu}\boldsymbol{p})^{\top}\boldsymbol{\Theta}\boldsymbol{\Sigma}_j\boldsymbol{\Theta}(\boldsymbol{y} - \boldsymbol{\mu}\boldsymbol{p}) \quad (\text{D.2}) $$

The Hessian, $\mathbf{H} \in \mathcal{M}_{J \times J}$, is given by Equation (D.3):

$$ \mathbf{H}_{i,i} = \frac{\partial^2 \ell}{\partial p_i^2} = -2 \operatorname{Tr}(\boldsymbol{\Theta}\boldsymbol{\Sigma}_i) + 4p_i^2 \operatorname{Tr}((\boldsymbol{\Theta}\boldsymbol{\Sigma}_i)^2) - 2p_i(\boldsymbol{y} - \boldsymbol{\mu}\boldsymbol{p})^{\top}\boldsymbol{\Theta}\boldsymbol{\Sigma}_i\boldsymbol{\Theta}\boldsymbol{\mu}_{.i} - \boldsymbol{\mu}_{.i}^{\top}\boldsymbol{\Theta}\boldsymbol{\mu}_{.i} - 2p_i(\boldsymbol{y} - \boldsymbol{\mu}\boldsymbol{p})^{\top}\boldsymbol{\Theta}\boldsymbol{\Sigma}_i\boldsymbol{\Theta}\boldsymbol{\mu}_{.i} - (\boldsymbol{y} - \boldsymbol{\mu}\boldsymbol{p})^{\top}\boldsymbol{\Theta}(4p_i^2\boldsymbol{\Sigma}_i\boldsymbol{\Theta}\boldsymbol{\Sigma}_i - \boldsymbol{\Sigma}_i)\boldsymbol{\Theta}(\boldsymbol{y} - \boldsymbol{\mu}\boldsymbol{p}), \quad i \in \tilde{J} \quad (\text{D.3}) $$

$$ \mathbf{H}_{i,j} = \frac{\partial^2 \ell}{\partial p_i \partial p_j} = 4p_j p_i \operatorname{Tr}(\boldsymbol{\Theta}\boldsymbol{\Sigma}_j \boldsymbol{\Theta}\boldsymbol{\Sigma}_i) - 2p_i(\boldsymbol{y} - \boldsymbol{\mu}\boldsymbol{p})^{\top}\boldsymbol{\Theta}\boldsymbol{\Sigma}_i\boldsymbol{\Theta}\boldsymbol{\mu}_{.j} - \boldsymbol{\mu}_{.i}^{\top}\boldsymbol{\Theta}\boldsymbol{\mu}_{.j} - 2p_j(\boldsymbol{y} - \boldsymbol{\mu}\boldsymbol{p})^{\top}\boldsymbol{\Theta}\boldsymbol{\Sigma}_j\boldsymbol{\Theta}\boldsymbol{\mu}_{.i} - 4p_i p_j(\boldsymbol{y} - \boldsymbol{\mu}\boldsymbol{p})^{\top}\boldsymbol{\Theta}\boldsymbol{\Sigma}_i\boldsymbol{\Theta}\boldsymbol{\Sigma}_j\boldsymbol{\Theta}(\boldsymbol{y} - \boldsymbol{\mu}\boldsymbol{p}), \quad (i,j) \in \tilde{J}^2, i \neq j $$
# D.1 Optimisation and calculus

## D.1.1 Multivariate distributions and basic algebra properties

### Definition D.1.1: Multivariate Gaussian distributions

If random vector $\mathbf{X}$ of size $G$ follows a random multivariate Gaussian distribution, $\mathbf{X} \sim N_G(\mu, \mathbf{\Sigma})$, then its distribution is given by:

$$ \mathrm{Det}(2\pi\Sigma)^{-\frac{1}{2}} \exp\left(-\frac{1}{2}(\mathbf{x} - \mu)^{\mathrm{T}}\Sigma^{-1}(\mathbf{x} - \mu)^{\mathrm{T}}\right) $$

in which:

* $\mu = \mathbf{X}$ is the $G$-dimensional mean vector

* $\mathbf{\Sigma}$ is a $G \times G$ positive-definite Definition D.1.2 covariance matrix, whose diagonal terms, $\mathrm{Diag}(\mathbf{\Sigma}) = [(\mathrm{Var}[X_{i,j}]), \forall (i, j) \in \tilde{G}^2, i = j]^\mathrm{T}$ are the individual variances of each purified gene transcript in population $j$ and off-diagonal terms, $\mathbf{\Sigma}_{i,j} = \mathrm{Cov}[X_i, X_j], \forall (i, j) \in \tilde{G}^2, i \neq j$ are the covariance between variables. We note $\Theta = \mathbf{\Sigma}^{-1}$, the inverse of the covariance matrix, called the precision matrix.

### Property D.1.1: Affine invariance property of multivariate GMMs

The two following properties hold for a multivariate Gaussian distribution:

* if $\mathbf{X} \sim N_G(\mu, \mathbf{\Sigma})$, then $p\mathbf{X}$, with $p$ a constant, follows itself a multivariate Gaussian distribution, given by: $p\mathbf{X} \sim N_G(p\mu, p^2\mathbf{\Sigma})$

* given two independent random vectors $\mathbf{X}_1 \sim N_G(\mu_1, \mathbf{\Sigma}_1)$ and $\mathbf{X}_2 \sim N_G(\mu_2, \mathbf{\Sigma}_2)$ following a multivariate Gaussian distribution, then the random variable $\mathbf{X}_1 + \mathbf{X}_2$ follows itself the multivariate Gaussian distribution:

$$ X + Y \sim N_G(\mu_1 + \mu_2, \mathbf{\Sigma}_1 + \mathbf{\Sigma}_2) $$

By induction, this property generalises to the sum of $J$ independent random vectors of same dimension $\mathbb{R}^G$.

Deriving the characteristic function of the multivariate GMM yields directly results reported in Property D.1.1.
Definition D.1.2: Definite matrix

A symmetric real matrix $\mathbf{A}$ of rank $G$ is positive-definite if:

$$ \mathbf{x}^\top \mathbf{A} \mathbf{x} > 0, \quad \mathbf{x} \in \mathbb{R}^G \qquad (D.4) $$

To gain a clearer grasp of the positive-definite constraint imposed on the covariance parameter of a multivariate Gaussian distribution, let's delve into the most straightforward scenario, in which we assume that any of the individual features exhibit pairwise independence. This particular setup is parametrised by a covariance matrix containing exclusively diagonal elements.

If the matrix is not strictly positive-definite, then some of the diagonal elements can display negative values, otherwise that the individual variances for some of the covariates are negative. It is not physically possible and leads to improper, degenerate probability distributions.
## D.1.2 Matrix calculus

Fundamental algebra calculus formulas used to derive first-order and second-order derivatives of the generative model of DeCovarT are reported in Property D.1.2 and Property D.1.3, respectively.

### Property D.1.2: First-order matrix calculus

Given two invertible matrices, $A = A(p)$ and $B = B(p)$, functions of a scalar variable $p$, the following matrix calculus hold:

$$ (a) \quad \frac{\partial \mathrm{Det}(A)}{\partial p} = \quad (b) \quad \frac{\partial UAV}{\partial p} = U \frac{\partial A}{\partial p}V \quad (c) \quad \frac{\partial A^{-1}}{\partial p} = -A^{-1} \frac{\partial A}{\partial p}A^{-1} \mathrm{Det}(A) \mathrm{Tr}(A^{-1} \frac{\partial A}{\partial p}) $$

From a) and fundamental linear algebra properties, we can readily compute applying the chain rule property on the logarithm:

$$ \frac{\partial \log(\mathrm{Det}(A))}{\partial p} = \mathrm{Tr}(A^{-1} \frac{\partial A}{\partial p}) \\ \frac{\partial \log(\mathrm{Det}(A^{-1}))}{\partial p} = -\mathrm{Tr}(A^{-1} \frac{\partial A}{\partial p}) $$

Finally, injecting these first-order matrix derivatives, we obtain:

$$ \frac{\partial(y - xp)^{\top}\Theta(y - xp)}{\partial p} = -2(y - xp)^{\top}\Theta x \\ = -2x^{\top}\Theta(y - xp) \\ \text{with } A = D = -x \in \mathbb{R}^{G}, \quad b = e = y, \quad C = \Theta \text{ symmetric} $$

### Property D.1.3: Second-order matrix calculus

Given an invertible matrix $A$ depending on a variable $p$, the following calculus formulas hold:

$$ (a) \quad \frac{\partial^2 A^{-1}}{\partial p_i \partial p_j} = A^{-1} \left( \frac{\partial A}{\partial p_i} A^{-1} \frac{\partial A}{\partial p_j} - \frac{\partial^2 A}{\partial p_i \partial p_j} + \frac{\partial A}{\partial p_j} A^{-1} \frac{\partial A}{\partial p_i} \right) A^{-1} \left( \begin{array}{c} (b) \\ \mathrm{Tr}\left(\frac{\partial \mathrm{Tr}(A)}{\partial p_i}\right) \end{array} \right) = $$

Combining Property D.1.2 with the linear property of the trace operator yields:

$$ \frac{\partial^2 \log(\mathrm{Det}(A^{-1}))}{\partial^2 p} = -\mathrm{Tr}\left[A^{-1} \frac{\partial^2 A}{\partial^2 p_i}\right] + \mathrm{Tr}\left[\left(A^{-1} \frac{\partial A}{\partial p_i}\right)^2\right] $$
### D.1.3 First and second-order derivation of constrained DeCovarT

To reparametrise the log-likelihood function (Equation (D.1)) in order to explicitly handling the unit simplex constraint (Equation (C.2)), we consider the following mapping function: $\psi: \theta \to p | \theta \in \mathbb{R}^{J-1}, p \in [0, 1]^J$ (Equation (D.5)):

1.

$$ p = \psi(\theta) = \begin{cases} p_j = \frac{e^{\theta_j}}{\sum_{k<j} e^{\theta_k} + 1}, & j < J \\ p_J = \frac{1}{\sum_{k<j} e^{\theta_j} + 1} \end{cases} \quad (D.5) $$

$$ 2. \theta = \psi^{-1}(p) = \left( \ln \left( \frac{p_j}{p_J} \right) \right)_{j \in \{1, \dots, J-1\}} $$

that is a $C^2$-diffeomorphism, since $\psi$ is a bijection between $p$ and $\theta$ twice differentiable. Its Jacobian, $J_{\psi} \in M_{J \times (J-1)}$ is given by Equation (D.6):

$$ J_{i,j} = \frac{\partial p_i}{\partial \theta_j} = \begin{cases} \frac{e^{\theta_i} B_i}{A^2}, & i = j, i < J \\ \frac{-e^{\theta_j} e^{\theta_i}}{A^2}, & i \neq j, i < J \\ \frac{-e^{\theta_j}}{A^2}, & i = J \end{cases} \quad (D.6) $$

with $i$ indexing vector-valued $p$ and $j$ indexing the first-order order partial derivatives of the mapping function, $A = \sum_{j'<J} e^{\theta_{j'}} + 1$ the sum over exponential (denominator of the mapping function) and $B = A - e^{\theta_i}$ the sum over ratios minus the exponential indexed with the currently considered index $i$.

The Hessian of the multi-dimensional mapping function $\psi(\theta)$ exhibits symmetry for each cell ratio component $j$, as anticipated in accordance with Schwarz's theorem. It is a third-order tensor of rank $(J-1)(J-1)J$, given by Equation (D.7):

$$ \frac{\partial^2 p_i}{\partial k \partial j} = \begin{cases} \frac{e^{\theta_i} e^{\theta_l} (-B_i + e^{\theta_l})}{A^3}, & (i < J) \land ((i \neq j) \oplus (i \neq k)) \quad (a) \\ \frac{2e^{\theta_i} e^{\theta_j} e^{\theta_k}}{A^3}, & (i < J) \land (i \neq j \neq k) \quad (b) \\ \frac{e^{\theta_i} e^{\theta_j} (-A + 2e^{\theta_j})}{A^3}, & (i < J) \land (j = k \neq i) \quad (c) \\ \frac{B_i e^{\theta_i} (B_i - e^{\theta_i})}{A^3}, & (i < J) \land (j = k = i) \quad (d) \\ \frac{e^{\theta_j} (-A + 2e^{\theta_j})}{A^3}, & (i = J) \land (j = k) \quad (e) \\ \frac{2e^{\theta_j} e^{\theta_k}}{A^3}, & (i = J) \land (j \neq k) \quad (f) \end{cases} \quad (D.7) $$

with $i$ indexing $p$, $j$ and $k$ respectively indexing the first-order and second-order partial derivatives of the mapping function with respect to $\theta$. In line (a), $\oplus$ refers to the Boolean XOR operator, $\land$ to the AND operator and $l = \{j, k\} \setminus i$.

To derive the log-likelihood function in Equation (D.2), we reparametrise $p$ to $\theta$, using a standard chain rule formula. Considering the original log-likelihood function, Equation (D.1), and the mapping function, Equation (D.5), the differential at the first order and at the second order is given by Equation (D.8) and Equation (D.9), respectively defined in $\mathbb{R}^{J-1}$ and $M_{(J-1) \times (J-1)}$:

$$ \left[ \frac{\partial \ell_{y|\zeta}}{\partial \theta_j} \right]_{j < J} = \sum_{i=1}^{J} \frac{\partial \ell_{y|\zeta}}{\partial p_i} \frac{\partial p_i}{\partial \theta_j} \quad (D.8) $$

$$ \left[ \frac{\partial \ell_{y|\zeta}^2}{\partial \theta_k \partial \theta_j} \right]_{j < J, k < J} = \sum_{i=1}^{J} \sum_{l=1}^{J} \left( \frac{\partial p_i}{\partial \theta_j} \frac{\partial^2 \ell_{y|\zeta}}{\partial p_i \partial p_l} \frac{\partial p_l}{\partial \theta_k} \right) + \sum_{i=1}^{J} \left( \frac{\partial \ell_{y|\zeta}}{\partial p_i} \frac{\partial^2 p_i}{\partial \theta_k \partial \theta_j} \right) \quad (d) \quad (D.9) $$
## D.2 A MCMC Algorithm for the Joint Distribution of Purified Profiles and Ratios

We introduce two variations of the MCMC algorithm, namely the *Metropolis-Hasting* and the *Gibbs sampling* algorithms. They are respectively tailored to approximate distributions for which no explicit form is known (Definition D.2.1) or streamline the optimisation of strongly dependent parameters (Definition D.2.2).
## D.2.1 An Introduction to Gibbs and Metropolis Hasting Samplers

### Definition D.2.1: Metropolis-Hasting algorithm

First of all, we introduce some key notations:

* Function $f(\theta|.)$ is the target distribution, usually the posterior distribution that encompasses both prior knowledge and new data.

* The distribution $q(\theta|\theta^{(q-1)})$ is the proposal distribution, alternatively known as the transition kernel, and the transient value sampled from it, $\theta^{(*)}$, is the proposal estimate.

* The probability of accepting the proposal is naturally called the acceptance probability.

Each iteration, indexed by $(q)$, of the MH algorithm includes the following steps:

1. Draw a proposal, $\theta^{(*)}$, from conditional distribution $q(\theta|\theta^{(q-1)})$.

2. Compute the acceptance probability parameter, noted $\alpha$:

$$K(\theta^{(*)}) = \min \left( \frac{f(\theta^{(*)}|.)}{f(\theta^{(q-1)}|.)} \frac{q(\theta|\theta^{(q-1)})}{q(\theta^{(q-1)}|\theta)}, 1 \right)$$

Typically, choices of the acceptance probability and the kernel distribution are tailored to satisfy the balance condition of the MCMC chains and ensure that the chain behaviour reproduces the sampling pattern of the desired distribution.

3. The decision of whether to accept or reject the new state proposal is determined by the previously computed acceptance probability parameter, denoted as $\alpha$. This parameter sets the threshold for accepting a value drawn from a standard Uniform distribution, $u \sim U[0, 1]$. This mathematical protocol is further described in Equation (D.10):

$$\theta^{(q)} = \begin{cases} \theta^{(*)}, & u \le K(\theta^{(*)}) \\ \theta^{(q-1)}, & u > K(\theta^{(*)}) \end{cases} \qquad (D.10)$$

In the first case, we say that the proposal is accepted, while in the second case, it is rejected.

One of the major advantages of the Metropolis-Hastings (MH) algorithm lies in the design of its acceptance probability function, denoted as $K(\theta^{(*)})$. Indeed, by involving the computation of a ratio between two density functions, the normalisation constant, which represents the value of the marginal likelihood and is usually intractable to compute is naturally cancelled out.

In our modelling framework, the acceptance function further simplifies with the choice of proposing a new proposal by adding an error term following a null-centred, multivariate and symmetric distribution (Property D.2.1):
### Property D.2.1: Random walk Metropolis-Hastings

Indeed, from the analytical properties of the kernel distribution (Equation (D.11) and Equation (D.12)):

$$q(\theta^{(*)}|\theta^{(q-1)}) = q(\epsilon) \qquad (D.11)$$

$$q(\theta^{(q-1)}|\theta^*) = q(-\epsilon) \qquad (D.12)$$

$$q(\epsilon) = q(-\epsilon) \qquad (D.12)$$

the acceptance probability function simplifies to Equation (D.13) ([Tab21]):

$$K(\theta^{(*)}) = \min \left( \frac{f(\theta^{(*)}|.)}{f(\theta^{(q-1)}|.)}, 1 \right) \qquad (D.13)$$

### Definition D.2.2: Gibbs sampling

The fundamental concept behind Gibbs sampling is to break down the joint posterior distribution of the parameters, into a product of conditional distributions of the parameters. To that end, it is generally assumed that there exists a natural partition of the hidden parameters allowing them to factorise in a meaningful way.

It is usually implemented when the joint conditional posterior distribution is intractable to compute, whereas the conditional distribution for a subset of the parameters, conditioned on all others, is rather straightforward to derive. This is especially the case when the set of hidden parameters is linked to each other, and that numerous numerical constraints linking them must be endorsed.

Using our notations, the following joint posterior distribution $f(p, X|D)$, with $D$ denoting the observed data, here $y$ is analytically complex to derive, while $f(p|X, D)$ and $f(X|p, D)$, the posterior cellular ratios and purified individual cell expression profiles, respectively, can be simply computed. In practice, like any MCMC framework, you start to initiate the values for all the parameters. Then, the iterated Gibbs process samples each parameter (or subset of parameters), one at a time, and updates its value conditioned on the other parameters at their current values. It can be proven that after a sufficient number of iterations, the corresponding Markov chain of parameters converges and approximates well the desired joint distribution.

## D.2.2 Pseudo-code Gibbs sampler

We detail in Algorithm 0 a potential pseudo-code to generate MCMC chains of the joint distribution of the parameters of interest, in which variable $q$ denotes the running index, $B$ the number of burn-in iterations to be discarded after sampling, and $Q$ the actual length of the resulting Markov chain.
**Algorithm 0:** Pseudo-code of the iterated optimisation method, to retrieve the parameters of DeCovarT's generative model.

**Input:**

* Prior estimates of the mean, $[\mu_{ji} \in \mathbb{R}^G]$, $j \in \{1, \dots, J\}$ and covariance, $[\Sigma_{ji} \in \mathcal{M}_{\mathbb{R}^{G \times G}}]$, $j \in \{1, \dots, J\}$ of each cell population.

* Initial estimates of cellular ratios, $p_{0i}$, and purified cell expression profiles, $X_{0i}$, for each individual. They should align with both the fundamental linear deconvolution assumption (Equation (C.1)), and the unit-simplex constraint (Equation (C.2)).

* Standard deviation, $\sigma_{0(\rho)}$ and $\sigma_{0(X)}$ of the additional residual term added to each cellular ratio and each individual cell profile of the transition kernel, respectively. [SK19], [And+18], [VK21] and [Mar+20] suggests tuning these hyper-parameters such that the acceptance rates in the long term are bounded between 0.234 and 0.574.

1. **for** $q \leftarrow 1$ **to** $(B+Q)$ **do**
2. &nbsp;&nbsp;&nbsp;&nbsp;**for** $i=1:N$ **do**
3. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ \epsilon_{i(\rho)} \sim \mathcal{N}_{J-1}(0, \sigma_{0(\rho)}^2) \mathbb{I}_{J-1} $$
4. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ \rho_i^{(*)} = \rho_i^{(q-1)} + \epsilon_{i(\rho)} $$
5. &nbsp;&nbsp;&nbsp;&nbsp;**for** $j=1:(J-1)$ **do**
6. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ u \sim \mathbb{U}(0, 1) $$
7. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $u < \min(1, K_{\rho}(\rho_j^{(q-1)} \rightarrow \rho_j^{(*)}))$ **then**
8. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ \rho_{ji}^{(q)} \leftarrow \rho_{ji}^{(*)} $$
9. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else**
10. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ \rho_{ji}^{(q)} \leftarrow \rho_{ji}^{(q-1)} $$
11. &nbsp;&nbsp;&nbsp;&nbsp;**end**
12. &nbsp;&nbsp;&nbsp;&nbsp;$$ p_i^{(q)} = \psi(\rho_i^{(q)}) \quad (i) $$
13. &nbsp;&nbsp;&nbsp;&nbsp;MCMC.$p_i \leftarrow p_i^{(q)} $$
14. &nbsp;&nbsp;&nbsp;&nbsp;**end**
15. &nbsp;&nbsp;&nbsp;&nbsp;**for** $j=1:(J-1)$ **do**
16. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ \epsilon_{ji(X)} \sim \mathcal{N}_G(0, \sigma_{0(X)}^2) \mathbb{I}_G $$
17. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ x_{ji}^{(*)} = x_{ji}^{(q-1)} + \epsilon_{ji(X)} $$
18. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ u \sim \mathbb{U}(0, 1) $$
19. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $u < \min(1, K_x(x_{ji}^{(q-1)} \rightarrow x_{ji}^{(*)}))$ **then**
20. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ x_{ji}^{(q)} \leftarrow x_{ji}^{(*)} $$
21. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else**
22. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$ x_{ji}^{(q)} \leftarrow x_{ji}^{(q-1)} $$
23. &nbsp;&nbsp;&nbsp;&nbsp;**end**
24. &nbsp;&nbsp;&nbsp;&nbsp;**end**
25. &nbsp;&nbsp;&nbsp;&nbsp;$$ x_{ji}^{(q)} \leftarrow \frac{y_i - \sum_{j=1}^{J-1} x_{ji}^{(q)} p_{ji}^{(q)}}{p_{ji}^{(q)}} \quad (ii) $$
26. &nbsp;&nbsp;&nbsp;&nbsp;MCMC.$X_i \leftarrow X_i^{(q)} $$
27. &nbsp;&nbsp;&nbsp;&nbsp;**end**
28. &nbsp;&nbsp;&nbsp;&nbsp;**end**
We implemented two reparametrisations at each global iteration (outer loop) to ensure that the kernel distribution generates proposals that fall within the “support” of the target distribution, and notably the adhesion to the fundamental linear deconvolution relation (Equation (C.1)) and the unit-simplex constraint over the cellular ratios (Equation (C.2)):

1. First, the mapping function, described in (i) and in eq. (D.5), enforces that the estimated cellular ratios adhere to the unit-simplex constraint (Equation (C.2)).

2. Second, the fundamental linearity of deconvolution, Equation (C.1), is endorsed by the update formula (ii). Formula Equation (C.1) implies that the last $J$ cellular expression profile, $x_{.j}$ is not a free parameter, and can be rewritten as a combination of the others, given by Equation (D.14):

$$ 
\begin{align} 
\mathbf{y}_i &= \sum_{j=1}^{J-1} p_{ji}^{(q)} \mathbf{x}_{ji}^{(q)} + p_{ji}^{(q)} \mathbf{x}_{ji}^{(q)} = \sum_{j=1}^{J} p_{ji}^{(q-1)} \mathbf{x}_{ji}^{(q-1)} \nonumber \\
&\iff \nonumber \\
\mathbf{x}_{ji}^{(q)} &= \frac{\mathbf{y}_i - \sum_{j=1}^{J-1} \mathbf{x}_{ji}^{(q)} p_{ji}^{(q)}}{p_{ji}^{(q)}} \tag{D.14}
\end{align} 
$$

### D.2.3 Derivation of the Acceptance Probability Function

By utilising a Random Walk MH approach, additionally cancelling out the normalisation constant (Equation (D.13)), the acceptance probability function to compute is simply the product of the prior distributions and the likelihood of the observed data, $f(\theta) \times f(D|\theta)$.

To simplify further this product of distributions, we preliminary suppose that the density distribution characterising the priors is improper, in other words, that $f(\theta)$ is always equal to 1 on $\Theta$. As we provide a closed form of the log-likelihood of our generative model (eq. (D.1)), and not the likelihood, we need to apply an exponential transformation to recover the desired acceptance probability (Equation (D.15)):

$$ K(\theta^{(*)}) = \min \left( 1, \exp \left( \ell(\theta^{(*)}|.) - \ell(\theta^{(q-1)}|.) \right) \right) \quad (D.15) $$

with:

* $\theta^{(*)}$, the current proposal (either $p^{(*)}$ or $X^{(*)}$)

* $\ell(\theta^{(*)}|.) = \log(f(Y|\theta^{(*)}, \zeta))$, the log-likelihood of the currently observed data, conditioned on the current values of the latent parameters to estimate and the user-defined parameters $\zeta$, given in our generative model by Equation (D.1).

* $\zeta = (\mu, \Sigma)$, $\mu = (\mu_{.j})_{j \in \tilde{J}} \in \mathcal{M}_{G \times J}$, $\Sigma = \sum_{j=1}^{J} p_j^2 \Sigma_j \in \mathcal{M}_{G \times G}$ denote the parameters provided by the user before conducting the study.