Recent advances in high-throughput sequencing, single-cell isolation, and bioinformatics have enabled multimodal, single-cell-level interrogation of biological systems (Wu et al., 2014; Hong et al., 2020; He et al., 2025). Among single-cell technologies, single-cell RNA sequencing (scRNA-seq) stands out for its ability to resolve transcriptional states at cellular resolution, permitting precise identification of cell types and subpopulations, reconstruction of developmental trajectories, and dissection of molecular mechanisms underlying health and disease (Lopez et al., 2018; Ranjan et al., 2021). As single-cell studies scale across tissues, conditions, and laboratories, robust computational methods (Kinker et al., 2020; Flores et al., 2022) that can extract meaningful biological signals from noisy, sparse measurements are increasingly essential for biological discovery and translational applications such as biomarker identification and therapeutic target prioritization (Haque et al., 2017; Su et al., 2022).

Clustering remains a foundational step in scRNA-seq analysis, but it is challenged by properties inherent to these data: extreme sparsity (many zero or near-zero counts), high dimensionality (tens of thousands of genes per cell), andmixed sources of variation, including true biological heterogeneity and technical noise (dropout events, varying capture efficiency, and limited sequencing depth) (Conesa et al., 2016; Qi et al., 2020; Zhao et al., 2021; Ghorbani et al., 2024; Gong et al., 2018). Sparsity and high dimensionality exacerbate the curse of dimensionality, reducing the discriminative power of conventional distance metrics and degrading the performance of classical clustering methods. Moreover, zero inflation and measurement noise obscure subtle but biologically important gene–gene relationships that are critical for accurate cell-type separation and downstream interpretation (Stegle et al., 2015; Camara, 2018).

Existing ZINB-based models primarily focus on modeling count distributions and dropout events, but often overlook explicit modeling of gene–gene dependencies. Conversely, recent deep clustering methods emphasize representation learning but typically rely on generic reconstruction objectives that are insensitive to biological sparsity patterns. As a result, existing approaches struggle to simultaneously address technical noise, zero inflation, and contextual gene relationships within a unified framework.

To address these challenges, we introduce scDMAC, a unified framework that couples principled probabilistic denoising with contextual masked reconstruction to produce compact, biologically informative embeddings for clustering. Specifically, scDMAC differs fundamentally from existing composite methods such as scDeepCluster and scziDesk in three aspects:

The ZINB-based denoising module in scDMAC is used as an explicit pre-denoising stage rather than being jointly optimized with clustering, which stabilizes subsequent representation learning;

Together, these design choices allow scDMAC to address both zero inflation and gene dependency learning in a coordinated manner, going beyond architectural variations of existing deep clustering approaches. We evaluate scDMAC on multiple widely used scRNA-seq benchmarks and demonstrate consistent improvements in clustering accuracy, stability, and robustness to dropout compared with state-of-the-art methods. scDMAC delivers clearer separation of canonical cell types, more reliable identification of rare populations, and improved reproducibility across noisy conditions. Collectively, these results show that combining a statistically grounded noise model with masked contextual learning is an effective strategy for extracting biologically meaningful embeddings from scRNA-seq data, thereby improving downstream tasks such as cell-type annotation, trajectory inference, and differential expression analysis.

Figure 1 shows the flowchart of the single-cell RNA sequencing cluster method based on denoising and masking learning. To optimize single-cell RNA sequencing (scRNA-seq) data for contrastive learning models, the data first undergoes normalization and log transformation, followed by gene filtering, and finally construction of a k-nearest neighbor (KNN) graph based on cosine distance.

scRNA-seq data often exhibits substantial variation in total gene expression per cell due to differences in sequencing depth, which compromises the comparability of expression values across cells. To mitigate this technical bias, expression values are normalized. Specifically, for the expression value X_ij of gene j in cell i, the normalized value is calculated in Equation 1. X i j n o r m = X i j ∑ j = 1 G X i j · s 0 (1) where s ₀ is a scaling factor (set to 10,000 in this study). This adjusts the total expression of each cell to a common scale, reducing the impact of sequencing depth.

Log1p normalization, which is defined as Equation 2, is applied after denoising and is only used for downstream masked representation learning and clustering. This separation ensures both statistical validity of the ZINB likelihood and numerical stability for deep representation learning. X i j log = ln X i j n o r m + 1 (2)

This helps approximate a normal distribution, making the data more suitable for contrastive learning. Subsequently, highly variable genes are selected using the Scanpy package to minimize the influence of uninformative features. The resulting preprocessed matrix X serves as input to the model. Finally, a KNN graph is constructed using cosine distance to represent cell neighborhoods.

At the same time, this study uses data augmentation to enhance model performance by generating variations of the original data. For scRNA-seq data, designed augmentation helps simulate real-world data distribution, improving imputation and downstream tasks such as clustering. This section applies three augmentation strategies:

Masking Gene Expressions: Randomly selected genes (10%) have their expressions set to zero. This mimics “dropout” events common in scRNA-seq data, encouraging the model to learn contextual relationships for predicting masked values.

Adding Gaussian Noise: To simulate technical variability from sequencing or handling, Gaussian noise with a variance of 0.6 is added, which improves model robustness to noise, especially useful in datasets with high technical variation.

Swapping Expressions Between Neighboring Cells: Based on the KNN graph, a cell’s expressions are randomly swapped with those of its neighbors at a ratio of 0.2. This promotes local structural variability and reduces over-reliance on fixed neighborhood patterns.

To denoise the gene expression matrix X and capture key characteristics of scRNA-seq data, such as high sparsity and overdispersion. A Zero-Inflated Negative Binomial (ZINB)-based autoencoder is employed, as shown in Figure 2. This model integrates an encoder with a denoising autoencoder architecture inspired by DCA, enhancing its ability to handle scRNA-seq noise and dropout effects. The ZINB module probabilistically models dropout events via the zero-inflation parameter, capturing technical zeros, while the masking strategy serves as a self-supervised regularization mechanism rather than an explicit zero generator. Masked values are only introduced during training and are not interpreted as biological zeros. By decoupling probabilistic dropout modeling from masking-induced perturbations, The preprocessed expression matrix X is input into a deep count autoencoder, which uses a ZINB-based loss to reconstruct a denoised expression matrix X _z . The ZINB distribution is parameterized by three components: mean (μ), dispersion (θ), and dropout probability (π). The probability of an observed count is defined as Equations 3, 4. NB X | μ , θ = Γ X + θ X ! Γ θ θ θ + μ θ μ μ + θ X (3) ZINB X | π , μ , θ = π δ X + 1 − π NB X | μ , θ (4) where Γ denotes the gamma function and δ(X)represents a point mass at zero. Unlike standard autoencoders, the ZINB model uses three separate fully connected output layers connected to the decoder’s final hidden layer to estimate the three parameters, as Equations 5–7. Θ = exp W θ D (5) M = diag s i × exp W μ D (6) Π = sigmoid W π D (7)

Here, Θ, M, and Π denote the matrices of dispersion, mean, and dropout probability, respectively. W _μ , W _θ and W _π are the weight matrices for each parameter head, and D is the output from the last hidden layer of the decoder. The exponential function ensures non-negativity for mean and dispersion, while sigmoid constrains the dropout probability to [0,1]. The scaling factor s _i , obtained during preprocessing, adjusts for cell-specific library size.

The loss function for denoising is the negative log-likelihood of the ZINB distribution as Equation 8. L z i n b = − log ZINB X | π , μ , θ (8)

Following the denoising step, variability and perturbation are introduced into the denoised gene expression matrix X _z through the following procedure:

First, the expression values of each gene are randomly shuffled within the matrix to preserve intra-gene correlations, resulting in a perturbed matrix X'.

Next, a masking matrix M is generated using a Bernoulli distribution Bernoulli (p _j ) for each gene as Equation 9, where p _j controls the masking probability for the jth gene: M i j ∼ Bernoulli p j (9)

Here, M _ij represents the element in the ith row and jth column of the mask.

Finally, the masked gene expression matrix X _M is obtained via element-wise operations as Equation 10. X M i j = X i j ′ · M i j + X z · 1 − M i j (10) where X M i j is an element of the masked matrix, and X i j ′ and X _z are elements from the shuffled and denoised matrices, respectively.

Importantly, the ZINB-based denoising autoencoder is trained on raw count data prior to masking, and masked entries are excluded from the ZINB likelihood. Masking is applied only to the denoised output for subsequent self-supervised representation learning.

Figure 3 shows the masking autoencoder, which consists of three main components: an encoder, a mask predictor, and a decoder. The encoder transforms the masked gene expression matrix X _M into a low-dimensional embedding Z. For an encoder with F layers, the output of the fth layer is computed as Equation 11. Z f = σ W f Z f − 1 + b f (11) where σ is the activation function. The final layer applies a linear transformation (i.e., identity activation), and its output Z _F serves as the embedding Z.

To address potential inaccuracies in the masked input, the model first uses a mask predictor to estimate which expression values have been modified, producing a predicted mask matrix M'. It is implemented as a linear layer trained with cross-entropy loss as Equation 12. L m = − ∑ i j M i j ⁡ log M i j ′ (12)

The decoder reconstructs the gene expression matrix using the embedding Zand the predicted mask M'. A weighted mean squared error (MSE) loss is applied to emphasize masked genes as Equation 13. And weight W_ij is defined as Equation 14. L r e c = 1 N ∑ i j W i j · X i j − X ∼ i j 2 (13) W i j = λ · M i j + 1 − λ · 1 − M i j (14)

Here, λ is a hyperparameter assigning a higher weight to masked genes. The total loss is a weighted combination of the two objectives as Equation 15. L m a s k = γ m · L m + 1 − γ m · L r e c (15) where γ _m balances the two terms. During clustering, only the embedding Z generated by the encoder is used.

To enhance the clustering performance of the model, a weighted soft clustering module is introduced. This module employs a weighted K-means approach to assign data points (cells) to cluster centers while preserving local similarity structures among cells with comparable gene expression profiles. The weighted K-means loss is defined as Equation 16. L k m e a n s = ∑ k = 1 n ∑ k = 1 C w i k z i − c k 2 (16) where w _ik is the weight for cell and cluster, z _i is the embedding of cell i, and c _k is the center of cluster k, updated by Equation 17. c k = ∑ k = 1 n w i k z i k ∑ k = 1 n w i k (17)

The weight is computed using Equation 18. w ¯ i k = exp − z i − c k 2 ∑ k = 1 C ⁡ exp − z i − c k 2 (18)

Subsequently, the weights are sharpened using a Markov-like inflation step by Equation 19. w i k = w ¯ i k α ∑ k = 1 C w ¯ i k α (19) where is a hyperparameter (default 1). To better capture similarity relationships, a student’s t-distribution is used to model pairwise cell similarities. The soft assignment probability q _ij is given by Equation 20: q i j = 1 + z j − z i 2 / t − t + 1 2 ∑ c ≠ i 1 + z j − z i 2 / t − t + 1 2 (20) with t set to 1. A target distribution p is derived from q to strengthen high-confidence assignments as Equation 21. p i j = q i j 2 / ∑ i ≠ j q i j ∑ c ≠ i q i j 2 / ∑ i ≠ j q i j (21)

The Equation 22 defines the clustering loss, the KL divergence between p and q: L c l u s t e r = ∑ i ∑ j p i j ⁡ log p i j q i j (22)

The overall training objective combines all losses as Equation 23. L = α L z i n b + β L m a s k + φ L k m e a n s + θ L c l u s t e r (23) where α, β, φ, θ are tunable hyperparameters.

This study evaluates the model’s clustering and classification performance using six annotated scRNA-seq datasets from both mouse and human. These datasets cover multiple biological systems, diverse cell types, and varying cluster sizes, and were generated using different RNA extraction protocols and sequencing platforms (e.g., Smart-seq, Drop-seq). The datasets—Adam, Deng, Muraro, Pollen, Chen, and Zeisel—are summarized in Table 1.

For fair comparison, all baseline methods were implemented using recommended settings from their original publications or official repositories. Highly variable genes (HVGs) were selected using Scanpy with default parameters unless otherwise specified. Latent dimensionality was set to 10–32 depending on method defaults.

Graph-based methods (Seurat, graph-sc) used k-nearest neighbor graphs with k = 15–30 and Leiden clustering with default resolution. Deep clustering baselines were run with identical train/validation splits and optimized using Adam. For stochastic methods, each experiment was repeated multiple times with different random seeds, and the best-performing configuration was reported.

The base environment of the experimental platform utilizes an Intel Xeon E5-2,630 v4 CPU, a Tesla 2080Ti GPU with 24 GB of VRAM, and 64 GB of memory. The operating system is Ubuntu 18.04.6, PyTorch version 1.13.1, and Python 3.8.16.

This paper proposes scDMAC, a clustering model for single-cell RNA sequencing data that integrates a denoising autoencoder with a masking autoencoder. The model first denoises scRNA-seq data using a ZINB-based denoising autoencoder to better approximate the underlying expression distribution. It then introduces variability by randomly shuffling expression values within genes and applies a Bernoulli-based masking strategy to generate perturbed gene expression profiles. These are encoded into low-dimensional embeddings through a masking autoencoder, which jointly optimizes feature reconstruction and mask prediction. Finally, a weighted soft clustering mechanism is applied to produce the clustering results.

Experimental results demonstrate that scDMAC achieves improved performance by effectively capturing gene-wise relationships and enhancing feature robustness. While simulation studies provide controlled ground truth, they often fail to capture the complex noise structure and biological heterogeneity of real scRNA-seq data. In this work, we prioritize evaluation on well-annotated benchmark datasets that are widely used in literature. Nevertheless, we acknowledge this limitation and plan to include simulation-based validation in future work.