
==== Front
Brief Bioinform
Brief Bioinform
bib
Briefings in Bioinformatics
1467-5463
1477-4054
Oxford University Press

10.1093/bib/bbae093
bbae093
Problem Solving Protocol
AcademicSubjects/SCI01060
Incorporating network diffusion and peak location information for better single-cell ATAC-seq data analysis
Yu Jiating School of Mathematics and Statistics, Nanjing University of Information Science & Technology, Nanjing 210044, China
IAM, MADIS, NCMIS, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China
School of Mathematical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China

https://orcid.org/0000-0002-0816-9690
Leng Jiacheng IAM, MADIS, NCMIS, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China
School of Mathematical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China
Zhejiang Lab, Hangzhou 311121, China

Hou Zhichao IAM, MADIS, NCMIS, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China
School of Mathematical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China

https://orcid.org/0000-0002-2802-6347
Sun Duanchen School of Mathematics, Shandong University, Jinan 250100, China

https://orcid.org/0000-0001-9487-0215
Wu Ling-Yun IAM, MADIS, NCMIS, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China
School of Mathematical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China

Corresponding authors: Ling-Yun Wu, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China. E-mail: lywu@amss.ac.cn; Duanchen Sun, School of Mathematics, Shandong University, Jinan 250100, China. E-mail: dcsun@sdu.edu.cn
3 2024
16 3 2024
16 3 2024
25 2 bbae09321 9 2023
22 12 2023
20 2 2024
© The Author(s) 2024. Published by Oxford University Press.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com

Abstract

Single-cell assay for transposase-accessible chromatin using sequencing (scATAC-seq) data provided new insights into the understanding of epigenetic heterogeneity and transcriptional regulation. With the increasing abundance of dataset resources, there is an urgent need to extract more useful information through high-quality data analysis methods specifically designed for scATAC-seq. However, analyzing scATAC-seq data poses challenges due to its near binarization, high sparsity and ultra-high dimensionality properties. Here, we proposed a novel network diffusion–based computational method to comprehensively analyze scATAC-seq data, named Single-Cell ATAC-seq Analysis via Network Refinement with Peaks Location Information (SCARP). SCARP formulates the Network Refinement diffusion method under the graph theory framework to aggregate information from different network orders, effectively compensating for missing signals in the scATAC-seq data. By incorporating distance information between adjacent peaks on the genome, SCARP also contributes to depicting the co-accessibility of peaks. These two innovations empower SCARP to obtain lower-dimensional representations for both cells and peaks more effectively. We have demonstrated through sufficient experiments that SCARP facilitated superior analyses of scATAC-seq data. Specifically, SCARP exhibited outstanding cell clustering performance, enabling better elucidation of cell heterogeneity and the discovery of new biologically significant cell subpopulations. Additionally, SCARP was also instrumental in portraying co-accessibility relationships of accessible regions and providing new insight into transcriptional regulation. Consequently, SCARP identified genes that were involved in key Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways related to diseases and predicted reliable cis-regulatory interactions. To sum up, our studies suggested that SCARP is a promising tool to comprehensively analyze the scATAC-seq data.

scATAC-seq
network diffusion
genomic distance
accessible regions
cellular heterogeneity
National Key Research and Development Program of China 10.13039/501100012166 2022YFA1004803 2020YFA0712402 National Natural Science Foundation of China 10.13039/501100001809 12231018 62202269
==== Body
pmcINTRODUCTION

Chromatin accessibility is closely linked to the occurrence of transcriptional regulation, as accessible chromosome fragments often contain transcription factor (TF) binding sites and other important cis-regulatory elements, such as enhancers and promoters [1, 2]. While single-cell transcriptome data can help us to construct gene regulatory networks from the level of transcripts and identify cellular diversity without bias, it is still challenging to reverse-engineer the mechanism of transcriptional regulation [2]. Fortunately, the development of single-cell epigenomic technology has brought new perspectives to solve this problem, such as single-cell assay for transposase-accessible chromatin using sequencing (scATAC-seq) [3] and single-cell combinatorial indexing assay for transposase-accessible chromatin with sequencing (sci-ATAC-seq) [4, 5], which can build chromatin accessibility profiles with single-cell resolution [6]. With the aid of scATAC-seq data, we were able to model the co-accessibility of chromosome segments from the genome level directly, which provides the basic conditions for transcriptional regulation to occur.

Compared with single-cell transcriptome data, scATAC-seq data exhibit more sparsity, higher dimension and near-binarization properties, which pose many challenges to scATAC-based analysis. Because for a diploid genome, only one or two reads can be captured at each accessible chromatin site [6], resulting in the lack of sequencing data, i.e. the curse of ‘missingness’ in scATAC-seq [1, 6]. In addition, accessible regions are usually annotated on the entire genome, so the features of scATAC-seq data can reach more than 100k dimensions or even higher. These properties make scRNA-seq-based algorithms less effective for analyzing scATAC-seq [7], although these algorithms are mathematically transferable since both genes and peaks can be understood as features to characterize cells.

Currently, various methods have been specially developed for scATAC-seq data to decipher their underlying information. For example, chromVAR [8] measures the accessibility gain or loss among peaks with the same motif or annotation and uses a bias-corrected deviation matrix for downstream analysis. However, the aggregation of similar peaks loses much information, making it less effective for cell clustering. SCALE [1] incorporates the variational autoencoder and Gaussian mixture model to extract latent features of cells. Nevertheless, it has numerous parameters and needs to be carefully tuned to get good performance [6]. Besides, cisTopic [2] uses the Latent Dirichlet Allocation model for co-optimal clustering of the cells and accessible regions, yet the collapsed Gibbs sampler makes it computationally slow. A recently proposed method, scAND [6], adopts the Katz index [9] diffusion method to alleviate the high sparsity in scATAC-seq data, and using the smoothed matrix as model input can achieve the best performance in their following analyses. However, scAND does not consider the differences between different chromosomes in the diffusion process and neglects the location associations of peaks.

In genetics studies such as recombination between two loci and linkage disequilibrium [10–13], the genetic map distance has been considered an essential factor. Genes on different chromosomes or distantly separated on the same chromosomes are assorted independently and described as physically unlinked [12], and the plausibility of modeling chromatin contacts as a decreasing function of genomic distances has been explained and verified [14]. Therefore, as scATAC-seq data give the accessibility of genome-wide segments in cells, incorporating genomic distance information of accessible fragments may be beneficial for the downstream data analyses, such as cell subpopulation identification and cis-regulatory network construction, which has not been considered in the state-of-the-art methods.

In addition, it is critical to exploit a superior network diffusion method (NDM) that can appropriately aggregate the neighborhood information to deal with the sparseness and missingness problems in scATAC-seq data. Our previously developed Network Refinement (NR) [15] diffusion method can improve the quality of noisy networks by aggregating path information of various lengths. NR employs a global diffusion mathematical operator, as defined by the random walk on a graph, to enhance the self-organization properties of complex networks. We have demonstrated that NR is a degree-normalized version of the Katz index and has exhibited superior performance in various complex network-related tasks, including network denoising [15] and link prediction [16]. We, therefore, expect to see that NR also handles scATAC-seq data well, with the rationale being that (1) the diffusion of accessibility relationships can compensate for missing information and depict similarities between nodes (i.e. cells and peaks) of the network and (2) peaks accessible in too many cells should reduce their impact on diffusion, as their ability to portray cell similarity is overestimated.

In this study, we proposed a novel network-based method to comprehensively analyze scATAC-seq data, named Single-Cell ATAC-seq Analysis via Network Refinement with Peaks Location Information (SCARP). Specifically, SCARP takes genetic location information of peaks into account and globally diffuses the accessibility relationships using our previously developed NR diffusion method. The output matrix derived from SCARP can be further processed by the dimension reduction method to obtain low-dimensional embeddings of cells and peaks, which can benefit the downstream analyses such as the cells clustering and cis-regulatory relationships prediction (Figure 1).

Figure 1 The workflow of SCARP. SCARP consists of two modules. The first module involves constructing the network from scATAC-seq data, with the nodes being cells or peaks, the edges between cells and peaks reflecting accessibility, and the edges between peaks reflecting prior information about their genome position. The second module employs the NR diffusion method to compensate for the high sparsity of scATAC-seq data as well as characterize the similarity between cells. The output matrix of SCARP with the low-dimensional embedding of cells and peaks can further benefit the downstream analyses, such as cell embedding visualization, cell subgroup identification, and cis-regulatory network construction. (A) Binarizing scATAC-seq data and building a bipartite network to model accessibility relationships between cells and peaks. (B) Calculating prior edge weights of adjacent peaks based on their genome distance. (C) Incorporating prior edges into the bipartite network. (D) Applying the NR diffusion method to subnetworks corresponding to different chromosomes individually. (E) Merging diffusion matrices obtained from different chromosomes and integrating cell–cell similarities. (F) Performing dimensionality reduction to obtain cell and peak embeddings.

We have demonstrated through sufficient experiments that SCARP facilitates superior analysis of scATAC-seq data, including improving cell clustering performance to better elucidate cell heterogeneity and reveal new biologically significant cell subpopulations and characterizing co-accessibility of genome segments to shed fresh light on transcriptional regulation. Our studies suggested that SCARP is a promising tool to comprehensively analyze the scATAC-seq data from a new perspective.

MATERIALS AND METHODS

Overview of SCARP

We proposed a novel network diffusion–based method, SCARP, to comprehensively analyze the scATAC-seq data. The workflow of SCARP is shown in Figure 1. First, SCARP constructed a bipartite network to model the accessible relationships between the cells and peaks (Figure 1A). Second, prior edge weights of adjacent peaks based on their genome distance were employed to better capture the co-accessibility of adjacent peaks (Figure 1B). Together, they constituted the input network for the next step (Figure 1C). Third, the NR diffusion method [15] was performed on the subnetworks corresponding to different chromosomes separately to obtain dense matrixes that can reflect accessibility similarities between the cells and peaks, as well as possible cell–peak accessible relationships (Figure 1D). After this, diffusion matrixes obtained from different chromosomes were spliced together and the cell–cell similarities were integrated (Figure 1E). Finally, dimensionality reduction technique was performed to map cells and peaks into low-dimensional space, and these representations were used for downstream analyses, such as cells clustering and construction of cis-regulatory networks (Figure 1F).

SCARP workflow

Given scATAC-seq data matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${A}_{c\times p}$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $c$\end{document} cells and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $p$\end{document} peaks, we first binarized \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $A$\end{document} to better cater to its biological interpretation of depicting accessibility. Specifically, we set all values greater than 1 to 1, making \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $A$\end{document} a Boolean matrix with entries from the Boolean domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\left\{0,1\right\}$\end{document}, reflecting whether a cell is accessible on a certain genome region. Then, a symmetric square root graph normalization \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\overset{\sim }{A}={D}_L^{-\frac{1}{2}}A{D}_R^{-\frac{1}{2}}$\end{document} was applied to eliminate the library size differences between cells and peaks, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${D}_L$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${D}_R$\end{document} are diagonal degree matrixes with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\left({D}_L\right)}_{ii}=\sum_j{A}_{ij}$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\left({D}_R\right)}_{jj}=\sum_i{A}_{ij}$\end{document}.

After the pre-processing step, SCARP constructed a bipartite network with an adjacency matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${B}_{N\times N}$\end{document} to model the accessible relationship between the cells and peaks, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $N=c+p$\end{document} is the number of nodes of the network, i.e. the peaks and cells are all treated as nodes of the network. To incorporate genetic location information for a better depiction of the co-accessibility of adjacent peaks, we computed the prior weights between adjacent peaks using the negative Haldane’s genetic map function. We can now represent the network as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ {B}_{N\times N}=\left[\begin{array}{@{}cc@{}}0& \overset{\sim }{A}\\{}{\overset{\sim }{A}}^T& H\end{array}\right] $$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $H$\end{document} is a sub-diagonal matrix of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $p\times p$\end{document}, whose element represents the prior weight between the adjacent peaks in the same chromosome.

The NR diffusion method was then applied to compute accessibility similarities between the cells and peaks, as well as predict possible cell–peak accessible relationships. Considering that those peaks on the same chromosome can better learn from each other, the NR diffusion process was calculated on the subnetwork of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $B$\end{document}, which is induced from nodes of cells and those peaks on the same chromosome and then spliced back together. Specifically, assuming there are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $M$\end{document} chromosomes, the graph \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $B$\end{document} was divided into \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $M$\end{document}subgraphs \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${B}_{N_1\times{N}_1}^{(1)},\dots, {B}_{N_M\times{N}_M}^{(M)}$\end{document}, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${N}_k=c+{p}_k$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${p}_k$\end{document} is the number of peaks in the subgraph \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${B}^{(k)}$\end{document}, and NR was performed on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${B}^{(k)}$\end{document} to get a diffused matrix. This diffusion process can be performed in parallel, which reduces the time cost. Notably, when splitting peaks, we divided them according to whether they were on the same chromosome or not. Given the sparsity defect of the scATAC-seq data, we did not recommend further subdividing one chromosome into smaller pieces, as this may lead to insufficient information and potential disruption of chromosomal information chains (Supplementary Figure 1). When there were fewer peaks, we merged the adjacent chromosomes and set the threshold parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\gamma$\end{document} (default as 3000) to determine whether we need to merge the adjacent chromosomes. Specifically, when the number of peaks in a chromosome is less than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\gamma$\end{document}, we combine them with the peaks in the next chromosome to form a group, which is then subjected to subsequent diffusion analysis steps. Conversely, if the number of peaks in a chromosome is more than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\gamma$\end{document}, they are treated as a separate group.

Notice that after parallel diffusion, we got \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $M$\end{document} dense matrixes, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $${S}^{(k)}=\left[\begin{array}{@{}cc@{}}{S}_{c\times c}^{(k)(11)}& {S}_{c\times{p}_k}^{(k)(12)}\\{}{S}_{p_k\times c}^{(k)(21)}& {S}_{p_k\times{p}_k}^{(k)(22)}\end{array}\right]\left(k=1,\dots, M\right)$$\end{document}, each of which has a submatrix of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $c\times c$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${S}_{c\times c}^{(k)(11)}$\end{document}, representing the cell similarity calculated based on the peaks of the current subgraph, and when splicing them back into matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${D}_{N\times N}$\end{document},

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ {D}_{N\times N}=\left[\begin{array}{@{}cc@{}}{D}_{c\times c}^{(11)}& {D}_{c\times p}^{(12)}\\{}{D}_{p\times c}^{(21)}& {D}_{p\times p}^{(22)}\end{array}\right]=\left[\begin{array}{@{}cccc@{}}{D}_{c\times c}^{(11)}& {S}_{c\times{p}_1}^{(1)(12)}& \cdots & {S}_{c\times{p}_M}^{(M)(12)}\\{}{S}_{p_1\times c}^{(1)(21)}& {S}_{p_1\times{p}_1}^{(1)(22)}&\ & \kern0.5em \\{}\vdots &\ & \ddots &\ \\{}{S}_{p_M\times c}^{(M)(21)}& \kern0.5em &\ & {S}_{p_M\times{p}_M}^{(M)(22)}\end{array}\right] $$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${D}_{c\times c}^{(11)}$\end{document} was the average over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $M$\end{document} matrixes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${S}_{c\times c}^{(k)(11)}$\end{document} (Figure 1).

Finally, to get cell and peak embeddings, the dimensional reduction method, such as UCPCA, was performed on the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\left[{D}_{c\times c}^{(11)},\kern0.5em {D}_{c\times p}^{(12)}\right]$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${D}_{p\times c}^{(21)}$\end{document} separately. The number of kept components was determined based on the proportion of variance that can be explained, as detailed below.

Prior weights between adjacent peaks

The prior edge weights between the adjacent peaks in this study were calculated using the negative Haldane’s genetic map function, motivated by the traditional use of gene map function [10] to estimate the relation connecting recombination fractions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $r$\end{document} and genomic distances \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $d$\end{document}:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ r={e}^{-d}+{e}^{-d}\frac{d^3}{3!}+\dots =\frac{1}{2}\left(1-{e}^{-2d}\right) $$\end{document}

We assumed that the probability of chromosome segments being co-accessible should be a decreasing function of its genomic distance, as opposed to the probability of recombination [17, 18]. Specifically, the distance between two peaks \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\mathrm{peak}}_1=\left[{s}_1,{e}_1\right]$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\mathrm{peak}}_2=\left[{s}_2,{e}_2\right]$\end{document} (any two peaks do not intersect with each other) was \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $d\left({\mathrm{peak}}_1,{\mathrm{peak}}_2\right)=\max \left({s}_2-{e}_1,{s}_1-{e}_2\right)$\end{document}, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${s}_i$\end{document} stands for the start position and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${e}_i$\end{document} stands for the end position on the genome, and the prior weight of the two peaks was computed as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ w\left({\mathrm{peak}}_1,{\mathrm{peak}}_2\right)={e}^{-\frac{2d\left({\mathrm{peak}}_1,\kern0.5em {\mathrm{peak}}_2\right)}{\beta }} $$\end{document}

Here, we set a parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\beta$\end{document} (default as 5000) to control the extent to which prior edge weight decays as peak distance increases.

NR diffusion method

We used our previously proposed method, NR [15], to cope with the sparsity of scATAC-seq data. Specifically, NR takes the adjacency matrix of a network as input, uses a diffusion process defined by a random walk on the graph to enhance the self-organization properties of complex networks and outputs a network with adjusted edge weights. When the NR diffusion method is applied to the original cell–peak graph, the resulting network exhibits the following characteristics: (1) edge weights between cells are imputed. The edges between cells weighted zero in the original bipartite graph. After NR diffusion, the corresponding edge weights in the diffused network can reflect the similarity between cells regarding chromatin accessibility across all peaks. (2) Edge weights between cells and peaks are adjusted. The edges between cells and peaks originally had weights representing the level of peak accessibility in the cells from scATAC-seq data. After NR diffusion, the edge weights can be interpreted as predictive scores for peak accessibility in cells based on the aggregation of multiple source information. (3) Edge weights between peaks are imputed. The edges between peaks had prior weights computed based on genome distance. After NR diffusion, the edge weights in the diffused network can reflect the similarity between peaks in terms of chromatin accessibility across all cells.

We used \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $G=\left(V,E,W\right)$\end{document} to denote a graph/network, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $V$\end{document} represents vertex set of the graph, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $E\subset \left(V,V\right)$\end{document} represents edge set of the graph and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $W$\end{document} captures the weights of the elements in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $E$\end{document}. NR takes network\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${G}_{in}\left(V,E,W\right)$\end{document} as input and outputs a diffused network represented as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${G}_{out}\left(V,E,{W}^{\prime}\right)$\end{document}. The core of the NR method is the graph diffusion operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${F}_m$\end{document}, which is defined based on the random walk diffusion operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${f}_m$\end{document}. Specifically, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${f}_m$\end{document} transforms a transition matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $P$\end{document} to another by adding the probability of all paths of different lengths joining two nodes, with a smaller weight coefficient \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $1/{m}^k$\end{document} for a longer path of length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $k$\end{document}:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ {f}_m:\kern0.75em \mathbf{\mathcal{P}}\to \mathbf{\mathcal{P}} $$\end{document}

(1) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*} {f}_m(P)=\frac{1}{\sum_{k=1}^{\infty}\left(1/{m}^k\right)}\sum_{k=1}^{\infty}\frac{P^k}{m^k}=\left(m-1\right)P{\left( mI-P\right)}^{-1} \end{equation*}\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\sum_k\left(1/{m}^k\right)$\end{document} is a normalization factor and\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $m>1$\end{document} to ensure that the series converges when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $k$\end{document} approaches infinity [15]. The parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $m$\end{document} controls the diffusion intensity. The smaller the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $m$\end{document}, the higher the degree of diffusion (default as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $m=1.5$\end{document}). Notably, the operator \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${f}_m$\end{document} has described a signal enhancement process by accumulating the transition probabilities of all lengths of paths between the two nodes, which can also be regarded as a process of signal diffusion.

In order to perform the diffusion process defined by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${f}_m$\end{document}, we should first define the random walk on the observed network as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $P={D}_W^{-1}W$\end{document} and recover the diffused random walk matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${f}_m(P)$\end{document}to the undirected graph. Specifically, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ {F}_m(W)=\frac{1}{\sum_{k=1}^{\infty}\left(1/{m}^k\right)}\sum_{k=1}^{\infty}\frac{D_W{\left({D}_W^{-1}W\right)}^k}{m^k} $$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ =\frac{m-1}{m}{D}_W\left({D}_W^{-1}W\right)+\frac{m-1}{m^2}{D}_W{\left({D}_W^{-1}W\right)}^2+\dots $$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ =\frac{m-1}{m}W+\frac{m-1}{m^2}W{D}_W^{-1}W+\frac{m-1}{m^3}W{D}_W^{-1}W{D}_W^{-1}W\dots $$\end{document}

As we described in our previous work [15], NR is a degree-normalized version of the Katz index, which makes the path with higher intermediate-node degrees less important because their information is dispersed through more adjacent edges.

Dimensionality reduction

We used the uncentred principal component analysis (UCPCA) [19] to obtain the low-dimensional embeddings of the corresponding cells and peaks in this study. The number of kept components was determined based on the proportion of variance that can be explained [20]. Specifically, we kept the first 50 (if peaks number less than 50 000) or 100 (if peaks number more than 50 000) principal components (PCs) and plotted the standard deviation (SD) of each PC. The number of retained PCs was determined based on the SD variance of the unselected PCs. If the SD of the remaining PCs did not change significantly, we did not retain them. Specifically, denoting the SD vector of the top 50 or 100 PCs as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${v}_{SD}$\end{document}, we use \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${v}_{SD}\left[k:\right]$\end{document} to represent the SD vector of the last \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $k$\end{document}unselected PCs. The number of retained PCs is chosen to be \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $K$\end{document} if the variance of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${v}_{SD}\left[K:\right]$\end{document} is greater than 5e-3, and the variance of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${v}_{SD}\left[\left(K+1\right):\right]$\end{document} is less than 5e-3.

We give in Supplementary Figure 2 the SD plots and the number of PCs retained on all scATAC-seq datasets in this study. Finally, we applied L2-normalization to each of the low-dimensional representations, which was a classical processing step [20, 21] after dimensionality reduction. We compared the impact of different dimensionality reduction methods and clustering methods on the results of SCARP on subsequent benchmark scATAC-seq datasets and found that UCPCA exhibited relatively stable performance (Supplementary Figure 3).

Evaluation metrics

We used the normalized mutual information (NMI) and the adjusted Rand index (ARI) to evaluate the clustering performance. Assuming that there are two partitions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $X=\left\{{X}_1,\dots{X}_r\right\}$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $Y=\left\{{Y}_1,\dots{Y}_s\right\}$\end{document} for\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $N$\end{document}vertexes, ARI is computed as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ \mathrm{ARI}=\frac{\sum_{ij}\left(\begin{array}{@{}c@{}}{n}_{ij}\\{}2\end{array}\right)-\frac{\sum_i\left(\begin{array}{@{}c@{}}{a}_i\ \\{}2\end{array}\right)\sum_j\left(\begin{array}{@{}c@{}}{b}_j\\{}2\end{array}\right)}{\left(\begin{array}{@{}c@{}}n\\{}2\end{array}\right)}}{\frac{1}{2}\left(\sum_i\left(\begin{array}{@{}c@{}}{a}_i\ \\{}2\end{array}\right)+\sum_j\left(\begin{array}{@{}c@{}}{b}_j\\{}2\end{array}\right)\right)-\frac{\sum_i\left(\begin{array}{@{}c@{}}{a}_i\ \\{}2\end{array}\right)\sum_j\left(\begin{array}{c}{b}_j\\{}2\end{array}\right)}{\left(\begin{array}{@{}c@{}}n\\{}2\end{array}\right)}} $$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${n}_{ij}$\end{document} is the number of vertexes in partitions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${X}_i$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${Y}_j$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${a}_i$\end{document} is the number of vertexes in partition \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${X}_i$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${b}_j$\end{document} is the number of vertexes in partition \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${Y}_j$\end{document}. The ARI takes values from −1 to 1, and a higher value of ARI indicates better performance. Besides, NMI is computed as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ NMI\left(X,Y\right)=\frac{MI\left(X,Y\right)}{\sqrt{H(X)H(Y)}} $$\end{document}

where the mutual information of two partitions and the entropy of a partition are computed as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ MI\left(X,Y\right)=\sum_{ij}\frac{n_{ij}}{N}\log \frac{n_{ij}/N}{a_i/N\times{b}_j/N} $$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ H(X)=\sum_i\frac{a_i}{N}\log \frac{a_i}{N},H(Y)=\sum_j\frac{b_j}{N}\log \frac{b_j}{N} $$\end{document}

NMI takes values from 0 to 1, and a higher value of NMI indicates better performance.

We also used the silhouette coefficient to evaluate the distance between clusters, given a partition for cells. Specifically, for cell \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $x$\end{document}:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $$ Silhouette\ (x)=\frac{b(x)-a(x)}{\max \left\{a(x),b(x)\right\}} $$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $a(x)$\end{document} is the average distance between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $x$\end{document} and other cells in the same cluster and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $b(x)$\end{document} is the average distance between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $x$\end{document} and cells in the closest different cluster [22]. The distance is measured by their Euclidean distance on a two-dimensional visualization (such as Uniform Manifold Approximation and Projection (UMAP) and t-Distributed Stochastic Neighbor Embedding (t-SNE)). The silhouette coefficient of the entire partition was the average over all cells. The silhouette coefficient takes values from −1 to 1, and a higher value indicates better performance.

RESULTS

SCARP exhibited superior cell clustering performance on benchmarking scATAC-seq datasets

We selected nine benchmarking scATAC-seq datasets with reference cell type annotations. For these datasets, the number of peaks ranges from 7000 to 140 000 with different sparsity levels (Figure 2A and Supplementary Table 1). We compared SCARP with six state-of-the-art methods (scOpen [22], DCA [23], MAGIC [24], PCA, cisTopic [2] and scAND [6]) as well as SCARP without using prior weights of adjacency peaks (denoted as SCARP*). Notably, DCA and MAGIC were especially designed for scRNA-seq analysis，and we wanted to verify whether they were also valid for scATAC-seq data (Supplementary Materials). The running time of SCARP is advantageous when compared with these methods (Figure 2B and Supplementary Table 3). Using blood2K dataset as an example, we found stronger consistency between the identified cell clusters using SCARP cell embeddings and the annotated cell types when compared with other candidate methods (Figure 2C and Supplementary Figures 4 and 5). The SCARP-derived low dimensional visualization of cells also exhibits more explicit cell type boundaries and a consistent developmental trajectory with known Fluorescence Activated Cell Sorting (FACS)-sorted population [6, 25] (Figure 2D and Supplementary Figures 6–9). For all benchmarking datasets, SCARP achieved the top-ranked clustering performances with different evaluation metrics (Figure 2E and Supplementary Figure 10). It is noteworthy that SCARP* ranked second, meaning that both NR diffusion and prior weights play essential roles in obtaining better low-dimensional cell representations. Together with the observations that SCARP had the best averaged clustering performances regarding different metrics (Figure 2F and G), all these pieces of evidence comprehensively demonstrated the superiority of SCARP over other methods in cell feature extraction.

Figure 2 Cell clustering performance on benchmarking scATAC-seq datasets. (A) Data presentation. Each dot represents a scATAC-seq dataset, with the x-axis (y-axis) showing its number of peaks (cells), the size of the dots showing the number of annotated cell types, and the color of the dots showing the data sparsity. (B) Running times (y-axis) of various methods on nine benchmarking datasets (x-axis). One method is shown by one line. The line labeled 'SCARP' represents the time taken to obtain the NR diffused matrix, while the line marked 'SCARP (+)' indicates the time required to obtain the final cell embeddings. (C) The confusion matrix of SCARP-derived cell embeddings clustered by Louvain (y-axis) and the annotated cell types (x-axis) on the blood2K dataset. (D) The UMAP plots of four methods on the blood2K dataset, colored by the annotated cell types. The rightest legend was the known FACS-sorted population of origin, which shows the process of cell development and differentiation. (E) The cell clustering results of nine methods (x-axis) on eleven datasets (y-axis) under ARI evaluation. Colors represent the row z-scores calculated from the ARI score of the corresponding dataset. (F) Cell clustering accuracy of various methods. Scores were calculated by averaging the ARI scores and NMI scores for a method over all datasets. (G) Cell distance accuracy of various methods. Scores were calculated by averaging the silhouette coefficient scores for a method over all datasets.

Since cell type annotations based on scRNA-seq data are more readily available, we further investigate whether the scATAC-derived cell clustering is consistent with the scRNA-based results, using the SNARE-seq data, which are paired multi-omics data that simultaneously profile gene expressions and chromatin accessibilities [26]. We used the cell clustering results on RNA-seq data processed by Seurat [20] and Signac [21] as a reference and compared their consistency with those of SCARP, scAND and cisTopic on scATAC-seq data. We observed that SCARP yielded higher cell clustering concordance between paired multi-omics data (Figure 3A).

Figure 3 Robustness of SCARP. (A) Four UMAP visualizations show the cell clustering results of Seurat & Signac based on scRNA-seq data, as well as cisTopic, scAND, and SCARP based on scATAC-seq data. Sil refers to silhouette score. The UMAP plots were colored by the annotated cell types obtained by Seurat. The right panel displays the clustering accuracy scores below the UMAP plot title. (B) The x-axis represents different peak filtering strategies, and the y-axis on the left panel represents the number of peaks in the blood2K dataset, the y-axis on the right panel represents the running time of SCARP. (C) Barplots of cell clustering performance for three methods under different strategies of peak filtering. Scores were calculated by averaging the ARI score and NMI score for one method. (D) The ARI scores (y-axis) of various methods on the Leukemia dataset under different downsampling ratios (x-axis). (E) The average of NMI and ARI scores (y-axis) when taking different values of parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\beta$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $m$\end{document}, and\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\gamma$\end{document} on various scATAC-seq datasets (x-axis).

We next explored the robustness of SCARP in dealing with different peak numbers, missing data and different model parameters. Since the blood2K dataset has the maximum number of peaks, we used it to test the impact of peak filtering to the performance of SCARP. We found that the number of peaks and the running time of SCARP decreased as the filtering level increased under both count-based and variance-based peaks filtering strategies (Figure 3B). In contrast, the accuracy of cell clustering of SCARP had not been much affected, and the cell clustering performances were still better than that of two other competitive methods for most filtering levels (Figure 3C). Besides, we randomly down-sampled the counts in Leukemia dataset, which has the minimum number of peaks, to simulate the missing values in scATAC-seq data. We observed that our method was not sensitive to missing data for different down-sampling ratios and got the best performance overall (Figure 3D and Supplementary Figure 11). Furthermore, SCARP has three model parameters (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\beta$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $m$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\gamma$\end{document}). We found that SCARP is robust to the parameter selection since the cell clustering results only slightly fluctuated (Figure 3E and Supplementary Figures 12–14).

Taking together, SCARP is a fast and robust method to extract informative low-dimensional cell embeddings and has an outstanding performance of cell clustering compared with all the other methods.

SCARP identified a new CD14 monocyte subpopulation characterized with high differentiation activities

We applied SCARP to a 10X multiome peripheral blood mononuclear cell (PBMC) dataset [27], where scRNA-seq and scATAC-seq data can be acquired simultaneously. The cell clusters generated using SCARP-derived cell embeddings were well separated between different cell types (Figure 4A). Interestingly, we observed that the cells originally annotated as CD14 monocytes were divided into two groups (labeled as cluster 1 and cluster 2) in UMAP visualization (Figure 4A and Supplementary Figure 15). After excluding the influence of data quality (Supplementary Figure 16), this indicates the existence of possible subpopulations of CD14 monocytes. Differential expression analysis showed that cluster 2 monocytes were characterized by the overexpression of genes BCL11B, CD247 and some interleukin-related genes (Figure 4B and D). Several marker peaks were also identified corresponding to each monocyte subpopulation (Figure 4C and E and Supplementary Figures 17–19). Additionally, since certain recognized marker genes (and their promoters) of known cell types were not significantly differentially expressed (accessible), we ruled out the potential that these two cell subpopulations are other recognized cell types, like T cells and macrophages (Figure 4B and C).

Fig 4 SCARP identified a new CD14 monocyte subpopulation. (A) UMAP visualization of SCARP -derived cell embeddings colored by the cell types annotated by Seurat. (B, C) Dot plots showing the expression levels of the top 10 marker genes (B) and accessibility levels of the top 10 marker peaks (C) for two cell subpopulations (‘clu’ means cluster), and the known marker genes for monocytes, T-cells, and macrophages. (D, E): The UMAP visualizations of cells colored by the expressions of marker genes (D) and the accessibility of marker peaks (E) for two clusters. (F–H): The UMAP visualizations of cells in external M-CSF scRNA-seq dataset, colored by gene signature scores for CD14 monocytes (donor 1) under M-CSF stimulation at day 0 (f), day 3 (G) and day 6 (H), along with the violin plots showing the difference between the two cell clusters. The p-value was calculated based on the Wilcoxon test. (I) Top biological processes that gene signatures are involved in. The y-axis represents the name of the KEGG pathway, and the x-axis represents the ratio of genes enriched in a pathway. The number of genes enriched in a pathway is related to the size of the circle (gene count less than 10) or triangle (gene count greater than 10), and the color represents the corresponding FDR-adjusted P-value.

To ascertain whether the newly identified CD14 monocyte subpopulation is biologically meaningful, we built a gene signature mono-c2 using the top 50 marker genes of cluster 2 monocytes. After this, we applied our mono-c2 signature to an external time-series scRNA-seq dataset of human CD14 monocytes, where CD14 monocytes were stimulated by macrophage colony-stimulating factor (M-CSF) at days 0, 3 and 6 [28]. We observed that only 15% of monocytes were clustered together, and 25% of them had relatively higher mono-c2 signature scores at day 0 (Figure 4F, Wilcoxon-test P = 1.42e-37). This finding became more apparent at day 3, with approximately half of the monocytes having significantly higher signature scores (Figure 4G, Wilcoxon-test P = 1.1e-74). With the continuous stimulation of M-CSF, almost all cells (86%) exhibited high signature scores, and the difference between the two clusters was less significant (Figure 4H, Wilcoxon-test P = 2.77e-08). These observations can be reproduced using other replicates (Supplementary Figures 20 and 21). Besides, we also identified the new CD14 monocyte subpopulation in another external CD14 monocyte scRNA-seq dataset with higher signature scores (Supplementary Figure 22).

Since monocytes are precursors to macrophages and can differentiate into dendritic cells, we hypothesized that the mono-c2 signature scores can serve as indicators of monocyte differentiation activities, and the monocytes with high signature scores are prone to differentiate into other cells. We thereby conducted GO functional enrichment analysis using the signature genes. We found that these genes were associated with mononuclear cell differentiation and its related processes, such as T-cell differentiation [29] and immune response-activating signal transduction [30] (Figure 4I). These results confirmed the ability of our gene signature to identify CD14 monocyte subtypes.

In conclusion, we demonstrated that the SCARP-derived low-dimensional representation of cells can refine cell types and reveal biologically meaningful cell subtypes.

SCARP discovered key factors of melanoma progression by portraying co-accessibility relationships of peaks

Besides for the low-dimensional cell embeddings, SCARP can also generate biologically meaningful peak embeddings. In this experiment, we want to investigate whether the SCARP-derived peak embeddings can uncover potential regulatory relationships in melanoma disease progression. To achieve this, we applied SCARP on a SOX10 knockdown time series scATAC-seq dataset (SOX10KD) [2], which captures peak accessibility at 0, 24, 48 and 72 h after SOX10 knockdown in two melanoma cultures (MM057 and MM087) derived from patient biopsies.

We first annotated two regions related to SOX10, i.e. the promoter (chr22:38380499-38380726) and 3’UTR (chr22:38364975-38365257) of SOX10 (Supplementary Figure 23). Since the promoter initiates gene transcription while 3’UTR can decrease the expression of the corresponding gene, it is reasonable to see that the accessibility of SOX10 promoter increases and 3’UTR decreases as knockdown time passes. By calculating the Spearman correlations between peaks’ low-dimensional representations, we obtained the top 10 co-accessible regions with the highest correlation to each of the SOX10 promoter and 3’UTR. These loci have the more similar time-varying trend as the SOX10 promoter (3’UTR) than that found by using original data, indicating that SCARP can well characterize the co-accessibility of accessible segments (Figure 5A). To better interpret these co-accessible regions, we obtained the corresponding genes using the promoter regions (Supplementary Figure 24) and selected the top 50 genes with the most co-accessible relationships with SOX10. By constructing a gene regulatory network (GRN) between these genes, we found abundant database-supported regulatory relationships (Figure 5C and Supplementary Table 4). Notably, SCARP owned the largest number of curated edges and the second number of verified genes in the constructed GRNs compared to other candidate methods, such as cisTopic and scAND (Figure 5B).

Figure 5 SCARP portrayed the co-accessibility of peaks on SOX10KD scATAC-seq dataset. (A) Changes in accessibility for those loci (x-axis) that were highly correlated with SOX10 promoter and 3’UTR as knockdown time walked by (y-axis). MM057 and MM087 are two melanoma cultures derived from patient biopsies. The ‘MM057_0hr’ cell group contains the cells in melanoma culture MM057, 0 h after SOX10 was knocked down. Each dot shows the mean accessibility across all cells within a cell group, corresponding to a specific peak. To quantitatively compare the time-varying trend similarity, the violin plot at the right bottom showed the trend consistency metric (Supplementary materials) for the top 10 peaks co-accessible with the promoter (3’UTR) of SOX10. (B) The transcriptional regulatory relationships are validated by the Reactome database. Solid lines indicate experimentally validated regulatory relationships, while dashed lines indicate database-predicted regulatory relationships. (C) For each of the differential accessibility gene sets obtained by four methods, we compared the number of genes (y-axis of the lower panel) that were verified as directly or indirectly related by the database, and the number of genes that were not verified; as well as the number of edges (y-axis of the upper panel) in the corresponding GRNs, including predicted edges, i.e. those with a functional interaction (FI) type of ‘predicted’, and database edges, i.e. those with an FI type other than ‘predicted’, such as ‘activate’ and ‘inhibited by’. (D) Results of KEGG pathway analysis. The legends are the same as in Figure 4I. (E) The results of survival analysis on genes TRIM69, MITF and STK17B.

Among the genes in SCARP-derived GRN, numerous genes are strongly associated with melanoma disease. For example, the upregulations of STMN1 and ROR1 could contribute to tumor migration and proliferation [31–33]. LTA4H regulates the cell cycle process and can lead to skin carcinogenesis by its overexpression [34]. Notably, we uniquely identified that HIF1A, which is known for promoting cancer progression and causing therapy resistance [35, 36], interacts with the key master regulator of melanocyte development and melanoma deterioration MITF [37, 38] (Figure 5C). Numerous studies have reported that hypoxia is closely related to the invasiveness, angiogenesis and therapy response in melanoma [39, 40]. In studies of its pathogenic mechanism, MITF has been shown to stimulate the transcriptional activity of HIF1A to supply oxygen to cancer cells [41, 42].

We further conducted the functional enrichment analysis using the genes in SCARP-derived GRN. As a result, we found many melanoma progression–related KEGG pathways, e.g. melanogenesis [43, 44], MAPK signaling pathways [45, 46] and adherens junction [47, 48] (Figure 5D and Supplementary Figure 26). Notably, the HIF-1 signaling pathway was again exclusively identified by SCARP (Figure 5D, Supplementary Figures 28–30 and Supplementary Table 5) [43, 49]. In addition, GO enrichment categories contained several melanoma-related terms, such as muscle tissue development [50] and membrane raft [51] (Supplementary Figure 27). Furthermore, we performed survival analysis using the genes in SCARP-derived GRN and found that many genes’ abnormal expressions were associated with the poor prognoses, such as TRIM69, MITF and STK17B (Figure 5E and Supplementary Figure 31). In particular, the low expression of gene STK17B, which has been shown in the literature to play an important role in immune infiltration, was associated with skin cancer, making it a key biomarker for the diagnosis of melanoma disease [52].

SCARP predicted reliable cis-regulatory interactions supported by external evidence

We next explored the potential of SCARP-derived peak embeddings and investigated their cis-regulatory relationships. We defined the co-accessibility score of two peaks as the cosine similarity of their low-dimensional peak embeddings, which can serve as the confidence of corresponding cis-regulatory interactions. We then used external evidence to assess the reliability of this score, such as promoter-capture Hi-C (PCHi-C) [17, 53] and ChIP-seq [54] data.

Specifically, we selected those peaks belonging to promoter regions of certain genes and focused on the regulatory relationships between these peaks (i.e. promoters) with other peaks. As PCHi-C data depicted whether there were physical interactions between the promoters (baits) and the entire genome (other ends), we can treat whether SCARP-derived cis-regulatory relationships were validated by PCHi-C data as a binary classification problem and measure its accuracy by calculating AUROC score. The result shows a significant difference between PCHi-C-validated and unvalidated SCARP-derived cis-regulatory scores (Figure 6A, Wilcoxon-test P < 1e-38), and the receiver operating characteristic (ROC) curve demonstrates the superiority of SCARP over other methods (Figure 6B).

Figure 6 SCARP-derived cis-regulatory interactions are meaningful and validated by external evidence. (A) The boxplot showing the difference between PCHi-C-validated and -unvalidated SCARP-derived transcriptional regulatory relationships and their co-accessibility scores (y-axis), followed by the Wilcoxon rank-sum P-value to test the difference. (B) ROC curve of four methods (SCARP, original, cisTopic, scAND), treating whether method-derived transcriptional regulatory relationships were validated by PCHi-C data as a binary classification problem. (C, D) Visualizations of co-accessibility scores (y-axis) of SCARP-derived cis-regulatory interactions, PCHi-C evidence, and CHIP-seq evidence around the promoters of genes (C) BCL3 and (D) ZBTB14. BS means binding site. Those regions with strong co-accessibility to the BCL3 promoter and ZBTB14 promoter were highly enriched in motifs MA1508 and MA1961, respectively.

Besides, we also used ChIP-seq data of human TF to validate SCARP-derived cis-regulatory relationships, which gives the experimentally verified interaction between proteins and DNA to help identify TF binding sites. As a demonstration example, we plotted cis-regulatory interactions between the promoter of gene BCL3, which is a human TF associated with human chronic lymphocytic leukemia and B-cell malignancies [55, 56], and nearby accessible regions based on multiple biological data (Figure 6C). We found that those regions with strong co-accessibility to the BCL3 promoter were highly enriched in motif MA1508, which preferred to be bound by TF IKZF1 (Figure 6C). The IKZF1 belongs to the IKAROS family and plays an important regulatory role in acute myeloid leukemia disease [57]. It has been shown that IKZF1 encodes a tumor suppressor that inhibits the expression of BCL family genes in B cells [58], leading to the proliferation and development of leukemia cells. In this example, we demonstrated that the peak co-accessibility obtained by SCARP can be partially verified by ChIP-seq or PCHi-C data, some of which were direct regulatory relationships (overlapping sites of BCL3-BS and ATAC peaks), and some were co-regulated by a certain gene (overlapping sites of IKZF1-BS and ATAC peaks).

In another example, we also visualized those highly supported interactions between the promoter of TF ZBTB14 (also known as ZFP161), which regulated the development of monocytes and macrophages [59] and those nearby regions (Figure 6D), and similar conclusions are confirmed. To sum up, we show that the peak co-accessibility portrayed by the low-dimensional representations obtained by SCARP can be supported by external evidence data.

CONCLUSION AND DISCUSSION

In this study, we proposed a novel network NDM-based computational approach to comprehensively analyze scATAC-seq data. We have demonstrated through sufficient experiments that SCARP promotes superior analyses of scATAC-seq data compared to other state-of-the-art methods. Of all methods evaluated, SCARP had the best overall clustering performances on benchmarking scATAC-seq data and was robust to different peak filter strategies and parameter selections. The running time of SCARP is also advantageous since we used some appropriate tricks to reduce the time cost. For example, to make those peaks on the same chromosome better learn from each other, we performed NR diffusion on subgraphs in parallel and spliced them back together. Apart from the benchmarking dataset, we have also designed three innovative case studies to explore what kind of analyses could be done with the scATAC-seq data to contribute to our further understanding of biological regulatory mechanisms.

Single-cell ATAC-seq data established genome-wide chromatin accessibility profiles with single-cell resolution that could complement single-cell transcriptome data, and together, they provided new insights into transcriptional regulatory mechanisms. The accessible regions of scATAC-seq data, which can be seen as cell features, exhibited many different properties from genes of scRNA-seq data, such as higher sparsity, higher dimension and near-binarization [1, 6]. Thus, even if both peaks and genes can be used to characterize cellular heterogeneity, and many scRNA-seq-based algorithms are mathematically transferable to analyze scATAC-seq, we should not be limited to previous attempts and ideas for dealing with scRNA-seq data given the different nature of the two features; instead, we should develop new approaches that focus on the specific properties and information of the scATAC-seq data itself and elaborate on cellular heterogeneity from new perspectives.

Although the severe signal deficiencies and high sparsity of scATAC-seq data are well known, the importance of filling in missing data and denoising was not given enough attention in some commonly used scATAC methods [22]. Considering the high sparsity of the data, it is very appropriate and natural to use the diffusion method for aggregating the neighborhood information between two nodes. Our previously developed NR diffusion method [15] has achieved outstanding success in many other applications, such as network denoising [15] and link prediction [16], and in this study, it also proved to have excellent performance in handling scATAC-seq data.

In addition, to mine and exploit the unique information of scATAC-seq data, we have borrowed ideas from many genomic studies and used the genetic map distance as an important priori factor in the analysis of epigenomic data. To further test the soundness of this hypothesis, we measured the co-accessibility likelihood between two peaks using the P-value of Fisher’s exact test and then investigated its relevance to genomic distance of peaks. The results shown in Supplementary Figure 32A confirm that peaks closer in genome are more likely to be co-accessible. The peaks with network weights higher (or lower) than expected after diffusion show a certain pattern on the genome (Supplementary Figure 32B–E), i.e. the posterior probabilities of co-accessibility in some regions are higher (or lower) than that predicted solely based on genomic distance. The enrichment of the different region types within the SCARP-derived low-dimensional peak embedding features also supported that there were clear differences in accessibility between the different region types, which were closely related to their regulatory functions (Supplementary Figure 25). Furthermore, the results of our experiments supported that the use of prior edge weights between adjacent peaks facilitates SCARP to mine more valuable information from scATAC-seq data and helped to obtain better downstream analysis performance. The sources of prior information extend beyond merely genomic distance. We could potentially integrate further information, including high-confidence co-accessibility scores derived from alternative methodologies, as well as the interactions between regulatory elements associated with peaks.

In the future, we expect to explore more unique properties of the scATAC-seq data to better handle it and consider doing joint analyses with other multi-omics data, such as data alignment and integration.

Key Points

SCARP is a novel network diffusion–based computational method to comprehensively analyze scATAC-seq data. It can extract meaningful cell and peak embeddings.

SCARP incorporates distance information of adjacent peaks on the genome to better capture the co-accessibility among peaks and the NR diffusion method to alleviate the high sparsity of scATAC-seq data.

SCARP has shown to be a promising tool for superior analysis of scATAC-seq data, including improving cell clustering performance to better elucidate cell heterogeneity, and characterizing co-accessibility of genome segments to shed fresh light on transcriptional regulation.

SCARP-derived low-dimensional representation of cells can refine cell types and reveal biologically meaningful cell subtypes.

FUNDING

This work has been supported by the National Key Research and Development Program of China (No. 2022YFA1004803 and No. 2020YFA0712402) and the National Natural Science Foundation of China (No. 12231018, 62202269).

DATA AVAILABILITY

All datasets analyzed in this study were publicly available. The detailed information for the scATAC-seq and the multi-omics datasets used in this study can be found in Supplementary Table 1. The gene expression and phenotype data of the Skin Cutaneous Melanoma (SKCM) project were downloaded from TCGA. For promoter-capture Hi-C data, it can be obtained from the original publication [17, 53]. For TF ChIP–seq data, it can be downloaded from the ENCODE project (https://www.encodeproject.org/) and further processed using GLUE’s script (https://github.com/gao-lab/GLUE). SCARP is available at https://github.com/Wu-Lab/SCARP and https://github.com/Wu-Lab/SCARP-reproduce.

Supplementary Material

Supplementary_Materials_Fig_1-32_and_Table_1-3_bbae093

Supplementary_Table_4_Cytoscape_Results_bbae093

Supplementary_Table_5_Gene_Enrichment_Analysis_bbae093

Author Biographies

Jiating Yu is an assistant professor at School of Mathematics and Statistics, Nanjing University of Information Science & Technology, and obtained her PhD from the Academy of Mathematics and System Sciences, Chinese Academy of Sciences, China. Her research interests include bioinformatics, complex networks and deep learning.

Jiacheng Leng received his PhD degree from the Academy of Mathematics and System Sciences, Chinese Academy of Sciences, China. He is currently a Senior Research Fellow with Zhejiang Lab, China. His research interests include bioinformatics, network inferences and deep learning.

Zhichao Hou is a graduate student at the Academy of Mathematics and System Sciences, Chinese Academy of Sciences, China.

Duanchen Sun is a professor at the School of Mathematics, Shandong University. His research interests are bioinformatics and operations research.

Ling-Yun Wu is a professor and the director of Bioinformatics Center at the Academy of Mathematics and System Sciences, Chinese Academy of Sciences, China. His research interests include bioinformatics, optimization and machine learning.
==== Refs
References

1. Xiong  L, Xu  K, Tian  K, et al.  SCALE method for single-cell ATAC-seq analysis via latent feature extraction. Nat Commun  2019;10 :1–10.30602773
2. Bravo González-Blas  C, Minnoye  L, Papasokrati  D, et al.  cisTopic: cis-regulatory topic modeling on single-cell ATAC-seq data. Nat Methods  2019;16 :397–400.30962623
3. Buenrostro  JD, Wu  B, Litzenburger  UM, et al.  Single-cell chromatin accessibility reveals principles of regulatory variation. Nature  2015;523 :486–90.26083756
4. Cusanovich  DA, Daza  R, Adey  A, et al.  Multiplex single-cell profiling of chromatin accessibility by combinatorial cellular indexing. Science  2015;348 :910–4.25953818
5. Cusanovich  DA, Reddington  JP, Garfield  DA, et al.  The cis-regulatory dynamics of embryonic development at single-cell resolution. Nature  2018;555 :538–42.29539636
6. Dong  K, Zhang  S. Network diffusion for scalable embedding of massive single-cell ATAC-seq data. Sci Bull  2021;66 :2271–6.
7. Liu  Y, Zhang  J, Wang  S, et al.  Are dropout imputation methods for scRNA-seq effective for scATAC-seq data?  Brief Bioinform  2022;23 :1–12.
8. Schep  AN, Wu  B, Buenrostro  JD, Greenleaf  WJ. ChromVAR: inferring transcription-factor-associated accessibility from single-cell epigenomic data. Nat Methods  2017;14 :975–8.28825706
9. Katz  L . A new status index derived from sociometric analysis. Psychometrika  1953;18 :39–43.
10. Speed  T.P.  Genetic Map Functions. In: Armitage P, Colton T (eds.) Encyclopedia of Biostatistics. John Wiley & Sons, Ltd, Chichester, 2005.
11. Altshuler  D, Daly  MJ, Lander  ES. Genetic mapping in human disease. Science  2008;322 :881–8.18988837
12. Robinson  M.A . Linkage disequilibrium. In: Delves PJ (ed.) Encyclopedia of Immunology. Elsevier, Oxford, 1998, 1586 –1588.
13. Teare  MD, Barrett  JH. Genetic linkage studies. The Lancet  2005;366 :1036–44.
14. Dekker  J, Marti-Renom  MA, Mirny  LA. Exploring the three-dimensional organization of genomes: interpreting chromatin interaction data. Nat Rev Genet  2013;14 :390–403.23657480
15. Yu  J, Leng  J, Sun  D, et al.  Network refinement: denoising complex networks for better community detection. Physica A Stat Mech Appl  2023;617 :128681.
16. Yu  J, Wu  LY. Multiple order local information model for link prediction in complex networks. Physica A Stat Mech Appl  2022;600 :127522.
17. Cao  ZJ, Gao  G. Multi-omics single-cell data integration and regulatory inference with graph-linked embedding. Nat Biotechnol  2022;40 :1458–66.35501393
18. Pliner  HA, Packer  JS, McFaline-Figueroa  JL, et al.  Cicero predicts cis-regulatory DNA interactions from single-cell chromatin accessibility data. Mol Cell  2018;71 :858–871.e8.30078726
19. Cadima  J, Jolliffe  I. On relationships between uncentred and column-centred principal component analysis. Pakistan Journal of Statistics  2009;25 :473–503.
20. Stuart  T, Butler  A, Hoffman  P, et al.  Comprehensive integration of single-cell data. Cell  2019;177 :1888–1902.e21.31178118
21. Stuart  T, Srivastava  A, Madad  S, et al.  Single-cell chromatin state analysis with Signac. Nat Methods  2021;18 :1333–41.34725479
22. Li  Z, Kuppe  C, Ziegler  S, et al.  Chromatin-accessibility estimation from single-cell ATAC-seq data with scOpen. Nat Commun  2021;12 :6386.34737275
23. Eraslan  G, Simon  LM, Mircea  M, et al.  Single-cell RNA-seq denoising using a deep count autoencoder. Nat Commun  2019;10 :1–14.30602773
24. van  Dijk  D, Sharma  R, Nainys  J, et al.  Recovering gene interactions from single-cell data using data diffusion. Cell  2018;174 :716–729.e27.29961576
25. Buenrostro  JD, Corces  MR, Lareau  CA, et al.  Integrated single-cell analysis maps the continuous regulatory landscape of human hematopoietic differentiation. Cell  2018;173 :1535–1548.e16.29706549
26. Chen  S, Lake  BB, Zhang  K. High-throughput sequencing of the transcriptome and chromatin accessibility in the same cell. Nat Biotechnol  2019;37 :1452–7.31611697
27. PBMC from a healthy donor, single cell multiome ATAC gene expression demonstration data by cell ranger ARC 1.0.0. 10X genomics. https://www.10xgenomics.com/resources/datasets/pbmc-from-a-healthy-donor-granulocytes-removed-through-cell-sorting-10-k-1-standard-1-0-0.
28. Hie  B, Bryson  B, Berger  B. Efficient integration of heterogeneous single-cell transcriptomes using Scanorama. Nat Biotechnol  2019;37 :685–91.31061482
29. Zimmermann  HW, Bruns  T, Weston  CJ, et al.  Bidirectional transendothelial migration of monocytes across hepatic sinusoidal endothelium shapes monocyte differentiation and regulates the balance between immunity and tolerance in liver. Hepatology  2016;63 :233–46.26473398
30. Wierda  RJ, Goedhart  M, van  Eggermond  MCJA, et al.  A role for KMT1c in monocyte to dendritic cell differentiation: epigenetic regulation of monocyte differentiation. Hum Immunol  2015;76 :431–7.25843229
31. Chen  J, Abi-Daoud  M, Wang  A, et al.  Stathmin 1 is a potential novel oncogene in melanoma. Oncogene  2012;32 :1330–7.22665054
32. Fernández  NB, Lorenzo  D, Picco  ME, et al.  ROR1 contributes to melanoma cell growth and migration by regulating N-cadherin expression via the PI3K/Akt pathway. Mol Carcinog  2016;55 :1772–85.26509654
33. Hojjat-Farsangi  M, Ghaemimanesh  F, Daneshmanesh  AH, et al.  Inhibition of the receptor tyrosine kinase ROR1 by anti-ROR1 monoclonal antibodies and siRNA induced apoptosis of melanoma cells. PloS One  2013;8 :e61167.23593420
34. Oi  N, Yamamoto  H, Langfald  A, et al.  LTA4H regulates cell cycle and skin carcinogenesis. Carcinogenesis  2017;38 :728–37.28575166
35. D’Aguanno  S, Mallone  F, Marenco  M, et al.  Hypoxia-dependent drivers of melanoma progression. J Exp Clin Cancer Res  2021;40 :159.33964953
36. Jing  X, Yang  F, Shao  C, et al.  Role of hypoxia in cancer therapy by regulating the tumor microenvironment. Mol Cancer  2019;18 :1–15.30609930
37. Levy  C, Khaled  M, Fisher  DE. MITF: master regulator of melanocyte development and melanoma oncogene. Trends Mol Med  2006;12 :406–14.16899407
38. Hartman  ML, Czyz  M. MITF in melanoma: mechanisms behind its expression and activity. Cell Mol Life Sci  2015;72 :1249–60.25433395
39. Hwang  HW, Baxter  LL, Loftus  SK, et al.  Distinct microRNA expression signatures are associated with melanoma subtypes and are regulated by HIF1A. Pigment Cell Melanoma Res  2014;27 :777–87.24767210
40. Lakhter  AJ, Lahm  T, Broxmeyer  HE, Naidu  SR. Golgi associated HIF1a serves as a reserve in melanoma cells. J Cell Biochem  2016;117 :853–9.26375488
41. Feige  E, Yokoyama  S, Levy  C, et al.  Hypoxia-induced transcriptional repression of the melanoma-associated oncogene MITF. Proc Natl Acad Sci U S A  2011;108 :E924–33.21949374
42. Buscà  R, Berra  E, Gaggioli  C, et al.  Hypoxia-inducible factor 1α is a new target of microphthalmia-associated transcription factor (MITF) in melanoma cells. J Cell Biol  2005;170 :49–59.15983061
43. Slominski  A, Kim  TK, Brozyna  AA, et al.  The role of melanogenesis in regulation of melanoma behavior: Melanogenesis leads to stimulation of HIF-1α expression and HIF-dependent attendant pathways. Arch Biochem Biophys  2014;563 :79–93.24997364
44. Slominski  A, Zbytek  B, Slominski  R. Inhibitors of melanogenesis increase toxicity of cyclophosphamide and lymphocytes against melanoma cells. Int J Cancer  2009;124 :1470–7.19085934
45. Sullivan  RJ, Flaherty  K. MAP kinase signaling and inhibition in melanoma. Oncogene  2012;32 :2373–9.22945644
46. Fecher  LA, Amaravadi  RK, Flaherty  KT. The MAPK pathway in melanoma. Curr Opin Oncol  2008;20 :183–9.18300768
47. Lee  DJ, Kang  DH, Choi  M, et al.  Peroxiredoxin-2 represses melanoma metastasis by increasing E-cadherin/β-catenin complexes in adherens junctions. Cancer Res  2013;73 :4744–57.23749642
48. Korla  PK, Chen  CC, Gracilla  DE, et al.  Somatic mutational landscapes of adherens junctions and their functional consequences in cutaneous melanoma development. Theranostics  2020;10 :12026–43.33204327
49. Malekan  M, Ebrahimzadeh  MA, Sheida  F. The role of hypoxia-inducible factor-1alpha and its signaling in melanoma. Biomed Pharmacother  2021;141 :111873.34225012
50. Moss  ALH, Rees  MJW. Metastatic malignant melanoma in muscle. Br J Plast Surg  1984;37 :250–2.6713164
51. Baruthio  F, Quadroni  M, Rüegg  C, Mariotti  A. Proteomic analysis of membrane rafts of melanoma cells identifies protein patterns characteristic of the tumor progression stage. Proteomics  2008;8 :4733–47.18942674
52. Shi  X, Zhou  Q, Huang  B, et al.  Prognostic and immune-related value of STK17B in skin cutaneous melanoma. PloS One  2022;17 :e0263311.35171924
53. Javierre  BM, Sewitz  S, Cairns  J, et al.  Lineage-specific genome architecture links enhancers and non-coding disease variants to target gene promoters. Cell  2016;167 :1369–1384.e19.27863249
54. de  Souza  N . The ENCODE project. Nat Methods  2012;9 :1046.23281567
55. Yabumoto  K, Ohno  H, Doi  S, et al.  Involvement of the BCL3 gene in two patients with chronic lymphocytic leukemia. Int J Hematol  1994;59 :211–8.8011990
56. Mckeithan  TW, Takimoto  GS, Ohno  H, et al.  BCL3 rearrangements and t(14;19) in chronic lymphocytic Leukemia and other B-cell malignancies: a molecular and cytogenetic study. Genes Chromosomes Cancer  1997;20 :64–72.9290956
57. Zhang  X, Zhang  X, Li  X, et al.  The specific distribution pattern of IKZF1 mutation in acute myeloid leukemia. J Hematol Oncol  2020;13 :140.33081843
58. Ge  Z, Zhou  X, Gu  Y, et al.  Ikaros regulation of the BCL6/BACH2 axis and its clinical relevance in acute lymphoblastic leukemia. Oncotarget  2017;8 :8022–34.28030830
59. Deng  Y, Wang  H, Liu  X, et al.  Zbtb14 regulates monocyte and macrophage development through inhibiting pu.1 expression in zebrafish. Elife  2022;11 :e80760.36205309
