1. 为什么单细胞表观基因组学也需要“主题建模”不少人第一次听到cisTopic这个名字第一反应是“又一个聚类工具”。它确实经常被归类在聚类教程里但它做的事情和常规的k-means、层次聚类很不一样。cisTopic的全称是“cis-regulatory topic”核心思路是把每个单细胞当成一篇文档、把调控区域当成单词用概率主题模型LDA家族去抽取出若干个共调控的调控主题。这样得到的不是简单的细胞分群而是“细胞状态”和“调控程序”之间的概率关联换句话说你不仅知道这些细胞是一类还能知道驱动它们成为一类的顺式调控线索是什么。这种思路很适合单细胞表观基因组学数据。以scATAC-seq为例单细胞分辨率下每个细胞只贡献很少的可用片段所以细胞-峰矩阵的稀疏度经常超过95%甚至超过99%。这种矩阵如果直接拿去做欧氏距离聚类结果往往被几个高覆盖峰主导而主题模型天然是为稀疏离散数据设计的它不假设每个细胞只能属于一个类而是允许一个细胞同时包含多个主题的活性。这种特性特别契合真实生物学一个细胞里往往同时活跃着基础调控程序、细胞周期程序、谱系特异性程序硬聚类一刀切很容易丢掉这种叠加信息。所以这篇教程会沿着cisTopic的完整使用链路走一遍从输入矩阵的准备到模型运行、主题数选择再到细胞聚类、区域得分和转录因子富集最后把cisTopic和其他常见聚类路线放在一起做对比。适合已经跑通过scATAC-seq基础流程、准备往调控机制方向深挖的读者如果你只想要一个快速的细胞分群不一定需要上cisTopic但如果你想找“哪些调控区域共同定义了某个细胞状态”那这部分功夫就是值得的。1.1 单细胞表观数据到底是什么先说数据形态。不管你是用10x Genomics的Cell Ranger ATAC输出还是用ArchR、SnapATAC2处理过的结果最终进入cisTopic前通常要整理成一张细胞-峰矩阵行是基因组上的候选顺式调控区域列是细胞矩阵里的值一般取0/1或归一化后的“该细胞在该区域检测到可访问片段”的次数。这张矩阵的规模很夸张一个典型的人类样本动辄几十万个峰、数万个细胞。但真正非零的条目很少因为单细胞ATAC本身的覆盖深度有限一个细胞只能覆盖几百到几千个可访问位点。高稀疏度意味着很多在批量ATAC-seq里显著的信号在单细胞里会被稀释得非常厉害。这也是为什么直接用普通聚类算法处理之前通常要先降维、先做某种隐变量提取cisTopic的思路更彻底一点它不降维后再聚类而是在原始二值矩阵上直接做概率分解。另外一个值得注意的点是单细胞表观矩阵不像转录组的UMI计数那样有比较厚的动态范围它更像是“开/关”二值信号。这种离散二项式结构的噪声模型和主题模型里的词袋假设非常接近。换句话说把每个细胞当成一篇文章、把每个峰当成一个词不是一个比喻而是数学模型层面的一种天然适配。1.2 为什么k-means这类聚类算法在这里不够用很多聚类教程都会用k-means作为入门算法而且也确实有人在scATAC-seq的LSI降维结果上跑k-means能得到相当不错的细胞分群。但k-means有一个核心限制它假设每个样本可以被划分到唯一一个簇样本之间是互斥关系。可生物学上一个细胞往往同时处在多个调控程序的“叠加态”里。举个例子一个处于G1/S期的T细胞它既延续T细胞谱系调控程序又启动细胞周期相关程序硬要让它在T细胞簇和增殖细胞簇之间二选一多少会丢失信息。层次聚类的问题类似它给出的是树状划分虽然可以看到不同层级的子结构但本质上还是硬聚类而且在大规模数据上计算距离矩阵的代价很高。欧氏聚类对稀疏高维二值矩阵尤为不适零膨胀会让距离度量失真一个峰覆盖度差异很大的细胞反而可能被分到一起。如果你想了解“哪些区域协同变化、共同定义一种细胞状态”上述这些聚类方法都无法直接回答。cisTopic则不同它用概率主题模型把每个细胞拆解成若干个主题的混合。每个主题对应一组峰这些峰倾向于在某些细胞里共同开放。最终你能拿到两类概率细胞-主题概率每个细胞在多大程度上依赖某个主题和主题-峰概率每个峰在不同主题里的权重。前者可以做传统聚类和可视化后者则可以直接用来定义顺式调控程序往下对接motif分析、转录因子富集。所以可以把cisTopic理解成“聚类的上游工具”它先把数据拆成可解释的调控模块模块层面的细胞分组远比你直接基于原始峰矩阵聚类要稳定得多。2. 模型直觉从LDA到顺式调控主题本质上是“词袋”用在基因组上cisTopic的核心模型脱胎于自然语言处理里的LDALatent Dirichlet Allocation。在LDA里一篇文档是单词的集合文档会以不同比例属于若干个主题而每个主题又有一个在所有单词上的概率分布。cisTopic把整套概念平移到了基因组上每个细胞是一篇“文档”每个开放的峰是一个“词”每个主题就是一个“顺式调控主题”。模型要做的就是反推出细胞里隐藏的主题混合比例以及每个主题下峰的贡献权重。这种“文档-主题-词”的三层结构可以用一个很直白的例子理解。假如你有十篇小说小说里反复出现“大海”“航船”“海盗”那么模型大概率会抽出一个“航海主题”另一批小说如果频繁出现“城堡”“骑士”“王后”就会抽出“宫廷主题”。有些小说二者都有那它在这两个主题上的得分都会偏高。细胞也一样一个细胞里如果同时高权重命中造血干细胞的调控区域和细胞周期调控区域那它就不是“单标签”的它在这两个主题上都会有可观的活性。2.1 细胞-主题矩阵和主题-峰矩阵cisTopic输出中有两个关键矩阵理解它们你才能顺利做下游分析。第一个是细胞-主题矩阵。假设你有10000个细胞、设置了8个主题那这个矩阵就是10000行、8列每一行代表一个细胞在8个主题上的活性分布。因为每行加起来等于1所以它本质上是个概率分布。你拿这8个维度去做UMAP、t-SNE或者直接算相关系数做层次聚类都会比在原始高维峰矩阵上运行稳定得多。第二个是主题-峰矩阵。假设有500000个峰、8个主题那这个矩阵就是8行、500000列每个主题给出了它在所有峰上的权重分布。权重高的峰就是该主题的“代表词”。把这些高权重峰映射回基因组坐标就能看到该主题对应的调控区域在哪里集中把这些区域拿来跑motif富集就能推测哪些转录因子可能主导这个调控程序。这两个矩阵是cisTopic一切下游分析的基础。所谓“顺式调控主题的概率建模”最终产出其实就是这两个概率分布。2.2 cisTopic和LDA的差异点cisTopic不是把LDA原封不动搬过来它在几处做了针对基因组数据的调整。比如它会考虑峰与峰之间的基因组距离和染色质可及性的分布特征而不是完全忽略词序另一个常见调整是它会用二项式/多项式的框架去适配0/1稀疏矩阵的特性让模型对极低覆盖度细胞更稳健。另外主题数量是用户给定的。这一点和k-means需要指定簇数一样没有免费的午餐。cisTopic给了一套基于模型对数似然和主题稳定性的选择机制后面第4节我会专门讲如何选主题数这也是整个流程里最容易翻车的环节。在实际实现上cisTopic通常采用折叠吉布斯采样collapsed Gibbs sampling来估计后验分布。你不需要手动实现采样R包封装好了但理解它的意义能帮你设置迭代次数采样需要足够多的迭代让马尔可夫链收敛否则不同主题之间会出现大量“复制粘贴式”的相似结果。3. 输入数据构造把散落片段变成细胞-峰矩阵的基本功cisTopic是个R包建议在R/RStudio环境里用。安装时依赖比较多如果你之前装过Seurat、Signac、ArchR等包通常不容易冲突因为cisTopic的核心依赖是Matrix、data.table、topicmodels这些老牌R包。建议直接用Bioconductor安装安装命令如下if (!requireNamespace(BiocManager, quietly TRUE)) install.packages(BiocManager) BiocManager::install(cisTopic)安装完先不要急着跑模型先确认一下版本。cisTopic在不同版本里的函数名有细微调整比如有些版本用runModels有些版本更推荐runCGSModels作为并行入口。我的经验是跑任何函数前先看一眼帮助文档library(cisTopic) ?runModels3.1 从fragments.tsv.gz构建细胞-峰矩阵如果你是10x平台的数据最原始的输入是fragments.tsv.gz和对应的singlecell.csv里面记录每个片段来自哪个细胞、在哪条染色体哪段区间。cisTopic直接接收的通常不是这种片段文件而是已经聚合好的峰矩阵所以你需要先完成两步转换获得全样本共享的峰集合。最简单的办法是从Cell Ranger的peaks.bed读取或者用MACS2在合并后的伪批量数据上重新call peak。统计每个细胞落在每个峰里的片段数生成稀疏矩阵。这个过程可以用Signac的FeatureMatrix()也可以直接用ArchR里的getMatrixFromProject(ArchRProj, useMatrix PeakMatrix)或者用SnapATAC2的snapToFragments。矩阵准备好后一般保存成RDS或者Matrix Market格式。如果你手头只有10x的filtered_peak_bc_matrix.h5可以先在R里读出来library(Matrix) library(Seurat) mat - Read10X_h5(filtered_peak_bc_matrix.h5)这里的mat是一个稀疏矩阵行名是“chr1-10000-20000”这种坐标列名是细胞Barcode。这个格式已经可以直接进cisTopic。3.2 使用createcisTopicObject构建对象拿到稀疏矩阵后创建cisTopic对象很简单ct - createcisTopicObject( count.matrix mat, project.name my_scatac_project )需要注意这个函数默认把输入当作“cell by peak”格式处理但在构造对象时有些版本要求行是峰、列是细胞所以如果发现细胞数和峰数反了先转置矩阵再传入。创建对象时还可以传入细胞元数据比如样本来源、批次、预期细胞类型注释。方法是在createcisTopicObject里使用cell.meta参数或者后续用addCellMetadata()补充。这一步不是必须的但我强烈建议一开始就带上样本ID否则等你跑完模型发现主题主要按批次分群时才想起来没有样本信息可以对照会很被动。数据质量过滤也建议在进入cisTopic前完成。常见的过滤标准包括每个细胞至少要有500个唯一片段、TSS富集分数大于某个阈值、线粒体片段比例要低。不同平台标准不一不必教条但一定要记录过滤前后细胞数量变化不然别人问起数据量你连基线都说不清。3.3 染色体命名和坐标格式是隐藏陷阱一个特别容易被忽略的坑是染色体命名。你的峰坐标可能来自hg19或hg38但如果你混用了不同版本或者某些区域名写成“1:10000-20000”、某些写成“chr1:10000-20000”cisTopic在后续计算主题得分和重叠区域时会报错或者更麻烦的是悄悄漏掉一部分区域导致结果偏差。我在处理公共数据时遇到过这种情况一个公共矩阵用了NCBI的1_10000_20000命名另一个用了UCSC的chr1-10000-20000命名合并后如果不统一后续所有分析都是错的。所以进模型前建议用gsub清洗一遍坐标命名至少要保证全部统一为chr前缀格式并用GenomicRanges验证一下排序。4. 核心运行与主题数选择跑模型只是开始难的是决定几个主题关键代码无法保证准确的函数名。用通用描述。构建好对象后运行模型大概是这样的set.seed(123) ct - runModels( ct, topic c(2, 4, 6, 8, 10, 12), seed 123, burnin 125, iterations 250, addProps 0.05 )参数说明topic是你要尝试的主题数量集合可以直接传一组整数burnin是吉布斯采样的预热期iterations是采样的迭代次数。我一般会把迭代次数设置到500以上虽然这会明显变慢但收敛后的主题可解释性会好很多。如果时间紧张可以先跑250次看看趋势再用小迭代跑一组、大迭代跑一组做对比。cisTopic会为topic里每个候选值单独训练一个模型保留在对象的models槽里。之后用selectModel挑选最终模型ct - selectModel(ct, select c(6, 8, 10)) modelSelPlot(ct)selectModel会基于对数似然和稳定性综合排序选出你指定候选里相对合适的那个。modelSelPlot会画出一条主题数和对数似然的关系曲线你可以像看PCA碎石图一样找出“拐点”。4.1 不要只看数值还要看主题是否可解释主题选择是cisTopic流程里最需要经验的地方。很多教程会让你直接看modelSelPlot但这个图只能告诉你哪个主题数在统计上更合理不能告诉你哪个主题数在生物学上更有用。我个人的筛选流程有三个步骤。第一步看稳定性。同一主题数下用不同随机种子各跑一遍然后计算主题之间的一致性。如果两次运行得到的主题内容高度相似说明这个主题数比较稳如果两次运行的主题五花八门说明模型没有收敛或者主题数设置得不自然。第二步看主题分离度。把细胞-主题矩阵画成热图检查是不是有几个主题在所有细胞上都是均匀分布。如果出现“万金油主题”几乎每个细胞都有很高的占比那可能是主题数太多模型把一个公共调控程序硬拆成了多个碎片。第三步也是最关键的把主题映射回生物学。选主题数不比选k值不是为了数学最优而是为了“每一个主题都能讲出一个故事”。当你设置8个主题时理想状态应该是每个主题对应一个明确的细胞类型或者一个明确的生物过程比如CD4 T细胞主题、CD8 T细胞主题、单核细胞主题、细胞周期主题等。如果你发现两个主题几乎只差在一两个峰上那你的主题数可能偏多如果某个主题把好几类细胞的信号混在一起那可能偏少。4.2 多模型对比怎么看runModels会为每个主题数各训练一个模型selectModel会帮你选一个但我建议不要只依赖自动选择。我会把所有候选模型的细胞-主题矩阵都提取出来分别做一遍UMAP然后用已知的细胞类型注释去看聚类结果。有些时候主题数为6的对数似然更高但主题数为8才能把你关心的稀有细胞类型分开。对数似然是对整体数据的拟合度量它不会主动在乎你关心的那个稀有亚群。所以正确做法是用modelSelPlot缩小候选范围再用生物学注释做最终决策。这也是概率建模和纯机器学习最大的不同——下游解释才是检验标准。4.3 迭代次数和收敛性吉布斯采样类方法有个问题迭代次数不足时不同主题之间会出现“重复化”。比如你设了8个主题结果其中3个主题的峰权重分布高度相似差异只来自随机噪声。这种时候先别急着改主题数试着把迭代次数从200翻到1000很多“伪主题”会自己合并。我在真实数据上观察过低迭代次数下细胞周期相关峰会被拆到好几个主题里提高到足够迭代次数后它们会稳定归到一个主题。这让我养成了一个习惯——小规模调参时用快速设置一旦确定了主题数范围最后一定用长迭代完整跑一遍并保存对象副本。5. 下游实操细胞聚类、可视化、motif与转录因子富集模型跑完只是“烧完钱”真正产生生物学价值的从这一步开始。cisTopic的下游分析可以分成四条线细胞聚类、主题区域可视化、主题归一到基因组、转录因子富集。5.1 用细胞-主题矩阵做聚类和可视化先提取细胞-主题概率cell_topic - getCellData(ct, select c(8), dims 2)这里select参数选择你最终确定的那个模型。getCellData会返回一个包含降维结果的对象你要的细胞-主题矩阵在这个对象的cell.data里。通常我们再对主题概率做一次归一化或转换然后跑UMAPlibrary(umap) umap_out - umap(cell_topic$cell.data[, grep(Topic, colnames(cell_topic$cell.data))])如果你更喜欢用Seurat生态也可以把这个细胞-主题矩阵直接塞进CreateSeuratObject()然后走RunUMAP()和FindClusters()。很多团队会对比直接基于LSI的聚类和基于cisTopic主题的聚类发现后者往往在稀有细胞类型的分离上更有优势。另一个常见操作是基于主题概率做层次聚类。因为主题数通常很少比如8个、12个你完全可以把每个细胞看作一个8维或12维的概率向量计算它们的相关性距离然后画层次聚类树并用pheatmap展示。library(pheatmap) hm - pheatmap(cor(t(cell_topic$cell.data[, topic_cols])))这种做法的好处是你可以看到细胞在主题层面的树状结构层次聚类能告诉你哪些主题更接近哪些主题分化更早。把聚类结果和已知注释对照能快速判断主题是否有生物学意义。5.2 主题活性在基因组上的映射cisTopic的一个独特能力是给每个峰打分找出每个主题的“代表峰”。做法大致是region_scores - getRegionsScores(ct, select c(8), method NormTop, scale TRUE)method通常有“Z-score”或“NormTop”等选项。得到的region_scores里每个主题对应一组区域分数分数越高说明该区域对这个主题越重要。拿到这个分数后你可以做几件事筛选每个主题的高分区域导出成BED文件去IGV里看具体信号。把这些区域和基因注释、增强子数据库如ENCODE、ENdb重叠判断主题涉及哪些功能元件。把这些区域作为输入跑转录因子富集。比如你的主题8里高分区大量落在某个基因的启动子和近端增强子附近那这个主题就可能与该基因所在细胞类型的功能分化有关。把高分区做motif富集可以直接推测调控该主题的核心转录因子。5.3 转录因子富集怎么做cisTopic官方教程里通常会把主题区域导出给RcisTarget或者HOMER做motif富集。RcisTarget的优势是不需要把区域严格二值化它可以用区域排序来做富集更平滑。简单流程是# 假设已经安装并载入RcisTarget library(RcisTarget) motifRankings - importRanks(hg38_500bp_up_and_100bp_down_tss.mc9nr.feather)然后按主题的高分峰列表做addMotifAnnotation和addSignificantGenes。最终你会得到一张表哪个motif在某个主题的调控区域里显著富集以及对应的转录因子和可能的靶基因。这里有个实操经验不要只看p值最小的motif要看“motif富集区域是否真的是该主题特有的”。比如有些基础转录因子如p300、CTCF结合基序几乎在每个主题里都富集因为它们本来就广泛参与染色质结构维持。更有信息量的是“差异化富集”也就是某个motif在主题A里显著富集但在其他主题里不富集。我建议把每个主题的高分区域批量跑一遍富集然后做对比矩阵而不是单独看一个主题。5.4 主题和差异可及性分析的结合另一个常见下游是差异可及性分析。你可以在UMAP聚类后把不同细胞簇的cisTopic主题活性做差异比较找出某个簇特异激活的主题。这种“先聚类再映射到主题”的路线比直接跑差异峰要稳因为主题已经把无数单峰信号压缩成了几个可解释的调控模块。例如你发现了两个与细胞周期相关的主题一个在G1/S期细胞里高一个在G2/M期细胞里高。两边高分区对应的motif可能完全不同这可以为你提供非常有针对性的转录因子线索。这类结果放到论文里比单纯给出一堆差异peak要有说服力得多。6. 和其他聚类路线对照k-means、层次聚类、LSI之后到底该用谁因为cisTopic经常被放在“聚类教程”这个检索场景里我觉得有必要把它和k-means、层次聚类、子空间聚类这些常见路线放在一起做一次明确对比免得大家盲目套用。方法输入聚类类型主要输出适用场景k-means低维向量如LSI降维结果硬聚类细胞簇标签快速分群、大型数据集层次聚类低维向量或距离矩阵硬聚类/树状图树、聚类标签小样本、想看层级关系子空间聚类高维矩阵软/硬子空间划分高维信号集中在局部维度cisTopic原始稀疏峰矩阵软主题聚类主题分布、峰权重找调控程序、机制解释这张表不是想说cisTopic比k-means高级而是想强调它们回答的问题不一样。k-means回答的是“每个细胞该贴哪个标签”cisTopic回答的是“这个细胞的状态由哪些调控程序叠加而成”。在很多分析流水线里两者可以串联使用先用cisTopic提取主题再用k-means对主题概率向量聚类得到细胞簇然后用层次聚类展示主题层级。如果你已经有Seurat/Signac的完整流程跑过LSI降维和FindClusters再去叠加cisTopic会有点重复但不冲突。我通常的做法是Signac先给我一个快速可靠的细胞注释然后我用注释去校验cisTopic主题的可解释性而不是反过来让cisTopic来做最基础的细胞类型鉴定。思路是细胞注释交给成熟的快速流程调控机制挖掘交给cisTopic。6.1 什么时候可以不上cisTopic如果你的目标只是“把细胞分出来”而且你的数据量特别大、时间很紧直接用LSI加k-means完全够用没必要上主题模型。主题模型训练慢调参空间大解释起来也费劲用它回答一个简单分群问题属于高射炮打蚊子。如果你的数据里有非常明显的结构比如几个大细胞类型各占一大块LSI聚类效果通常不会差。但如果你关心的是细胞类型内部的连续状态比如分化轨迹上的中间态或者同一个细胞类型里不同调控程序的协同变化k-means硬聚类就有点不够看了。这时候cisTopic的软主题结构就体现出价值它不会强迫中间态细胞选边站而是给出一组“高不成低不就”的主题组合恰恰是这种组合量化了过渡状态。6.2 关于“子空间聚类概述及python完整实现”热词的一点延伸有朋友搜索“子空间聚类”这类方法在高维数据里寻找低维子空间然后再聚类与主题建模的哲学有一些相通之处。但子空间聚类对混合模型假设有很强的依赖应用于单细胞表观组学这种极度稀疏数据时需要对缺失机制做额外假设。cisTopic走的是更直接的生成模型路线它通过对“峰-主题-细胞”的联合概率建模让缺失和噪声都统一在概率框架下处理所以应用起来更顺手。如果你本身就是Python栈不想切到R也可以考虑用gensim或scikit-learn里的NMF来做类似的主题分解。NMF在非负稀疏矩阵上也能提取“主题”式的因子结构但没有Dirichlet先验结果的可解释性和稳定性通常没有cisTopic好。如果你用Python做最常见的路线是先做TF-IDF或LSI再用KMeans聚类注意这里因为二值数据没有“词频”概念TF-IDF要改造成二值加权形式否则效果会偏移。6.3 聚类热图、趋势图、富集条目的串联还有一种常见需求是想把聚类结果跟富集条目串成一张“展示图”。比如你已经把细胞分成了几簇每个簇有一组标记峰想看看这些峰的motif富集条目和趋势。一般流程是用主题的高分峰取并集或交集画一个主题×峰的聚类热图再在旁边标注富集到的转录因子条目。cisTopic已经帮你把“峰降维成主题”了所以这一步做起来比直接从几万个峰里挑特征要干净很多。7. 我踩过的坑与调参经验用cisTopic这两年我踩过的坑不少。很多坑不是模型本身的问题而是使用习惯问题。挑几个印象最深的分享一下。7.1 主题太多反而丢失稀有细胞类型第一次跑正式数据时我参考某些教程设了20个主题结果日志似然很好看但稀有细胞类型对应的主题被切割得七零八落几乎没有一个是干净的。后来把主题数降到10左右稀有亚群的信号才集中起来。这不是说10比20好而是说“主题越多每个主题的平均丰度越低”稀有信号更容易被拆碎。如果你特别关注稀有细胞建议主题数取偏小的区间比如5到10然后配合已知marker基因去验证。7.2 批次效应会伪装成主题多批次数据尤其容易出这个问题。如果主题A在批次1的细胞里高、主题B在批次2的细胞里高但两个批次注明的细胞类型完全相同那要么是数据质量差异太大要么是批次效应被模型捕捉了。遇到这种情况先别急着下生物学结论。可以试两个方向一是检查峰矩阵是否做过批次间的归一化比如用chromVAR的偏差校正或者简单的分位数归一化二是在聚类前用Harmony或Seurat的整合功能在主题概率矩阵上做批次整合。主题概率矩阵因为维度很低整合速度也很快。7.3 染色体坐标版本混用前面提过一次但值得再强调。hg19和hg38坐标混用反映在结果里很可能就是一个主题的区域几乎全部来自某个坐标版本。这个问题排查起来很隐蔽因为模型不会报错只有等下游重叠数据库时才会发现大量区域重叠不上。建议在创建对象之前用GenomicRanges::seqinfo()统一检查一下染色体长度和命名并保留一份坐标版本记录。7.4 内存和时间优化几十万峰、几万细胞的数据跑多主题数、多迭代的模型确实会慢。我一般在正式跑之前先做一个“下采样版”随机抽5000个细胞、50000个峰快速跑一遍确定主题数范围和迭代次数再用全量数据正式跑。这样能节省大量试错时间。另外虽然runModels支持多主题数同时跑但我还是建议分多次运行每次隔一段时间保存一个RDS。一旦某个主题数的模型因为内存问题崩溃你不需要全部重来。7.5 保存一切包括随机种子概率模型的随机性不可忽视。同一份数据同样参数两次运行的主题顺序可能完全不同甚至内容有一定波动。所以我在每次runModels前都固定种子并把种子写进项目元数据每个阶段都保存一个独立RDS文件比如ct_raw.rds、ct_model.rds、ct_selected.rds。这样即使后面分析叉了也可以回到任意一个节点不用从头跑。最后再分享一个小技巧选定主题后把你的主题编号、对应细胞类型注释、每个主题的代表峰个数、代表性转录因子做成一张总表。这张表既是和合作者讨论的基础也是后续写论文时的素材库。做分析的时候你可能觉得这些记录很冗余但几个月后回头看你会感谢当时的自己。
1. 为什么单细胞表观基因组学也需要“主题建模”不少人第一次听到cisTopic这个名字第一反应是“又一个聚类工具”。它确实经常被归类在聚类教程里但它做的事情和常规的k-means、层次聚类很不一样。cisTopic的全称是“cis-regulatory topic”核心思路是把每个单细胞当成一篇文档、把调控区域当成单词用概率主题模型LDA家族去抽取出若干个共调控的调控主题。这样得到的不是简单的细胞分群而是“细胞状态”和“调控程序”之间的概率关联换句话说你不仅知道这些细胞是一类还能知道驱动它们成为一类的顺式调控线索是什么。这种思路很适合单细胞表观基因组学数据。以scATAC-seq为例单细胞分辨率下每个细胞只贡献很少的可用片段所以细胞-峰矩阵的稀疏度经常超过95%甚至超过99%。这种矩阵如果直接拿去做欧氏距离聚类结果往往被几个高覆盖峰主导而主题模型天然是为稀疏离散数据设计的它不假设每个细胞只能属于一个类而是允许一个细胞同时包含多个主题的活性。这种特性特别契合真实生物学一个细胞里往往同时活跃着基础调控程序、细胞周期程序、谱系特异性程序硬聚类一刀切很容易丢掉这种叠加信息。所以这篇教程会沿着cisTopic的完整使用链路走一遍从输入矩阵的准备到模型运行、主题数选择再到细胞聚类、区域得分和转录因子富集最后把cisTopic和其他常见聚类路线放在一起做对比。适合已经跑通过scATAC-seq基础流程、准备往调控机制方向深挖的读者如果你只想要一个快速的细胞分群不一定需要上cisTopic但如果你想找“哪些调控区域共同定义了某个细胞状态”那这部分功夫就是值得的。1.1 单细胞表观数据到底是什么先说数据形态。不管你是用10x Genomics的Cell Ranger ATAC输出还是用ArchR、SnapATAC2处理过的结果最终进入cisTopic前通常要整理成一张细胞-峰矩阵行是基因组上的候选顺式调控区域列是细胞矩阵里的值一般取0/1或归一化后的“该细胞在该区域检测到可访问片段”的次数。这张矩阵的规模很夸张一个典型的人类样本动辄几十万个峰、数万个细胞。但真正非零的条目很少因为单细胞ATAC本身的覆盖深度有限一个细胞只能覆盖几百到几千个可访问位点。高稀疏度意味着很多在批量ATAC-seq里显著的信号在单细胞里会被稀释得非常厉害。这也是为什么直接用普通聚类算法处理之前通常要先降维、先做某种隐变量提取cisTopic的思路更彻底一点它不降维后再聚类而是在原始二值矩阵上直接做概率分解。另外一个值得注意的点是单细胞表观矩阵不像转录组的UMI计数那样有比较厚的动态范围它更像是“开/关”二值信号。这种离散二项式结构的噪声模型和主题模型里的词袋假设非常接近。换句话说把每个细胞当成一篇文章、把每个峰当成一个词不是一个比喻而是数学模型层面的一种天然适配。1.2 为什么k-means这类聚类算法在这里不够用很多聚类教程都会用k-means作为入门算法而且也确实有人在scATAC-seq的LSI降维结果上跑k-means能得到相当不错的细胞分群。但k-means有一个核心限制它假设每个样本可以被划分到唯一一个簇样本之间是互斥关系。可生物学上一个细胞往往同时处在多个调控程序的“叠加态”里。举个例子一个处于G1/S期的T细胞它既延续T细胞谱系调控程序又启动细胞周期相关程序硬要让它在T细胞簇和增殖细胞簇之间二选一多少会丢失信息。层次聚类的问题类似它给出的是树状划分虽然可以看到不同层级的子结构但本质上还是硬聚类而且在大规模数据上计算距离矩阵的代价很高。欧氏聚类对稀疏高维二值矩阵尤为不适零膨胀会让距离度量失真一个峰覆盖度差异很大的细胞反而可能被分到一起。如果你想了解“哪些区域协同变化、共同定义一种细胞状态”上述这些聚类方法都无法直接回答。cisTopic则不同它用概率主题模型把每个细胞拆解成若干个主题的混合。每个主题对应一组峰这些峰倾向于在某些细胞里共同开放。最终你能拿到两类概率细胞-主题概率每个细胞在多大程度上依赖某个主题和主题-峰概率每个峰在不同主题里的权重。前者可以做传统聚类和可视化后者则可以直接用来定义顺式调控程序往下对接motif分析、转录因子富集。所以可以把cisTopic理解成“聚类的上游工具”它先把数据拆成可解释的调控模块模块层面的细胞分组远比你直接基于原始峰矩阵聚类要稳定得多。2. 模型直觉从LDA到顺式调控主题本质上是“词袋”用在基因组上cisTopic的核心模型脱胎于自然语言处理里的LDALatent Dirichlet Allocation。在LDA里一篇文档是单词的集合文档会以不同比例属于若干个主题而每个主题又有一个在所有单词上的概率分布。cisTopic把整套概念平移到了基因组上每个细胞是一篇“文档”每个开放的峰是一个“词”每个主题就是一个“顺式调控主题”。模型要做的就是反推出细胞里隐藏的主题混合比例以及每个主题下峰的贡献权重。这种“文档-主题-词”的三层结构可以用一个很直白的例子理解。假如你有十篇小说小说里反复出现“大海”“航船”“海盗”那么模型大概率会抽出一个“航海主题”另一批小说如果频繁出现“城堡”“骑士”“王后”就会抽出“宫廷主题”。有些小说二者都有那它在这两个主题上的得分都会偏高。细胞也一样一个细胞里如果同时高权重命中造血干细胞的调控区域和细胞周期调控区域那它就不是“单标签”的它在这两个主题上都会有可观的活性。2.1 细胞-主题矩阵和主题-峰矩阵cisTopic输出中有两个关键矩阵理解它们你才能顺利做下游分析。第一个是细胞-主题矩阵。假设你有10000个细胞、设置了8个主题那这个矩阵就是10000行、8列每一行代表一个细胞在8个主题上的活性分布。因为每行加起来等于1所以它本质上是个概率分布。你拿这8个维度去做UMAP、t-SNE或者直接算相关系数做层次聚类都会比在原始高维峰矩阵上运行稳定得多。第二个是主题-峰矩阵。假设有500000个峰、8个主题那这个矩阵就是8行、500000列每个主题给出了它在所有峰上的权重分布。权重高的峰就是该主题的“代表词”。把这些高权重峰映射回基因组坐标就能看到该主题对应的调控区域在哪里集中把这些区域拿来跑motif富集就能推测哪些转录因子可能主导这个调控程序。这两个矩阵是cisTopic一切下游分析的基础。所谓“顺式调控主题的概率建模”最终产出其实就是这两个概率分布。2.2 cisTopic和LDA的差异点cisTopic不是把LDA原封不动搬过来它在几处做了针对基因组数据的调整。比如它会考虑峰与峰之间的基因组距离和染色质可及性的分布特征而不是完全忽略词序另一个常见调整是它会用二项式/多项式的框架去适配0/1稀疏矩阵的特性让模型对极低覆盖度细胞更稳健。另外主题数量是用户给定的。这一点和k-means需要指定簇数一样没有免费的午餐。cisTopic给了一套基于模型对数似然和主题稳定性的选择机制后面第4节我会专门讲如何选主题数这也是整个流程里最容易翻车的环节。在实际实现上cisTopic通常采用折叠吉布斯采样collapsed Gibbs sampling来估计后验分布。你不需要手动实现采样R包封装好了但理解它的意义能帮你设置迭代次数采样需要足够多的迭代让马尔可夫链收敛否则不同主题之间会出现大量“复制粘贴式”的相似结果。3. 输入数据构造把散落片段变成细胞-峰矩阵的基本功cisTopic是个R包建议在R/RStudio环境里用。安装时依赖比较多如果你之前装过Seurat、Signac、ArchR等包通常不容易冲突因为cisTopic的核心依赖是Matrix、data.table、topicmodels这些老牌R包。建议直接用Bioconductor安装安装命令如下if (!requireNamespace(BiocManager, quietly TRUE)) install.packages(BiocManager) BiocManager::install(cisTopic)安装完先不要急着跑模型先确认一下版本。cisTopic在不同版本里的函数名有细微调整比如有些版本用runModels有些版本更推荐runCGSModels作为并行入口。我的经验是跑任何函数前先看一眼帮助文档library(cisTopic) ?runModels3.1 从fragments.tsv.gz构建细胞-峰矩阵如果你是10x平台的数据最原始的输入是fragments.tsv.gz和对应的singlecell.csv里面记录每个片段来自哪个细胞、在哪条染色体哪段区间。cisTopic直接接收的通常不是这种片段文件而是已经聚合好的峰矩阵所以你需要先完成两步转换获得全样本共享的峰集合。最简单的办法是从Cell Ranger的peaks.bed读取或者用MACS2在合并后的伪批量数据上重新call peak。统计每个细胞落在每个峰里的片段数生成稀疏矩阵。这个过程可以用Signac的FeatureMatrix()也可以直接用ArchR里的getMatrixFromProject(ArchRProj, useMatrix PeakMatrix)或者用SnapATAC2的snapToFragments。矩阵准备好后一般保存成RDS或者Matrix Market格式。如果你手头只有10x的filtered_peak_bc_matrix.h5可以先在R里读出来library(Matrix) library(Seurat) mat - Read10X_h5(filtered_peak_bc_matrix.h5)这里的mat是一个稀疏矩阵行名是“chr1-10000-20000”这种坐标列名是细胞Barcode。这个格式已经可以直接进cisTopic。3.2 使用createcisTopicObject构建对象拿到稀疏矩阵后创建cisTopic对象很简单ct - createcisTopicObject( count.matrix mat, project.name my_scatac_project )需要注意这个函数默认把输入当作“cell by peak”格式处理但在构造对象时有些版本要求行是峰、列是细胞所以如果发现细胞数和峰数反了先转置矩阵再传入。创建对象时还可以传入细胞元数据比如样本来源、批次、预期细胞类型注释。方法是在createcisTopicObject里使用cell.meta参数或者后续用addCellMetadata()补充。这一步不是必须的但我强烈建议一开始就带上样本ID否则等你跑完模型发现主题主要按批次分群时才想起来没有样本信息可以对照会很被动。数据质量过滤也建议在进入cisTopic前完成。常见的过滤标准包括每个细胞至少要有500个唯一片段、TSS富集分数大于某个阈值、线粒体片段比例要低。不同平台标准不一不必教条但一定要记录过滤前后细胞数量变化不然别人问起数据量你连基线都说不清。3.3 染色体命名和坐标格式是隐藏陷阱一个特别容易被忽略的坑是染色体命名。你的峰坐标可能来自hg19或hg38但如果你混用了不同版本或者某些区域名写成“1:10000-20000”、某些写成“chr1:10000-20000”cisTopic在后续计算主题得分和重叠区域时会报错或者更麻烦的是悄悄漏掉一部分区域导致结果偏差。我在处理公共数据时遇到过这种情况一个公共矩阵用了NCBI的1_10000_20000命名另一个用了UCSC的chr1-10000-20000命名合并后如果不统一后续所有分析都是错的。所以进模型前建议用gsub清洗一遍坐标命名至少要保证全部统一为chr前缀格式并用GenomicRanges验证一下排序。4. 核心运行与主题数选择跑模型只是开始难的是决定几个主题关键代码无法保证准确的函数名。用通用描述。构建好对象后运行模型大概是这样的set.seed(123) ct - runModels( ct, topic c(2, 4, 6, 8, 10, 12), seed 123, burnin 125, iterations 250, addProps 0.05 )参数说明topic是你要尝试的主题数量集合可以直接传一组整数burnin是吉布斯采样的预热期iterations是采样的迭代次数。我一般会把迭代次数设置到500以上虽然这会明显变慢但收敛后的主题可解释性会好很多。如果时间紧张可以先跑250次看看趋势再用小迭代跑一组、大迭代跑一组做对比。cisTopic会为topic里每个候选值单独训练一个模型保留在对象的models槽里。之后用selectModel挑选最终模型ct - selectModel(ct, select c(6, 8, 10)) modelSelPlot(ct)selectModel会基于对数似然和稳定性综合排序选出你指定候选里相对合适的那个。modelSelPlot会画出一条主题数和对数似然的关系曲线你可以像看PCA碎石图一样找出“拐点”。4.1 不要只看数值还要看主题是否可解释主题选择是cisTopic流程里最需要经验的地方。很多教程会让你直接看modelSelPlot但这个图只能告诉你哪个主题数在统计上更合理不能告诉你哪个主题数在生物学上更有用。我个人的筛选流程有三个步骤。第一步看稳定性。同一主题数下用不同随机种子各跑一遍然后计算主题之间的一致性。如果两次运行得到的主题内容高度相似说明这个主题数比较稳如果两次运行的主题五花八门说明模型没有收敛或者主题数设置得不自然。第二步看主题分离度。把细胞-主题矩阵画成热图检查是不是有几个主题在所有细胞上都是均匀分布。如果出现“万金油主题”几乎每个细胞都有很高的占比那可能是主题数太多模型把一个公共调控程序硬拆成了多个碎片。第三步也是最关键的把主题映射回生物学。选主题数不比选k值不是为了数学最优而是为了“每一个主题都能讲出一个故事”。当你设置8个主题时理想状态应该是每个主题对应一个明确的细胞类型或者一个明确的生物过程比如CD4 T细胞主题、CD8 T细胞主题、单核细胞主题、细胞周期主题等。如果你发现两个主题几乎只差在一两个峰上那你的主题数可能偏多如果某个主题把好几类细胞的信号混在一起那可能偏少。4.2 多模型对比怎么看runModels会为每个主题数各训练一个模型selectModel会帮你选一个但我建议不要只依赖自动选择。我会把所有候选模型的细胞-主题矩阵都提取出来分别做一遍UMAP然后用已知的细胞类型注释去看聚类结果。有些时候主题数为6的对数似然更高但主题数为8才能把你关心的稀有细胞类型分开。对数似然是对整体数据的拟合度量它不会主动在乎你关心的那个稀有亚群。所以正确做法是用modelSelPlot缩小候选范围再用生物学注释做最终决策。这也是概率建模和纯机器学习最大的不同——下游解释才是检验标准。4.3 迭代次数和收敛性吉布斯采样类方法有个问题迭代次数不足时不同主题之间会出现“重复化”。比如你设了8个主题结果其中3个主题的峰权重分布高度相似差异只来自随机噪声。这种时候先别急着改主题数试着把迭代次数从200翻到1000很多“伪主题”会自己合并。我在真实数据上观察过低迭代次数下细胞周期相关峰会被拆到好几个主题里提高到足够迭代次数后它们会稳定归到一个主题。这让我养成了一个习惯——小规模调参时用快速设置一旦确定了主题数范围最后一定用长迭代完整跑一遍并保存对象副本。5. 下游实操细胞聚类、可视化、motif与转录因子富集模型跑完只是“烧完钱”真正产生生物学价值的从这一步开始。cisTopic的下游分析可以分成四条线细胞聚类、主题区域可视化、主题归一到基因组、转录因子富集。5.1 用细胞-主题矩阵做聚类和可视化先提取细胞-主题概率cell_topic - getCellData(ct, select c(8), dims 2)这里select参数选择你最终确定的那个模型。getCellData会返回一个包含降维结果的对象你要的细胞-主题矩阵在这个对象的cell.data里。通常我们再对主题概率做一次归一化或转换然后跑UMAPlibrary(umap) umap_out - umap(cell_topic$cell.data[, grep(Topic, colnames(cell_topic$cell.data))])如果你更喜欢用Seurat生态也可以把这个细胞-主题矩阵直接塞进CreateSeuratObject()然后走RunUMAP()和FindClusters()。很多团队会对比直接基于LSI的聚类和基于cisTopic主题的聚类发现后者往往在稀有细胞类型的分离上更有优势。另一个常见操作是基于主题概率做层次聚类。因为主题数通常很少比如8个、12个你完全可以把每个细胞看作一个8维或12维的概率向量计算它们的相关性距离然后画层次聚类树并用pheatmap展示。library(pheatmap) hm - pheatmap(cor(t(cell_topic$cell.data[, topic_cols])))这种做法的好处是你可以看到细胞在主题层面的树状结构层次聚类能告诉你哪些主题更接近哪些主题分化更早。把聚类结果和已知注释对照能快速判断主题是否有生物学意义。5.2 主题活性在基因组上的映射cisTopic的一个独特能力是给每个峰打分找出每个主题的“代表峰”。做法大致是region_scores - getRegionsScores(ct, select c(8), method NormTop, scale TRUE)method通常有“Z-score”或“NormTop”等选项。得到的region_scores里每个主题对应一组区域分数分数越高说明该区域对这个主题越重要。拿到这个分数后你可以做几件事筛选每个主题的高分区域导出成BED文件去IGV里看具体信号。把这些区域和基因注释、增强子数据库如ENCODE、ENdb重叠判断主题涉及哪些功能元件。把这些区域作为输入跑转录因子富集。比如你的主题8里高分区大量落在某个基因的启动子和近端增强子附近那这个主题就可能与该基因所在细胞类型的功能分化有关。把高分区做motif富集可以直接推测调控该主题的核心转录因子。5.3 转录因子富集怎么做cisTopic官方教程里通常会把主题区域导出给RcisTarget或者HOMER做motif富集。RcisTarget的优势是不需要把区域严格二值化它可以用区域排序来做富集更平滑。简单流程是# 假设已经安装并载入RcisTarget library(RcisTarget) motifRankings - importRanks(hg38_500bp_up_and_100bp_down_tss.mc9nr.feather)然后按主题的高分峰列表做addMotifAnnotation和addSignificantGenes。最终你会得到一张表哪个motif在某个主题的调控区域里显著富集以及对应的转录因子和可能的靶基因。这里有个实操经验不要只看p值最小的motif要看“motif富集区域是否真的是该主题特有的”。比如有些基础转录因子如p300、CTCF结合基序几乎在每个主题里都富集因为它们本来就广泛参与染色质结构维持。更有信息量的是“差异化富集”也就是某个motif在主题A里显著富集但在其他主题里不富集。我建议把每个主题的高分区域批量跑一遍富集然后做对比矩阵而不是单独看一个主题。5.4 主题和差异可及性分析的结合另一个常见下游是差异可及性分析。你可以在UMAP聚类后把不同细胞簇的cisTopic主题活性做差异比较找出某个簇特异激活的主题。这种“先聚类再映射到主题”的路线比直接跑差异峰要稳因为主题已经把无数单峰信号压缩成了几个可解释的调控模块。例如你发现了两个与细胞周期相关的主题一个在G1/S期细胞里高一个在G2/M期细胞里高。两边高分区对应的motif可能完全不同这可以为你提供非常有针对性的转录因子线索。这类结果放到论文里比单纯给出一堆差异peak要有说服力得多。6. 和其他聚类路线对照k-means、层次聚类、LSI之后到底该用谁因为cisTopic经常被放在“聚类教程”这个检索场景里我觉得有必要把它和k-means、层次聚类、子空间聚类这些常见路线放在一起做一次明确对比免得大家盲目套用。方法输入聚类类型主要输出适用场景k-means低维向量如LSI降维结果硬聚类细胞簇标签快速分群、大型数据集层次聚类低维向量或距离矩阵硬聚类/树状图树、聚类标签小样本、想看层级关系子空间聚类高维矩阵软/硬子空间划分高维信号集中在局部维度cisTopic原始稀疏峰矩阵软主题聚类主题分布、峰权重找调控程序、机制解释这张表不是想说cisTopic比k-means高级而是想强调它们回答的问题不一样。k-means回答的是“每个细胞该贴哪个标签”cisTopic回答的是“这个细胞的状态由哪些调控程序叠加而成”。在很多分析流水线里两者可以串联使用先用cisTopic提取主题再用k-means对主题概率向量聚类得到细胞簇然后用层次聚类展示主题层级。如果你已经有Seurat/Signac的完整流程跑过LSI降维和FindClusters再去叠加cisTopic会有点重复但不冲突。我通常的做法是Signac先给我一个快速可靠的细胞注释然后我用注释去校验cisTopic主题的可解释性而不是反过来让cisTopic来做最基础的细胞类型鉴定。思路是细胞注释交给成熟的快速流程调控机制挖掘交给cisTopic。6.1 什么时候可以不上cisTopic如果你的目标只是“把细胞分出来”而且你的数据量特别大、时间很紧直接用LSI加k-means完全够用没必要上主题模型。主题模型训练慢调参空间大解释起来也费劲用它回答一个简单分群问题属于高射炮打蚊子。如果你的数据里有非常明显的结构比如几个大细胞类型各占一大块LSI聚类效果通常不会差。但如果你关心的是细胞类型内部的连续状态比如分化轨迹上的中间态或者同一个细胞类型里不同调控程序的协同变化k-means硬聚类就有点不够看了。这时候cisTopic的软主题结构就体现出价值它不会强迫中间态细胞选边站而是给出一组“高不成低不就”的主题组合恰恰是这种组合量化了过渡状态。6.2 关于“子空间聚类概述及python完整实现”热词的一点延伸有朋友搜索“子空间聚类”这类方法在高维数据里寻找低维子空间然后再聚类与主题建模的哲学有一些相通之处。但子空间聚类对混合模型假设有很强的依赖应用于单细胞表观组学这种极度稀疏数据时需要对缺失机制做额外假设。cisTopic走的是更直接的生成模型路线它通过对“峰-主题-细胞”的联合概率建模让缺失和噪声都统一在概率框架下处理所以应用起来更顺手。如果你本身就是Python栈不想切到R也可以考虑用gensim或scikit-learn里的NMF来做类似的主题分解。NMF在非负稀疏矩阵上也能提取“主题”式的因子结构但没有Dirichlet先验结果的可解释性和稳定性通常没有cisTopic好。如果你用Python做最常见的路线是先做TF-IDF或LSI再用KMeans聚类注意这里因为二值数据没有“词频”概念TF-IDF要改造成二值加权形式否则效果会偏移。6.3 聚类热图、趋势图、富集条目的串联还有一种常见需求是想把聚类结果跟富集条目串成一张“展示图”。比如你已经把细胞分成了几簇每个簇有一组标记峰想看看这些峰的motif富集条目和趋势。一般流程是用主题的高分峰取并集或交集画一个主题×峰的聚类热图再在旁边标注富集到的转录因子条目。cisTopic已经帮你把“峰降维成主题”了所以这一步做起来比直接从几万个峰里挑特征要干净很多。7. 我踩过的坑与调参经验用cisTopic这两年我踩过的坑不少。很多坑不是模型本身的问题而是使用习惯问题。挑几个印象最深的分享一下。7.1 主题太多反而丢失稀有细胞类型第一次跑正式数据时我参考某些教程设了20个主题结果日志似然很好看但稀有细胞类型对应的主题被切割得七零八落几乎没有一个是干净的。后来把主题数降到10左右稀有亚群的信号才集中起来。这不是说10比20好而是说“主题越多每个主题的平均丰度越低”稀有信号更容易被拆碎。如果你特别关注稀有细胞建议主题数取偏小的区间比如5到10然后配合已知marker基因去验证。7.2 批次效应会伪装成主题多批次数据尤其容易出这个问题。如果主题A在批次1的细胞里高、主题B在批次2的细胞里高但两个批次注明的细胞类型完全相同那要么是数据质量差异太大要么是批次效应被模型捕捉了。遇到这种情况先别急着下生物学结论。可以试两个方向一是检查峰矩阵是否做过批次间的归一化比如用chromVAR的偏差校正或者简单的分位数归一化二是在聚类前用Harmony或Seurat的整合功能在主题概率矩阵上做批次整合。主题概率矩阵因为维度很低整合速度也很快。7.3 染色体坐标版本混用前面提过一次但值得再强调。hg19和hg38坐标混用反映在结果里很可能就是一个主题的区域几乎全部来自某个坐标版本。这个问题排查起来很隐蔽因为模型不会报错只有等下游重叠数据库时才会发现大量区域重叠不上。建议在创建对象之前用GenomicRanges::seqinfo()统一检查一下染色体长度和命名并保留一份坐标版本记录。7.4 内存和时间优化几十万峰、几万细胞的数据跑多主题数、多迭代的模型确实会慢。我一般在正式跑之前先做一个“下采样版”随机抽5000个细胞、50000个峰快速跑一遍确定主题数范围和迭代次数再用全量数据正式跑。这样能节省大量试错时间。另外虽然runModels支持多主题数同时跑但我还是建议分多次运行每次隔一段时间保存一个RDS。一旦某个主题数的模型因为内存问题崩溃你不需要全部重来。7.5 保存一切包括随机种子概率模型的随机性不可忽视。同一份数据同样参数两次运行的主题顺序可能完全不同甚至内容有一定波动。所以我在每次runModels前都固定种子并把种子写进项目元数据每个阶段都保存一个独立RDS文件比如ct_raw.rds、ct_model.rds、ct_selected.rds。这样即使后面分析叉了也可以回到任意一个节点不用从头跑。最后再分享一个小技巧选定主题后把你的主题编号、对应细胞类型注释、每个主题的代表峰个数、代表性转录因子做成一张总表。这张表既是和合作者讨论的基础也是后续写论文时的素材库。做分析的时候你可能觉得这些记录很冗余但几个月后回头看你会感谢当时的自己。
