直接上干货。手里刚拿到一批scATAC-seq的原始下机数据第一件事不是急着跑peak calling而是先把“一个细胞对应哪些片段”这件事彻底搞清楚。scATAC-seq数据分析里barcode拆分是真正意义上的第一步这一步如果糊弄过去后面的细胞聚类、peak calling、motif富集全部都会带着噪声跑偏。这篇实验记录就聚焦在最核心的问题上如何把scATAC-seq数据按barcode拆分到单个细胞以及在解决barcode问题时我踩过的坑、试过的方法和最终落地的方案。适合刚接触单细胞ATAC数据分析、被barcode处理绕晕的初学者也适合已经跑通流程但想搞清楚底层逻辑的同行。1. 先搞清楚barcode在scATAC-seq里到底扮演什么角色1.1 一张图看懂scATAC-seq的数据结构很多人第一次接触scATAC-seq数据时会下意识地把它当成“单细胞版的ATAC-seq”然后用bulk ATAC-seq的思路去处理结果越走越偏。这里有一个本质区别bulk ATAC-seq的建库过程没有细胞身份信息所有细胞的染色质开放信息混在一起而scATAC-seq在建库时每个细胞会被分配到一个独一无二的barcode序列这个barcode就是细胞的“身份证号”。10X Genomics平台是目前最主流的scATAC-seq平台它的GemCode技术把每个细胞包裹在单独的凝胶微滴里然后在这个微滴内完成转座和加barcode的步骤。下机数据里每个测序片段的两端都带有这个barcode信息要么在Read 1里要么在Index序列里。所以拆分barcode的过程本质上就是“按身份证号把混在一起的所有片段重新归类到各自的主人名下”。1.2 barcode拆分的三个层次一开始我以为barcode拆分就是把序列比对后按tag一筛就完事了实际做下来才发现这个问题至少有三个层次处理深度不同下游分析质量差异非常大。第一个层次是原始FASTQ层面的拆分。也就是说在比对之前就根据barcode序列把reads分到不同的文件里。这种方式适合自己手动建库、barcode序列已知的实验效率高但对barcode的容错性要求很高一个碱基的错误就会导致大量reads被错分或丢弃。第二个层次是BAM层面的拆分。比对完成后barcode信息以标签tag的形式存在每条read的BAM记录里比如10X数据最终会输出CBcell barcode标签。这种方式的好处是可以结合比对信息来过滤低质量比对然后在拆分的同时做去重、过滤等操作是目前比较主流的处理方式。第三个层次是矩阵层面的拆分也就是最终生成细胞-峰矩阵或细胞-片段矩阵。这个层次已经不单纯是“拆分”了而是将barcode与生物学信息关联起来确认哪些barcode真的对应一个细胞哪些只是背景噪声或双细胞。我在这篇记录里重点讲前两个层次因为这是解决barcode问题的核心战场。第三个层次和细胞数鉴定有关会顺带提一部分。1.3 拆分方案选型为什么我没有一上来就用cellranger-atac先说结论如果你用的是10X平台的scATAC-seq数据最省心、最推荐的方式确实是用Cell Ranger ATAC现在叫Cell Ranger ARC或cellranger-atac来做拆分和比对它对barcode的处理逻辑是封装好的、经过大规模数据验证的。但我这次为什么没有直接无脑跑cellranger-atac因为我的这批数据来源比较杂有一部分是老搭档用改进的微流控平台做的实验barcode序列不是标准10X的16bp结构还混入了一些接头污染。这种数据直接喂给cellranger-atac会在barcode识别环节就丢掉大量有效reads甚至可能出现barcode误判。所以我的策略是先理解barcode的物理结构然后自己做一轮FASTQ层面的barcode提取和校正最后再决定是走cellranger-atac的完整流程还是用sinto等工具做轻量化的拆分。后面我会详细展开这两种路线的实操过程。2. barcode处理前必须搞清楚的四个关键问题2.1 barcode到底长什么样、在哪条read里很多人拿到数据后连barcode在哪都不知道就开始处理这是最要命的。10X scATAC-seq的barcode序列在不同版本的试剂盒里位置并不完全一样。在10X Chromium Single Cell ATAC v1试剂盒中barcode序列位于Read 1的起始部分长度是16bp紧接着的是转座酶识别位点序列。Read 2则是测序得到的基因组DNA插入片段序列。i7 index用于样本 multiplexing不是细胞barcode。但要注意的是如果你用的是10X Next GEM ATAC v2试剂盒结构上会有一些调整。所以动手之前第一步一定是要去看文库构建说明书或者直接看下机数据的fastq文件开头序列对照一下barcode的位置。这一步不需要什么高级工具直接zcat看几条reads就能确认。这里分享一个我常用的快速确认方法把R1的fastq文件随机抽几千条reads统计前30bp每个位置的碱基分布。如果前16bp位置的碱基分布比较均匀、没有明显的偏好性而后面的位置出现明显的接头序列特征那基本可以确定barcode就在前16bp。2.2 一管数据里到底“藏着”多少个细胞这是barcode拆分前最容易被忽略、但又最影响拆分策略的问题。标准的10X单细胞实验通常上样目标是回收5000到10000个细胞但由于细胞悬液质量、微流控芯片上样效率等因素实际回收的细胞数和目标值会有出入。我在处理数据时习惯先跑一轮非常快速的barcode计数也就是把所有reads按barcode序列做极简的统计看一眼barcode种类的数量分布。正常情况下你会看到少数barcode的reads数极高几万到几十万中间大量的barcode reads数从中等到偏低然后就是一条长长的尾巴——这些尾巴就是背景噪声、环境DNA或测序错误导致的“假barcode”。这个分布曲线直接决定了你的拆分阈值怎么设定。如果细胞数在预期范围内曲线呈典型的长尾分布说明实验质量不错。如果高reads数的barcode特别多可能说明上样量偏大或存在细胞团如果几乎全是低reads数的barcode说明细胞回收效率极低后面分析的意义也不大了。2.3 一个关键的误区barcode拆分不等于直接按reads分组这是我在实际操作中最想纠正的一个误区。很多人拿到数据后第一反应是写一个脚本根据barcode序列把fastq里的reads分到不同文件然后就认为拆分完成了。并不是。这样做虽然能实现“分组”但完全忽略了单细胞数据的一个基本事实测序过程中会有PCR重复同一个DNA片段会被扩增出很多拷贝这些拷贝的barcode完全相同、比对位置也完全相同。如果不做去重一个细胞的数据量会被高估几倍甚至几十倍下游的可变性分析、peak calling都会受到影响。标准的处理流程应该是barcode提取、barcode校正、比对、按barcode分组、标记/去除PCR重复、过滤低质量比对、生成细胞-峰矩阵。每一个环节都是独立的不能简单地用“按reads分组”替代。2.4 barcode有错配怎么办校正与白名单机制16bp的barcode理论上有4的16次方种组合空间非常大。但在实际测序过程中barcode序列本身也会引入测序错误比如一个碱基被读错。如果不做校正同一个细胞的一部分reads会被错误地分到另一个不存在的“假barcode”下。10X的官方流程里barcode白名单whitelist机制就是为了解决这个问题首先找到数据中丰度最高的barcode集合然后用这些barcode作为参考把与白名单barcode相差一个碱基的barcode合并到白名单里。在我自己处理非标准数据时也沿用了这个思路。我不会把所有barcode都当作真实细胞而是先统计barcode丰度取reads数较高的前N个barcode作为“种子”然后计算所有低丰度barcode与种子barcode之间的编辑距离。如果某个低丰度barcode与某个高丰度barcode只差一个碱基大概率就是测序错误直接合并。这个逻辑自己用Python实现并不难后面我会给出代码示例。3. 实操过程从fastq到“单细胞”的完整拆分流程3.1 环境准备与工具选型我的处理环境是Linux服务器conda管理环境。核心工具包括fastp用于FASTQ质控、接头去除、barcode提取bwa mem或bowtie2用于序列比对scATAC-seq的数据我更习惯用bwa mem速度更快对转座酶随机插入产生的短片段比对效果也更好samtoolsBAM文件处理sinto专门用于按barcode拆分BAM的工具非常好用自己写的Python脚本用于barcode统计、校正、可视化SnapATAC2或ArchRR包用于后续的单细胞分析这里只涉及拆分不展开如果走cellranger-atac路线还需要单独安装cellranger-atac目前官方推荐的是cellranger-atac 2.1.0版本或者更新的cellranger ARC。但我要强调即使最终走cellranger-atac流程理解上面的基础逻辑也很有必要因为很多报错和异常结果都需要从barcode的角度去排查。3.2 方案A用cellranger-atac fastq文件直接跑标准流程如果你的数据完全是10X标准建库、barcode位置标准、无严重污染我的建议还是直接用cellranger-atac count。这是最简单、最稳定、后面做可视化也最方便的方式。cellranger-atac count \ --idsample_01 \ --reference/path/to/refdata-cellranger-arc-GRCh38-2020-A-2.0.0 \ --fastqs/path/to/fastq \ --sampleSampleName \ --localcores32 \ --localmem128这个流程会依次完成FASTQ质控与barcode处理、比对、过滤、去重、peak calling、生成单细胞矩阵。最终输出在outs/目录下核心文件是filtered_feature_bc_matrix用于下游分析和singlecell.csv包含每个barcode的质控指标。但这里有个需要注意的细节cellranger-atac对barcode的识别是严格按照10X白名单来的如果你的数据barcode本身就有问题比如接头污染严重它在第一步就会报“Low fraction of valid barcodes”的警告。我遇到过一份数据valid barcode fraction只有40%左右这种情况下继续跑下去出来的结果基本没有生物学意义。这时候就要考虑走方案B。3.3 方案B手工拆分自己控制每一步当你确定cellranger-atac不能直接用的时候手工拆分的路线就派上用场了。第一步从fastq里提取barcode并统计丰度。我习惯用Python脚本处理读取R1的前16bp作为barcode如果数据是双barcode结构P5和P7侧均有barcode则把两个barcode拼接成一个更长的barcode。统计每个barcode出现的reads数输出分布直方图。import gzip from collections import Counter def extract_barcodes(r1_file, barcode_len16, max_reads2000000): barcode_counter Counter() with gzip.open(r1_file, rt) as f: for idx, line in enumerate(f): if idx % 4 1: # 序列行 barcode line.strip()[:barcode_len] barcode_counter[barcode] 1 if idx max_reads * 4: break return barcode_counter这里限制只统计前200万条reads目的是快速拿到barcode分布概况不需要全量统计等到正式拆分时再全量处理。第二步确定有效barcode集合。我通常做一个双维度筛选第一维是reads数阈值在UMI水平上单个细胞的reads数通常不会低于1000基于我的经验不同建库方式会有差异第二维是累积占比取累积占比达到90%以上的barcode作为有效集合。这两步筛完后再做一个基于编辑距离的校正把距离有效barcode一个碱基的barcode合并进去。第三步用sinto按barcode拆分BAM。先在比对后对BAM文件添加CB标签然后用sinto的filterbarcodes或barcodeprofiling功能。# 添加barcode标签到BAM用sinto的annotation功能 sinto annotate \ --bam sample.sorted.bam \ --fragment barcodes.tsv \ --barcodefile barcodes.txt \ -o sample.annotated.bam # 按单细胞拆分BAM sinto filterbarcodes \ --bam sample.annotated.bam \ --barcodes barcodes.txt \ --outdir single_cells/sinto的filterbarcodes会根据barcode标签自动把BAM拆分成每个细胞一个BAM文件。但我在实际使用时发现当细胞数量超过5000个时这种拆分方式会产生大量小文件对inode和磁盘IO的压力很大。所以我更多时候是直接把CB标签保留在BAM里后续用Python按标签分组处理避免小文件灾难。第四步生成片段文件。scATAC-seq的核心分析对象是“片段”而不是“reads”。转座酶会把开放的染色质区域切出来一个片段的两端会分别比对上在BAM里表现为一对reads。片段文件格式是bed格式每行包含染色体、起始、结束、CB标签、重复数。sinto也有专门的fragments命令来做这件事sinto fragments \ --bam sample.annotated.bam \ --barcodes barcodes.txt \ --outdir fragments_out/生成的fragments文件才是下游分析真正需要的东西。SnapATAC2、ArchR这些工具也都是在片段文件的基础上进行稀疏矩阵构建的。3.4 自己写脚本拆分的完整演示实际上在日常工作中如果数据量不大几千个细胞用sinto工具就足够了。但如果你的需求比较特殊比如要按自定义规则拆分、或者需要把拆分和过滤逻辑整合在一起我建议直接用pysam写一个拆分脚本。下面是一个可以直接取用和修改的示例脚本import pysam from collections import defaultdict def split_bam_by_barcode(input_bam, output_prefix, valid_barcodes): # 打开BAM文件 bamfile pysam.AlignmentFile(input_bam, rb) # 按barcode缓存read cells defaultdict(list) for read in bamfile.fetch(): if read.is_unmapped: continue try: bc read.get_tag(CB) except KeyError: continue if bc in valid_barcodes: cells[bc].append(read) # 写入拆分后的BAM文件 for bc, reads in cells.items(): out_path f{output_prefix}_{bc}.bam with pysam.AlignmentFile(out_path, wb, templatebamfile) as out_bam: for read in reads: out_bam.write(read) print(fWritten {len(reads)} reads for {bc}) bamfile.close() # 使用方式 if __name__ __main__: valid_bc set(open(valid_barcodes.txt).read().strip().split()) split_bam_by_barcode(sample.sorted.bam, single_cell, valid_bc)这段脚本的逻辑其实也就是把“按barcode把BAM拆分”这件事用代码实现了一遍。写出来的好处是你可以在写的过程中逐步加入自己的过滤逻辑比如过滤掉比对上线粒体基因组的reads或者过滤掉比对质量低于某个阈值的reads。我当时处理那批非标准数据时就是先用这种脚本跑了一遍拆分确认每个barcode的片段分布、TSS富集分数都正常之后再决定要不要用cellranger-atac重新跑一遍标准流程做结果交叉验证。两者结论一致我才放心把结果交给下游分析。4. 常见问题与排查技巧实录4.1 valid barcode比例偏低是怎么回事如果你跑cellranger-atac时发现valid barcode fraction低于70%大概率不是细胞活性问题而是文库结构问题。我遇到过的几种情况第一barcode位置判断错误。有些版本试剂盒的barcode序列不在Read 1的最前16bp而是在Read 2的起始处或Index序列里。如果你没细看文库说明书直接按Read 1前16bp处理自然提取不到对的barcode。解决办法是用MultiQC结合fastq的碱基质量分布来判断或者直接查看Space Ranger/Cell Ranger的文库配置。第二接头污染。建库过程中如果转座酶的接头序列没洗干净测序时会读到大量接头序列这些read的barcode区域实际上是接头序列的一部分会形成大量错误的barcode。解决办法是在QC阶段用fastp加--detect_adapter_for_pe选项先做一轮严格的接头去除。第三样本混合测序时Index hopping。多通道测序时如果index序列本身不兼容或出现错误会有一部分reads被错误分配。这个相对少见但如果数据集中某一个样本的valid barcode特别低而其他样本正常需要考虑这个因素。4.2 拆分后单个细胞的比对率和片段数差异巨大这是一个非常常见的情况。拆分后你会发现有的细胞有几十万条片段有的细胞只有几百条。不要急着把所有低片段数的细胞丢掉先关注两个指标一是细胞核的完整性和健康状况。死细胞、破碎细胞核的染色质开放性会显著降低转座酶可进入的区域变少碎片数自然少。二是测序深度。同一个文库中不同细胞的reads数本身就呈负二项分布少量细胞数据量极高是正常的。我在处理时通常不会在线性尺度上过滤而是看片段数的分布直方图寻找一个明显的“拐点”。在ArchR里可以用createArrowFile时的minFrags参数来设置阈值推荐先设为1000如果细胞数还是明显高于预期再逐步提高到3000到5000。但在调整阈值时一定要结合TSS富集分数来看有些细胞片段数很少但TSS富集很高说明这些细胞是真实且高质量的只是测序深度低这类细胞可能代表了稀有细胞类型过滤时要谨慎。4.3 barcode校正后突然少了很多细胞这种情况通常是因为校正逻辑写得过于激进。如果白名单里的barcode数量过大或者编辑距离阈值设置得过高比如允许两个以上的错配会把两个真实细胞的barcode错误地合并到一个barcode下。10X官方用的原则是“一个碱基错配且丰度低于真实barcode的1%”这个经验值可以直接借用。如果你用自定义的校正规则一定要先可视化barcode丰度分布看看不同丰度区间barcode的数量比例再调整合并阈值。另外要注意的是barcode校正对低深度细胞的影响更大。一个细胞如果只有几百条reads某一条read带着一个碱基的错误校正后可能就把这条read并到了错误barcode下。所以对低深度barcode做校正时建议只用比对质量高的reads来计算校正候选避免噪声干扰。4.4 拆分后做Peak Calling发现信号极其分散这是拆分过程中非常容易踩的坑。如果你在比对前就按barcode把fastq拆分了然后对每个细胞单独做bwa比对比对率通常会很低。为什么因为bwa mem这类比对工具是对短read进行全局比对单细胞ATAC的数据里大量片段长度在50bp以下单独比对时很多reads会因比对位置不唯一或多重匹配而被丢弃或标记为低质量。正确的做法是先对所有reads做统一的批量比对比对完成后再根据barcode信息做分组和拆分。不要先拆分后比对这是一个反模式。我在初次处理这类数据时就吃过这个亏拆分后比对率只有30%后来重新按“先比对后拆分”的流程走比对率恢复到了95%以上。4.5 barcode标签在BAM里丢了用sinto或自定义脚本处理时有时候会发现BAM文件里的CB标签丢失了。这个问题通常出在BAM排序环节samtools sort默认只保留standard tags如果你没有显式指定-t CB参数自定义标签会在sort过程中被丢弃。# 排序并保留CB标签 samtools sort -t CB -o sample.sorted.bam sample.unsorted.bam在写pipeline脚本时这一点特别容易疏忽。建议在每一步BAM处理命令里都检查一下确保CB标签存活。5. barcode拆分之后还需要做什么拆分解决的是“每个细胞有哪些片段”的问题。但拿到每个细胞的片段之后数据还不能直接用来做聚类中间还有几个环节我用简短篇幅过一下方便你把拆分放到整个流程里理解。第一个环节是生成细胞-峰矩阵或细胞-bin矩阵。把基因组划分成固定大小的窗口比如500bp的bin统计每个细胞的片段在每个窗口的覆盖情况得到一个稀疏矩阵。这个矩阵和单细胞RNA-seq的表达矩阵在形式上非常类似但数值是二进制的一个bin是否有片段覆盖而不是连续的表达量。第二个环节是维度归约和聚类。scATAC-seq的数据稀疏性极高一个细胞通常只有几千个可检测的开放区域所以直接做PCA效果很差。SnapATAC2、ArchR等工具会先用LSILatent Semantic Indexing或TF-IDF变换再降维到低维空间最后做聚类。第三个环节是细胞类型注释。和scRNA-seq不同scATAC-seq注释细胞类型通常需要用到染色质开放性区域的差异、motif富集分数、以及与参考图谱的比对。这也是为什么拆分质量直接决定了注释的准确度——如果barcode拆分不干净会有大量“合成细胞”出现在聚类结果中导致某一群细胞数虚高最后注释出的细胞类型也会是错的。所以在实际项目中拆分完成后我不会立刻冲向下游分析而是会做至少两个质控检查一是随机抽几个barcode用IGV看一下对应区域的片段分布是否正常二是检查所有barcode在TSS区域的富集分数如果整体偏低可能是文库质量问题也可能是拆分逻辑有问题需要回头排查。写在最后一点个人的体会做scATAC-seq barcode拆分说难也难说容易也容易。难就难在大部分初学者面对一堆BAM、fastq、barcode时根本没有一个清晰的“处理地图”一个环节处理错了后面全部的质控指标都看不懂。容易的地方在于只要你理解barcode就是细胞的身份证号、拆分的核心就是把身份证号对应的所有reads/fragments完整准确地归位再复杂的流程也能拆解成一个个可独立验证的小步骤。我建议大家第一次处理这类数据时不要只跑一条cellranger-atac命令拿到输出就交差。花半天时间对照本文写几个小脚本把barcode的分布、白名单的构建、拆分的逻辑都手动走一遍。这个过程能帮你建立非常扎实的直觉以后再遇到数据异常排查起来比别人快好几倍。我自己就是在这条路上一点一点摸过来的走弯路不可怕怕的是走完了路还不知道自己走过弯路。