前面在写ATACseq学习笔记的时候，后台有很多粉丝询问scATAC-seq有没有相关资料，现在就来学习，开个新专辑《scATAC-seq分析》~
scATAC-seq用的最多的包就是 Signac 包了，这个包来自Seurat包的扩展，分析步骤等跟Seurat非常像，下面来看看。学习官网为：
https:///signac/
单细胞的ATAC-seq学习还是需要了解一些背景知识的 ，不然就是茫然的纯跑代码：ATAC-Seq 数据分析2025
这个包可以分析的测序数据类型包括 scATAC-seq ， scCUT&Tag ， scNTT-seq ，功能包括： 与 Seurat 、 SeuratWrappers 、 SeuratDisk 以及 SeuratData （Seurat数据集分发）的功能实现无缝对接。 安装 ## 使用西湖大学的 Bioconductor镜像 options(BioC_mirror= "https://mirrors./bioconductor" ) options( "repos" =c(CRAN= "https://mirrors./CRAN/" )) setRepositories(ind=1:3) install.packages( "Signac" ) rm(list=ls()) ## 加载包 library(Signac) library(Seurat) library(GenomicRanges) library(ggplot2) library(patchwork) 先来学习一波官网的案例，与Seurat软件一样，官方给了很多示例，今天学习：https:///signac/articles/pbmc_vignette
0、示例数据 这个数据来自10x Genomics，为组织类型为： human peripheral blood mononuclear cells (PBMCs) 。
下载地址为：
rawdata ：https://cf./samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_filtered_peak_bc_matrix.h5metadata ：https://cf./samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_singlecell.csvfragment文件 ：https://cf./samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_fragments.tsv.gzfragement index文件 ：https://cf./samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_fragments.tsv.gz.tbi下载：
# bash命令，或者自己选择用浏览器，迅雷都可以 wget https://cf./samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_filtered_peak_bc_matrix.h5 wget https://cf./samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_singlecell.csv wget https://cf./samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_fragments.tsv.gz wget https://cf./samples/cell-atac/2.1.0/10k_pbmc_ATACv2_nextgem_Chromium_Controller/10k_pbmc_ATACv2_nextgem_Chromium_Controller_fragments.tsv.gz.tbi 1、数据预处理 上游 CellRanger的输出数据为：
Peak/Cell matrix ：矩阵的每一行代表基因组的一个区域（一个peak），该区域被预测为开放染色质的区域。矩阵中的每个值表示每个barcode（即一个细胞）映射到每个峰值内的Tn5整合位点的数量。
更详细的说明：https://support./single-cell-atac/software/pipelines/latest/output/matrices
这里下载的是h5文件，如果是文件夹，则是标准的三个文件：
filtered_peak_bc_matrix ├── barcodes.tsv ├── peaks.bed └── matrix.mtx peaks.bed的内容：
head filtered_peak_bc_matrix/peaks.bed chr1 237588 237917 chr1 564444 565537 chr1 567478 568248 chr1 569021 569641 chr1 713461 715293 chr1 752379 753032 chr1 762073 763379 chr1 773651 774064 chr1 779547 780286 chr1 793345 794375 读取数据并创建seurat对象： # 1.读取数据： counts <- Read10X_h5(filename = "PBMCs/10k_pbmc_ATACv2_nextgem_Chromium_Controller_filtered_peak_bc_matrix.h5" ) metadata <- read.csv( file = "PBMCs/10k_pbmc_ATACv2_nextgem_Chromium_Controller_singlecell.csv" , header = TRUE , row.names = 1 ) # 简单的探索一下数据 head(metadata) dim(counts) # 每一行为一个peak区间，每一列为一个细胞 counts[ 1 : 5 , 1 : 5 ] # 创建对象 chrom_assay <- CreateChromatinAssay( counts = counts, sep = c( ":" , "-" ), fragments = "PBMCs/10k_pbmc_ATACv2_nextgem_Chromium_Controller_fragments.tsv.gz" , min.cells = 10 , min.features = 200 ) # 生成seurat对象 pbmc <- CreateSeuratObject( counts = chrom_assay, assay = "peaks" , meta.data = metadata ) pbmc ## An object of class Seurat ## 165434 features across 10246 samples within 1 assay ## Active assay: peaks (165434 features, 0 variable features) ## 2 layers present: counts, data ATACseq数据存储在 ChromatinAssay 里面：
pbmc[[ 'peaks' ]] # ChromatinAssay data with 165434 features for 10246 cells # Variable features: 0 # Genome: # Annotation present: FALSE # Motifs present: FALSE # Fragment files: 1 # 查看基因组范围 granges(pbmc) # GRanges object with 165434 ranges and 0 metadata columns: # seqnames ranges strand # <Rle> <IRanges> <Rle> # [1] chr1 9772-10660 * # [2] chr1 180712-181178 * # [3] chr1 181200-181607 * # [4] chr1 191183-192084 * # [5] chr1 267576-268461 * # ... ... ... ... # [165430] KI270713.1 13054-13909 * # [165431] KI270713.1 15212-15933 * # [165432] KI270713.1 21459-22358 * # [165433] KI270713.1 29676-30535 * # [165434] KI270713.1 36913-37813 * # ------- # seqinfo: 35 sequences from an unspecified genome; no seqlengths 去除非染色体序列的部分： # 去除非染色体的部分 standardChromosomes(granges(pbmc)) seqnames(granges(pbmc)) peaks.keep <- seqnames(granges(pbmc)) % in % standardChromosomes(granges(pbmc)) table(peaks.keep) pbmc <- pbmc[as.vector(peaks.keep), ] pbmc 添加基因信息 这里需要注意与上游分析使用的参考基因需要是同一个版本 ，这个数据的10x Genomics用的版本为： GRCh38-2020-A ，对应 Ensembl v98 patch release ，他们的对应关系可以从这里获取：
https://www./info/website/archives/assembly.html
# 添加基因注释 library(AnnotationHub) ah <- AnnotationHub() save(ah, file = "ah.RData" ) # 保存以便下次使用 # Search for the Ensembl 98 EnsDb for Homo sapiens on AnnotationHub query(ah, "EnsDb.Hsapiens.v98" ) ensdb_v98 <- ah[[ "AH75011" ]] # extract gene annotations from EnsDb annotations <- GetGRangesFromEnsDb(ensdb = ensdb_v98) # change to UCSC style since the data was mapped to hg38 seqlevels(annotations) <- paste0( 'chr' , seqlevels(annotations)) genome(annotations) <- "hg38" # 将基因信息添加到对象中 Annotation(pbmc) <- annotations 2、计算QC指标 scATAC-seq的数据指标与bulk ATAC-seq的原理一致，有下面这些指标：
Nucleosome banding pattern ：DNA片段大小的直方图应该显示出一个强烈的核小体条带模式，这对应于缠绕在单个核小体周围的DNA长度。我们针对每个单细胞计算这一模式，并定量估算单核小体片段与无核小体片段的近似比例（存储为 nucleosome_signal ）； TSS富集打分 ：ENCODE项目定义了一个基于转录起始位点（TSS）中心的片段与TSS侧翼区域片段比例的ATAC-seq靶向评分（详见 https://www./data-standards/terms/ ）。通常，质量不佳的ATAC-seq实验会有较低的TSS富集评分。 我们可以使用 TSSEnrichment() 函数为每个细胞计算这一指标，结果存储在元数据的 TSS.enrichment 列中； peaks区域内的片段总数 ：用于衡量细胞测序深度/复杂度的指标。测序深度较低的细胞可能由于read数过少而需要被排除。而具有极高片段数的细胞可能代表双细胞、细胞核聚集或其他异常情况 ； peaks区域内的片段比例 ：表示所有片段中落入ATAC-seq的peaks区域的比例。具有低值（例如<15%-20%）的细胞通常代表低质量细胞或技术假象，应该被去除。 需要注意的是，这一比例可能对所使用的peaks集合较为敏感； 基因组黑名单区域的reads比例 ：ENCODE项目提供了一个黑名单区域列表，这些区域中的reads通常与假象信号相关联。 如果细胞中映射到这些区域的reads比例较高（与映射到peaks区域的reads相比），则这些细胞通常代表技术假象，应该被去除。Signac软件包中包含了人类（hg19和hg38）、小鼠（mm9和mm10）、果蝇（dm3和dm6）以及秀丽隐杆线虫（ce10和ce11）的ENCODE黑名单区域。可以使用 FractionCountsInRegion() 函数来计算每个细胞在给定区域集合内的所有计数的比例。我们可以利用这个函数和黑名单区域来计算每个细胞的黑名单计数比例。 这三个指标： peaks区域内的片段总数 ， peaks区域内的片段比例 ， 基因组黑名单区域的reads比例 在CellRanger 的标准输出结果中：
设置 quantiles=TRUE ，可以快速的看出合适的过滤cutoff。
## 2.计算QC指标 # compute nucleosome signal score per cell pbmc <- NucleosomeSignal(object = pbmc) # compute TSS enrichment score per cell pbmc <- TSSEnrichment(object = pbmc) # add fraction of reads in peaks pbmc $pct_reads_in_peaks <- pbmc $peak_region_fragments / pbmc $passed_filters * 100 # add blacklist ratio pbmc $blacklist_ratio <- FractionCountsInRegion( object = pbmc, assay = 'peaks' , regions = blacklist_hg38_unified ) head(pbmc@meta.data) colnames(pbmc@meta.data) # 可视化 DensityScatter(pbmc, x = 'nCount_peaks' , y = 'TSS.enrichment' , log_x = TRUE, quantiles = TRUE) 我们还可以观察所有细胞的片段长度周期性，并根据核小体信号强度的高低对细胞进行分组。从前面的图表可以看出，单核小体/无核小体比例异常的细胞具有不同的核小体条带模式。其余细胞表现出典型的成功ATAC-seq实验的模式。
pbmc $nucleosome_group <- ifelse(pbmc $nucleosome_signal > 4, 'NS > 4' , 'NS < 4' ) FragmentHistogram(object = pbmc, group.by = 'nucleosome_group' ) 还可以分别使用小提琴图来绘制每个质量控制指标的分布：
VlnPlot( object = pbmc, features = c( 'nCount_peaks' , 'TSS.enrichment' , 'blacklist_ratio' , 'nucleosome_signal' , 'pct_reads_in_peaks' ), pt.size = 0.1, ncol = 5 ) 根据上面的数据分布，过滤条件如下，如果是自己的数据需要根据项目背景进行过滤cutoff调整：
pbmc <- subset( x = pbmc, subset = nCount_peaks > 9000 & nCount_peaks < 100000 & pct_reads_in_peaks > 40 & blacklist_ratio < 0.01 & nucleosome_signal < 4 & TSS.enrichment > 4 ) pbmc # An object of class Seurat # 165376 features across 9651 samples within 1 assay # Active assay: peaks (165376 features, 0 variable features) # 2 layers present: counts, data 还剩下9k多个细胞。
3、标准化和线性降维 分为三个小步骤：
标准化 ：Signac进行 TF-IDF 标准化 ，额还不知道是个什么标准化，mark一下后面再来学习这个细节特征选择 ：scATAC-seq数据的低动态范围使得我们难以像处理scRNA-seq数据那样进行可变特征选择。但可以选择仅使用排名top n%的特征（peaks）来进行降维，或者使用 FindTopFeatures() 函数移除出现在少于n个细胞中的特征 降维 ：在上述选择的特征（峰值）所构成的TF-IDF矩阵上运行奇异值分解（SVD） ## 3.标准化降维 pbmc <- RunTFIDF(pbmc) pbmc <- FindTopFeatures(pbmc, min.cutoff = 'q0' ) pbmc <- RunSVD(pbmc) DepthCor(pbmc) 第一个 LSI 成分通常捕获测序深度（技术变异）而非生物学变异 。如果情况如此，那么该成分应从下游分析中移除。可以使用 DepthCor() 函数来评估每个LSI成分与测序深度之间的相关性：
在这里，我们看到第一个 LSI 成分与细胞的总计数数量之间存在非常强的相关性（图中看到是就是前10个LSI，除了第一个基本上都在0.3以内）。我们将在后续步骤中不使用这个成分，因为我们不想根据细胞的总测序深度将它们分组，而是根据细胞类型特异性峰值的可及性模式对细胞进行分组。
4、非线性标准化和降维 现在细胞已经被嵌入到低维空间中，我们可以使用通常用于scRNA-seq数据分析的方法来进行基于图的聚类和非线性降维以用于可视化。 RunUMAP() 、 FindNeighbors() 和 FindClusters() 这些函数都来自Seurat包。
## 4.非线性标准化聚类 pbmc <- RunUMAP(object = pbmc, reduction = 'lsi' , dims = 2:30 ) pbmc <- FindNeighbors(object = pbmc, reduction = 'lsi' , dims = 2:30 ) pbmc <- FindClusters(object = pbmc, verbose = FALSE, algorithm = 3) DimPlot(object = pbmc, label = TRUE) + NoLegend() 聚类结果如下：
5、构建基因活性矩阵 UMAP图展示了人类血液中存在多个细胞群体。如果你熟悉PBMC的scRNA-seq分析，你可能会在scATAC-seq数据中识别出某些髓系和淋巴系细胞群体的存在。
然而，在scATAC-seq数据中注释和解释细胞簇更具挑战性，因为人们对非编码基因组区域的功能角色的了解远少于对编码蛋白的基因区域的了解 。
我们可以通过评估与基因相关的染色质可及性来量化基因组中每个基因的活性，并从scATAC-seq数据中创建一个新的基因活性检测方法。在这里，我们将使用一种简单的方法，即汇总与基因主体和启动子区域相交的片段（我们还建议探索Cicero工具，它可以实现类似的目标，我们在这里提供了一个展示如何在Signac工作流中运行Cicero的教程：
https:///signac/articles/cicero）。
为了创建基因活性矩阵，我们提取基因坐标，并将其扩展以包含上游2kb的区域（因为启动子的可及性通常与基因表达相关）。然后，我们使用 FeatureMatrix() 函数统计每个细胞映射到这些区域的片段数量。这些步骤可以通过 GeneActivity() 函数自动完成：
## 5.构建基因活性矩阵 gene.activities <- GeneActivity(pbmc) dim(gene.activities) gene.activities[1:5,1:5] # add the gene activity matrix to the Seurat object as a new assay and normalize it pbmc[[ 'RNA' ]] <- CreateAssayObject(counts = gene.activities) pbmc <- NormalizeData( object = pbmc, assay = 'RNA' , normalization.method = 'LogNormalize' , scale.factor = median(pbmc $nCount_RNA ) ) 现在可以可视化典型标记基因的活性，以帮助解释ATAC-seq聚类结果。以下是一些基因所对应的细胞类型：
B细胞
MS4A1（CD20） ：是B细胞的一个经典标记基因，主要在成熟的B细胞上表达。 T细胞
CD3D ：是T细胞表面的CD3复合物的一部分，主要在T细胞上表达。 LEF1 ：虽然LEF1在多种细胞类型中都有表达，但在T细胞中也有较高表达，且与T细胞的某些功能相关。 NK细胞
NKG7 ：是NK细胞的一个标记基因，主要在NK细胞中表达。 单核细胞/巨噬细胞
TREM1 ：主要在单核细胞/巨噬细胞等髓系细胞上表达。 LYZ ：编码溶菌酶，是单核细胞/巨噬细胞等髓系细胞的一个标记基因 需要注意的是，这些活性的测量值会比单细胞RNA-seq的测量值更加嘈杂。 这是因为它们代表了稀疏染色质数据的测量，并且它们假定基因主体/启动子的可及性与基因表达之间存在一般对应关系，而这种关系并不总是成立。尽管如此，我们仍然可以根据这些基因活性特征初步区分单核细胞、B细胞、T细胞和NK细胞的群体。然而，仅依靠监督分析来进一步细分这些细胞类型是具有挑战性的。
# 可视化 DefaultAssay(pbmc) <- 'RNA' FeaturePlot( object = pbmc, features = c( 'MS4A1' , 'CD3D' , 'LEF1' , 'NKG7' , 'TREM1' , 'LYZ' ), pt.size = 0.1, max.cutoff = 'q95' , ncol = 3 ) 6、与单细胞数据进行整合 为了帮助解释scATAC-seq数据，我们可以基于同一组织类型（人类PBMC）的scRNA-seq实验来对细胞进行分类。
我们利用跨模态整合和标签转移的方法：https:///10.1016/j.cell.2019.05.031，旨在识别基因活性矩阵和scRNA-seq数据集中共享的相关性模式，以识别两种模态之间的匹配生物学状态。这一过程将为每个细胞返回每个基于单细胞RNA-seq定义的聚类标签的分类得分。
单细胞数据为10x Genomics提供的预处理过的人类PBM scRNA-seq数据集，下载：https://signac-objects.s3./pbmc_10k_v3.rds
## 6.与单细胞数据整合 # Load the pre-processed scRNA-seq data for PBMCs pbmc_rna <- readRDS( "PBMCs/pbmc_10k_v3.rds" ) pbmc_rna <- UpdateSeuratObject(pbmc_rna) DimPlot(object = pbmc_rna, label = TRUE) + NoLegend() 已经注释好的单细胞：
找到 anchors：
transfer.anchors <- FindTransferAnchors( reference = pbmc_rna, query = pbmc, reduction = 'cca' ) transfer.anchors predicted.labels <- TransferData( anchorset = transfer.anchors, refdata = pbmc_rna $celltype , weight.reduction = pbmc[[ 'lsi' ]], dims = 2:30 ) predicted.labels 的内容如下：
可视化预测结果：
# 添加到对象中 pbmc <- AddMetaData(object = pbmc, metadata = predicted.labels) plot1 <- DimPlot( object = pbmc_rna, group.by = 'celltype' , label = TRUE, repel = TRUE) + NoLegend() + ggtitle( 'scRNA-seq' ) plot2 <- DimPlot( object = pbmc, group.by = 'predicted.id' , label = TRUE, repel = TRUE) + NoLegend() + ggtitle( 'scATAC-seq' ) plot1 + plot2 现在这个图就是很多文献中scRNA-seq+scATAC-seq联合分析的结果啦：
可以看到，基于scRNA-seq的分类结果与使用scATAC-seq数据计算得到的UMAP可视化结果是一致的。然而，需要注意的是，在scATAC-seq数据集中，有一小部分细胞被预测为血小板。这是出乎意料的，因为血小板是无核的，本不应被scATAC-seq检测到。 有可能被预测为血小板的细胞实际上是血小板的前体细胞——巨核细胞，它们主要存在于骨髓中，但在健康患者的外周血液中很少见，例如这些PBMC所来自的个体。鉴于巨核细胞在正常骨髓中已经极为罕见（< 0.1%），这种可能性似乎不大。
那是为什么？ 我猜测是因为scRNA-seq中存在小部分的血小板，而scATAC-seq中的细胞类型预测标签都来自于scRNA-seq。
比如我们前面学习的github文献： Github带有全套代码分享的文献复现2025
探索一下：platelets 绘制被分配到每个标签的细胞的预测得分，结果显示“血小板”细胞的得分相对较低（<0.8），这表明对其分配的细胞身份的可信度较低。在大多数情况下，这些细胞最可能的细胞身份是“CD4 naive”。
VlnPlot(pbmc, 'prediction.score.max' , group.by = 'predicted.id' ) 提取每种细胞的预测打分：
# Identify the metadata columns that start with "prediction.score." metadata_attributes <- colnames(pbmc[[]]) prediction_score_attributes <- grep( "^prediction.score." , metadata_attributes, value = TRUE) prediction_score_attributes <- setdiff(prediction_score_attributes, "prediction.score.max" ) prediction_score_attributes # [1] "prediction.score.CD4.Memory" "prediction.score.CD14..Monocytes" "prediction.score.NK.dim" # [4] "prediction.score.pre.B.cell" "prediction.score.NK.bright" "prediction.score.CD4.Naive" # [7] "prediction.score.CD8.Naive" "prediction.score.pDC" "prediction.score.Double.negative.T.cell" # [10] "prediction.score.CD16..Monocytes" "prediction.score.Platelet" "prediction.score.CD8.effector" # [13] "prediction.score.B.cell.progenitor" "prediction.score.Dendritic.cell" # Extract the prediction score attributes for these cells predicted_platelets <- which (pbmc $predicted .id == "Platelet" ) predicted_platelets platelet_scores <- pbmc[[]][predicted_platelets, prediction_score_attributes] platelet_scores # Order the columns by their average values in descending order ordered_columns <- names(sort(colMeans(platelet_scores, na.rm = TRUE), decreasing = TRUE)) ordered_platelet_scores_df <- platelet_scores[, ordered_columns] print (ordered_platelet_scores_df) 由于被分类为“血小板”的细胞数量极少（少于20个），很难确定它们的确切细胞身份。需要更大的数据集才能有信心识别这一细胞群体的特定峰值，并进行进一步分析以正确注释它们。因此，在下游分析中，我们将移除这些被预测的极为罕见的细胞状态，仅保留总细胞数超过20个的细胞注释。
# 移除血小板 predicted_id_counts <- table(pbmc $predicted .id) predicted_id_counts # Identify the predicted.id values that have more than 20 cells major_predicted_ids <- names(predicted_id_counts[predicted_id_counts > 20]) pbmc <- pbmc[, pbmc $predicted .id % in % major_predicted_ids] # change cell identities to the per-cell predicted labels Idents(pbmc) <- pbmc $predicted .id 7、细胞类型差异peaks分析 为了找到细胞亚群之间的差异可及性区域，我们可以进行差异可及性（DA）检验。一种简单的方法是执行Wilcoxon秩和检验，而 presto 包已经实现了一种极其快速的Wilcoxon检验，可以直接在Seurat对象上运行。
筛选 Naive CD4 cells and CD14 monocytes 之间的差异peaks：
## 7.差异peaks分析 # remotes::install_github('immunogenomics/presto') # change back to working with peaks instead of gene activities DefaultAssay(pbmc) <- 'peaks' # wilcox is the default option for test.use da_peaks <- FindMarkers( object = pbmc, ident.1 = "CD4 Naive" , ident.2 = "CD14+ Monocytes" , test.use = 'wilcox' , min.pct = 0.1 ) head(da_peaks) 选一个peaks可视化：
# 可视化差异peaks rownames(da_peaks)[1] plot1 <- VlnPlot( object = pbmc, features = rownames(da_peaks)[1], pt.size = 0.1, idents = c( "CD4 Naive" , "CD14+ Monocytes" ) ) plot2 <- FeaturePlot( object = pbmc, features = rownames(da_peaks)[1], pt.size = 0.1 ) plot1 | plot2 对peaks进行注释，与基因关联：
# 注释 open_cd4naive <- rownames(da_peaks[da_peaks $avg_log2FC > 3, ]) open_cd14mono <- rownames(da_peaks[da_peaks $avg_log2FC < -3, ]) closest_genes_cd4naive <- ClosestFeature(pbmc, regions = open_cd4naive) closest_genes_cd14mono <- ClosestFeature(pbmc, regions = open_cd14mono) head(closest_genes_cd4naive) 8、绘制 genomic regions 我们可以使用 CoveragePlot() 函数，针对pbmc中存储的任何基因组区域，按cluster、细胞类型或其他 metadata 对细胞进行分组，绘制Tn5整合在整个基因组区域的频率。这些表示 pseudo-bulk可及性轨迹，其中来自一个组内所有细胞的信号被平均在一起，以可视化一个区域内的DNA可及性（感谢Andrew Hill在他的优秀博客文章中为这个函数提供了灵感）。除了这些可及性轨迹外，我们还可以可视化其他重要信息，包括基因注释、峰值坐标和基因组链接（如果它们存储在对象中）。更多信息请参阅可视化教程：https:///signac/articles/visualization。
为了绘图的需要，将相关细胞类型放在一起会更好。我们可以通过运行 SortIdents() 函数，根据注释细胞类型之间的相似性自动对绘图顺序进行排序：
# 可视化 pbmc <- SortIdents(pbmc) # find DA peaks overlapping gene of interest regions_highlight <- subsetByOverlaps(StringToGRanges(open_cd4naive), LookupGeneCoords(pbmc, "CD4" )) CoveragePlot( object = pbmc, region = "CD4" , region.highlight = regions_highlight, extend.upstream = 1000, extend.downstream = 1000 ) 文末友情宣传 
