A versatile informative diffusion model for single-cell ATAC-seq data generation and analysis

Authors: lei  huang
lei  xiong
na  sun
zunpeng  liu
ka-chun  wong
manolis  kellis
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-02 17:24:33

Abstract: The rapid advancement of single-cell ATAC sequencing (scATAC-seq) technologies holds great promise for investigating the heterogeneity of epigenetic landscapes at the cellular level. The amplification process in scATAC-seq experiments often introduces noise due to dropout events, which results in extreme sparsity that hinders accurate analysis. Consequently, there is a significant demand for the generation of high-quality scATAC-seq data in silico. Furthermore, current methodologies are typically task-specific, lacking a versatile framework capable of handling multiple tasks within a single model. In this work, we propose ATAC-Diff, a versatile framework, which is based on a diffusion model conditioned on the latent auxiliary variables to adapt for various tasks. ATAC-Diff is the first diffusion model for the scATAC-seq data generation and analysis, composed of auxiliary modules encoding the latent high-level variables to enable the model to learn the semantic information to sample high-quality data. Gaussian Mixture Model (GMM) as the latent prior and auxiliary decoder, the yield variables reserve the refined genomic information beneficial for downstream analyses. Another innovation is the incorporation of mutual information between observed and hidden variables as a regularization term to prevent the model from decoupling from latent variables. Through extensive experiments, we demonstrate that ATAC-Diff achieves high performance in both generation and analysis tasks, outperforming state-of-the-art models.

Main paper: 
 introduction: 
 Assay for Transposase Accessible Chromatin with sequencing (ATAC-seq) data has shed light on genome research to dissect gene regulatory landscapes and cellular heterogeneity. Particularly, singlecell ATAC-seq (scATAC-seq) can be harnessed to probe the chromatin accessibility profile at singlecell level to reserve the diversity of cell types in a heterogeneous tissue, which can reveal important genomic elements for transcription factor binding and regulating downstream gene expression, particularly at distal non-coding regulatory regions [1, 2] . However, the analysis of such data poses challenges due to the high levels of noise, sparsity, and data scale encountered. Furthermore, the complexity of biological systems, with their intricate intercellular communications and numerous through the time step setting t = 1...T . The transform distribution is formulated as a fixed Markov chain which gradually adds Gaussian noise to the data with a user-defined variance schedule:EQUATIONwhere x t is the state of x at the time step t which is obtained by adding Gaussian noise to the former state x t-1 and β t controls the degree of the noise. The distribution delineates a Gaussian distribution centered on the incrementally corrupted data state x t . We could obtain a closed form from the input data x 0 to x T in a tractable way. The posterior distribution could be factorized as:EQUATIONLet ᾱt = t s=1 1 -β s , we could sample x t at any arbitrary time step t in a closed form by utilizing the reparameterization trick [9] :EQUATIONIf the time step is sufficiently large, the final distribution will approach that of a standard Gaussian distribution. The reverse process aims to reconstruct the original data x 0 from the diffused data x T , sampled from Gaussian distribution N (0, I), which is accomplished through the diffusion process.The reverse process can also be factorized as a Markov chain:EQUATIONWe aim to investigate the iterative reverse process p(x t-1 |x). However, Estimating the distribution p (x t-1 | x t ) is hard to estimate unless the gap between t -1 and t is infinitesimally small (T = ∞) [10] . Therefore, a learned Gaussian transitions p θ (x t-1 | x t ) is devised to approximate p (x t-1 | x t ) at every time step:EQUATIONFollowing previous work ( [9] ), µ θ (x t , t) can be modeled as:EQUATIONwhere ϵ θ is a neural network w.r.t trainable parameters θ. Having formulated the reverse process, we could maximize the likelihood of the training data as our object. Since directly calculating the likelihood is intractable, we adopt the evidence lower bound (ELBO) [9] to optimize.EQUATION3 MethodIn this section, we formally present our proposed model ATAC-Diff for scATAC-seq data analysis and generation. As illustrated in Figure 1 , ATAC-Diff is based on the conditional latent diffusion model equipped with the auxiliary module which provides the additional latent variables. Specifically, we use fragment counts to represent the scATAC-seq data. Then we compress the scATAC-seq data (fragment counts) into lower dimensional latent space. Here, we adopt the autoencoder (AE) framework, where the encoder Enc AE encodes raw data x raw into latent domain x 0 = Enc AE (x raw ) and the decoder Dec AE learns to decode x 0 back to raw data x raw . This framework can be trained by minimizing the reconstruction objective. The auxiliary encoder z = Enc ϕ (x 0 ) learns to map an input ATAC-seq data to a semantically meaningful representation z. By incorporating a latent variable of this nature, the diffusion model can prioritize high-level semantic information during the generation process, thereby producing data of superior quality. Additionally, this approach enables the utilization of meaningful representations at the hidden layer for downstream task analysis. Our work is inspired by the recent success of Diffusion Autoencoders [19] , but learning latent variables for the sparse scATAC-seq data is however challenging. We address this by introducing GMM to extract the distinct features from the latent space and utilize mutual information to maximize the shared information between latent variables and real scATAC-seq data. We elaborate on the auxiliary module in Section 3.1 and the conditional diffusion model in Section 3.2. Finally, we briefly describe the training and sampling scheme in section 3.3.X !"# X ! X $ X %Auxiliary Decoder 𝒛Maximization 𝐼(𝑿 $ , 𝑧)Conditional LDM𝑝 & 𝑿 !"# 𝑿 ! , 𝑧, 𝑦 𝑞 𝑿 ! 𝑿 !"# , 𝑿 $Sample 𝒛F igure 1: Overview of ATAC-Diff framework. ATAC-Diff is a versatile informative diffusion model of scATAC-seq data analysis and generation. It leverages the auxiliary model to yield low-dimensional meaningful semantic variables as the conditional information to help generate high-quality data while the reserved potential biological information could be applied for downstream tasks like cell type identification. 
semantic encoder: 
 The goal of the informative auxiliary module is two-fold. The first is to summarize an input scATACseq data into a semantic representation z = Enc ϕ (x 0 ) which contains fine-grained information to assist the diffusion decoder p θ (x t-1 |x t , z) denoise and predict the output scATAC-seq data. The second one is that the meaningful representation yielded by the auxiliary module could be utilised for downstream analysis such as cell relationships visualization and cell heterogeneity identification.The latent variables z, unlike the latent variables in the diffusion process, are flexible and can represent a low-dimensional vector of latent factors of variation which could be either continuous or discrete. In order to infer z, we design the approximate posterior q θ (z|x 0 ) which could be any architecture for the encoder. In our experiments, the encoder shares the same architecture in the diffusion backbone as the transformer encoder, which is formulated as:EQUATION)where FFL is the feed-forward layer, MHA is the multi-head self-attention layer.To learn the latent distribution of the input data and enable the unconditional sampling of x 0 in the diffusion model, we model the semantic prior p(z) over data features. Usually, the prior over the latent variables is commonly an isotropic Gaussian [22] . However, utilizing a Gaussian distribution as the prior probability distribution p(z) restricts the latent representation to a single mode while the raw features of the scATAC-seq data are multi-modal with different cell types [23] . Thus, we apply GMM as the posterior to learn multi-modal latent representations [8, 24] . We predefined categorical c ∈ {1, ...K} where K is a predefined number of components in the mixture, e.g., the number of cell types. The joint probability p(x 0 , z, c) could be modeled as:EQUATIONwhere p(z|c) is mixture of Gaussian distribution parametrized by a series of µ c and σ c . Then we define each factorized probability as:EQUATIONWe employ Kullback-Leibeler (KL) divergence to calculate the regularization term to enforce the latent variable z to the GMM manifold: R GMM = D KL (q ϕ (z, c | x 0 )∥p(z, c))Diffusion models can be considered as a special realisation of hierarchical VAE [25] which the encoder contains no learnable parameters. The iterative decoder model a complex decoding distribution p θ (x 0 |x 1:T , z) but suffers from the ignoring of low-dimensional latent variables z [20] since highdimensional x 0 does not depend on z, leading to degrading to unconditional generation. We tackle this issue by maximizing the mutual information (MI) between z and x 0 , assuming meaningful semantic variables reserve high MI with the real data point. We define MI of x 0 and z as:EQUATIONwhere H(z) and H(z|x 0 ) are the marginal information entropy and conditional information entropy, and q ϕ (z) is the parameterized posterior of mixture of Gaussian distribution. Maximization of the mutual information enables the model to generate x 0 which can infer z, thus avoiding the ignorance issue.Finally, we employ an auxiliary decoder to reconstruct the latent variables to recover the real data point. The auxiliary decoder could help the auxiliary module to consistently generate semantic variables z by maximizing the likelihood of p(x 0 |z), which can be considered as complementary information to mutual information I(x 0 ; z). Furthermore, this could prevent the latent variable z from generating into the pure mixture of Gaussians, which in turn interferes with the generation of x 0 . We define the reconstruction loss as:EQUATIONHaving formulated the latent auxiliary module, we conditioned the diffusion models on the latent variables z and other conditional information y such as cell types, tissue, and other omics data (e.g. scRNA-seq data). Then the learning objectivity is converted to p θ (x t-1 | x t , z, y). Hence, Eq. 4 and Eq. 5 becomes the conditional format:EQUATIONEQUATIONSince the diffusion model predefined the diffusion process by the user-defined variance schedules, we could encode x 0 to the stochastic x t by running the deterministic process in Eq. 3:EQUATION)where ϵ is sampled from N(0, I). Similar to the auxiliary module, we also utilize transformer to model µ θ (x t , t, z, y). For the conditional information y embedding, we also utilize MLPs to convert the raw representation of the features to the final conditional set: ϕ(y) = {MLP i (y)}. where i ∈ 1, ...L.Then we utilize the cross-attention mechanism to enable the prediction of x 0 to be conditioned on the specific attribute information y and z.EQUATIONwhere W Q , W K , W V are parameterized linear layers. For the categorical condition information such as cell type, we train a classifier on the latent embeddings z, and then select the corresponding the latent embeddings for certain cell type.Having formulated the informative auxiliary and the conditional diffusion decoder, the training of the reverse process is performed by optimizing the evidence lower bound (ELBO) on negative loglikelihood (conditional form of Eq. 7) with the auxiliary module regulation. We train ATAC-Diff by the following form:EQUATIONThe full derivation can be found in the Appendix. The third term and fourth term can be considered as the ELBO objectivity of VAE model, which are feasible to calculate. However, the last term is intractable to calculate since we cannot evaluate q ϕ (z). Following [20] , we could sample z ∼ q ϕ (z|x 0 ) which x 0 ∼ q(x 0 ) and then optimize it by likelihood free optimization techniques [26] .Empirically, the ELBO of diffusion models could be simplified as:L simple t = E x0,xt ∥x t -x θ (x t , t)∥ 2 = E x0,xt ∥x t -x θ ( √ ᾱt x 0 + √ 1 -ᾱt x t , t, z, y)∥ 2In fact, we can directly model µ θ by utilizing x θ instead of ϵ θ since Eq.6 can be rewritten as:EQUATIONIn practice, we find that predicting x 0 will enable the model to converge faster. We speculate that the reason for this phenomenon is that the scATAC-seq data is highly sparse with noise, predicting ϵ from x t which is close to the white noise is very hard, leading to a decrease in sampling quality. For the reconstruction loss of the AE, we adopt L2 norm loss.ATAC-Diff is different from vanilla Diffusion models which is conditioned on the latent variables. In order to sample from ATAC-Diff model, we need to design a sampling strategy to sample z from the latent distribution. Since the prior distribution of latent variables is the mixture of Gaussians, it is hard to sample z for unconditional generation. Therefore, we could train any arbitrary generative model to sample z. In this study, we employ a vanilla VAE to sample the latent variables. For data denoising and imputation, we just use the same data as the inputs for both diffusion model and the auxiliary module. In general, we could calculate the mean of the reverse Gaussian transitions u θ by Eq. 18. To sample the scATAC-seq data, we first sample the chaotic state x T from N (0, I) or use the desired noised data as inputs. z is sampled by the generative model or has the same value of x T . The next less chaotic state x t-1 is generated by Eq. 5. The final state x 0 is iteratively sampled x t-1 for T times. 
experiments: 
 In this section, we evaluate our proposed model ATAC-Diff across a range of experiments on three benchmark datasets, utilizing metrics which span generation quality, denoising and imputation effects, and latent space representation analysis. The extensive experimental results suggest that ATAC-Diff achieves higher or competitive performances compared with SOTA models which are designed for individual tasks.We adopt three datasets to benchmark our model and baseline models: Forebrain [27] , Hematopoiesis [28] , and PBMC10k3 . Forebrain dataset is derived from P56 mouse forebrain. Hematopoiesis dataset includes 2,000 cells during hematopoietic differentiation through FACS. PBMC10k dataset comprises peripheral blood mononuclear cells isolated from a healthy donor, with granulocytes selectively removed through cell sorting, resulting in a cell population of approximately 10,000 cells for detailed analysis.For the latent representation analysis, we extract the latent representation of cells. We evaluate the clustering results of latent representations by Normalised Mutual Information (NMI), Adjusted Rand Index (ARI), Homogeneity score (Homo) and Average Silhouette Width(ASW). To evaluate the similarities between the generated and real cells, we adopt Spearman Correlation Coefficient (SCC) and Pearson Correlation Coefficient (PCC) for performance comparison.We explore the quality of the latent representations by clustering to examine if it could be used to identify cell types. For performance comparison, we apply Highly Variable Peak (HVP) method from SCANPY [29] , cisTopic [7] , SCALE [8] , and PeakVI [5] to obtain the corresponding dimensionality reduction features. Specifically, we leverage K-Means to obtain the clusters. We evaluate the clustering performance based on NMI, ARI, Homo, and ASW. The clustering results of ATAC-Diff and the baseline models on 3 benchmark datasets are presented in Table 8 . The first, second and third highest values are colored by red, green and purple. The higher the value if the metrics, the better the clustering performance performance. Among the 3 datasets based on the four metrics, ATAC-Diff yields 8 best scores and comparable scores compared with the baseline models. For instance, ATAC-Diff outperforms all the baseline models across all the metrics on Hematopoiesis dataset, defeats SOTA models by at least 1.25% in NMI, 7.66% in ARI, 0.91% in Homo, and 10.16% in ASW. Furthermore, we observe that ATAC-Diff and SCALE achieve similar performance on Forebrain dataset but ATAC-Diff surpasses it and achieves competitive performances on PBMC10k dataset. To determine whether the methods effectively separate the cell types in the latent space, we employ UMAP to visualize the extracted features (Figure 2 ). Our analysis reveals that ATAC-Diff exhibits superior performance in distinguishing the cell types, whereas HVP fails to differentiate between the cell types. Furthermore, some of the other methods show overlapping results for certain cell types. In addition, our findings demonstrate that ATAC-Diff has the ability to unveil the distance between distinct cell subpopulations and their developmental trajectories. For instance, we observed that the three excitatory neuron cell types (EC1, EC2, EC3) from the Forebrain dataset exhibit close proximity to each other in the latent space of ATAC-Diff. Conversely, the distances between different cell types are significantly larger, indicating distinct cellular identities. In the Hematopoiesis dataset, we observed that the Multipotent Progenitor Cells (MPPs), which are downstream progenitor cells derived from Hematopoietic Stem Cells (HSCs), exhibit close proximity to each other in the latent space of ATAC-Diff. This finding aligns with their known differentiation relationship, further supporting the accuracy and biological relevance of our analysis.In this section, we have devised two strategies to evaluate the generation capabilities of ATAC-Diff. The first strategy involves unconditional generation, where we assess whether the scATAC-seq data generated by ATAC-Diff exhibits realistic characteristics. If our model excels in capturing the data distribution, it can be leveraged to create simulated data for augmentation purposes rather than sequencing additional cells, thus conserving time and resources. Furthermore, visualizing the synthetic data generated can enhance researchers' comprehension of the distribution patterns and fundamental structure of scATAC-seq data. In this case, we anticipate that ATAC-Diff will be able to generate cells that possess the specific subpopulation features associated with the provided conditional information. Since we are the first to generate scATAC-seq data from scratch, there are no baseline models for comparison. Therefore, we have made modifications to SCALE and PeakVI, which are VAE-based models, to enable unconditional generation by utilizing random noise as inputs.For conditional generation, we have replaced the VAE module of SCALE with a conditional VAE module 4 . For the unconditional generation, we average all single cells as the ground truth. For the conditional generation, we average all single cells of the same biological cell type as the ground truth.To address the requirement of ATAC-Diff for latent variables as auxiliary information, we employ a vanilla VAE to sample the latent variables. We evaluate the generation performance based on the SCC and PCC between the mean values of generated samples and the two ground truth data.We report the results of unconditional generation in Table 6 . We observe that ATAC-Diff outperforms other methods in generating scATAC-seq data from scratch, yielding more realistic results. Specifically, ATAC-Diff achieves the highest SCC of 0.892 and PCC of 0.969 on Forebrain dataset, SCC of 0.822 and PCC of 0.973 on Hematopoiesis dataset, and PCC of 0.983 on PBMC10k dataset. Table 6 also displays the mean correlation values of each cell type in the context of conditional generation. We observe that ATAC-Diff outperforms conditional SCALE among all the metrics. In addition, the performance of ATAC-Diff without conditional information degraded, indicating that ATAC-Diff can fuse the conditional cell type information to generate high-quality scATAC-seq data of specific cell types. Overall, we prove that ATAC-Diff is capable of generating realistic scATAC-seq data. The scATAC-seq data always contains both noised and a large number of missing values due to dropout events [30] . We design two approaches to test the ability of the method for denoising and recovering missing values (imputation) to address the issues under real scenarios. For scATACseq data denoising, analyzing real datasets poses a challenge due to the lack of ground truth data. Following previous work [8] , we mitigate this challenge to average all single cells within the same cell type, resulting in a meta-cell that serves as a good approximation of those individual cells.For scATAC-seq data imputation, previous works chose to create a corrupted matrix by randomly dropping out some non-zero entries. In our case, we follow the setting used in DCA [31] to simulate dropout events by masking 10% of the non-zero counts and setting them to zero. The probability of masking a particular non-zero count follows an exponential distribution which peaks with lower expression levels have a higher likelihood of dropout compared to genes with higher expression levels [32] . For comparison, we choose SCALE and peakVI as baseline models and SCC and PCC as metrics.The results of scATAC-seq data denoising and imputation of ATAC-Diff and baseline models are presented in Table 7 where the best result is highlighted in bold. For the data denoising evaluation, ATAC-Diff achieves the highest correlation of the denoised single cells with the corresponding meta-cells of each cell type on PBMC10k dataset and competitive results on the other two datasets.Although the other generative models PeakVI and SCALE perform well on some of the datasets, their performances are not robust. One potential reason is that the auxiliary latent variables enable ATAC-Diff to capture the high-level semantics information. Similarly, ATAC-Diff performs stable on scATAC-seq data imputation, where the generated data is closely related to the meta-cells. 
conclusion: 
 In this work, we propose a versatile framework ATAC-Diff which is based on a diffusion model conditioned on the latent variables for scATAC-seq data analysis and generation. ATAC-Diff could utilize the proposed auxiliary module to yield the latent variables which contain high-level meaningful information; it could be incorporated to assist the data sampling during diffusion. Additionally, these variables could be examined for downstream analyses such as cell type annotation. Based on the formulation of the mutual information, we prevent the model from ignoring the low-dimensional latent variables. Furthermore, we employ a Gaussian Mixture Model (GMM) as the latent prior and auxiliary decoder to enhance the recovery of real data, allowing the latent variables to learn the genomic semantics.We conducted extensive experiments to evaluate the performance of ATAC-Diff on three datasets and compared it with SOTA models, indicating that the ATAC-Diff model can generate high-quality data in silico scATAC-seq data while effectively disentangling the cell embeddings. An intriguing potential direction for ATAC-Diff involves exploring various conditional generation scenarios for genomic discovery, such as identifying the connections between scRNA-seq data and perturbation prediction. We anticipate that ATAC-Diff can contribute to the advancement of genomic analysis in the near future of personalized medicine at the single-cell level.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [No] Justification: We will release the full source code soon. The link is https://github.com/Layne-Huang/ATAC-Diff.Guidelines:The answer NA means that paper does not include experiments requiring code. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
acknowledgment: 
 The  
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [Yes] Justification: This paper poses no such risk.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes]Justification: [TODO]Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.• If this information is not available online, the authors are encouraged to reach out to the asset's creators. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: We provide proofs for the derivation of several properties in ATAC-Diff.In the diffusion process, we have the marginal distribution of the data at any arbitrary time step t in a closed form:q (x t | x 0 ) = N (x t ; √ ᾱt x 0 , (1 -ᾱt ) I) .Proof: recall the posterior q (x t | x t-1 ) in Eq. 3, we can obtain x t using the reparameterization trick. A property of the Gaussian distribution is that if we add N (0, σ 2 1 I) and N (0, σ 2 2 I), the new distribution is N (0, (σ 2 1 + σ 2 2 )I)EQUATIONwhere α t = 1 -β t , ϵ and ε are sampled from independent standard Gaussian distributions. Then we could derive Eq. 3A learned Gaussian transitions p θ (x t-1 | x t ) is devised to approximate the q (x t-1 | x t ) of every time step:p θ (x t-1 | x t ) = N x t-1 ; µ θ (x t , t) , σ 2 t I . µ θ is parameterized as follows: µ θ (x t , t) = 1 √ α t x t - β t √ 1 -ᾱt ϵ θ (x t , t) .Proof: the distribution q (x t-1 | x t ) can be expanded by Bayes' rule:EQUATIONEQUATIONwhere C (x t , x 0 ) is a constant. We can find that q (x t-1 | x t ) is also a Gaussian distribution. We assume that:EQUATIONwhere βt = 1/ αt βt +1 1-ᾱt-1 = 1-ᾱt-1 1-ᾱt • β t and μt (x t , x 0 ) = √ αt βt x t + √ ᾱt-1 1-ᾱt-1 x 0 / αt βt + 1 1-ᾱt-1 = √ αt(1-ᾱt-1) 1-ᾱt x t + √ ᾱt-1βt 1-ᾱt x 0 .Then we could parameterize µ θ (x t , t)= √ αt(1-ᾱt-1) 1-ᾱt x t + √ ᾱt-1βt1-ᾱt x θ , which is presentned in Eq. 18. From Eq. 19, we have x t == √ ᾱt x 0 + √ 1 -ᾱt ε. We take this into μ:μt = √ α t (1 -ᾱt-1 ) 1 -ᾱt x t + √ ᾱt-1 β t 1 -ᾱt 1 √ ᾱt x t - √ 1 -ᾱt ϵ t = 1 √ α t x t - β t √ 1 -ᾱt ϵ tµ θ is designed to model μ. Therefore, µ θ has the same formulation as μ but parameterizes ϵ:µ θ (x t , t) = 1 √ αt x t -βt √ 1-ᾱt ϵ θ (x t , t) .It is hard to directly calculate the conditional log-likelihood of the data. Instead, we can derive its ELBO objective for optimizing. For simplicity, the c in GMM is omittedE [log p θ (x 0 )] = E q(x0) log E q(z) p θ (x 0:T |z) dx 1:T = E q(x0) log p θ (x 0:T , z) dx 1:T dz = E q(x0) log q (x 1:T | x 0 ) q ϕ (z|x 0 ) p θ (x 0:T , z) q (x 1:T | x 0 ) q ϕ (z|x 0 ) dx 1:T dz = E q(x0) log E q(x 1:T |x0) E q ϕ (z|x0) p θ (x 0:T , z) q (x 1:T | x 0 ))q ϕ (z|x 0 ) ≥ E q(x 0:T ) log E q ϕ (z|x0) p θ (x 0:T , z) q (x 1:T | x 0 )q ϕ (z|x 0 )) (22)Then we further derive the conditional ELBO objective:EQUATIONFor simplicity, we utilize the last term log p θ (z) q ϕ (z|x0) and the regular term of GMM and mutual information to obtain the final form of the optimizing objectivity in Eq. 7 (We also omit c in GMM for simplicity).E q log p θ (z) q ϕ (z|x 0 ) + αI(x 0 ; z) + λR GMM = E q log p θ (z) q ϕ (z|x 0 ) + α log q ϕ (z) q ϕ (z|x 0 ) + λ log q ϕ (z|x 0 ) p(z) = E q log q ϕ (z|x 0 ) λ-α-1 q ϕ (z) α p(z) λ-1 = E q log q ϕ (z|x 0 ) λ-α-1 q ϕ (z) α p(z) λ-α-1 p(z) α = E q log q ϕ (z|x 0 ) λ-α-1 p(z) λ-α-1 + log q ϕ (z) α p(z) α = (λ -α -1)D KL (q ϕ (z|x 0 )||p(z)) + αD KL (q ϕ (z)||p(z)) (24)Finally, we could derive the objectivity as:EQUATIONWe summarize the statistic information of all processed datasets in Table 4 4 : The statistic summary of datasets For all three datastes, ATAC-Diff is trained by AdamW [33] optimizer for 10K iterations with a batch size of 256 and a learning rate of 0.0001 on one NVIDIA V100 GPU card. We set α as 0.1 and λ as 1.1001.In this study, we have removed the GMM and MI modules to conduct the ablation study for each task. We present a similarity matrix of the latent embedding based on Euclidean distances. Specifically, we average the latent embeddings within each cell type population and calculate the Euclidean distances across cell types. Our results show that ATAC-Diff can effectively capture the relationships between distinct cell subpopulations and developmental trajectories. For example, in the Forebrain dataset, the three excitatory neuron subtypes (EC1, EC2, EC3), which are biologically similar, cluster closely in ATAC-Diff's latent space based on the correlation of their embeddings. In contrast, more biologically distant cell types display greater separation. Similarly, in the Hematopoiesis dataset, Multipotent Progenitor Cells (MPPs) and Hematopoietic Stem Cells (HSCs) exhibit proximity, reflecting their known differentiation path and reinforcing the biological relevance of our findings. We have visualized the extracted features by different methods on PBMC10k dataset through UMAP, which is shown in Figure 3 HVP PCA peakVI cisTopic SCALE ATAC-DiffFigure 3 : UMAP visualization of the highly variable peak values and extracted features from PCA, cisTopic, SCALE, PeakVI, and ATAC-Diff on PBMC10k dataset.BackgroundThe diffusion model is formulated as two Markov chains: diffusion process and reverse process. Given a variance preserve schedule, the diffusion process transforms the real data x 0 to a latent variable distribution p(x 0:T ) and finally the data is diffused into predefined Gaussian noise x TDownloaded PBMC10k dataset from https://support.10xgenomics.com/single-cell-multiome-atac-gex/datasets/1.0.0/pbmc_ granulocyte_sorted_10k, generated from a healthy donor after removing granulocytes through cell sortingUnfortunately, due to the high integration of PeakVI into the scvi-tools package, we were unable to make modifications to it. 


Decision: The reviewers name the theoretical motivation and analysis of the method as well as good empirical performance as strengths. Writing is also generally speaking good, although some unclear parts were pointed out and should be clarified. Among weaknesses, some relevant comparisons were missing from the original manuscript and should be included for a more comprehensive evaluation with SOTA. Moreover, there are concerns about the evaluation metric (due to possibly misleading Euclidean distances in UMAP space).

Overall, it seems that the authors were able to address the most significant critical comments in their rebuttal (with the exception of the UMAP problem which remains an issue but, in my opinion, not a fatal one), and no further concerns were raised by the reviewers.

Reviews: review: paper_summary: paper_summary: Generating simulated scATAC-seq data is important for developing new methods and gaining a deeper understanding of the data. However, the simulation is challenging due to dropout and high noise in the data. Authors proposed a diffusion + VAE type of method to solve the problem. The general idea is first to use VAE to project original scATAC data to a lower embedding space. Then impose a diffusion process in the embedding space. The lower embedding space is a GMM rather than the classical isotropic Normal, which is the novelty part of the method. This configuration makes biological sense as cells can be grouped as different cell types. Due to the introduction of GMM distribution in the hidden GMM space, there exists complications in generalising the diffusion loss function. The authors have shown nice and solid derivations in the appendix. The method is then applied to three datasets on three different tasks and achieved comparable  performance as SOTA methods. It seems that authors have provided a convincing solution to the research question. While this paper is clearly written and the general idea is relatively easy to follow, it will be nice if the author could help to answer the following questions.
1. In Eq. 11, what is the parametric form of q_\phi(z|x_0)
2. In Eq. 14, what is the actual meaning of conditional information y, could you list the parameters
3. Could you provide a concrete network architecture of your network in a supplementary figure, i.e. including the tensors and their dimensions?
4. In sec 4.4.1, what is the dropout rate distribution in your simulated data, are they similar to the real data?
Depending on the answers, I may change my ratings in the future.
strength_weakness: strength_weakness: Generating simulated scATAC-seq data is important for developing new methods and gaining a deeper understanding of the data. However, the simulation is challenging due to dropout and high noise in the data. Authors proposed a diffusion + VAE type of method to solve the problem. The general idea is first to use VAE to project original scATAC data to a lower embedding space. Then impose a diffusion process in the embedding space. The lower embedding space is a GMM rather than the classical isotropic Normal, which is the novelty part of the method. This configuration makes biological sense as cells can be grouped as different cell types. Due to the introduction of GMM distribution in the hidden GMM space, there exists complications in generalising the diffusion loss function. The authors have shown nice and solid derivations in the appendix. The method is then applied to three datasets on three different tasks and achieved comparable  performance as SOTA methods. It seems that authors have provided a convincing solution to the research question. While this paper is clearly written and the general idea is relatively easy to follow, it will be nice if the author could help to answer the following questions.
1. In Eq. 11, what is the parametric form of q_\phi(z|x_0)
2. In Eq. 14, what is the actual meaning of conditional information y, could you list the parameters
3. Could you provide a concrete network architecture of your network in a supplementary figure, i.e. including the tensors and their dimensions?
4. In sec 4.4.1, what is the dropout rate distribution in your simulated data, are they similar to the real data?
Depending on the answers, I may change my ratings in the future.
questions: questions: Generating simulated scATAC-seq data is important for developing new methods and gaining a deeper understanding of the data. However, the simulation is challenging due to dropout and high noise in the data. Authors proposed a diffusion + VAE type of method to solve the problem. The general idea is first to use VAE to project original scATAC data to a lower embedding space. Then impose a diffusion process in the embedding space. The lower embedding space is a GMM rather than the classical isotropic Normal, which is the novelty part of the method. This configuration makes biological sense as cells can be grouped as different cell types. Due to the introduction of GMM distribution in the hidden GMM space, there exists complications in generalising the diffusion loss function. The authors have shown nice and solid derivations in the appendix. The method is then applied to three datasets on three different tasks and achieved comparable  performance as SOTA methods. It seems that authors have provided a convincing solution to the research question. While this paper is clearly written and the general idea is relatively easy to follow, it will be nice if the author could help to answer the following questions.
1. In Eq. 11, what is the parametric form of q_\phi(z|x_0)
2. In Eq. 14, what is the actual meaning of conditional information y, could you list the parameters
3. Could you provide a concrete network architecture of your network in a supplementary figure, i.e. including the tensors and their dimensions?
4. In sec 4.4.1, what is the dropout rate distribution in your simulated data, are they similar to the real data?
Depending on the answers, I may change my ratings in the future.
limitations: limitations: Generating simulated scATAC-seq data is important for developing new methods and gaining a deeper understanding of the data. However, the simulation is challenging due to dropout and high noise in the data. Authors proposed a diffusion + VAE type of method to solve the problem. The general idea is first to use VAE to project original scATAC data to a lower embedding space. Then impose a diffusion process in the embedding space. The lower embedding space is a GMM rather than the classical isotropic Normal, which is the novelty part of the method. This configuration makes biological sense as cells can be grouped as different cell types. Due to the introduction of GMM distribution in the hidden GMM space, there exists complications in generalising the diffusion loss function. The authors have shown nice and solid derivations in the appendix. The method is then applied to three datasets on three different tasks and achieved comparable  performance as SOTA methods. It seems that authors have provided a convincing solution to the research question. While this paper is clearly written and the general idea is relatively easy to follow, it will be nice if the author could help to answer the following questions.
1. In Eq. 11, what is the parametric form of q_\phi(z|x_0)
2. In Eq. 14, what is the actual meaning of conditional information y, could you list the parameters
3. Could you provide a concrete network architecture of your network in a supplementary figure, i.e. including the tensors and their dimensions?
4. In sec 4.4.1, what is the dropout rate distribution in your simulated data, are they similar to the real data?
Depending on the answers, I may change my ratings in the future.
paper_summary: paper_summary: This submission introduces ATAC-Diff, a conditional latent diffusion model for scATAC-seq data. ATAC-Diff incorporates a few components to a basic diffusion model, including a Gaussian Mixture Model as a semantic prior over the latent variables and mutual information as a regularizer. ATAC-Diff is benchmarked on several tasks: latent representation clustering, (contitional) data generation, and denoising/imputation. ATAC-Diff is shown to do better than PeakVI and SCALE (among other methods for clustering).


[UPDATE 8/9]: Adjusting score from Reject to Accept due to clarifications.
strength_weakness: strength_weakness: - presentation of ATAC-seq data is not clear. The cell x peak representation is a processed version that doesn't necessarily encapsulate the full ATAC-seq data. Clarity that they are not generating ATAC-seq data but rather processed ATAC-seq data in the form of cell x binary peaks is necessary.
- Example: other approaches to analyze ATAC-seq data incorporate sequences, eg. chromBPnet and AI-TAC.
- Weak evaluations. The evaluations of cell clustering are not clear. It is not clear how ground truth was determined, given it is likely given by another computational method. The metrics are not clear what they are and their strengths and limitations. The umap visualization is not reliable. The authors state that they observe EC1, EC2, EC3 from the Forebrain dataset are close proximity in latent space. But latent space can be highly warped, rendering any Euclidean distances in umap space not meaningful. 
- Generation quality task is confusing. 
	- When is unconditional generation needed? Why is this a meaningful task?
	- Not clear how the metrics are calculated. What is the ground truth here?
	- Why are generated cells averaged? Shouldn't generation be assessed at the single-cell level? 
	- The ATAC-Diff w.o con seems quite high. This affirms my concern that the evaluations might be weak. 
- The proposed components GMM and MI regularizer are not evaluated via an ablation study. Are they even needed?
- This study treats ATAC-seq data as a cell x (binary) peaks representation, instead of fragment counts, which improves scATAC-seq analysis (see Martens et al, Nature Methods, 2023).
questions: questions: Questions are integrated within Weaknesses (above).
limitations: limitations: This study seems interesting and their method ATAC-Diff may be an advance, but it is difficult to assess given the poor presentation. There are several weaknesses such as questionable ground truth within the benchmarks. Moreover, the task itself seems to be dated, focusing on cells x binary peaks. Moreover, the components of the proposed approach are not tested in any ablation study, bringing their purpose questionable. This approach seems to be a new method, but one of many in the growing universe of scATAC-seq data and it's not clear what impact this work will have due to these limitations.
paper_summary: paper_summary: The paper introduces a new diffusion model ATAC-diff for scATAC-seq data generation and analysis. ATAC-diff using a latent diffusion model which is conditioned on latent auxiliary models to encode latent variables, and integrate GMM as the latent prior to capture genetic information. This paper introduces a mutual information regularization to maintain the connection between observed and latent variables. Extensive experiments show that ATAC-diff outperforms state-of-the-art models in both data generation and analysis tasks.
strength_weakness: strength_weakness: Multiple baseline performance on the benchmark dataset does not match what was reported in previous publications. For example, [1] shows much higher PCA performance on PBMC10k dataset. This discrepancy makes me concern about the reliability of the experiment results. 
If would be great to include a few benchmark dataset to help cross validate the model performance, for example Buenrostro2018. 
There are top baselines that are not included in the experiment, please consider add the comparison on the latest baseline as well, such as scBasset. 
[1] Yuan H, Kelley D R. scBasset: sequence-based modeling of single-cell ATAC-seq using convolutional neural networks[J]. Nature Methods, 2022, 19(9): 1088-1096.
questions: questions: 1. Explain the discrepancy on baseline performance and cross validation with previous publications. 
2. Include more benchmark dataset for better cross validation of the model performance, such as Buenrostro2018.
3. Add top baselines, such as scBasset into comparison.
limitations: limitations: yes

Comments: title: Inquiry about the feedback from the reviewers
satisfactory response.
comment: Dear Reviewer,

We would like to express our gratitude for your valuable comments again. Please do not hesitate to reach out if you have any further concerns or require additional information.

Thanks for your attention. We are looking for your reply!
Thank you very much for the comments, and we deeply appreciate your valuable suggestions. Those comments are of great value for improving the quality of the manuscript.
Thank you for taking the time to review our response for your comments. We would like to express our sincere gratitude for your thorough review and valuable feedback throughout the entire review process. Your insights and suggestions have been instrumental in helping us improve the quality and clarity of our work.
I have read the rebuttal and it is ok to accept the paper.
Dear Authors,
You are right. Most likely your rebuttal is sufficient to resolve any unclarities at this point, but if there are any further questions, they'll hopefully be posted asap so that you may still react.
thanks for your patience,
Area Chair
Dear ACs,

I am sorry to bother you. Could you kindly assist in verifying if the reviewers have received the rebuttal information? It has been six days since the rebuttal period concluded, and we only have two days left to discuss my response to their comments and their feedback on this manuscript. 

Your prompt assistance is greatly appreciated. I am looking forward to your response. Thank you for your help.
The authors have addressed my concerns.  I still don't agree with UMAP analysis as it can also warp local spaces. I will adjust my scores accordingly.