1. 从单细胞染色质数据到转录因子调控解析做单细胞ATAC-seq分析的朋友们都知道找到差异开放区域只是第一步真正有意思的是理解这些区域背后的调控机制。最近我在分析小鼠大脑神经元数据时发现Pvalb和Sst两种抑制性神经元之间存在大量染色质可及性差异这让我特别好奇到底是哪些转录因子在调控这些差异Signac这个R包真是帮了大忙。它整合了JASPAR数据库的motif信息可以直接在单细胞染色质数据中寻找富集的转录因子结合位点。我花了三周时间反复测试这个流程现在把最实用的操作方法和避坑指南分享给大家。先说最重要的两点经验第一一定要在集群上跑chromVAR计算本地机器很容易内存溢出第二差异peak分析时记得调低min.pct参数scATAC-seq数据比scRNA-seq稀疏得多。2. 环境配置与数据准备2.1 软件包安装的注意事项第一次安装Signac相关工具链时我被依赖包搞得头大。后来发现用Bioconductor管理最省事。下面这个安装组合是我测试过最稳定的版本if (!requireNamespace(BiocManager, quietly TRUE)) install.packages(BiocManager) BiocManager::install(c(JASPAR2020, TFBSTools, BSgenome.Mmusculus.UCSC.mm10, motifmatchr, chromVAR)) install.packages(c(Signac, Seurat, ggseqlogo, patchwork))特别提醒两点1) BSgenome要选对参考基因组版本小鼠用mm10人用hg382) JASPAR2020包含脊椎动物的核心motif集如果研究其他物种需要换数据库。我在测试时发现用最新版JASPAR2022反而会出现兼容性问题所以建议先用2020版。2.2 数据加载与检查官方提供的小鼠大脑数据集是个很好的起点。加载数据后一定要检查对象结构library(Signac) brain - readRDS(adult_mouse_brain.rds) brain # 输出应显示两个assaypeaks和RNA # 确认有LSI降维和UMAP坐标这个数据集包含3517个细胞分为15个亚群。我重点关注的是Pvalb和Sst神经元它们在UMAP图上分别位于cluster 6和7。如果用自己的数据记得先用RunTFIDF和RunSVD做好基础分析。3. Motif信息整合与富集分析3.1 从JASPAR获取motif信息获取motif数据时tax_group参数特别重要。研究小鼠大脑时我用的是vertebrates脊椎动物但如果做肿瘤样本可能需要改为human。下面代码获取位置频率矩阵(PFM)pfm - getMatrixSet( x JASPAR2020, opts list(collection CORE, tax_group vertebrates, all_versions FALSE) )这里有个小技巧用listMotifs(pfm)可以查看所有motif名称。我通常会先过滤掉低质量的motif保留那些有明确转录因子注释的。3.2 将motif添加到Seurat对象AddMotifs函数需要指定参考基因组这是容易出错的地方brain - AddMotifs( object brain, genome BSgenome.Mmusculus.UCSC.mm10, pfm pfm )成功运行后可以用GetMotifData(brain)检查motif信息。我遇到过基因组版本不匹配报错这时需要确认BSgenome包和原始数据的参考基因组一致。4. 差异分析与motif富集4.1 寻找差异开放区域在Pvalb和Sst神经元间找差异peak时关键是要调整min.pct参数da_peaks - FindMarkers( object brain, ident.1 Pvalb, ident.2 Sst, only.pos TRUE, min.pct 0.05, # 比scRNA-seq的默认值0.1低 test.use LR, latent.vars nCount_peaks ) top_peaks - rownames(da_peaks[da_peaks$p_val 0.005, ])实测发现用LR检验比wilcox更适合scATAC-seq数据。latent.vars参数控制技术噪音我一般用nCount_peaks或nucleosome_signal。4.2 Motif富集分析实战FindMotifs会自动匹配GC含量相似的背景区域这是它的聪明之处enriched_motifs - FindMotifs( object brain, features top_peaks ) head(enriched_motifs[order(enriched_motifs$p.adjust), ])输出表格中fold.enrichment 2且p.adjust 0.05的motif值得关注。我发现在神经元数据中NeuroD1、Mef2c等神经发育相关TF经常富集。4.3 Motif可视化技巧用ggseqlogo可以自定义motif图library(ggseqlogo) motif_plot - MotifPlot( object brain, motifs c(MA0497.1, MA0501.1) ) theme_classic()如果想保存矢量图建议用ggsaveggsave(motif_plot.pdf, motif_plot, width 8, height 4)5. 单细胞水平的motif活性分析5.1 chromVAR计算注意事项RunChromVAR是计算量最大的步骤我的经验是一定要在集群上运行至少80G内存可以先取子集测试参数用future包实现并行加速library(future) plan(multicore, workers 8) brain - RunChromVAR( object brain, genome BSgenome.Mmusculus.UCSC.mm10 )5.2 活性可视化方法FeaturePlot展示motif活性时调整cutoff很重要DefaultAssay(brain) - chromvar FeaturePlot(brain, features MA0497.1, min.cutoff q10, max.cutoff q90, pt.size 0.5)我发现Pvalb神经元中Mef2c活性明显高于Sst神经元这与之前的富集结果相互印证。5.3 差异活性分析FindMarkers也可以用于比较motif活性diff_activity - FindMarkers( object brain, ident.1 Pvalb, ident.2 Sst, assay chromvar, mean.fxn rowMeans )分析结果显示Pvalb中Mef2c活性是Sst的2.3倍p1e-6这可能是两种神经元功能差异的关键调控因子。6. 分析结果解读与验证在实际项目中我会用三种方式验证motif分析结果与已发表的ChIP-seq数据交叉验证检查差异motif对应的TF基因是否在scRNA-seq中差异表达用GREAT工具分析motif所在peak的基因功能例如在小脑数据中发现Olig2 motif在少突胶质细胞前体细胞中富集这与已知的发育生物学知识一致。这种多组学交叉验证能大大提高结果的可信度。最后提醒大家motif分析只是提出假设的工具真正的机制还需要实验验证。我通常会把排名前5的motif作为候选设计后续的CRISPR干扰实验。保存分析结果时建议用write.csv保存富集表格用saveRDS保存整个Seurat对象方便后续复查。