# Inferring cellular heterogeneity with mixture models for DNA methylation rates

Hugo Barbot, David Causeur, Yuna Blum, Magali Richard

<figure><img src="image_1.png" /></figure>

To cite this version:

Hugo Barbot, David Causeur, Yuna Blum, Magali Richard. Inferring cellular heterogeneity with mixture models for DNA methylation rates. JdS, Jun 2025, Marseille, France. hal-05065027

HAL Id: hal-05065027

https://hal.science/hal-05065027v1

Submitted on 13 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# INFERRING CELLULAR HETEROGENEITY WITH MIXTURE MODELS FOR DNA METHYLATION RATES

Hugo Barbot¹ & David Causeur² & Yuna Blum³ & Magali Richard⁴

¹ IRMAR - UMR CNRS 6625, France, hugo.barbot@institut-agro.fr

² IRMAR - UMR CNRS 6625, France, david.causeur@institut-agro.fr

³ IGDR - UMR CNRS 6290, France, yuna.blum@univ-rennes.fr

⁴ TIMC - UMR CNRS 5525, France, magali.richard@univ-grenoble-alpes.fr

**Résumé.** L'hétérogénéité de la composition en types cellulaires d'échantillons biologiques est un marqueur important de la progression d'une maladie, utile pour son diagnostic. Cette composition cellulaire est cependant difficile à évaluer à partir de profils moléculaires d'un échantillon composite, la contribution de chaque type cellulaire aux signaux observés étant inconnue. La déconvolution cellulaire vise à estimer les proportions des différents types cellulaires à partir de ces profils moléculaires. Deux types de données omiques sont fréquemment employés dans cet objectif: l'expression des gènes et le taux de méthylation de l'ADN. La déconvolution cellulaire s'appuie sur l'hypothèse que le profil moléculaire de l'échantillon composite peut être approché par une somme pondérée de profils moléculaires spécifiques des mêmes gènes pour chaque type cellulaire considéré, les poids étant les proportions inconnues de ces types cellulaires. La plupart des méthodes statistiques utilisées pour la déconvolution cellulaire sont basées sur des extensions de l'algorithme des moindres carrés ordinaires, sous les contraintes de positivité et de somme à un sur les coefficients du mélange. L'utilisation de cet algorithme suppose implicitement l'indépendance, l'homoscédasticité et la normalité des erreurs résiduelles, conditions sous lesquelles il offre des garanties d'optimalité. Dans le cas présent, chacune de ces trois hypothèses est discutable. D'une part, la nature intrinsèque des données omiques requiert des modèles mieux adaptés à leur sur-dispersion. D'autre part, la structure de dépendance induite par le réseau de régulation des gènes est très forte.

Le but de ce travail est de proposer un cadre statistique respectant les caractéristiques inhérentes des données de méthylation de l'ADN. L'analyse exploratoire des distributions des taux de méthylation pour chaque type cellulaire à l'échelle du génome révèle en effet l'existence d'une structure latente en groupes des gènes, chaque composante de cette structure pouvant être associée à un modèle spécifique de déconvolution cellulaire. Nous proposons un mélange de modèles de régression Bêta estimé par maximisation sous contraintes de la vraisemblance grâce à l'utilisation d'un algorithme EM. L'application à des données *bench-mark in vitro*, générées à partir de populations cellulaires isolées présentes dans le cancer du pancréas [1], révèle que les modèles de déconvolution cellulaire spécifiques à chaque composante du mélange sont en général très différents. De plus, la sélection des gènes constitutifs de la meilleure composante conduit en général à des améliorations importantes de l'estimation des proportions des types cellulaires. La présentation abordera la question de la sélection des gènes par le choix judicieux d'une composante du mélange pour la déconvolution cellulaire.

**Mots-clés.** Déconvolution cellulaire, Modèle de mélange, régression Bêta.

**Abstract.** Cellular heterogeneity in biological tissues reflects progression of disease state
and is therefore useful for improved diagnostic and prognosis. Cellular composition of tissues
is however difficult to assess from bulk molecular profiles, with all cells present in the tissue
contributing to the recorded signals. Cell deconvolution is a common approach to unravel
the heterogeneous molecular profiles observed in bulk tissues, by inferring the underlying
relative abundance of individual cell types. Two omics data types are frequently used for cell
deconvolution: RNA-seq gene expressions and DNA methylation rates. Cellular deconvolu-
tion models assume that bulk omic profiles result from weighted sums of so-called signature
cell-specific omic profiles, weights being the unknown proportions of those cell types. Consis-
tently, most statistical methods used for cellular deconvolution are based on extensions of the
Ordinary Least Squares (OLS) optimization algorithm, under nonnegativity and sum-to-one
constraints on those unknown mixing coefficients. Using OLS implicitly assumes indepen-
dence, homoscedasticity and normality of the residual errors, conditions under which OLS
optimization guarantees optimal estimation. In cellular deconvolution applied to bulk molec-
ular profile, all three assumptions are highly questionable. Indeed, strong violations of those
assumptions may be due to the intrinsic nature of omics data, DNA methylation rates being
over-dispersed percentages, or to the dependence structure induced by the gene regulatory
network, some key genes being more influent on deconvolution accuracy than others.

The goal of this work is to provide a well defined statistical framework that respects the inherent characteristics of biological data for deconvolution, using DNA methylation rates. Whole-genome cell-type specific distributions of DNA methylation rates actually shows a latent group structure, that can explain poor estimation accuracy when fitting deconvolution models on the whole genome. Therefore, we propose a mixture of constrained beta regression models estimated using an EM algorithm. When applied to a *in vitro* benchmark dataset generated from isolated cell populations of Pancreatic cancer(see [1]), results confirms the marked sensitivity of cell deconvolution to the choice of a latent component in the mixture. Interestingly, deconvoluting using the geneset within the best latent component generally leads to a strong improvement with respect to the whole-genome deconvolution. The presentation will compare different gene selection strategies based on the choice of a latent component for cell deconvolution.

**Keywords.** Cell deconvolution, Mixture models, Beta regression.

# 1 A brief introduction to single-omic cell deconvolution

Cell deconvolution data produced for the study of Pancreatic ductal adenocarcinoma (PDAC) are used below for an illustrative purpose. This dataset is generated with the aim of serving as a benchmark reference for comparison of cell deconvolution algorithms (see [1]). Therefore, it is obtained under a strict control of the true proportions of each of the 9 cell types in the composition of each bulk. These data consist of DNA methylation rates at more than 800,000 CpG sites over the genome for each of 30 bulks and correspondingly for *p* = 9 cell types. In order to favor the future simultaneous use of DNA methylation and gene expression data in the cell deconvolution task, those methylation rates are aggregated into gene-level
measurements by averaging over all values at CpG sites in the promoter region of each gene. A profile of $m = 18735$ gene methylation rates is available on $n = 30$ independent bulks. The signature matrix contains the methylation rates of those $m$ genes in $p = 9$ cell types.

The basic principles of cell deconvolution are introduced hereafter, based on a single genetic profile of DNA methylation rates. For $1 \le i \le n$ and $1 \le j \le m$, let $y_{ij}$ denote the methylation rates of gene $j$ for bulk $i$. Let $y_i = (y_{i1}, \dots, y_{im})'$ denote the complete measurements profile for bulk $i$. The aggregated methylation rates $x_j = (x_{j1}, \dots, x_{jp})'$, $j = 1, \dots, m$ of all genes for $p$ cell types of interest is also available. In standard so-called supervised cell deconvolution models, the methylation profiles in bulks are assumed to result from a linear combination of the signature cell types, the mixing coefficients of this linear decomposition being the unknown proportions $\beta_i = (\beta_{i1}, \dots, \beta_{ip})'$ of each cell type within each bulk, up to an additive error term. Consistently, most cell deconvolution approaches are variants of the following constrained least-squares minimization issue:

$$ (\hat{\beta}_{0i}, \hat{\beta}'_i) = \operatorname{argmin}_{(\beta_{0i}, \beta'_i)} \sum_{j=1}^{m} (y_{ij} - \beta_{0i} - x'_j \beta_i)^2, $$

where $\beta_{0i}$ is an intercept and the coefficients $\beta_i$ are constrained to lie within the $p$-simplex $S_p = \{\beta = (\beta_1, \dots, \beta_p), 0 \le \beta_r \le 1, \sum_{r=1}^p \beta_r = 1\}$. The nonnegativity and sum-to-one constraints on $\beta$ makes the above constrained optimization issue more challenging, which explains the variety of algorithmic solutions available for this task [2].

Unconstrained least-squares optimization is well-studied and guarantees desirable properties, such as unbiasedness and minimum variance, under the standard assumptions of the linear regression model. However, independence, homoscedasticity and normality of the residual error terms are required to ensure optimality of least squares estimation. In the present situation, all three assumptions are highly questionable. Indeed, whereas in the standard approach of omic data analysis, genes are usually considered as features measured on independent statistical units being different biological samples, in cell deconvolution methods, statistical units are genes and features are bulks. Yet, the distribution of DNA methylation rates notoriously shows overdispersion and is driven by a gene regulatory network that induces a graph-structured stochastic dependence pattern across genes.

## 2 Beta-regression based cell deconvolution

For gene $j$, given the $p$-profile $x_j = (x_{j1}, \dots, x_{jp})'$ of signature DNA methylation rates in the $p$ cell types of interest, it is thereafter assumed that the methylation rates $y_i, i = 1, \dots, n$, in bulk $i$ are independent and distributed according to a Beta distribution with density [3, 4]:

$$ \varphi(y_{ij} | x_j) = \frac{\Gamma(\phi_i)}{\Gamma(\mu_i(x_j)\phi_i)\Gamma((1-\mu_i(x_j))\phi_i)} y_{ij}^{\mu_i(x_j)\phi_i - 1} (1 - y_{ij})^{(1-\mu_i(x_j))\phi_i - 1}, \quad (1) $$

where $\phi_i > 0$, $\mu_i(x_j) = \beta_{0i} + \beta_{1i}x_{j1} + \dots + \beta_{pi}x_{jp} > 0$, $\beta_{0i}$ is an intercept parameter and $\beta_i = (\beta_{i1}, \dots, \beta_{ip})' \in S_p$ is the vector of proportions of each cell type.
A Cyclic Coordinate Descent (CCD) algorithm is implemented to optimize the log-likelihood of model (1). In order to ensure nonnegativity of estimation, at each update of an estimated proportion parameter, if the marginal maximization of the log-likelihood provides a negative update, then the current value of the proportion parameter is set to zero.

Maximum-Likelihood estimations of the proportions of the nine reference cell types in the 30 bulks of the PDAC study are now calculated, using the above CCD algorithm under the assumption of the beta regression model (1). In the present situation where the true proportions of cell types in each bulk are controlled and therefore can be assumed to be known, cell deconvolution performance is evaluated using Mean-Squared-Errors (MSEs) of estimation, for each bulk over the nine cell types.

Figure 1 (left plot) displays boxplots of MSEs over the 30 bulks for NNLS, a standard deconvolution algorithm, and the beta regression estimation, referred to as BR. It shows that BR is slightly outperformed by NNLS, both providing poorly accurate estimates of the reference cell type proportions on average over the 30 bulks. Indeed, the right plot in Figure 1 illustrates this lack of accuracy for bulk 18, for which the MSE is close to the median MSE over the 30 bulks (0.006).

<figure><img src="image_3.png" /><figcaption>Figure 1: PDAC study. Left plot: estimation accuracy for cell deconvolution using NNLS and BR - Distribution of MSEs over the 30 bulks. Right plot: comparison of true and estimated proportions, using BR, for bulk 18.</figcaption></figure>

Kernel-based estimations of the density functions of the cell type-specific distributions of DNA methylation rates (see left plot in Figure 2) actually questions the assumption of a whole-genome unique deconvolution model. It can be hypothesized that the DNA methylation rates in the signature cell types are not homogeneously distributed according to a single whole-genome Beta distribution, but according to a mixture of distributions whose components are related to latent classes of genes.

Consistently, a mixture model of constrained Beta regression models is assumed thereafter.
<figure><img src="image_4.png" /><figcaption>Figure 2: PDAC study. Left plot: marginal distributions of methylation rates in each cell type. Right plot: EM estimation of a 4-component mixture of Beta distributions for DNA methylation rates in Bulk 1, which is the number of component minimizing BIC between 1 and 6.</figcaption></figure>

Let $\mathcal{C} = (\mathcal{C}_1, \dots, \mathcal{C}_K)$ denote a K-latent class variable, distributed according to a multino- mial distribution $\mathcal{M}(1, \pi_1, \dots, \pi_K)$, with $K \le 1$, $\pi_k = \mathbb{P}(C_k = 1) > 0$ and $\sum_{k=1}^K \pi_k = 1$. Conditionally on $\mathcal{C}$, the bulk DNA methylation rates $Y_{ij}, j = 1, \dots, m$ are assumed to be dis- tributed according to a Beta distribution, conditionally on the p-profile of DNA methylation rates $x_{jr}, j = 1, \dots, m, r = 1, \dots, p$ in the reference cell types:

$$Y_{ij} | \{C_{jk} = 1\} \sim \mathcal{B}(\mu_i^{(k)}(x_j) = \beta_{0i}^{(k)} + \beta_{1i}^{(k)} x_{j1} + \dots + \beta_{pi}^{(k)} x_{jp}, \phi_i^{(k)}).$$

where $\beta_{0i}^{(k)}, k = 1, \dots, K$ and $\beta^{(k)} = (\beta_{1i}^{(k)}, \dots, \beta_{pi}^{(k)}) \in \mathcal{S}_p$ are respectively component-specific intercepts and vectors of proportions for the p cell types in bulk $i$. The log-likelihood of the above mixture model is maximized using an EM algorithm, in which the M step consists of weighted maximization of componentwise constrained Beta regression models using the CCD algorithm introduced above.

Interestingly, componentwise estimated proportions in early steps of the algorithm offer markedly different viewpoints of the bulk cellular heterogeneity. Both estimated proportions and the corresponding MSEs at the initialisation step are reproduced in Table 1 for Bulk 1 chosen arbitrarily for illustration. Here, the first two components show a markedly better accuracy than the others compared to the whole genome estimation. More generally, for the 30 bulks, it turns out that at least one of the latent component improves accuracy up to 45% on average.

The presentation will address the dual question of the best stopping criterion for the EM algorithm and the selection of the latent component favoring optimal estimation accuracy.
<table><caption>Table 1: PDAC study. Bulk 1. Whole-genome and componentwise cell deconvolution using BR. The true proportions of the 9 reference cell types are provided in the last column. The last row gives the MSE of estimation.</caption><thead><tr><th rowspan="2">Cell types</th><th rowspan="2">Whole genome</th><th colspan="4">Components</th><th rowspan="2">True proportions</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><td>Macrophages</td><td>0.00</td><td>0.17</td><td>0.03</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>Neutrophils</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>B cells</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>CD4 T cells</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.01</td><td>0.01</td></tr><tr><td>CD8 T cells</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.01</td><td>0.04</td></tr><tr><td>Cancer classical</td><td>0.11</td><td>0.12</td><td>0.12</td><td>0.11</td><td>0.11</td><td>0.05</td></tr><tr><td>Cancer basal</td><td>0.27</td><td>0.12</td><td>0.21</td><td>0.26</td><td>0.14</td><td>0.21</td></tr><tr><td>Endothelial</td><td>0.48</td><td>0.25</td><td>0.40</td><td>0.48</td><td>0.53</td><td>0.30</td></tr><tr><td>Fibroblasts</td><td>0.15</td><td>0.34</td><td>0.23</td><td>0.15</td><td>0.20</td><td>0.38</td></tr><tr><td>MSE (×10<sup>-2</sup>)</td><td>1.08</td><td>0.51</td><td>0.49</td><td>1.07</td><td>1.07</td><td></td></tr></tbody></table>

References

[1] Elise Amblard, Vadim Bertrand, Luis Martin Pena, et al. A robust workflow to benchmark deconvolution of multi-omic data. *bioRxiv*, 2024.

[2] Clémentine Decamps, Alexis Arnaud, Florent Petitprez, et al. DECONbench: a benchmarking platform dedicated to deconvolution methods for tumor heterogeneity quantification. *BMC Bioinformatics*, 22(1):473, October 2021.

[3] Timothy J Triche Jr, Peter W Laird, and Kimberly D Siegmund. Beta regression improves the detection of differential dna methylation for epigenetic epidemiology. *BioRxiv*, page 054643, 2016.

[4] Silvia Ferrari and Francisco Cribari-Neto. Beta regression for modelling rates and proportions. *Journal of applied statistics*, 31(7):799-815, 2004.