做CUTTAG实验的同行应该都有同感这种技术从原理到操作都让人舒服细胞量只要几千到几万个信噪比高建库还快。但真正让数据分析头疼的问题之一就是样品之间的标准化。不同批次的细胞状态、抗体浓度、Tn5酶活性、PCR循环数任何一个环节稍微抖一下最后比对出来的reads数就会有明显差异。如果只按总的mapped reads去均一化很容易把真实的生物学差异给淹没掉。所以现在越来越多的CUTTAG实验会在体系里加入spike-in——最常见的就是E. coli DNA——用来作为系统效率的内参把批次间、样品间的技术波动校正回来。这篇就把加了spike-in的CUTTAG样品从下机数据到标准化BigWig、再到peak calling的完整分析流程捋一遍里面不少坑都是我自己踩过的。1. 为什么CUTTAG需要spike-in1.1 CUTTAG的原理回顾与标准化的特殊性CUTTagCleavage Under Targets and Tagmentation本质上是用Protein A/G-Tn5融合蛋白做“靶向转座”。抗体结合到目标蛋白比如H3K27ac、H3K4me3这些组蛋白修饰上之后pA-Tn5被招募到目标位置然后加入Mg2激活Tn5在目标DNA附近进行切割同时直接把测序接头连接到DNA片段上。整个流程比ChIP-seq少了很多步骤背景低重复性好还能在单细胞级别展开。但正是因为它把“片段化”和“建库”合并成一步引入了很多变量。比如Tn5酶的活性在不同批次试剂之间会有波动同一批实验中不同管的细胞数量可能偏差10%到20%抗体结合效率也受样品质量影响。如果只按比对到参考基因组的总reads数去做文库间均一化问题就来了样品A的抗体效率稍低背景reads比例高样品B的抗体效率正常目标区域reads占比高。两者总reads差不多但真正的信号强度可能差了1倍以上。这就是为什么CUTTAG需要外源spike-in做内参——它能告诉你每一管反应里“有效CUTTAG信号”到底有多少而不是单纯看测序深度。需要特别提醒的一点是spike-in标准化的本质是用外源已知量DNA的信号强度去反映整个实验系统的“传递效率”。这个效率包括Tn5切割效率、接头连接效率、PCR扩增效率、测序质量等等。和RNA-seq里用ERCC spike-in校正转录本绝对丰度不同CUTTAG的E. coli spike-in主要用于相对校正让不同样品之间可比较。它校正的不是绝对的分子数而是batch effect和文库间系统性差异。1.2 spike-in方案选型与加样量估算目前CUTTAG常见的spike-in方案主要有三种我直接做个表对比一下。spike-in类型原理优点缺点适用场景E. coli基因组DNA片段外源裸DNA与染色质一起经历Tn5切割和建库便宜、易操作、不依赖物种交叉反应没有核小体结构只能反映Tn5活性和回收效率无法反映抗体效率常规CUTTAG样品标准化含Flag标签的外源组蛋白通过抗Flag抗体捕获外源核小体可以同时评估抗体和Tn5效率需要构建或购买工具细胞系成本高抗体批次质控、多组学联合实验Drosophila S2细胞物种内参有核小体结构提供完整核小体底物接近真实染色质需要抗体有物种交叉反应细胞培养成本高跨物种比较或复杂实验设计其中E. coli DNA方案用的人最多原因很简单操作门槛低结果稳定。E. coli K12 MG1655基因组大约4.6 Mb序列和人类基因组没有明显同源性比对时几乎不会串扰。不过这个方案有一个局限E. coli DNA是裸露DNA没有核小体所以pA-Tn5在抗体引导下不会特异性地在它上面切割它能覆盖的是“Tn5转座活性建库回收效率”。换句话说它无法校正抗体结合这一步的效率。如果你的实验目的纯粹是比较同一种抗体在不同样品之间的组蛋白修饰信号E. coli spike-in已经足够如果你要评估不同抗体之间的效率差异那就要考虑组蛋白级别的spike-in了。关于加样量很多人第一次做都会纠结。一个经验值是这样的加E. coli DNA的量大约占样品总DNA量的1%左右。以10万个HeLa细胞为例一个二倍体人类细胞基因组约6 pg10万细胞总DNA约60 ng那就加600 pg的E. coli DNA。如果按测序量算假设一个样品最终测20M mapped readsE. coli spike-in比例在1%左右也就是大约200K reads用于E. coli比对这个数量足够计算稳定的标准化因子了。加太少了比如0.05%E. coli reads只有1万左右波动很大统计上不稳定加太多了比如5%以上白白浪费大量测序量目标基因组的数据量会不够。一般建议控制在0.5%到2%之间。操作上还有几个小细节值得注意E. coli DNA需要用超声打断到200到600 bp左右的片段接近染色质片段长度。直接用完整基因组DNA或过长片段会导致Tn5切割不均和PCR扩增偏好最后spike-in reads的片段分布会偏离目标片段。另外E. coli DNA要在实验最开始就加入细胞悬液充分吹打混匀。这一步极其关键如果混合不均匀同一个样品不同技术重复之间的E. coli reads差异可能超过2倍后面标准化算出来的因子就是错的。2. 下机数据的预处理与双重比对2.1 质控与接头去除拿到下机数据之后第一步永远是看一下原始质量。CUTTAG下机数据一般是双端150 bp的FASTQ文件建议先跑FastQC检查一下base quality、GC含量、接头比例、duplication level这些指标。这里有个需要注意的地方CUTTAG的文库在Tn5转座之后插入片段大多集中在核小体区域加上两端的测序接头实际测序的时候经常会遇到“短插入片段”的情况。也就是说两条reads有可能会重叠测到一半就已经读到接头上去了。所以FastQC里overrepresented sequences那一栏如果出现比较高的接头序列不要太慌这在CUTTAG里常见。接着用cutadapt去掉接头。Tn5转座建库使用的是标准的Illumina TruSeq接头命令可以这样写cutadapt \ -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCA \ -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT \ -o sample_trim_R1.fastq.gz \ -p sample_trim_R2.fastq.gz \ -m 20 \ -q 20 \ sample_R1.fastq.gz sample_R2.fastq.gz这里-m 20表示过滤掉修剪后长度小于20 bp的reads-q 20做末端质量修剪。CUTTAG的片段比较短阈值设得不要太激进否则会损失大量有效数据。之后我通常还会用FastQC再确认一下修剪效果重点看adapter content是否降下去了。2.2 目标基因组与spike-in基因组的比对CUTTAG的比对推荐用bowtie2而且建议加--local --very-sensitive参数。原因在于Tn5转座本身有序列偏好切割位点附近的碱基不一定能和参考基因组完美匹配加上CUTTAG的插入片段短局部比对--local能容忍reads两端的soft clipping避免因为末端一个碱基的错配导致整条read被丢弃。比对策略有两种一种是只比对到目标基因组比如hg38然后把unmapped的reads再比对到E. coli基因组另一种是同时比对到将两个基因组拼接起来的混合参考索引。我更推荐前者分开比对、分开统计逻辑清楚也能避免E. coli序列与人类基因组中的线粒体、核糖体DNA区域产生潜在的多重比对干扰。实际操作流程是这样的# 建立索引人类基因组和E. coli K12 MG1655基因组 bowtie2-build hg38.fa hg38_index bowtie2-build ecoli_K12_MG1655.fa ecoli_index # 比对到人类基因组 bowtie2 --local --very-sensitive -p 16 \ -x hg38_index \ -1 sample_trim_R1.fastq.gz -2 sample_trim_R2.fastq.gz \ | samtools view -bS -F 4 - \ | samtools sort - 8 -o sample.hg38.sorted.bam # 同样的reads比对到E. coli基因组 bowtie2 --local --very-sensitive -p 16 \ -x ecoli_index \ -1 sample_trim_R1.fastq.gz -2 sample_trim_R2.fastq.gz \ | samtools view -bS -F 4 - \ | samtools sort - 8 -o sample.ecoli.sorted.bam比对完成后统计E. coli比对reads数samtools view -c sample.ecoli.sorted.bam同时也统计人类基因组的比对情况samtools view -c sample.hg38.sorted.bam这里要补充一个过滤步骤。建议在统计spike-in和下游分析之前先过滤掉低质量比对和多重比对reads。CUTTAG中推荐过滤flag 1804也就是排除未比对、次要比对、QC失败、重复、配对不完整这些readssamtools view -b -F 1804 -q 30 sample.hg38.sorted.bam sample.hg38.filtered.bam不过在实际操作中有人会省略重复去除这一步因为CUTTAG的Tn5转座反应本质上是“低起始量”的建库duplicate rate高并不一定代表文库质量差。我个人的习惯是先保留所有reads做完spike-in标准化和可视化如果做peak calling或定量分析时发现重复率非常高再用Picard MarkDuplicates评估一下。CUTTAG文库中较高等的duplicate通常来自Tn5在同一个位置上的连续切割不是真正的PCR重复所以果断去除反而会损失信号。比对统计出来后通常可以快速判断实验质量。人类基因组比对率正常在70%到90%之间E. coli spike-in比对率根据加样量不同而变化一般在0.5%到2%。如果E. coli比对率异常低或者重复间差异特别大就回到第1节提到的那些操作细节去排查。3. spike-in标准化的完整实操3.1 标准化因子的计算逻辑与R实操spike-in标准化的核心逻辑一句话就能讲清楚相同投入量的E. coli DNA在哪个样品里被“做出来”的reads多说明哪个样品的实验系统效率高那么真正目标区域的信号在这个样品里也会被同比率放大。所以要用每个样品的E. coli reads数反推一个scale factor把所有样品归一化到同一个基准上。实际操作中基准值一般取所有样品E. coli reads的中位数或最小值然后每个样品的scale factor 基准值 / 该样品E. coli reads数。因子大于1说明这个样品系统效率偏低需要放大信号因子小于1说明系统效率偏高需要压缩信号。下面用R做一次完整的计算示例library(dplyr) sample_info - data.frame( sample c(WT_rep1, WT_rep2, KO_rep1, KO_rep2), ecoli_reads c(120000, 80000, 60000, 150000), hg38_reads c(18e6, 20e6, 16e6, 17e6) ) # 取中位数作为基准 ref - median(sample_info$ecoli_reads) # 计算标准化因子 sample_info - sample_info %% mutate(scale_factor ref / ecoli_reads) sample_info结果会是sampleecoli_readshg38_readsscale_factorWT_rep112000018e60.75WT_rep28000020e61.125KO_rep16000016e61.5KO_rep215000017e60.6这里中位数是100000。KO_rep1的E. coli reads只有60000说明它的系统效率最低scale factor为1.5也就是要把它的信号放大1.5倍才能和其他样品站在同一条起跑线上。WT_rep1的E. coli reads是120000系统效率偏高因子0.75需要压缩信号。这里有一个很多新手会犯的错误在算完spike-in scale factor之后又跑去叠加RPKM或者CPM标准化。其实这是多余的甚至是有害的。spike-in标准化已经校正了文库间系统差异这时候再除以文库大小相当于又引入了一遍“总reads数一致”假设反而会把真实差异抹掉。CUTTAG中目标区域reads占总reads的比例本来就很低且不稳定RPKM这种基于总文库大小的标准化方式会让信号被背景稀释。所以正确的做法是用了spike-in就直接用scale factor生成标准化后的信号文件不要再叠其他标准化。3.2 从BAM到spike-in标准化BigWig拿到scale factor之后下一步就是把过滤后的BAM文件生成spike-in标准化的BigWig用来做IGV可视化、计算相关性、差异分析、以及作为部分peak caller的输入。最常用的是deepTools里的bamCoveragebamCoverage \ --bam KO_rep1.hg38.filtered.sorted.bam \ --outFileName KO_rep1.spikein.bigWig \ --scaleFactor 1.5 \ --binSize 10 \ --smoothLength 30 \ --effectiveGenomeSize 2913022398 \ -p 8这里解释几个关键参数。--scaleFactor就是我们上一步算出来的spike-in因子。--binSize因为是CUTTAG这种定点富集的信号我习惯用10 bp的窗口。--smoothLength是平滑窗口设成30 bp相当于3个bin能在不损失太多分辨率的前提下降低噪音。--effectiveGenomeSize是对应物种的可比对基因组大小人类hg38的常规值就是2913022398这个值主要用于depth归一化但注意加上它之后bamCoverage默认会做RPGCreads per genome coverage标准化如果不想叠加就不要加这个参数或者显式设置--normalizeUsing None。我们这里已经用scaleFactor做spike-in标准化了所以不需要再指定--normalizeUsing只传scaleFactor即可。生成之后可以用deepTools的bigWigCompare做一下处理组和控制组的log2 ratio这样能直观看到差异区域bigWigCompare \ --bigwig1 KO_rep1.spikein.bigWig \ --bigwig2 WT_rep1.spikein.bigWig \ --operation log2 \ --outFileName KO_vs_WT_log2ratio.bw \ --binSize 10 -p 8在IGV里加载标准化前后的BigWig你会发现标准化后的信号量级更能反映真实的生物学差异。比如KO组如果一开始因为系统效率低导致总体信号很低标准化后会和WT组处于同一个量级这时再比较哪个位点真正上调或下调才有意义。3.3 Peak calling时的标准化衔接Peak calling这步也要注意不是所有工具都会自动读取你的scale factor。CUTTAG最常用的两个peak caller是SEACR和MACS2两者在spike-in标准化上的处理方式不同。SEACR是专为CUTRUN/CUTTAG这种稀疏富集信号开发的peak caller它的输入是bedgraph格式而且SEACR不会自己做文库大小归一化所以你应该把spike-in标准化的bedgraph直接喂给它。生成bedgraph可以用bedtools genomecovbedtools genomecov -bg -ibam KO_rep1.hg38.filtered.sorted.bam KO_rep1.fragments.bedgraph然后运行SEACRSEACR_1.3.sh KO_rep1.fragments.bedgraph norm stringent KO_rep1_seacr这里的norm参数表示使用标准化信号stringent和relaxed是两种阈值模式。我一般两个都会跑然后挑stringent结果做后续分析用relaxed结果做补充。如果你用MACS2情况稍微麻烦一点。MACS2默认假设Poisson分布会使用整个基因组的背景信号来估计阈值这跟CUTTAG这种低背景高富集的模式不太匹配容易产生很多假阳性。这时可以在callpeak时加上--scale-to参数或者通过--SPMR模式结合外部factor来缩放但总体体验不如SEACR顺手。我自己现在基本只用SEACR做CUTTAG的peak callingMACS2仅用于一些需要和ChIP-seq结果直接比较的场景。如果是多组比较比如WT vs KO每个组几个重复我建议的流程是每个重复样本分别做spike-in标准化、分别call peaks然后取重复间的confidence peaks比如至少2个重复共有再对每个样本用featureCounts统计peak区域的reads count最后做DESeq2差异分析。注意这里DESeq2要加spike-in的scale factor作为normalization factor而不是用它默认的median-of-ratios。具体操作是先根据3.1计算的scale factor构造一个matrix传给DESeq2的normalizationFactors参数。这一步很多人会漏掉导致spike-in标准化做了等于白做。4. 常见问题与排查技巧4.1 spike-in读数异常的排查方向做了这么多次CUTTAG最常碰到的异常情况就是E. coli spike-in的reads比例不在预期范围内。我整理了一个问题速查表方便大家对着排查。现象可能原因排查方向与建议E. coli比对率远低于0.5%spike-in加样量不足或加样后未混匀检查加样记录下次操作时提前将E. coli DNA稀释到合适浓度加大体积加入细胞悬液并充分吹打混匀E. coli比对率过高超过5%spike-in加得太多或目标细胞数远低于预期计算一下细胞实际数量必要时降低spike-in加量过高的spike-in会挤占目标数据的测序量技术重复间E. coli reads差异超过2倍spike-in混合不均匀、移液误差、Tween浓度不均这个最常见。建议将E. coli DNA先稀释到每个样本取用的体积在5到10 μL之间再加到细胞悬液中减少移液误差E. coli reads高但目标基因组比对低细胞降解、抗体失效、Tn5活性低检查目标基因组比对率是否同时偏低如果只是E. coli高大概率是细胞起始量不足E. coli DNA没有对应的染色质“竞争”导致相对比例偏高加样后所有样本的E. coli reads都极低Tn5酶失活、转座反应步骤漏加Mg2这个属于实验事故需要回看实验记录建议做正对照如H3K4me3/H3K27ac确认整个流程4.2 标准化后质量评估的几个关键指标费劲做完了spike-in标准化怎么知道结果到底靠不靠谱我有几个百试百灵的检查手段推荐大家养成习惯。第一看spike-in标准化前后的重复样本相关性。在标准化之前同组两个重复的pearson相关性如果只有0.85左右标准化之后应该能升到0.9以上如果标准化后相关性反而下降了大概率是scale factor不对回查E. coli reads的统计是否准确。第二看已知阳性区域的信号富集程度。比如H3K4me3的CUTTAGIGV里应该看到非常典型的启动子附近峰或者用deepTools的plotProfile看所有基因转录起始位点TSS附近是否出现明显的“峰谷峰”模式。如果标准化后这些已知模式没有出现而总信号还特别高可能是因为背景区域的信号也被同步放大了——这种情况说明样品本身质量就有问题单纯标准化救不回来。第三检查FRiPFraction of Reads in Peaks评分。CUTTAG的FRiP通常比ChIP-seq高得多在0.3到0.6甚至更高都很正常。如果FRiP低于0.1要么是peak calling参数不对要么是抗体选择有问题。这里可以用featureCounts计算reads数后自己除一下也可以直接用deeptools的plotEnrichment做一个快速评估。第四还有一个容易被忽视的点加了spike-in之后总reads的组成里有一部分会被E. coli占用。比如E. coli占1%你最终目标数据量其实只有99%。所以在设计测序量时要把这部分预留出来每样本需要20M有效reads的话实际测序量要到20.2M以上。如果E. coli占5%还不自知目标数据读不出来后面会非常尴尬。最后再分享一个我自己常踩的坑spike-in标准化因子不需要做得过于精细。有的代码教程会基于每个位点的E. coli read depth做分bin回归或者用复杂的multidimensional scaling去估计因子。对于绝大多数CUTTAG项目直接用每个样本E. coli比对reads数算一个scalar factor就够了。复杂模型在重复数很少比如每组只有2个重复的情况下反而容易过拟合让你误判样品之间的真实差异。老老实实把混匀做好、把比对统计做对比花哨的算法管用得多。