自己做了一批加了 E. coli spike-in 的 CUTTag 数据本来以为只是“多了几步比对的活”结果处理到一半发现事情没那么简单。最开始我图省事把所有的 reads 直接一股脑比对到人类参考基因组上后面出来一堆分布在 gene body 和基因间区的“宽峰”细看很多其实是 E. coli 来源的序列被错误比对过来的。更麻烦的是不同批次样本的 signal 完全没法横向比因为根本不知道每个样本里有效细胞数差了多少。这个项目做完以后我最大的感受是加了 spike-in 的 CUTTag 分析不是“普通 CUTTag 分析 一个额外步骤”而是整套标准化逻辑都要换掉。从双轮比对的顺序到 scale factor 的计算再到 peak 定量和差异分析每一步都得把 spike-in 考虑进去。这篇就把我最后跑通、也验证过结果的完整流程捋一遍重点讲清楚为什么这样做、参数怎么定、以及我在里面踩过的几个坑。适合手里已经有加了 spike-in 的 CUTTag 数据的分析人员也适合准备做实验、但不太确定数据分析该走什么流程的朋友。1. 为什么 CUTTag 需要 spike-in没有 input 的实验拿什么当基准1.1 CUTTag 和 ChIP-seq 在定量本质上的区别做 ChIP-seq 的都知道样品组总要配一个 input用来反映“染色质本身的可及性”和“建库过程的偏好”。但 CUTTag 不存在这样的 input。它的核心是让 pA/G-Tn5 融合蛋白在抗体引导下定位到目标蛋白结合区域然后加入 Mg2 激活 Tn5直接在原位切 DNA 加上接头。这个反应发生在细胞核里没有超声打断也没有通用的“未富集对照”。所以如果你试图用普通 CUTTag 的结果除以一个 input会发现根本没有合适的 input 可用。IgG 或未加抗体的对照组确实可以拿来评估非特异背景但它只反映抗体乱结合的水平反映不了 Tn5 本身在完整染色质上的随机切刻强度更反映不了不同样本之间“有效反应概率”的差异。1.2 spike-in 细胞的原理内置的“每细胞背景标尺”我们实验室在样本制备阶段会在同一个管子里按比例掺入 E. coli 细胞然后让样本细胞和 E. coli 细胞一起经历后续所有的步骤透化、抗体孵育、二抗、pA/G-Tn5 结合、转座反应、建库扩增。E. coli 基因组上并没有你加的抗体对应的靶蛋白它产生的 reads 基本上都来自非特异性的结合和 Tn5 的随机背景切刻。这几千个 E. coli reads 是什么它们就是这条管子里“背景反应强度”的一个非常清晰可见的代理指标。同一批实验里如果某个样本细胞数少、或者破膜效果差、或者转座酶活性被杂质抑制它的背景切刻水平也会成比例地下降映射出来的就是 E. coli reads 数量变少。反过来也一样。所以 spike-in 本质上给了你一条可以通过测序直接观测的“实验效率内参线”。1.3 不同种类 spike-in 的适用场景现在市面上的 CUTTag 试剂盒里常见的 spike-in 主要有这么几类类型常见来源适用场景加入形式细菌染色质E. coli如 K12 MG1655人和小鼠等哺乳动物样本细胞混合酵母染色质S. cerevisiae植物、真菌等参考基因组比较复杂的样本细胞混合商品化 DNA spike-inEpiCypher 等用于绝对定量和批次质控游离 DNA人、小鼠样本现在最普及的还是 E. coli 染色质因为参考基因组小、序列差异大和哺乳动物基因组比对时非常容易区分干净。植物样本因为细胞壁和蛋白结构的差异很多人会换酵母。如果用的是商品化 DNA spike-in比如 SNAP-CUTANA 甲基化 panel那它是纯 DNA 形式加入不走抗体和转座酶反应所以它反映的其实是“建库与测序”的技术偏差而不是“转座酶背景”。分析逻辑上类似但不能简单拿同一个标准化公式硬套这点我会在第 6 章单独展开。1.4 spike-in 加入比例怎么定我们实验室默认按样本细胞数的 1% 加入也就是 1:100。这个比例用下来对大多数组蛋白修饰H3K27ac、H3K4me3、H3K27me3 这种都比较合适最后测序数据里 E. coli reads 占比大概在 0.5% 到 5% 之间既不会影响参考基因组上的信号又能提供足够的背景计数来做归一化。如果你做的是低丰度转录因子比如某些只在特定细胞类型表达的 TF它的 signal 占比会小很多此时可以用到 1:50 甚至 1:20 来提高 spike-in counts 的可靠性。但反过来spike-in 比例太高会挤占有效测序资源本身就是成本浪费。有一个基本判断标准每个样本最终能比对到 spike-in 的 unique reads 不要低于 1000 条否则声标准化因子非常不稳定。我之前因为细胞数量估算失误有一个样本 E. coli 只比对了 400 多条 reads那个样本后续和同组样本比较的结果明显抖动后来又补做了才救回来。2. 分析前先别急着比对数据质量检查和 spike-in 信号初判2.1 用 MultiQC 把所有 FastQC 报告汇总看一遍CUTTag 建库其实非常简单所以数据质量上最容易出问题的反而不是接头而是 GC 偏差和插入片段分布。建议拿到 demux 后的 fastq 后先用 fastqc 跑一遍再顺手 MultiQC 汇总。主要看三件事平均质量分数Q30 占比低于 75% 的样本要留意看是不是整个 run 的问题还是单个样本降解了。Adapter 残留因为 CUTTag 片段普遍较短如果测序平台没有把 adapter 完全切干净下游比对会出现大量 reads 比不过去的情况。平台保证 clean data 的情况下一般还好。GC 分布如果 GC 曲线出现明显双峰可能是混入了大量 E. coliGC 约 50.8%和人类基因组GC 约 41%的数据这恰好也是 spike-in 存在的一个直观提醒。2.2 快速估算 spike-in 读数比例正式流程前我习惯先随机抽 10 万条 reads 快速比对到 spike-in 基因组上看比例是否和预期一致。方法很简单用 seqkit sample 抽子集然后 bowtie2 比对到 E. coli 索引再数一下比对上的 reads 数。# 抽 10 万条 reads seqkit sample -2 -n 100000 -o sub.R1.fastq.gz -o sub.R2.fastq.gz sample_R1.fastq.gz sample_R2.fastq.gz # 快速比对 bowtie2 --very-sensitive --end-to-end \ -x ecoli_k12_mg1655 \ -1 sub.R1.fastq.gz -2 sub.R2.fastq.gz \ -S sub_vs_ecoli.sam 2 align_summary.txt cat align_summary.txt | grep overall alignment rate如果整体比对率在 0.5%~5% 之间说明 spike-in 加得没问题可以安心往下走。如果 10 万条里一条都没有先怀疑索引文件是不是搞错了其次怀疑实验步骤里 spike-in 是否真的加入了。这一步成本很低能省掉后面所有基于错误前提的分析。2.3 insert size 分布是判断建库质量的硬指标CUTTag 有一个特点因为转座反应发生在核小体两侧所以插入片段会呈现明显的 mono-nucleosome 模式主峰在 150~250 bp 附近另外还有一部分 sub-nucleosomal 的小片段 120 bp在双端测序的插入片段分布图上能看到一个“小峰 大峰”的组合形态。这个分布可以直接用 fastp 或者比对后 deeptools 的 bamPEFragmentSize 看。我遇到过 insert size 主峰跑到 500 bp 以上的样本后来排查发现是细胞裂解过度、染色质降解成了大片断污染。这种样本即便测序量达标下游信号也会非常碎建议重新处理。CUTTag 分析中不要一上来就按 ChIP-seq 的习惯做去重Tn5 切刻本身就有一定的插入偏好生物学重复之间同样的 insert 位置重复出现是正常现象去重会砍掉大量真实信号。2.4 留意线粒体 reads 比例CUTTag 的常见背景问题之一就是线粒体基因组占据大量 reads。线粒体染色质上没有核小体Tn5 更容易切所以背景高的样品质控时线粒体 reads 比例可能到 20% 甚至更高。这个问题和 spike-in 没有直接关系但会影响你对整体背景水平的判断。如果一批样本里线粒体 reads 比例差异很大建议在比对参考基因组时直接把线粒体染色体过滤掉比如 hg38 的 chrM或者用 --un 参数单独分出来统计。3. 双轮比对流程把 spike-in reads 先“捞干净”再比对参考基因组3.1 为什么不用合并索引一次性比对有些人会把 E. coli 基因组和参考基因组拼成一个组合索引所有 reads 一次比对再根据比对到的染色体名分开。这条路技术上可行但我个人在实操层面不太推荐。原因有两点容易误分配如果参考基因组上有和 E. coli 相似的低复杂度序列bowtie2 在最佳匹配的情况下会把 reads 全部压到一边导致 spike-in 比例失真。不方便检查问题双轮比对每一轮都会产生清晰的统计日志哪一步出问题一目了然。合并索引一旦出现比例异常很难快速定位是建库问题还是比对策略问题。所以这边采用“先比对 E. coli提取未比对 reads再比对参考基因组”的双轮策略。步骤多一点但能确保每一步都心中有数。3.2 第一轮比对把 reads 匹配到 spike-in 基因组bowtie2 的索引我用的是 E. coli K-12 MG1655NCBI 的 GCF_000005845.2 版本这个版本在绝大多数 CUTTag 试剂盒说明里都有对应描述。如果没有现成索引直接下载 fasta 后用 bowtie2-build 构建即可。THREADS16 REFecoli_k12_mg1655 # 构建索引如果没有的话 bowtie2-build GCF_000005845.2_ASM584v2_genomic.fna $REF # 第一轮比对到 E. coli bowtie2 --very-sensitive --end-to-end \ --no-mixed --no-discordant \ -I 10 -X 700 \ -p $THREADS \ -x $REF \ -1 sample_R1.fastq.gz \ -2 sample_R2.fastq.gz \ --un-conc-gz not_ecoli_R%.fastq.gz \ -S vs_ecoli.sam 2 vs_ecoli.log参数解释--very-sensitive让比对更严格减少错误匹配到 E. coli 的几率。对于已经确定要分开的两个基因组这个参数比较稳妥。--end-to-endCUTTag 片段短不需要局部比对整体比对更干净。--no-mixed --no-discordant只保留两个 mate 都正确配对到同一条染色体上的 readsproper pair。CUTTag 的核心信号都来自插入片段不用管嵌合比对。-I 10 -X 700插入片段长度范围。CUTTag 正常分布是 40~700 bp这个范围已经覆盖得比较充分。如果发现有很多 reads 因为 insert 过大被丢弃可以放宽到 1000。--un-conc-gz是这一步最关键的部分它会把没有比对到 E. coli 的一致双端 reads 输出为not_ecoli_R1.fastq.gz和not_ecoli_R2.fastq.gz。这里有个小坑%占位符只会出现在输出文件名中bowtie2 会把%替换成1和2这里的 1/2 对应 read1/read2 的前缀数字但不同版本可能长度不一样简单起见直接用_R%.fastq.gz这样的命名方式即可后面再利用通配符引用。3.3 从比对到 E. coli 的 SAM 中统计 spike-in 比例--un-conc-gz会帮你把未比对的 reads 提取出来但比对到 E. coli 的 reads 信息仍然保留在vs_ecoli.sam里。这一步顺手统计 spike-in 的真实 reads 数。samtools view - $THREADS -bS vs_ecoli.sam vs_ecoli.bam samtools sort - $THREADS -o vs_ecoli.sorted.bam vs_ecoli.bam samtools index vs_ecoli.sorted.bam # spike-in 比对 reads 数这里数的是双端 read 条数不是 fragment samtools view -c -F 0x904 vs_ecoli.sorted.bam-F 0x904的作用是排除 unmapped 和 secondary alignment 等不需要的记录这样统计出来的数字才是真正有效的唯一比对 reads。后续标准化就是基于这个数字来的。3.4 第二轮比对提取出的 reads 进入参考基因组用第一轮提取出的not_ecoli_R1/R2.fastq.gz进第二轮比对REF_GENOMEhg38 bowtie2 --very-sensitive --end-to-end \ --no-mixed --no-discordant \ -I 10 -X 700 \ -p $THREADS \ -x $REF_GENOME \ -1 not_ecoli_R1.fastq.gz \ -2 not_ecoli_R2.fastq.gz \ -S vs_ref.sam 2 vs_ref.log samtools view - $THREADS -bS vs_ref.sam vs_ref.bam samtools sort - $THREADS -o sample_ref.sorted.bam vs_ref.bam samtools index sample_ref.sorted.bam这里很多人会问为什么第二轮不继续用--un-conc-gz其实没有必要。能比对到 E. coli 的 reads 已经在前一轮排除掉了剩下 reads 比对到参考基因组后肯定还会有极少一部分两个 mate 都没有匹配上这些 reads 对下游分析毫无贡献丢就丢了。真正需要关注的是比对到参考基因组的 reads 中是不是还有一小部分其实是 E. coli 来源的残留。如果你前面第一轮比对用的非常严格的参数这部分残留一般不会超过 0.01%可以忽略。3.5 合并染色质质量过滤得到sample_ref.sorted.bam后下游分析前我一般会再压一道过滤去掉比对上线粒体如果是 human 的话 chrM、去掉 Q30 以下、去掉 secondary alignment。这些步骤看起来繁琐但对后续标准化影响挺大的。命令可以参考FILTchrM samtools view - $THREADS -b \ -f 2 -F 0x904 \ -q 30 \ sample_ref.sorted.bam \ -U sample_ref.chrM.bam \ -o sample_ref.filt.sorted.bam chr1 chr2 ... chr22 chrX chrY samtools index sample_ref.filt.sorted.bam-U可以把剔除的 reads 单独输出到一个文件里方便后续检查线粒体占比。检验过滤效果的一个粗暴方式是过完这步如果 bam 文件大小缩水一半甚至更多问题不算大但如果缩水到原来的 10%说明建库质量存在更大问题建议回头看一眼 FastQC 报告。4. spike-in 标准化从原始 coverage 到能直接比较的 bigWig4.1 为什么 total reads 标准化在这里不够用很多第一次接触 spike-in 数据的朋友最直接的困惑就是我直接用 CPM 或者 RPKM 不就行了吗CPM 的思路是把每个样本的 reads 总数拉平假设的是“测序量一致即代表样本可比”。但 CUTTag 的坐标系里reads 总量由两个变量决定实验效率和测序深度。spike-in 真正解决的问题是实验效率也就是不同样本之间有效细胞数、转座酶活性、抗体结合效率的综合差异。举个例子样本 A 的转座酶活性是样本 B 的两倍同时 A 的背景切刻也会是 B 的两倍。CPM 只把 reads 总数拉平并不会纠正“A 的每个真实结合位点都自带 2 倍信号放大”这个偏差。而 spike-in reads 数会跟实验效率正相关拿 spike-in counts 来算 scaling factor才能把这种系统性偏差去掉。4.2 计算每个样本的 spike-in scaling factor核心逻辑是让不同样本的 spike-in counts 对齐到同一个尺度。我用的方法是取所有样本 spike-in mapped reads 的中位数作为参考值然后每个样本的 scale factor 等于这个参考值除以它自己的 spike-in counts。以一个三样本的批次为例假设样本spike-in mapped readstotal raw reads参考值中位数scale factorWT850012,000,00065006500/8500 0.7647KO1650015,000,00065001.0KO242009,500,00065006500/4200 1.5476KO2 的 spike-in counts 只有中位数的 65%说明这个样本的实验效率更低它的信号应该被放大处理以补偿scale factor 1。WT 则相反spike-in counts 偏高说明它本身实验效率好所有值乘一个 1 的系数往下压。计算代码随手写一个 Python/R 都行懒得写就直接手算样本超过 30 个的时候才建议自动化脚本处理。4.3 生成标准化 bigWigbamCoverage 实操deeptools 的bamCoverage可以直接应用这个 scale factor 生成 bigWig。关键点是在 coverage 计算完成后只保留 scale factor不再额外做 CPM。bamCoverage -b sample_ref.filt.sorted.bam \ -o sample_spikenorm.bw \ --scaleFactor 0.7647 \ --binSize 10 \ --smoothLength 30 \ --numberOfProcessors 8 \ --extendReads--extendReads是比较重要的一项CUTTag 双端测序本质上是成对比对不需要 extend reads 到期望片段长度的步骤因为双端 reads 本身已经告诉你片段的边界。严格来说是应该去掉--extendReads的bamCoverage 在对双端 BAM 处理时本来就按 fragment 来算所以不要画蛇添足加这一项。我这边删掉bamCoverage -b sample_ref.filt.sorted.bam \ -o sample_spikenorm.bw \ --scaleFactor 0.7647 \ --binSize 10 \ --smoothLength 30 \ --numberOfProcessors 8--binSize 10和--smoothLength 30是我个人比较习惯的配置。10 bp 的 bin 足够细保留峰的形状30 bp 的平滑可以有效弱化背景噪音又不至于把窄峰磨平。想看基因组浏览器展示的话这个配置出来的图很干净。4.4 怎么验证标准化是不是生效了标准化完成后不能直接开跑下游先检查一下标准化是否真的让样本之间的背景水平对齐了。最直接的做法是选一个已知恒定表达的基因位点比如 ACTB 启动子区域或者一个“预期信号基本恒定”的区域看看不同样本在 IGV 里标准化后的 track 是否符合预期。也可以用 deeptools 的multiBigwigSummary做全局相关性检查。如果两个生物学重复在 spike-in 标准化前相关性 0.85标准化后应该涨到 0.95 左右如果降了重点看看是不是某个样本本身实验质量太差。multiBigwigSummary bins -b rep1.bw rep2.bw -p 8 -o correlation.npz plotCorrelation -in correlation.npz --method spearman --plotNumbers -o corr.png对照检查、相关性变化这些都是成本很低的验证手段建议每次跑完标准化都做一遍。数据放到审稿阶段的时候审稿人基本都会期待看到这部分“质量验证”的证据。4.5 什么时候不适合做 spike-in 标准化spike-in 标准化也不是万能的。有一种情况我遇到过如果某个样本的 spike-in counts 占总比对 reads 的比例非常低低于 0.05%这说明 E. coli 细胞几乎没有参与到反应里此时算出来的 scale factor 会异常大几十甚至上百把真实信号放大得不正常结果完全不可信。这种样本我会直接标记为失败不会硬做标准化。反过来spike-in counts 比例大于 15% 的样本也要警惕。通常是样本细胞数估算错误比如细胞大量裂解或者抗体和 E. coli 有非特异交叉反应。这些异常样本在标准化前就应该被筛选掉。5. peak calling 与信号定量有 spike-in 之后差异分析怎么做才对5.1 选 SEACR 还是 MACS2CUTTag 的峰 calling 现在比较主流的选择就两个SEACR 和 MACS2。SEACRSparse Enrichment Analysis for CUTRUN最早是给 CUTRUN 写的后来 CUTTag 社区大量使用之后也逐渐成为标配。它不依赖模型直接基于 signal 和 background 的分布来找峰特别适合 CUTTag 这种信噪比高、但背景不是 Poisson 分布的数据。如果样本里加了 spike-inSEACR 甚至可以直接用不含 spike-in 的 bam 来跑不必依赖 input。MACS2 的优势是社区大、参数熟悉、后续和很多数据库能对接。但 MACS2 的模型本来就是从 ChIP-seq 数据里学出来的直接套到 CUTTag 上如果不调整参数经常会把你引到那些宽而平的区域去。如果一定要用 MACS2建议关掉模型估计走--nomodel --shift -50 --extsize 200同时--keep-dup all不主动去重这样可以跑出比较符合 CUTTag 形态的峰来。我的默认选择是 SEACR步骤简单结果也更符合 Tn5 信号的特征。5.2 SEACR 使用要点从 BAM 转 BEDGRAPHSEACR 的输入是 BEDGRAPH或 BED格式而且必须提供 unstranded 的信号值。deeptools 的bamCoverage能产 bedgraphSEACR 官方建议用bedtools genomecov来生成bedtools genomecov -bg -trackline \ -ibam sample_ref.filt.sorted.bam \ sample.bedgraph # 跑 SEACR推荐使用 non 模式stringent 模式在低深度样本上太激进 SEACR_1.3.sh sample.bedgraph \ 0.05 \ non \ sample_seacr_peaks第二个参数0.05是 FDR 的阈值你会发现自己调得比 ChIP-seq 常见阈值更严一些才合理。non表示使用非固定阈值模式。这个阈值需要根据样本的深度和信号强度来调整我一般先跑 0.05视觉检查一个已知阳性位点如果峰太多太碎就降到 0.01如果漏峰就放宽到 0.1。5.3 差异分析用 spike-in 替代 DESeq2 自带的 size factor差异 peak 分析的标准工具链是“每个峰一个 count 矩阵 DESeq2”。问题在于DESeq2 默认的 size factor 基于全样本所有基因的中位数比值这在 CUTTag 上会严重受到有信号区域不均衡的干扰。加了 spike-in 之后正确的姿势是把自己算好的 spike-in size factor 直接喂给 DESeq2。第一步用 bedtools 对 SEACR 得到的 merged peak set 做 coverage 计数# 先合并所有样本的 peak 成一个统一区间集合 cat sample1_seacr_peaks.bed sample2_seacr_peaks.bed sample3_seacr_peaks.bed \ | sort -k1,1 -k2,2n \ | bedtools merge all_peaks.bed # 逐个样本计数 for bam in *.filt.sorted.bam; do base$(basename $bam .filt.sorted.bam) bedtools coverage -counts -a all_peaks.bed -b $bam \ ${base}.counts.txt done第二步把这些 count 合并成矩阵后进入 R# 伪代码输入 counts_matrix 和 spikein_factor 向量 library(DESeq2) counts_matrix - read.table(counts_matrix.txt, headerT, row.names1) spikein_size - c(0.7647, 1.0, 1.5476) # 前面算好的 scale factor # 直接覆盖默认 size factor dds - DESeqDataSetFromMatrix(countData counts_matrix, colData data.frame(condition c(WT,KO1,KO2)), design ~ condition) sizeFactors(dds) - spikein_size dds - DESeq(dds)三组以上的比较、配对设计DESeq2 都可以继续沿用这套逻辑只要把 size factor 自定义即可。5.4 验证峰信号的可靠性做完差异峰分析之后还有一个很重要的验证步骤把显著差异的峰列表导回到 IGV 里随机挑 5~10 个看真的峰形是否干净。SEACR 有时在低信号区域会产生一个长条形的“伪峰”如果这些边缘区域大量出现在差异列表里需要重新审视阈值或考虑增加去除黑名单区域blacklist的步骤hg38 的 blacklist BED 可以从公开资源下载。这一步没什么高深技术但极其容易发现“p value 显著但生物学荒谬”的问题强烈建议不要跳过。6. 实战中的坑清单从实验设计到生信分析的翻车点6.1 spike-in 基因组索引版本不一致E. coli 基因组有好几种常见版本K-12 MG1655、DH10B、BL21 各自序列不完全相同。你实验时用的 E. coli 菌株是哪个索引就必须对应哪个。用错一个版本可能会导致 5%~20% 的 spike-in reads 比不上而造成比例假性偏低。最稳妥的办法是直接问试剂盒厂家要 spike-in 的 reference fasta可以省掉很多猜版本的时间。6.2 线粒体 reads 被当成 spike-in 的乌龙有一回我发现一个样本的 “spike-in 比例” 高达 20%第一反应以为是实验出了问题。后来查明是参考基因组比对时没有过滤 chrM大量线粒体 reads 虽然没有比对到 E. coli但也没被及时清走等算完比例才发现是计数方式把 chrM 也算进去了。所以规范化的流程里参考基因组比对后必须把线粒体染色体及其他不需要的 contig 过滤掉。6.3 去掉 spike-in reads 后 bam 大小骤减是正常的如果 spike-in 占比 2%去掉之后 bam 文件减小 2%再正常不过。但如果你的 spike-in 比例达到 8% 以上那过滤完 bam 变小得就更夸张这不是数据丢失而是你一开始的预期就不对。spike-in counts 不是“污染”它们是特意加进去的定量标尺丢掉的时候就要明确知道自己在丢掉什么。6.4 多批次数据的批次效应在较大规模的 CUTTag 项目里几十个样本不同批次之间 spike-in 比例本身就可能漂移比如一批早期实验用的 E. coli 细胞制备批次和后面一批不同。spike-in 标准化能处理一部分技术差异但不能完全替代实验设计上的批间平衡。理想状态是每个批次的各个条件都要有生物学重复然后用批次作为协变量加入差异分析设计。如果你样本量比较大建议做一下主成分分析看看第一批和第二批样本是否明显分离成了两个 cluster如果是spike-in 标准化解决不了所有问题适当加批次信息进 model 会更稳。6.5 DNA spike-in 与染色质 spike-in 分析上的差别如果试剂盒用的是游离 DNA spike-in比如 EpiCypher SNAP-CUTANA它在建库前某一步加入不经过抗体-Tn5-转座过程所以它的 reads 量不反映“转座酶背景”只反映“建库扩增和测序效率”。这种情况下用它的 counts 做 size factor 仍能校正建库层面的偏差但不能校正实验层面的“细胞有效量”。操作流程上还是双轮比对只是解释结果时要注意DNA spike-in 标准化后的信号跨抗体、跨样本比较时依然要警惕实验处理本身带来的改变。6.6 千万别把 spike-in reads 直接塞进 peak caller这个错误我见过不止一次。有人把比对到参考基因组和 spike-in 基因组的 bam 合并在一起去 call peak结果在 E. coli 染色体位置也 call 出了一堆“峰”然后下游分析彻底混乱。spike-in reads 的唯一归宿就是标准化计数进了 bamCoverage 之后就不要再参与任何 peak calling 和差异分析了。做完全部流程之后我自己养成了一个习惯每个样本最终交付前把“总 reads → 比对 spike-in → 过滤后参考基因组比对 → peak 数量 → 显著差异峰数量”这几层数字全部留档记录。这也让我在写论文时能很清楚地展示每一步的质控结果。加了 spike-in 的 CUTTag 整个分析思路并不轻松但它换来的是你跨样本比较时可以站在一个更扎实的定量基础上。以后再拿到带 spike-in 的数据按这套流程走一遍分析结论基本能经得起审稿人和自己实验的双重检验。 CUTTag spike-in数据分析全流程：从双轮比对到标准化及避坑指南
自己做了一批加了 E. coli spike-in 的 CUTTag 数据本来以为只是“多了几步比对的活”结果处理到一半发现事情没那么简单。最开始我图省事把所有的 reads 直接一股脑比对到人类参考基因组上后面出来一堆分布在 gene body 和基因间区的“宽峰”细看很多其实是 E. coli 来源的序列被错误比对过来的。更麻烦的是不同批次样本的 signal 完全没法横向比因为根本不知道每个样本里有效细胞数差了多少。这个项目做完以后我最大的感受是加了 spike-in 的 CUTTag 分析不是“普通 CUTTag 分析 一个额外步骤”而是整套标准化逻辑都要换掉。从双轮比对的顺序到 scale factor 的计算再到 peak 定量和差异分析每一步都得把 spike-in 考虑进去。这篇就把我最后跑通、也验证过结果的完整流程捋一遍重点讲清楚为什么这样做、参数怎么定、以及我在里面踩过的几个坑。适合手里已经有加了 spike-in 的 CUTTag 数据的分析人员也适合准备做实验、但不太确定数据分析该走什么流程的朋友。1. 为什么 CUTTag 需要 spike-in没有 input 的实验拿什么当基准1.1 CUTTag 和 ChIP-seq 在定量本质上的区别做 ChIP-seq 的都知道样品组总要配一个 input用来反映“染色质本身的可及性”和“建库过程的偏好”。但 CUTTag 不存在这样的 input。它的核心是让 pA/G-Tn5 融合蛋白在抗体引导下定位到目标蛋白结合区域然后加入 Mg2 激活 Tn5直接在原位切 DNA 加上接头。这个反应发生在细胞核里没有超声打断也没有通用的“未富集对照”。所以如果你试图用普通 CUTTag 的结果除以一个 input会发现根本没有合适的 input 可用。IgG 或未加抗体的对照组确实可以拿来评估非特异背景但它只反映抗体乱结合的水平反映不了 Tn5 本身在完整染色质上的随机切刻强度更反映不了不同样本之间“有效反应概率”的差异。1.2 spike-in 细胞的原理内置的“每细胞背景标尺”我们实验室在样本制备阶段会在同一个管子里按比例掺入 E. coli 细胞然后让样本细胞和 E. coli 细胞一起经历后续所有的步骤透化、抗体孵育、二抗、pA/G-Tn5 结合、转座反应、建库扩增。E. coli 基因组上并没有你加的抗体对应的靶蛋白它产生的 reads 基本上都来自非特异性的结合和 Tn5 的随机背景切刻。这几千个 E. coli reads 是什么它们就是这条管子里“背景反应强度”的一个非常清晰可见的代理指标。同一批实验里如果某个样本细胞数少、或者破膜效果差、或者转座酶活性被杂质抑制它的背景切刻水平也会成比例地下降映射出来的就是 E. coli reads 数量变少。反过来也一样。所以 spike-in 本质上给了你一条可以通过测序直接观测的“实验效率内参线”。1.3 不同种类 spike-in 的适用场景现在市面上的 CUTTag 试剂盒里常见的 spike-in 主要有这么几类类型常见来源适用场景加入形式细菌染色质E. coli如 K12 MG1655人和小鼠等哺乳动物样本细胞混合酵母染色质S. cerevisiae植物、真菌等参考基因组比较复杂的样本细胞混合商品化 DNA spike-inEpiCypher 等用于绝对定量和批次质控游离 DNA人、小鼠样本现在最普及的还是 E. coli 染色质因为参考基因组小、序列差异大和哺乳动物基因组比对时非常容易区分干净。植物样本因为细胞壁和蛋白结构的差异很多人会换酵母。如果用的是商品化 DNA spike-in比如 SNAP-CUTANA 甲基化 panel那它是纯 DNA 形式加入不走抗体和转座酶反应所以它反映的其实是“建库与测序”的技术偏差而不是“转座酶背景”。分析逻辑上类似但不能简单拿同一个标准化公式硬套这点我会在第 6 章单独展开。1.4 spike-in 加入比例怎么定我们实验室默认按样本细胞数的 1% 加入也就是 1:100。这个比例用下来对大多数组蛋白修饰H3K27ac、H3K4me3、H3K27me3 这种都比较合适最后测序数据里 E. coli reads 占比大概在 0.5% 到 5% 之间既不会影响参考基因组上的信号又能提供足够的背景计数来做归一化。如果你做的是低丰度转录因子比如某些只在特定细胞类型表达的 TF它的 signal 占比会小很多此时可以用到 1:50 甚至 1:20 来提高 spike-in counts 的可靠性。但反过来spike-in 比例太高会挤占有效测序资源本身就是成本浪费。有一个基本判断标准每个样本最终能比对到 spike-in 的 unique reads 不要低于 1000 条否则声标准化因子非常不稳定。我之前因为细胞数量估算失误有一个样本 E. coli 只比对了 400 多条 reads那个样本后续和同组样本比较的结果明显抖动后来又补做了才救回来。2. 分析前先别急着比对数据质量检查和 spike-in 信号初判2.1 用 MultiQC 把所有 FastQC 报告汇总看一遍CUTTag 建库其实非常简单所以数据质量上最容易出问题的反而不是接头而是 GC 偏差和插入片段分布。建议拿到 demux 后的 fastq 后先用 fastqc 跑一遍再顺手 MultiQC 汇总。主要看三件事平均质量分数Q30 占比低于 75% 的样本要留意看是不是整个 run 的问题还是单个样本降解了。Adapter 残留因为 CUTTag 片段普遍较短如果测序平台没有把 adapter 完全切干净下游比对会出现大量 reads 比不过去的情况。平台保证 clean data 的情况下一般还好。GC 分布如果 GC 曲线出现明显双峰可能是混入了大量 E. coliGC 约 50.8%和人类基因组GC 约 41%的数据这恰好也是 spike-in 存在的一个直观提醒。2.2 快速估算 spike-in 读数比例正式流程前我习惯先随机抽 10 万条 reads 快速比对到 spike-in 基因组上看比例是否和预期一致。方法很简单用 seqkit sample 抽子集然后 bowtie2 比对到 E. coli 索引再数一下比对上的 reads 数。# 抽 10 万条 reads seqkit sample -2 -n 100000 -o sub.R1.fastq.gz -o sub.R2.fastq.gz sample_R1.fastq.gz sample_R2.fastq.gz # 快速比对 bowtie2 --very-sensitive --end-to-end \ -x ecoli_k12_mg1655 \ -1 sub.R1.fastq.gz -2 sub.R2.fastq.gz \ -S sub_vs_ecoli.sam 2 align_summary.txt cat align_summary.txt | grep overall alignment rate如果整体比对率在 0.5%~5% 之间说明 spike-in 加得没问题可以安心往下走。如果 10 万条里一条都没有先怀疑索引文件是不是搞错了其次怀疑实验步骤里 spike-in 是否真的加入了。这一步成本很低能省掉后面所有基于错误前提的分析。2.3 insert size 分布是判断建库质量的硬指标CUTTag 有一个特点因为转座反应发生在核小体两侧所以插入片段会呈现明显的 mono-nucleosome 模式主峰在 150~250 bp 附近另外还有一部分 sub-nucleosomal 的小片段 120 bp在双端测序的插入片段分布图上能看到一个“小峰 大峰”的组合形态。这个分布可以直接用 fastp 或者比对后 deeptools 的 bamPEFragmentSize 看。我遇到过 insert size 主峰跑到 500 bp 以上的样本后来排查发现是细胞裂解过度、染色质降解成了大片断污染。这种样本即便测序量达标下游信号也会非常碎建议重新处理。CUTTag 分析中不要一上来就按 ChIP-seq 的习惯做去重Tn5 切刻本身就有一定的插入偏好生物学重复之间同样的 insert 位置重复出现是正常现象去重会砍掉大量真实信号。2.4 留意线粒体 reads 比例CUTTag 的常见背景问题之一就是线粒体基因组占据大量 reads。线粒体染色质上没有核小体Tn5 更容易切所以背景高的样品质控时线粒体 reads 比例可能到 20% 甚至更高。这个问题和 spike-in 没有直接关系但会影响你对整体背景水平的判断。如果一批样本里线粒体 reads 比例差异很大建议在比对参考基因组时直接把线粒体染色体过滤掉比如 hg38 的 chrM或者用 --un 参数单独分出来统计。3. 双轮比对流程把 spike-in reads 先“捞干净”再比对参考基因组3.1 为什么不用合并索引一次性比对有些人会把 E. coli 基因组和参考基因组拼成一个组合索引所有 reads 一次比对再根据比对到的染色体名分开。这条路技术上可行但我个人在实操层面不太推荐。原因有两点容易误分配如果参考基因组上有和 E. coli 相似的低复杂度序列bowtie2 在最佳匹配的情况下会把 reads 全部压到一边导致 spike-in 比例失真。不方便检查问题双轮比对每一轮都会产生清晰的统计日志哪一步出问题一目了然。合并索引一旦出现比例异常很难快速定位是建库问题还是比对策略问题。所以这边采用“先比对 E. coli提取未比对 reads再比对参考基因组”的双轮策略。步骤多一点但能确保每一步都心中有数。3.2 第一轮比对把 reads 匹配到 spike-in 基因组bowtie2 的索引我用的是 E. coli K-12 MG1655NCBI 的 GCF_000005845.2 版本这个版本在绝大多数 CUTTag 试剂盒说明里都有对应描述。如果没有现成索引直接下载 fasta 后用 bowtie2-build 构建即可。THREADS16 REFecoli_k12_mg1655 # 构建索引如果没有的话 bowtie2-build GCF_000005845.2_ASM584v2_genomic.fna $REF # 第一轮比对到 E. coli bowtie2 --very-sensitive --end-to-end \ --no-mixed --no-discordant \ -I 10 -X 700 \ -p $THREADS \ -x $REF \ -1 sample_R1.fastq.gz \ -2 sample_R2.fastq.gz \ --un-conc-gz not_ecoli_R%.fastq.gz \ -S vs_ecoli.sam 2 vs_ecoli.log参数解释--very-sensitive让比对更严格减少错误匹配到 E. coli 的几率。对于已经确定要分开的两个基因组这个参数比较稳妥。--end-to-endCUTTag 片段短不需要局部比对整体比对更干净。--no-mixed --no-discordant只保留两个 mate 都正确配对到同一条染色体上的 readsproper pair。CUTTag 的核心信号都来自插入片段不用管嵌合比对。-I 10 -X 700插入片段长度范围。CUTTag 正常分布是 40~700 bp这个范围已经覆盖得比较充分。如果发现有很多 reads 因为 insert 过大被丢弃可以放宽到 1000。--un-conc-gz是这一步最关键的部分它会把没有比对到 E. coli 的一致双端 reads 输出为not_ecoli_R1.fastq.gz和not_ecoli_R2.fastq.gz。这里有个小坑%占位符只会出现在输出文件名中bowtie2 会把%替换成1和2这里的 1/2 对应 read1/read2 的前缀数字但不同版本可能长度不一样简单起见直接用_R%.fastq.gz这样的命名方式即可后面再利用通配符引用。3.3 从比对到 E. coli 的 SAM 中统计 spike-in 比例--un-conc-gz会帮你把未比对的 reads 提取出来但比对到 E. coli 的 reads 信息仍然保留在vs_ecoli.sam里。这一步顺手统计 spike-in 的真实 reads 数。samtools view - $THREADS -bS vs_ecoli.sam vs_ecoli.bam samtools sort - $THREADS -o vs_ecoli.sorted.bam vs_ecoli.bam samtools index vs_ecoli.sorted.bam # spike-in 比对 reads 数这里数的是双端 read 条数不是 fragment samtools view -c -F 0x904 vs_ecoli.sorted.bam-F 0x904的作用是排除 unmapped 和 secondary alignment 等不需要的记录这样统计出来的数字才是真正有效的唯一比对 reads。后续标准化就是基于这个数字来的。3.4 第二轮比对提取出的 reads 进入参考基因组用第一轮提取出的not_ecoli_R1/R2.fastq.gz进第二轮比对REF_GENOMEhg38 bowtie2 --very-sensitive --end-to-end \ --no-mixed --no-discordant \ -I 10 -X 700 \ -p $THREADS \ -x $REF_GENOME \ -1 not_ecoli_R1.fastq.gz \ -2 not_ecoli_R2.fastq.gz \ -S vs_ref.sam 2 vs_ref.log samtools view - $THREADS -bS vs_ref.sam vs_ref.bam samtools sort - $THREADS -o sample_ref.sorted.bam vs_ref.bam samtools index sample_ref.sorted.bam这里很多人会问为什么第二轮不继续用--un-conc-gz其实没有必要。能比对到 E. coli 的 reads 已经在前一轮排除掉了剩下 reads 比对到参考基因组后肯定还会有极少一部分两个 mate 都没有匹配上这些 reads 对下游分析毫无贡献丢就丢了。真正需要关注的是比对到参考基因组的 reads 中是不是还有一小部分其实是 E. coli 来源的残留。如果你前面第一轮比对用的非常严格的参数这部分残留一般不会超过 0.01%可以忽略。3.5 合并染色质质量过滤得到sample_ref.sorted.bam后下游分析前我一般会再压一道过滤去掉比对上线粒体如果是 human 的话 chrM、去掉 Q30 以下、去掉 secondary alignment。这些步骤看起来繁琐但对后续标准化影响挺大的。命令可以参考FILTchrM samtools view - $THREADS -b \ -f 2 -F 0x904 \ -q 30 \ sample_ref.sorted.bam \ -U sample_ref.chrM.bam \ -o sample_ref.filt.sorted.bam chr1 chr2 ... chr22 chrX chrY samtools index sample_ref.filt.sorted.bam-U可以把剔除的 reads 单独输出到一个文件里方便后续检查线粒体占比。检验过滤效果的一个粗暴方式是过完这步如果 bam 文件大小缩水一半甚至更多问题不算大但如果缩水到原来的 10%说明建库质量存在更大问题建议回头看一眼 FastQC 报告。4. spike-in 标准化从原始 coverage 到能直接比较的 bigWig4.1 为什么 total reads 标准化在这里不够用很多第一次接触 spike-in 数据的朋友最直接的困惑就是我直接用 CPM 或者 RPKM 不就行了吗CPM 的思路是把每个样本的 reads 总数拉平假设的是“测序量一致即代表样本可比”。但 CUTTag 的坐标系里reads 总量由两个变量决定实验效率和测序深度。spike-in 真正解决的问题是实验效率也就是不同样本之间有效细胞数、转座酶活性、抗体结合效率的综合差异。举个例子样本 A 的转座酶活性是样本 B 的两倍同时 A 的背景切刻也会是 B 的两倍。CPM 只把 reads 总数拉平并不会纠正“A 的每个真实结合位点都自带 2 倍信号放大”这个偏差。而 spike-in reads 数会跟实验效率正相关拿 spike-in counts 来算 scaling factor才能把这种系统性偏差去掉。4.2 计算每个样本的 spike-in scaling factor核心逻辑是让不同样本的 spike-in counts 对齐到同一个尺度。我用的方法是取所有样本 spike-in mapped reads 的中位数作为参考值然后每个样本的 scale factor 等于这个参考值除以它自己的 spike-in counts。以一个三样本的批次为例假设样本spike-in mapped readstotal raw reads参考值中位数scale factorWT850012,000,00065006500/8500 0.7647KO1650015,000,00065001.0KO242009,500,00065006500/4200 1.5476KO2 的 spike-in counts 只有中位数的 65%说明这个样本的实验效率更低它的信号应该被放大处理以补偿scale factor 1。WT 则相反spike-in counts 偏高说明它本身实验效率好所有值乘一个 1 的系数往下压。计算代码随手写一个 Python/R 都行懒得写就直接手算样本超过 30 个的时候才建议自动化脚本处理。4.3 生成标准化 bigWigbamCoverage 实操deeptools 的bamCoverage可以直接应用这个 scale factor 生成 bigWig。关键点是在 coverage 计算完成后只保留 scale factor不再额外做 CPM。bamCoverage -b sample_ref.filt.sorted.bam \ -o sample_spikenorm.bw \ --scaleFactor 0.7647 \ --binSize 10 \ --smoothLength 30 \ --numberOfProcessors 8 \ --extendReads--extendReads是比较重要的一项CUTTag 双端测序本质上是成对比对不需要 extend reads 到期望片段长度的步骤因为双端 reads 本身已经告诉你片段的边界。严格来说是应该去掉--extendReads的bamCoverage 在对双端 BAM 处理时本来就按 fragment 来算所以不要画蛇添足加这一项。我这边删掉bamCoverage -b sample_ref.filt.sorted.bam \ -o sample_spikenorm.bw \ --scaleFactor 0.7647 \ --binSize 10 \ --smoothLength 30 \ --numberOfProcessors 8--binSize 10和--smoothLength 30是我个人比较习惯的配置。10 bp 的 bin 足够细保留峰的形状30 bp 的平滑可以有效弱化背景噪音又不至于把窄峰磨平。想看基因组浏览器展示的话这个配置出来的图很干净。4.4 怎么验证标准化是不是生效了标准化完成后不能直接开跑下游先检查一下标准化是否真的让样本之间的背景水平对齐了。最直接的做法是选一个已知恒定表达的基因位点比如 ACTB 启动子区域或者一个“预期信号基本恒定”的区域看看不同样本在 IGV 里标准化后的 track 是否符合预期。也可以用 deeptools 的multiBigwigSummary做全局相关性检查。如果两个生物学重复在 spike-in 标准化前相关性 0.85标准化后应该涨到 0.95 左右如果降了重点看看是不是某个样本本身实验质量太差。multiBigwigSummary bins -b rep1.bw rep2.bw -p 8 -o correlation.npz plotCorrelation -in correlation.npz --method spearman --plotNumbers -o corr.png对照检查、相关性变化这些都是成本很低的验证手段建议每次跑完标准化都做一遍。数据放到审稿阶段的时候审稿人基本都会期待看到这部分“质量验证”的证据。4.5 什么时候不适合做 spike-in 标准化spike-in 标准化也不是万能的。有一种情况我遇到过如果某个样本的 spike-in counts 占总比对 reads 的比例非常低低于 0.05%这说明 E. coli 细胞几乎没有参与到反应里此时算出来的 scale factor 会异常大几十甚至上百把真实信号放大得不正常结果完全不可信。这种样本我会直接标记为失败不会硬做标准化。反过来spike-in counts 比例大于 15% 的样本也要警惕。通常是样本细胞数估算错误比如细胞大量裂解或者抗体和 E. coli 有非特异交叉反应。这些异常样本在标准化前就应该被筛选掉。5. peak calling 与信号定量有 spike-in 之后差异分析怎么做才对5.1 选 SEACR 还是 MACS2CUTTag 的峰 calling 现在比较主流的选择就两个SEACR 和 MACS2。SEACRSparse Enrichment Analysis for CUTRUN最早是给 CUTRUN 写的后来 CUTTag 社区大量使用之后也逐渐成为标配。它不依赖模型直接基于 signal 和 background 的分布来找峰特别适合 CUTTag 这种信噪比高、但背景不是 Poisson 分布的数据。如果样本里加了 spike-inSEACR 甚至可以直接用不含 spike-in 的 bam 来跑不必依赖 input。MACS2 的优势是社区大、参数熟悉、后续和很多数据库能对接。但 MACS2 的模型本来就是从 ChIP-seq 数据里学出来的直接套到 CUTTag 上如果不调整参数经常会把你引到那些宽而平的区域去。如果一定要用 MACS2建议关掉模型估计走--nomodel --shift -50 --extsize 200同时--keep-dup all不主动去重这样可以跑出比较符合 CUTTag 形态的峰来。我的默认选择是 SEACR步骤简单结果也更符合 Tn5 信号的特征。5.2 SEACR 使用要点从 BAM 转 BEDGRAPHSEACR 的输入是 BEDGRAPH或 BED格式而且必须提供 unstranded 的信号值。deeptools 的bamCoverage能产 bedgraphSEACR 官方建议用bedtools genomecov来生成bedtools genomecov -bg -trackline \ -ibam sample_ref.filt.sorted.bam \ sample.bedgraph # 跑 SEACR推荐使用 non 模式stringent 模式在低深度样本上太激进 SEACR_1.3.sh sample.bedgraph \ 0.05 \ non \ sample_seacr_peaks第二个参数0.05是 FDR 的阈值你会发现自己调得比 ChIP-seq 常见阈值更严一些才合理。non表示使用非固定阈值模式。这个阈值需要根据样本的深度和信号强度来调整我一般先跑 0.05视觉检查一个已知阳性位点如果峰太多太碎就降到 0.01如果漏峰就放宽到 0.1。5.3 差异分析用 spike-in 替代 DESeq2 自带的 size factor差异 peak 分析的标准工具链是“每个峰一个 count 矩阵 DESeq2”。问题在于DESeq2 默认的 size factor 基于全样本所有基因的中位数比值这在 CUTTag 上会严重受到有信号区域不均衡的干扰。加了 spike-in 之后正确的姿势是把自己算好的 spike-in size factor 直接喂给 DESeq2。第一步用 bedtools 对 SEACR 得到的 merged peak set 做 coverage 计数# 先合并所有样本的 peak 成一个统一区间集合 cat sample1_seacr_peaks.bed sample2_seacr_peaks.bed sample3_seacr_peaks.bed \ | sort -k1,1 -k2,2n \ | bedtools merge all_peaks.bed # 逐个样本计数 for bam in *.filt.sorted.bam; do base$(basename $bam .filt.sorted.bam) bedtools coverage -counts -a all_peaks.bed -b $bam \ ${base}.counts.txt done第二步把这些 count 合并成矩阵后进入 R# 伪代码输入 counts_matrix 和 spikein_factor 向量 library(DESeq2) counts_matrix - read.table(counts_matrix.txt, headerT, row.names1) spikein_size - c(0.7647, 1.0, 1.5476) # 前面算好的 scale factor # 直接覆盖默认 size factor dds - DESeqDataSetFromMatrix(countData counts_matrix, colData data.frame(condition c(WT,KO1,KO2)), design ~ condition) sizeFactors(dds) - spikein_size dds - DESeq(dds)三组以上的比较、配对设计DESeq2 都可以继续沿用这套逻辑只要把 size factor 自定义即可。5.4 验证峰信号的可靠性做完差异峰分析之后还有一个很重要的验证步骤把显著差异的峰列表导回到 IGV 里随机挑 5~10 个看真的峰形是否干净。SEACR 有时在低信号区域会产生一个长条形的“伪峰”如果这些边缘区域大量出现在差异列表里需要重新审视阈值或考虑增加去除黑名单区域blacklist的步骤hg38 的 blacklist BED 可以从公开资源下载。这一步没什么高深技术但极其容易发现“p value 显著但生物学荒谬”的问题强烈建议不要跳过。6. 实战中的坑清单从实验设计到生信分析的翻车点6.1 spike-in 基因组索引版本不一致E. coli 基因组有好几种常见版本K-12 MG1655、DH10B、BL21 各自序列不完全相同。你实验时用的 E. coli 菌株是哪个索引就必须对应哪个。用错一个版本可能会导致 5%~20% 的 spike-in reads 比不上而造成比例假性偏低。最稳妥的办法是直接问试剂盒厂家要 spike-in 的 reference fasta可以省掉很多猜版本的时间。6.2 线粒体 reads 被当成 spike-in 的乌龙有一回我发现一个样本的 “spike-in 比例” 高达 20%第一反应以为是实验出了问题。后来查明是参考基因组比对时没有过滤 chrM大量线粒体 reads 虽然没有比对到 E. coli但也没被及时清走等算完比例才发现是计数方式把 chrM 也算进去了。所以规范化的流程里参考基因组比对后必须把线粒体染色体及其他不需要的 contig 过滤掉。6.3 去掉 spike-in reads 后 bam 大小骤减是正常的如果 spike-in 占比 2%去掉之后 bam 文件减小 2%再正常不过。但如果你的 spike-in 比例达到 8% 以上那过滤完 bam 变小得就更夸张这不是数据丢失而是你一开始的预期就不对。spike-in counts 不是“污染”它们是特意加进去的定量标尺丢掉的时候就要明确知道自己在丢掉什么。6.4 多批次数据的批次效应在较大规模的 CUTTag 项目里几十个样本不同批次之间 spike-in 比例本身就可能漂移比如一批早期实验用的 E. coli 细胞制备批次和后面一批不同。spike-in 标准化能处理一部分技术差异但不能完全替代实验设计上的批间平衡。理想状态是每个批次的各个条件都要有生物学重复然后用批次作为协变量加入差异分析设计。如果你样本量比较大建议做一下主成分分析看看第一批和第二批样本是否明显分离成了两个 cluster如果是spike-in 标准化解决不了所有问题适当加批次信息进 model 会更稳。6.5 DNA spike-in 与染色质 spike-in 分析上的差别如果试剂盒用的是游离 DNA spike-in比如 EpiCypher SNAP-CUTANA它在建库前某一步加入不经过抗体-Tn5-转座过程所以它的 reads 量不反映“转座酶背景”只反映“建库扩增和测序效率”。这种情况下用它的 counts 做 size factor 仍能校正建库层面的偏差但不能校正实验层面的“细胞有效量”。操作流程上还是双轮比对只是解释结果时要注意DNA spike-in 标准化后的信号跨抗体、跨样本比较时依然要警惕实验处理本身带来的改变。6.6 千万别把 spike-in reads 直接塞进 peak caller这个错误我见过不止一次。有人把比对到参考基因组和 spike-in 基因组的 bam 合并在一起去 call peak结果在 E. coli 染色体位置也 call 出了一堆“峰”然后下游分析彻底混乱。spike-in reads 的唯一归宿就是标准化计数进了 bamCoverage 之后就不要再参与任何 peak calling 和差异分析了。做完全部流程之后我自己养成了一个习惯每个样本最终交付前把“总 reads → 比对 spike-in → 过滤后参考基因组比对 → peak 数量 → 显著差异峰数量”这几层数字全部留档记录。这也让我在写论文时能很清楚地展示每一步的质控结果。加了 spike-in 的 CUTTag 整个分析思路并不轻松但它换来的是你跨样本比较时可以站在一个更扎实的定量基础上。以后再拿到带 spike-in 的数据按这套流程走一遍分析结论基本能经得起审稿人和自己实验的双重检验。