染色质可及性、单分子、新生哺乳动物这三个词并排放在一个项目标题里外行看着像一句绕口令懂行的人却知道分量不轻。所谓染色质可及性通俗讲就是基因组上某些区域对调控蛋白开门还是关门门开着的区域转录因子能上去增强子和启动子才能干活基因才有机会被转录门关着的区域再强的调控信号也进不去。普通ATAC-seq做的是把几千、几万个细胞的信号叠在一起告诉你某个区域整体看开了一半还是关着。而这个项目要往前走一步在单细胞、甚至单条染色质分子的水平上把每个区域的开与关拆开来看再拼出一张属于新生哺乳动物的开门地图。为什么盯住新生动物因为出生前后恰好是所有器官大规模换挡的窗口——肝脏要接管糖代谢肺要准备呼吸肠道要适应奶水免疫系统要认识一群新来的微生物。这些转变背后的司令部就是染色质可及性的重编程。这篇文章我把这个项目从实验设计到数据分析的完整思路、关键参数和踩过的坑一条条讲清楚。适合正准备开单细胞染色质可及性项目、或者刚接触表观调控数据分析的团队和同学参考。1. 项目思路拆解新生窗口与单分子视角的碰撞1.1 新生期本质上是一次全基因组范围的换挡哺乳动物从子宫内切换到子宫外不是简单挪个地方而是全身各个器官在一两天内完成功能接管。在子宫里营养靠胎盘氧气靠母体循环肠道基本处于待机状态出生那一刻起肺必须立即开始气体交换肝脏要在几小时内启动糖异生保证血糖不掉下去小肠上皮要改变通透性迎接母乳刺激皮肤要开始面对干燥环境。这些生理切换的底层是大量转录因子在特定增强子上重新部署而增强子在基因组上有没有被搭理最直接的观测指标就是染色质可及性。所以时间点设计上我们选了三个节点出生前1天、出生当天P0、出生后7天左右。这三个点拉开的不是日历上的几天而是宫内准备期—出生应激期—稳态适应期三个生物学状态的切换。在这个窗口里可及性信号的变动极其剧烈peak数量、峰强度、细胞组成都会大幅变化。反过来说如果只拿成年个体做这些扰动早就平复干净了什么都看不出来。做这类发育转型期的项目最怕的就是时间点拉得太松错过了真正有信息量的转折。还有一个很容易被忽略的事实新生动物的组织非常小一只新生鼠的肺、胰腺拿出来能提的细胞数量相当有限。这直接否定了按成人组织建库的常规做法逼着我们在样本合并、低起始量建库上花功夫。这个细节后面单独说。1.2 群体平均遮蔽了单分子的开关差异传统ATAC-seq给出的结果是某个区间内reads的累积信号。假设某个增强子在10万个细胞里有50%是打开的批量ATAC会告诉你这儿有个中等强度的peak但它没法回答是全部细胞都开了50%的染色质还是50%的细胞完全打开、50%的细胞完全关闭这两种情况生物学含义完全不同——前者可能是动态振荡或者受细胞周期影响后者则意味着存在两个真正不同的细胞状态只是被混匀了。这就是景观两个字的真正含义。要画的是不同细胞类型、不同发育阶段里每个调控元件的打开比例和打开方式而不是一张模糊的平均值热图。单细胞层面的scATAC-seq把这个问题拆开了对每个细胞我们拿到一份二值化的可及性profile对一群同类型细胞某个位点的可及性就变成一个0到1之间的比例值。这个比例比信号强度更接近开关语义。严格意义上讲单分子层面还需要再进一步。如果配合长读长测序平台用甲基转移酶在可及的位点上打标记我们就能在单条DNA分子上同时看到几个增强子的协同打开状态它们是在同一条纤维上一并开放还是各自独立我们项目把两条线都做了——主线用单细胞ATAC高通量画整个可及性景观支线用长读长在关键区间做正交验证。这样既有统计火力又有分子层面的物理证据文章说服力完全不一样。1.3 三个问题把整个项目框死动工之前我们和课题负责人把目标收敛成三个问题免得后面被海量数据带偏新生期各器官包含哪些细胞类型它们在出生前后各自经历怎样的可及性重编程哪些转录因子的结合位点在这个窗口集中获得或失去可及性构成器官成熟的核心调控模块关键增强子与启动子之间的连接co-accessibility是否随时间点重构能否解释已知的出生后表达变化这三个问题分别对应分析里的聚类注释、motif富集和共可及性网络每一步分析都围着它们转。我强烈建议所有做类似单细胞研究的团队动手前先写这样一张问题清单。否则数据一出来你只会花大量时间盯着UMAP图发呆然后被各种漂亮的降维图带着跑最后一篇文章一个明确结论都立不住。2. 样本与实验设计从组织块到可分析的文库2.1 多器官、多时间点的取样策略器官选择上我们最终定了肝、肺、小肠、脑四个。理由很简单这四个器官覆盖了新生期三类典型转变——代谢接管肝、外部环境直接接触肺、肠道、以及相对保守的神经发育窗口脑。多器官设计的好处是后面分析可以做跨器官对比能清晰区分出生本身带来的通用重编程和器官特异的成熟程序这是单器官项目给不了的视角。样本分组是门学问。我们每个时间点用4到6窝鼠每窝取两到三只雌雄尽量均衡然后按器官收集合并成一个样本。这里有个关键原则永远不要在单只动物上做生物学重复。新生动物个体间差异很大满肚子奶和刚出生两小时的鼠同一个基因的表达能差一倍不止。所以正确的做法是每盒做一个池子池子与池子之间才算真正独立的生物学重复。我们每组做了三个池后面差异分析才不至于被个体噪声淹没。组织解离这一步各器官脾气完全不同。肝脏松散好解但红细胞、肝实质细胞占大头要注意去除碎片肺是个难题胶原含量高需要酶消化加机械解离配合时间控制不好细胞活性掉得很快小肠上皮容易在解离过程中丢失得提前用EDTA把上皮层揭下来再单独处理。核提取阶段的细节直接决定文库质量这块多花一周优化是值得的。2.2 细胞核提取与转座反应一个参数能毁掉整个文库scATAC-seq的核心是Tn5转座酶。这个酶在结合了转座子末端序列的复合物里会一边切断DNA一边把测序接头插进去这个过程叫tagmentation。Tn5只对可进入的染色质敏感进入不了的区域它下不去刀。正因如此Tn5反应就成了整套实验的命门。我们用的核提取方案是低渗裂解加去垢剂尽可能少用机械剪切避免染色质结构被破坏。具体配方大概围绕10 mM Tris-HCl、10 mM NaCl、3 mM MgCl2调整渗透压加0.1% NP-40破膜全程在冰上操作离完心里默念着别让核膜碎掉。每步都要用DAPI染色看核的完整性和纯度——这也是我后来特别坚持的一个习惯建库前不看到清晰的核影像绝不停下。转座反应的时间必须做滴定。我们拿同样一批核分别设定5、10、20、30分钟四个条件看最终文库的插入片段分布和TSS富集。时间太短转座不充分文库复杂度低时间太长Tn5会开始在一些半开放区域随机落刀背景噪声暴涨。我们最终锁定了15分钟、37度、转座酶浓度固定的方案但不同批次的Tn5酶活性会有差异每换一个批号都要重新滴定一次这点尤其容易在赶进度时被省略。2.3 质控硬指标TSS富集、片段分布与线粒体占比单细胞ATAC建库完成测序之前有几项硬指标必须先看。这里我直接给一张我们内部一直用的质控表指标建议阈值不合格时最常见的坑TSS富集分数大于8细胞核纯度差、转座过度、核被切碎每个细胞有效片段数大于1000测序深度不足、细胞数估错、接头连接效率低片段在peak内的比例FRiP大于20%背景噪声高、比对到错误参考基因组插入片段呈核小体周期分布存在147bp左右的周期性非特异转座过多、反应过度线粒体reads占比小于20%核膜破损、线粒体DNA被当作开放染色质转座TSS富集分数说白了就是把所有基因的转录起始位点对齐看中间200bp区域的信号是不是明显高于两侧2kb的背景。这个分数如果小于8说明你的开放染色质信号很多是噪声后面所有峰值分析都会带上脏东西。我见过一些项目怀着先测了再说后面用算法洗的心态跳过这个检查最后只能花更多时间在数据里抠污染。线粒体占比也是一个常被忽略的信号。线粒体基因组是裸露的环状DNATn5对它的转座效率天然很高在完整的细胞核提取流程里占比应当在5%-10%以下。如果某个样本线粒体reads到了30%以上基本可以认定大量细胞核破裂细胞质里的线粒体DNA混了进来。这个指标在普通bulk ATAC里不太看但在单细胞水平它能把一批样本直接判死刑一定要提前卡住。3. 分析流程从FASTQ到可及性景观的核心环节3.1 比对、去重和fragment文件生成干实验流程的第一站是把测序原始数据变成一张哪个细胞的哪条片段来源于哪个基因组位置的表格。这个文件通常叫fragments.tsv.gz每一行记录一个Tn5插入事件染色体、起始、终止、细胞标签、唯一分子标签的重复计数。比对这块我们直接用快速比对工具配上默认参数然后按标准约定把插入位置做偏移校正Tn5转座后两条链的插入信号实际落在靶位点上游4bp和下游5bp的位置比对后要在正链加上4bp、负链减去5bp才能精确还原真正被切开的那个碱基。实际操作里PCR重复的去除也很关键。scATAC每个细胞的有效片段数本来就珍贵如果重复率过高相当于测序钱花在了同一个模板上。我们会用sambamba按坐标和UMI去重再按细胞标签拆分统计。下面这段是把fragment文件整理后做peak调用的常用操作# 已按细胞拆分后的片段文件合并去重 sort -k1,1 -k2,2n sample.all_fragments.bed | \ bedtools merge -i stdin sample.merged_fragments.bed # 基于合并片段调用peak标准的ATAC参数 macs2 callpeak \ -t sample.merged_fragments.bed \ -f BED \ -g mm \ -n sample \ -q 0.05 \ --shift -100 --extsize 200 --nomodel \ --keep-dup all这里--shift -100 --extsize 200的组合是ATAC-seq时代沉淀下来的经验值目的是把Tn5的切口位置扩展成完整的开放区间让peak更接近真实的开放染色质区域。--nomodel是因为Tn5的插入信号本身就是锐利的不需要像ChIP-seq那样估计片段长度分布。3.2 单细胞质控、降维、聚类与注释拿到fragment文件之后单细胞层面的分析我们用的主力框架是ArchR。这个工具比早期那套老流程稳不少尤其是对TSS富集过滤、doublet识别和基因活性打分都有内置实现。核心步骤可以浓缩成下面这段R代码library(ArchR) addArchRThreads(threads 12) set.seed(1) ArrowFiles - createArrowFiles( inputFiles c(SampleA A.fragments.tsv.gz, SampleB B.fragments.tsv.gz), sampleNames c(SampleA, SampleB), minTSS 8, minFrags 1000, addTileMat TRUE, addGeneScoreMat TRUE ) proj - ArchRProject(ArrowFiles) proj - filterDoublets(proj) proj - addIterativeLSI(proj, useMatrix TileMatrix) proj - addClusters(proj, resolution 0.8) proj - addUMAP(proj) # 用基因活性打分辅助注释 markers - c(Alb, Cdh5, Ptprc, Epcam) p - plotGenes(proj, markerGenes markers, plotAs umap)为什么ATAC的数据降维不直接用PCA而是用LSI原因很实在单细胞ATAC矩阵是超级稀疏的0/1矩阵绝大多数位点在绝大多数细胞里都是0直接做PCA会被这些零值主导。LSI来自文本挖掘的思路先做TF-IDF转换term frequency-inverse document frequency把在所有细胞里都开放的位点降权把只在少数细胞里开放的位点加权这样降维出来的空间才真正反映细胞身份差异。这个道理我们组刚上手的新同学经常想不通但只要你对比过PCA和LSI出来的聚类结果回不去了。细胞类型注释我们主要靠已知marker基因在UMAP上的基因活性得分分布。由于ATAC数据里没有直接的表达量基因活性得分是把该基因TSS上下游约100kb范围内的可及性按距离衰减加权求和得到的。这个分数和真实表达量有一定对应关系但要注意它不是万能的有些基因表达受转录后调控基因活性高不代表蛋白水平高。注释阶段一定要多个marker交叉验证并且回到原始reads分布上看峰形别光看一坨颜色。3.3 可及性景观的量化与呈现聚类注释完接下来就是把这个项目的重头戏——量化可及性景观。我们做的第一件事是把每个细胞类型在每个时间点的数据按细胞合并成伪bulk然后在伪bulk上重新调用peak。这一步非常关键单细胞层面的peak调用噪声大伪bulk能给出更干净的peak集再用这个peak集反过来对每个细胞统计是否至少有一条fragment覆盖。这样每个细胞类型和时间点的组合就得到一组开放频率值某个增强子在某个时间点、某类细胞里有30%的细胞开放另一个时间点可能涨到70%。这个百分比比reads强度更直观也更贴近开关的生物学语义。画图的时候我们习惯把这组比例值做成热图行是peak列是细胞类型乘以时间点。如果你看到一大片峰在出生当天颜色骤变那基本可以锁定向这个窗口内关键的调控元件。差异可及性分析我们是在伪bulk层面用常规的RNA-seq差异分析思路来做的。先把每个细胞类型的peak矩阵汇总成样本级别的count再用edgeR或者limma做配对比较。这里有个小讲究不同生物学重复之间的文库大小差异必须显式处理否则某些群体可及性升高的结论很可能是文库大小漂移造成的假象。我们通常会加一个总count归一化再加trimmed mean of M-valuesTMM校正双保险。3.4 调控线索挖掘motif富集、足迹分析与共可及性网络有了峰集和差异峰下面要回答第二个问题哪些转录因子在这个窗口接管了调控权。最常用的两招是motif富集和足迹分析。motif富集用HOMER跑一遍findMotifsGenome.pl输入差异获得的peak区域背景用所有peak做这样富集出来的是相对于全基因组可及区域更侧重的因子而不是相对于全基因组。这一步参数要是选错会给你报出一堆与生物学毫无关系的管家因子很常见。足迹分析则更进一步。Tn5在结合了转录因子的DNA区域是进不去的因此把大量片段在某个因子结合位点周围堆叠后你会发现结合位点正中心有一个小小的信号凹陷像踩在泥里的脚印。这就是足迹。我们用了TOBIAS做偏置校正——Tn5自身有序列偏好性不校正的话会把偏好性误判成足迹。足迹分析能告诉我们某个因子的结合位点不仅序列上存在而且在该细胞类型里很可能真的结合了因为那里确实被保护住了。第三个问题关于增强子-启动子连接我们用了Cicero从单细胞可及性共变里推断顺式调控连接。基本原理是同一条DNA链上两个位点如果经常在同一批细胞里一起开放它们之间就可能存在调控连接。这当然只是统计推断但在新生物种里没有大量表达数据做支撑时它依然能给出不错的候选增强子-启动子配对后面用长读长验证哪几对是真正共变的。4. 常见问题与排查技巧实录4.1 线粒体reads高到离谱第一轮测序结果回来我们有个肺样本线粒体占比接近35%。当时第一反应是怀疑转座反应过度但复检发现是核提取那步出了问题——肺组织解离后碎片多杂质里大量线粒体被带进了核悬液。我们后来在纯化环节多加了两次蔗糖垫离心把线粒体占比压回到8%以内。排查这类问题记住一个原则线粒体占比已经摆在那就别再花力气试图用过滤把它洗掉。你可以过滤掉高线粒体的细胞但剩下的细胞如果普遍质量差整个样本的统计效力都受损。早期在核提取后加一步流式分选比后期任何计算补救都管用。4.2 TSS富集不足小肠样本的TSS富集连续两次都只有4.5左右。一开始我们怪建库试剂后来发现是取样环节的RNA酶污染。新生小鼠小肠含大量消化酶组织解离时如果不加足量的RNA酶抑制剂和蛋白酶抑制剂混合物核悬液里满是消化酶残渣Tn5的活性被严重干扰。我们把解离缓冲液里的抑制剂体系换了一整套TSS富集立刻回到了9以上。这件事给我的教训是实验侧的质控要跟分析侧的数据一起看不能只看出了多少峰。TSS富集是核完整性和酶活性的综合反映它低一定说明上游有问题别指望用分析端的硬过滤解决。4.3 批次效应伪装成生物学差异新生小鼠出生当天的时间点样本我们分别在两个星期里各做了一轮实验。第一轮聚类时发现同一个时间点的两个池子居然各自聚成一团分得比不同器官之间还清楚。排查起来发现是Tn5酶批号换了新批号的酶活比旧批号低导致第二轮整体的转座效率下降开放区域的覆盖均匀度完全不同。我们最终做了两件事一是所有时间点重新补做一轮完整实验保证每个时间点的重复都在同一批内完成二是在分析里用harmony做批次整合。更重要的教训是实验一开始就该把同组重复标本打散到不同批次里让批次效应可以被统计模型估计掉。这个操作叫block randomization做单细胞实验少这一步后面再强的算法也很难洗干净。4.4 过渡态细胞群的注释陷阱新生期组织里充满了正在分化路上的细胞最容易出问题的是注释阶段硬把一个过渡群塞进已知类型里。我们的肝样本里有一个细胞群Epcam和Alb的基因活性都有点但都不到典型值。你要是按肝细胞硬归后面所有差异分析都会被污染。后来我们用了一个笨但有效的办法给每个候选细胞类型做成熟度打分——挑一批已知的终末分化marker把它们的基因活性加权求和再按分数把过渡态的细胞单独拎出来作为一类。这样既不影响成熟群体的注释又能保留出生后正在发生的分化轨迹信息。做单细胞分析宁可多留几个未知群也不要为了凑注释率强分硬归。5. 工具选型与参数速查5.1 单分子可及性分析工具对照表这个项目从头到尾我们用了一套相对稳定的工具组合这里整理成表供参考环节工具关键参数/说明比对bwa-mem默认参数注意对线粒体reads单独统计片段去重sambamba markdup按坐标加UMI去重统计重复率片段合并管理bedtools merge用于伪bulk的peak调用Peak调用MACS2--shift -100 --extsize 200 --nomodel单细胞质控/聚类ArchRminTSS8minFrags1000filterDoublets开可及性矩阵分析SnapATAC2备选对超大样本更快适合集群环境Motif富集HOMER背景必须选全peak集足迹分析TOBIAS先跑bias correction再跑footprint共可及性网络Cicero需先构建CDS对象计算量较大差异分析edgeR / limma伪bulk层面TMM归一化5.2 分析资源与时间估算很多人低估了单细胞ATAC数据在计算量上的胃口。我们的经验是一个10万细胞规模、覆盖4个器官3个时间点的项目每个样本的比对和fragment生成大约需要4到6个CPU核心并行半小时到一小时ArchR的Arrow生成和LSI迭代会吃比较多内存16GB内存的机器跑5万细胞没问题但到15万细胞就需要32GB以上。这一步建议直接上服务器或者云的配置型实例别在笔记本上硬扛。测序深度方面我们每个细胞的目标有效片段数定在3万到5万之间。如果经费紧张2万也可以但下游的双重过滤和低质量细胞剔除会把可用细胞数削掉一大截实际能用的细胞可能只有入池的一半。预算上建议多留20%-30%给补测和重复批次宁可首轮少测几个细胞也别把整个项目压在一批测序上。时间线上实验加建库大概一个月首轮测序两周分析到拿到完整可复现的景观图至少要一个月。如果算上验证实验和论文整理整个周期半年到八个月是常态。做课题规划的时候别按数据到手就能写文章来排时间表。最后再分享一点个人经验。这个项目做下来我最大的体会是单细胞可及性数据不缺工具缺的是对每个参数背后是什么生物学过程的理解。LSI为什么要做TF-IDF、TSS富集为什么会受酶活性影响、共可及性推断为什么需要长读长验证——这些问题的答案不在软件文档里而在你对染色质生物学和实验原理的掌握程度上。项目收尾阶段我们留了一小管P0和P7的核悬液冻存在液氮里本来是为了应付突发补测后来恰好用上做了关键区域的长读长验证。这算是我给所有做类似项目的人一个偏冷门的建议永远留一份高质量的核或者文库金样本它能在项目最困难的时候救你一命。