企业官网建设流程全解析

做单细胞ATAC测序（scATAC-seq）分析，很多人以为最难的环节是后续的peak calling和motif富集。但以我处理几十批项目的经验来看，真正容易翻车的是第一步：把混在一起的测序数据按细胞拆开。测序机下机就是一堆FASTQ，你面对的是成千上万个细胞混合在一起的插入片段，还有游离DNA、背景噪音和建库错误混在当中。要从这堆数据里还原出“哪条片段来自哪个细胞”，全靠barcode——一段16个碱基左右的短序列。barcode问题解决了，单细胞拆分才算迈过第一道门槛。这篇文章是实验记录式的复盘，面向两类人：刚拿到下机数据、想搞懂Cell Ranger背后逻辑的新手，以及打算自己写流程、不依赖官方管道的分析者。读完之后，你应该能独立从原始FASTQ中提取barcode、做白名单匹配和纠错，并避开那些我踩过无数次的坑。

1. 为什么拆分单细胞的第一步是barcode

1.1 scATAC-seq到底在测什么

ATAC-seq的原理，是Tn5转座酶在染色质开放区域插入测序接头，把“哪些地方是开放的、哪些地方是被紧紧包住的”给测出来。单细胞版本的差别在于，每一个细胞核被分隔进一个油包水的GEM微滴里，在同一个微滴内完成转座、打碎和标记，然后再把所有微滴的内容混在一起建库、上机测序。

所以，一开始拿到手里的下机数据，本质上是来自不同细胞的DNA片段混在一起。切开来看，每一条read都带着两样关键信息：一是片段本身在基因组上的位置，二是这个片段来自哪个细胞。前者靠比对参考基因组解决，后者就只能靠barcode。没有barcode，所有片段都在一个大池子里，根本分不清谁是谁。

这里有个容易混淆的地方：ATAC-seq和普通ChIP-seq不一样，它没有“抗体特异地抓某段DNA”这层富集，信号天然稀疏，一个细胞里能测到的有效片段非常有限。因此，barcode拆分的准确性直接关系到你后面能看到多少真正的染色质开放信号。拆分不准，轻则多出一堆低质量细胞，重则整个样本的信号都被噪声淹没。

1.2 barcode是数据里的身份证

在10x Chromium平台上，scATAC-seq的barcode通常由16个碱基组成，位置固定在测序Read 1的5'端。每种版本的试剂盒都附带一份官方barcode白名单，里面包含了建库时可能使用的候选序列，数量通常是几十万个级别。绝大部分情况下，一条read上的16个碱基一定能在白名单里找到对应序列，这样才能知道这段DNA当初被分配给哪个细胞。

但这套机制在实际数据里并不总是那么完美。测序过程会引入错误，16个碱基中间可能发生替换、插入或者缺失；barcode周围还混着接头序列、转座酶序列，如果不知道结构，很容易把不该算进来的碱基当成barcode。我在早期处理数据时，就犯过把R1末尾的辅助序列也一起提取进barcode的错，结果白名单匹配率低得吓人。

所以，所谓“解决barcode问题”，不只是把16个碱基读出来，而是要做三件事：第一，确认barcode在原始read里的准确位置；第二，与白名单做匹配；第三，对测序错误做合理的纠错。这三步做完，拆分才算有了可靠的地基。

1.3 拆分的目标到底是什么

把单细胞数据拆开，最终是要得到一张“细胞×特征”矩阵。放到ATAC-seq里，特征就是基因组区间或峰。你拿到的除非是构建好的矩阵，否则都得先经过“原始FASTQ → barcode拆分 → 比对 → 去重 → 区间计数”这条路。

很多新手会忽略一点：barcode拆分并不是一个孤立的“取前16个碱基”操作。它决定了下游所有步骤的数据质量。如果你在barcode阶段把一些片段错误归给某个细胞，后面的峰识别就会在该细胞里引入假阳性；如果你把大量片段直接丢弃，原本含有真实细胞信息的数据就白白损失了。我见过有项目barcode错误率高，导致最后数量矩阵里一半细胞只有个位数的片段，这种数据后期就算用再高级的聚类算法也救不回来。

因此，沿着“拆分”这个目标往前推，barcode其实是整个分析流程里第一个真正的决策点。理解了这一点，再看后面的技术细节就不会觉得琐碎。

2. 从序列层面确认barcode结构

2.1 10x文库构建的隐藏逻辑

10x的scATAC文库结构并不是随意设计的，它的核心目的，是把barcode放在一个方便读取、又不会干扰基因组比对的位置上。建库时，Tn5转座体在开放染色质的DNA双链上切出缺口，同时把包含测序接头的片段连接上去。经过PCR扩增以后，每个文库片段的两端分别带有P5和P7接头，而barcode则被放在靠近P5接头那一侧。

这就是为什么测序时R1往往很短，只有二十几个碱基，而R2和R3才是真正用于基因组比对的读段。因为R1从头到尾的主要任务就是把barcode读出来。可以说，R1是“身份识别通道”，R2/R3才是“基因组定位通道”。

理解了这层逻辑，就不会再犯“把R1也拿去做全基因组比对”的错误。R1的前面十几个碱基根本不属于基因组，直接比对会降低比对率，还会污染后续的峰信号。

2.2 R1、R2、R3各自的位置价值

具体到10x标准建库，不同版本的读段组成略有差异，但总体的角色分工是固定的：

• R1：长度通常在24到26 bp左右，前16 bp是barcode，后面是转座酶序列或样本索引的一部分。它的作用是身份识别，不应该参与基因组比对。
• R2：插入片段一侧的基因组序列，是真正的比对主力之一。
• R3：插入片段另一侧的基因组序列，与R2构成配对关系，帮助确定插入片段的真实长度和位置。

实际处理时，我习惯用seqkit stats先看文件基本信息，再用seqkit sample抽个两三百对reads，肉眼检查R1的长度分布和碱基组成。如果前16 bp的GC含量和大致多样性都正常，说明数据没被建库或者上机过程搞乱，可以放心往下走。

2.3 白名单与“有效barcode”定义

白名单的获取很简单，直接去10x官方支持页面，按试剂盒版本下载对应的barcode whitelist文件。文件格式是一行一个序列，没有任何表头。使用时，把它加载成一个Python集合或者哈希表，就能在极短时间内完成海量barcode的检索。

不过，白名单匹配只是第一步。它只能告诉你“这个16 bp序列是不是官方建库时可能用过的序列”，不能直接告诉你“这个序列对应哪个细胞”。在多数情况下，一个barcode就对应一个细胞，但也有少数情况是同一个细胞被多个barcode标记，或者一个barcode对应到了两个细胞。后者在细胞鉴定环节还需要更复杂的处理，但在barcode拆分这一步，我们只需要把每条read归类到某个barcode下。

这里有一个经验值：正常情况下，原始数据里能与白名单精确匹配的reads比例应该很高。如果这一比例明显低于50%，就要回头检查是不是提取了错误的碱基窗口，或者测序质量太差，而不是贸然增加容错度。

3. 实操：把barcode从原始FASTQ中分离出来

3.1 开始前先确认read结构

实际操作时，我从来不会直接上来就写脚本。先花三分钟检查数据，能省下后面一整天的返工时间。

先看文件命名，确认哪个是R1、哪个是R2、哪个是R3。再看每个FASTQ文件里reads的长度分布。如果R1的长度不是24或26，而是和R2相同、都是50 bp，那就要怀疑是不是测序平台或建库协议做过调整，这时不能想当然地取前16 bp作为barcode。

我会用一条命令快速抽查：

seqkit stats *.fastq.gz seqkit sample -n 100 -s 42 R1.fastq.gz | seqkit fx2tab | head -20

拿到的结果里，R1序列开头16个碱基应该呈现出比较丰富的多样性，而后面的碱基则可能有某种偏向性。如果连前16 bp看起来都很整齐、缺乏多样性，那八成是文库或者文件对应关系出了问题。

3.2 用cutadapt快速去除barcode

如果你想保留R2/R3去比对，只是不想要R1里的barcode干扰，可以用cutadapt把R1的前16个碱基裁掉，然后只对R2/R3做比对。命令很简单：

cutadapt -j 8 \ -u 16 \ -o R1_noBC.fastq.gz \ R1.fastq.gz

-u 16表示从每条read的5'端切除前16个碱基。这个操作不会动R2/R3，也不会影响后续比对。但要注意，这样处理之后，你已经把barcode从R1上永久删掉了。如果后续想再回溯某个片段来自哪个细胞，就会很麻烦。所以我更推荐保留barcode信息的做法，也就是下面要写的脚本。

3.3 基于Python的barcode提取与白名单匹配

我这几年在项目里反复用的一套逻辑是这样的：同时读R1、R2、R3三个FASTQ，取R1的前16个碱基作为raw barcode，用白名单做精确匹配和纠错，通过后把barcode写到输出read的header里，同时输出R2和R3。这样既保留了身份信息，又不影响下游比对。

一个可供参考的示例脚本如下：

import gzip def read_fastq(path): with gzip.open(path, 'rt') as f: while True: name = f.readline().strip() if not name: break seq = f.readline().strip() sep = f.readline().strip() qual = f.readline().strip() yield name, seq, qual def hamming(s1, s2): return sum(c1 != c2 for c1, c2 in zip(s1, s2)) def correct_barcode(raw, whitelist): if raw in whitelist: return raw candidates = [] for wb in whitelist: if hamming(raw, wb) == 1: candidates.append(wb) if len(candidates) == 1: return candidates[0] return None # 加载白名单 with open('barcode_whitelist.txt') as f: whitelist = {line.strip() for line in f} out_r2 = gzip.open('R2_withBC.fastq.gz', 'wt') out_r3 = gzip.open('R3_withBC.fastq.gz', 'wt') count = 0 kept = 0 for r1, r2, r3 in zip(read_fastq('R1.fastq.gz'), read_fastq('R2.fastq.gz'), read_fastq('R3.fastq.gz')): raw_bc = r1[1][:16] cb = correct_barcode(raw_bc, whitelist) count += 1 if cb is None: continue kept += 1 header_r2 = f"{r2[0]} BC:Z:{cb}" header_r3 = f"{r3[0]} BC:Z:{cb}" out_r2.write(f"{header_r2}\n{r2[1]}\n+\n{r2[2]}\n") out_r3.write(f"{header_r3}\n{r3[1]}\n+\n{r3[2]}\n") out_r2.close() out_r3.close() print(f"total: {count}, kept: {kept}, ratio: {kept/count:.2%}")

这段代码有几个关键点。第一，白名单用集合存放，匹配速度是O(1)级别。第二，纠错逻辑只允许唯一候选，如果一条raw barcode同时和两个白名单序列的距离都是1，就选择丢弃，而不是二选一。第三，整个脚本可以流式处理，内存占用很低，非常适合动辄几百G的原始数据。

3.4 读段是否要过滤低质量barcode

这里我踩过几次坑，想单独拿出来说。barcode只有16个碱基，其中任何一个碱基质量差，都可能导致白名单匹配失败。但如果在提取阶段就把整条read扔掉，又会损失后面R2/R3的有效信息。

我的建议是，不要在barcode阶段做太严格的质量过滤。更稳妥的方法是，先看raw barcode是否能精确匹配白名单，不行再允许1个碱基替换；如果替换后也无法匹配，说明可能是低质量碱基导致的，这时再看barcode对应的质量值，如果质量值确实很低，可以把它当作无法纠错的读段丢弃。这样相当于给数据多一次纠错机会，能多挽回不少有效片段。

4. 从barcode到单细胞矩阵的完整流程

4.1 比对与去重

拿到带barcode标签的R2/R3后，下一步是比对到参考基因组。我用得比较多的是bwa-mem2或者minimap2，两者对ATAC-seq这类短片段数据都处理得不错。比对命令示例：

bwa-mem2 mem -t 16 -M \ reference.fa \ R2_withBC.fastq.gz \ R3_withBC.fastq.gz \ | samtools view -bS - \ | samtools sort -@ 16 -o aligned.sorted.bam

这里有个重要的点：ATAC-seq的片段通常比较短，比对时要允许软剪裁，否则Tn5插入位点附近的少量错配会导致比对失败。-M参数是为了兼容下游标记重复，建议保留。

去重步骤，我通常用picard MarkDuplicates或者samtools markdup。ATAC-seq文库经过PCR扩增，同一个插入片段会产生多个拷贝，不去重的话，一个细胞在某个峰上的read数会被严重放大。结合barcode信息去重的方式，和普通WGS去重类似，但必须按barcode分组进行，否则会把不同细胞里天然相同的插入片段误判为重复。

4.2 生成细胞×峰矩阵

比对和去重完成之后，就可以生成矩阵了。我不建议对每一个barcode单独call peak，那样既慢又不稳定。更常用的策略是：把所有barcode的reads合在一起，用MACS2在全局上识别一组峰，然后把每个barcode在每条峰里的插入片段数统计出来。

这里有一个细节：统计时不应该只看read数，而应该以Tn5的切割位点为准。Tn5在插入位点会留下5'端偏移，处理时通常会把每条read转换成两个切割位点，再用一个固定窗口扩展来计算覆盖度。这样得到的计数更能反映真实的开放染色质信号。

实际工作中，也可以直接用ArchR或SnapATAC2这类R包，它们自带从BAM到矩阵的转换流程，只要你在BAM里保留好了CB标签，下游处理会很顺畅。

4.3 鉴别真细胞和空液滴

生成矩阵后，还有一道关卡：区分空液滴和真实细胞。10x平台上，一个建库反应会产生大量没有细胞核的GEM空微滴，这些空微滴也会带barcode，也可能产生少量reads。如果不剔除，最终矩阵里会出现一大批“幽灵细胞”。

通常的做法是画一个barcode的排序图：横轴是barcode按总reads数降序排列的序号，纵轴是总reads数。真实细胞对应的barcode会形成一个明显的拐点，拐点右边平台期或者悬崖式下降的部分，就是空液滴。手动定阈值也好，用DropletUtils这类工具也好，这一步一定要做，但不能在barcode拆分阶段就动手，否则会把低质量的真实细胞也误删。

我自己的经验是：宁可多留一些可疑barcode，也不要一开始就把阈值卡得太死。因为下游的聚类和双细胞鉴定还能进一步清理，但如果在源头把细胞弄丢了，就再也找不回来了。

5. 常见问题与踩坑记录

5.1 barcode在R1还是R2里？

这个问题听起来很基础，但我真的见过不止一次。10x标准建库，barcode一定在R1；但某些第三方平台或者自建流程，会把barcode放在index读段里，甚至放在R2的开头。拿到数据后先看实验记录或平台说明，不要想当然地套10x规则。

我遇到过一位同事，直接把非10x的某平台数据拿来做白名单比对，结果匹配率不到3%。原因就是该平台把barcode放在了index读段里，而R1只是普通基因组读段。后来重新按平台文档提取barcode，匹配率立刻恢复正常。

5.2 明明白名单匹配率很高，但拆分后细胞数很少

这种情况往往不是barcode的问题，而是后面的细胞筛选阈值卡得太狠。比如总reads数很低的barcode也被保留，但到细胞鉴定时被当作空液滴给剔掉了。另外一个可能原因是文库复杂度低、PCR重复率很高，导致去重后每个细胞的unique fragment太少。

排查思路是先看几个质控指标：reads比对率、去重后fragment数、barcode多样性。如果去重后fragment数中位数本来就不到几百，那问题多半在建库质量，而不是barcode拆分流程本身。

5.3 允许的编辑距离设成多少合适

我的默认值是1，也就是只允许1个碱基的替换。只有当数据质量非常差、且白名单匹配率明显偏低，并且有实验记录证明建库过程正常时，才会考虑放宽到2。但一定要记住，编辑距离越大，错误分配的风险越高。一个16 bp的barcode，如果允许3个碱基的错误，两个不同barcode之间距离可能只有4到5个碱基，那就会出现严重的串扰。

为了验证纠错是否合理，我通常会随机抽取几十条被纠错的read，人工看raw barcode和纠正后的barcode，确认只差一个碱基，而且该碱基的质量分通常较低。这能帮助判断纠错到底是在修测序错误，还是在乱分数据。

5.4 同一份数据两次分析结果不稳定

如果你用自己脚本处理数据，两次运行结果却不一样，先检查两点：一是白名单加载时是否用了set，而set的迭代顺序不固定，导致纠错时“唯一候选”的判断受影响；二是是否有随机抽样或者随机种子没有固定。前者看似无关紧要，但在极端情况下会让个别barcode流向不同细胞。

解决方法是，在所有可能影响结果的地方都固定随机种子，并且在输出文件名或者日志里记录参数哈希值。这样即使后面对比时发现差异，也知道是参数变了还是代码变了。

6. 自己写脚本还是用现成工具

我经常被问到一个问题：是不是非要自己写脚本？我的回答是，看目标是什么。

如果你只需要一个可靠的结果，直接用Cell Ranger ATAC或者新版Cell Ranger的ATAC流程就好。它会自动完成barcode提取、白名单匹配、纠错、比对、去重、矩阵生成，而且官方经过大量样本验证，稳定性很高。你只需要准备好参考基因组和原始FASTQ，跑一条count命令即可。

但如果你是做方法学的，或者想对数据分析有完全的控制权，那就值得自己写一遍。自己写流程的最大好处，是你能明确知道每一步发生了什么，遇到异常时可以快速定位是barcode问题、比对问题还是计数问题。坏处是工作量和维护成本都不小，尤其是barcode纠错这步，要考虑性能优化和边界情况。

无论选哪条路，我都建议先做一个小样本测试。抽取两三万对reads，跑通整条流程，检查每一步的输出是否符合预期，然后再全量运行。这个小习惯帮我省下了无数次全量重跑的麻烦。

最后说一个我自己的实操体验：拆分完成后，不要急着往下游冲。随手写一段小脚本，从每个barcode里抽出少量reads，去比对结果里看它们的比对位置是否合理，以及不同barcode的片段是否在基因组上呈现出互不干扰的分布。这一步只需要几分钟，却能提前发现大量潜在问题。毕竟，barcode拆分是整个单细胞数据分析的承重墙，地基稳了，后面才敢放心盖楼。