Clustering of scATAC-seq data

We used the peak matrixes and fragment files to create chromatin assay by “CreateChromatinAssay” function with default parameters in Signac (v1.8.0) (Stuart et al., 2021). Seurat object was created by “CreateSeuratObject” function for each library. Subsequently, all library datasets were merged together by “merge” function. Then, the data were normalized with the Term Frequency-Inverse Document Frequency methods by “RunTFIDF” function. Top variable features were calculated by the “FindTopFeatures” function with “min.cutoff = ‘q70’” and then the singular value decomposition was computed by “RunSVD” function. Next, Uniform Manifold Approximation and Projection analysis was performed by “RunUMAP” function with parameter “dims = 2:30, reduction = ‘lsi’, n.neighbors = 50, min.dist = 0.4.” “FindNeighbors” and “FindClusters” functions with parameters “algorithm = 3, resolution = 3” were performed to produce cell-type clusters. Finally, the gene activity scores for each gene in cells were calculated by “addGeneScoreMatrix” function, and the cluster-specific genes based on gene scores were calculated by “getMarkerFeatures” function with “useMatrix = ‘GeneScoreMatrix’” (FDR ≤ 0.01 and log2FC ≥ 1) in ArchR. Peak calling for each cell type was performed using “addReproduciblePeakSet” function in MACS2. The cell-type specific peaks were calculated by “getMarkerFeatures” function with “useMatrix = ‘PeakMatrix’” (FDR ≤ 0.01 and log2FC ≥ 1). Motif enrichment analysis was performed by “peakAnnoEnrichment” function with “cutOff = ‘FDR ≤ 0.01 & Log2FC ≥ 0.5.’”