上周隔壁课题组找我帮忙他们的变异检测结果已经拿到了VCF文件但问题很实际每个样本到底有多少个SNP听起来就是一行统计的事可真上手去做不少人会在这里卡上一整天。用Excel打开VCF几十万行直接卡死用bcftools stats跑完发现它默认按位点整体统计不按样本拆分自己写Python脚本又搞不清楚FORMAT里GT、AD、DP这些字段到底哪个该用来判断“这个样本在这个位点有没有变异”。这篇我就把平时处理这类需求的完整流程写出来从bcftools安装到纯Python解析再到生产环境用的cyvcf2方案附带一堆实测踩坑记录适合刚拿到VCF文件、需要快速给课题组交统计结果的生信同学参考。1. 拿到VCF别急着写代码先弄清楚要统计什么、口径怎么定1.1 VCF这张表的结构前8列是位点后面才是样本VCFVariant Call Format是存储变异位点的标准文本格式每一行代表一个变异位点而不是一个样本的检测结果。文件最前面的##行是元信息真正的工作区从#CHROM这一行开始。固定列是8列CHROM、POS、ID、REF、ALT、QUAL、FILTER、INFO从第9列开始每一列就是一个样本在该位点上的测序与比对结果。我经常用一个具体例子跟人讲比如这样一行经典记录#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001 NA00002 NA00003 20 14370 rs6054257 G A 29.0 PASS DP14 GT:GQ:DP 0|0:48:4 0|1:48:4 1|1:43:5这行的意思是20号染色体14370位置上参考碱基是G检测到A这个变异NA00001是0|0没有变异NA00002是0|1杂合NA00003是1|1纯合。FORMAT这一列用冒号分隔了多个字段其中GTGenotype一定排在最前面是判断样本有没有变异的根本依据。很多新手一上来就盯着INFO里的DP或AD类字段看那是测序深度和等位基因支持数不等于基因型判定结果。1.2 三种常见统计口径先跟提需求的人对齐“每个样本的SNP统计”这句话有歧义我在实际工作中至少遇到三种不同理解如果不先对齐后面代码白写口径含义统计代码关注点每条记录都算不管该样本GT是什么只要这个位点在VCF里就算一个SNP不需要看GT直接按位点加1只算变异型的样本GT为0/1、1/1这类携带变异的样本才计数0/0和./.不计必须解析每个样本的GT字段只算PASS且变异型在口径2基础上额外要求该位点FILTER列是PASS既要解析GT又要过滤FILTER大多数课题组问“每个样本多少SNP”默认其实是第二种——携带变异的样本才计数。但也有PI只要一个粗略总数哪种口径都能接受。我的经验是先按口径2写一个版本再额外输出一个PASS-only版本两个数字一起交付让选择权交回去。千万不要替对方做决定否则数字对不上报告解释成本远高于多跑一遍脚本的时间。1.3 准备一份样例文件用来验证脚本无论走哪条实现路线都需要一份带已知结果的VCF用于验证。实战中我通常先拿公开的千人基因组VCF或自己测序项目的一个小片段比如只保留1号染色体前5万行来做。小文件不仅跑得快还能人工核对数字确认脚本逻辑没问题后再丢到全基因组规模的VCF上跑。个人建议按下面的方式准备一份样例# 从完整VCF中截取一部分染色体区间 bcftools view -r chr1:1000000-2000000 big.vcf.gz test_region.vcf # 再随机保留一部分位点方便肉眼核对 bcftools view -H test_region.vcf | head -50 test_region_head.txt如果原文件没有按染色体排序先做一次归一化和排序再截取否则bcftools view -r会报区间不连续或者定位错误。2. bcftools装不上的一半人都是卡在同样的地方安装与验证全记录2.1 conda路线一条命令配好环境bcftools是生信环境里最常见的工具之一但对新手来说apt install bcftools装出来的版本可能老到连bcftools query的-f输出格式都跟你网上搜到的教程对不上。我自己更推荐用conda管理生信软件好处是版本锁定、环境隔离不会把系统Python或系统库搞乱。conda create -n bio -c conda-forge -c bioconda bcftools1.19 conda activate bio如果觉得conda装软件时解析依赖太慢建议换成mambaconda install -n base -c conda-forge mamba mamba create -n bio -c conda-forge -c bioconda bcftools1.19装完后执行bcftools --version会看到类似bcftools 1.19, using htslib 1.19的输出。这里有个小细节conda的channel优先级会影响版本选择conda-forge和bioconda两个channel同时使用时建议把conda-forge写在前面不然可能出现htslib依赖版本与bcftools不匹配的问题。2.2 apt/yum路线快但版本控制要留个心眼在Ubuntu或Debian服务器上很多人习惯直接用系统包管理器sudo apt update sudo apt install -y bcftools优点是快缺点也很明显Ubuntu 20.04自带的bcftools是1.10左右Ubuntu 22.04是1.13CentOS则可能停留在1.9。老版本不是不能用而是部分新功能缺失比如bcftools plugin相关命令和新的-Ou压缩输出行为。如果你只是做简单统计系统包够用如果你后续要跑比较复杂的过滤、合并或分割功能建议还是上conda或源码编译避免因为版本差异浪费时间。2.3 源码编译路线依赖处理与经典报错在没有root权限的服务器上源码编译是唯一可行方案。bcftools的源码包可以从GitHub官方仓库下载推荐下载release版的tar.bz2wget https://github.com/samtools/bcftools/releases/download/1.19/bcftools-1.19.tar.bz2 tar -xjf bcftools-1.19.tar.bz2 cd bcftools-1.19 ./configure --prefix$HOME/bcftools make -j4 make install编译过程最大的坑是系统缺少基础依赖库我列几个常见报错及对应解决方式报错信息缺失依赖解决方式zlib.h: No such file or directoryzlib开发包sudo apt install zlib1g-devliblzma not foundlzma开发包sudo apt install liblzma-devUnable to find htslibhtslib库bcftools 1.12自带htslib源码会自动编译1.9及更老版本需手动编译htslib并设置HTSDIRlibcurl not foundlibcurl开发包sudo apt install libcurl4-openssl-dev或configure时加--disable-libcurl在没有root权限、也没有conda的机器上如果依赖库实在装不全可以尝试./configure --disable-bz2 --disable-lzma --disable-libcurl跳过部分可选依赖的编译。bcftools最核心的VCF读取、GT统计功能不依赖这些可选库压缩读取用的是zlib这个是必须装的。2.4 装完之后先验证功能装完环境后我建议先用一个简单的query命令验证安装是否正常工作bcftools query -f %CHROM\t%POS\t%REF\t%ALT\n test_region.vcf | headbcftools query的-f参数是格式化输出字段%CHROM、%POS、%REF、%ALT分别对应染色体、位置、参考碱基、变异碱基。能正确输出这几列说明核心功能没问题。如果这一步就报错多半是环境变量没配好或者conda环境没有真正激活。3. 用bcftools query先拿到“标准答案”再谈优化3.1 一行命令看全貌把每个位点的每个样本GT铺开bcftools stats不适合按样本统计但bcftools query非常合适。它支持按样本维度输出GT把每个位点上每个样本的基因型横向展开bcftools query -f %CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE%GT]\n test_region.vcf | less -S[...]括起来的部分是对每个样本循环展开的语法。输出结果类似chr1 100001 A G sampleA0/1 sampleB0/0 sampleC1/1 chr1 100078 C T sampleA0/0 sampleB./. sampleC0/1这个输出对人类阅读非常友好三个样本、每个位点的基因型一目了然。这也是我后面验证Python脚本是否正确的重要参考标准。3.2 完整统计命令query加awk在管道里完成接下来把上面的输出通过管道交给awk按样本名累加“至少携带一个变异等位基因”的次数。GT规范化是重点VCF里的phased基因型会用竖线分隔比如0|1本质和0/1一样统计前必须先统一符号bcftools query -f [\t%SAMPLE%GT]\n test_region.vcf | \ awk -F\t { for (i1; iNF; i) { split($i, a, ) gt a[2] gsub(/\|/, /, gt) n split(gt, g, /) has_var 0 for (j1; jn; j) { if (g[j] ! 0 g[j] ! .) has_var 1 } if (has_var) count[a[1]] } } END { for (s in count) print s, count[s] } | sort解释一下这段awk的判定逻辑GT字段如果只有0和.表示该样本没有携带变异0/0是纯合参考./.是缺失只要出现任何数字大于0的等位基因索引比如1、2就说明该样本在这个位点携带了变异。最后的输出就是每个样本的SNP携带数量。这里我要提醒一个容易踩的坑上面的统计没有过滤SNP与INDEL也没有过滤FILTER列。%CHROM\t%POS\t%REF\t%ALT这几列如果不用可以直接不输出awk里只处理样本列即可。但如果你要按“只统计SNP”的口径走就要在query阶段加上判断awk里检查REF和ALT的长度是否为1。我实际跑的完整命令一般是这样的bcftools query -f %REF\t%ALT[\t%SAMPLE%GT]\n test_region.vcf | \ awk -F\t { ref $1; alt $2 if (length(ref) ! 1) next n_alt split(alt, a, ,) is_snp 1 for (k1; kn_alt; k) if (length(a[k]) ! 1) is_snp 0 if (!is_snp) next for (i3; iNF; i) { split($i, b, ) gt b[2] gsub(/\|/, /, gt) n split(gt, g, /) has_var 0 for (j1; jn; j) if (g[j] ! 0 g[j] ! .) has_var 1 if (has_var) count[b[1]] } } END { for (s in count) print s, count[s] } | sort -k2,2nr这么一步步判断能保证得到的数字是严格执行“SNP、按样本、携带变异”这三个条件的。3.3 用人工可读方式交叉验证结果为了方便后续核对我还会顺手输出一个可读性强的小摘要bcftools query -f %CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE%GT]\n test_region.vcf | \ awk -F\t NR10 {print}挑前10个位点人工数一下每个样本的变异型数再和3.2节的统计结果比对。这个小动作看似笨但能把“命令写错了但看着像那么回事”的风险降到最低。我在带新人时反复强调管道的每一环都可能有隐性问题不要信任一条没验证过的统计命令的输出。4. Python逐行解析VCF不依赖第三方库的统计脚本4.1 脚本设计三步走虽然bcftools已经能解决统计问题但很多场景下大家还是希望有一个不依赖外部工具的Python脚本原因不外乎几点一是最终要集成到已有的Python分析流程中二是需要对统计结果做更多自定义加工三是团队环境里不一定有权限装bcftools。既然标题是以Python实战为主这一步才是重头戏。纯Python解析VCF的思路可以拆成三步逐行读取文件跳过##开头的元信息遇到#CHROM行时解析样本名列表对每个非表头行先判断是否为SNP位点再解析每个样本的GT字段并计数。4.2 完整代码及各段逻辑注释我用gzip模块读取压缩的VCF这样无需手动解压也不会占用大量磁盘空间。判断SNP位点的方法是看REF和ALT是否为单碱基REF必须长度1ALT可能有多等位基因用逗号分隔所有ALT都必须长度1。GT的解析则要注意phased的竖线符号|统一替换成斜杠/再判断。import gzip from collections import defaultdict def count_snp_per_sample(vcf_path, only_passTrue): sample_counts defaultdict(int) samples [] with gzip.open(vcf_path, rt) as fin: for line in fin: if line.startswith(##): continue # 遇到#CHROM行解析样本名列表 if line.startswith(#CHROM): header line.strip().split(\t) samples header[9:] continue cols line.strip().split(\t) # 可选只统计FILTER列为PASS的位点 # 很多VCF把未过滤位点写成.PASS写成PASS if only_pass and cols[6] not in (PASS, .): continue ref cols[3].upper() alt cols[4].upper() # 判断SNPREF必须是单碱基所有ALT也必须是单碱基 if len(ref) ! 1: continue alt_set alt.split(,) if any(len(a) ! 1 for a in alt_set): continue # 解析FORMAT找到GT字段的位置 fmt_keys cols[8].split(:) try: gt_pos fmt_keys.index(GT) except ValueError: continue # 逐个样本解析GT并计数 for idx, sample in enumerate(samples): sample_fields cols[9 idx].split(:) if gt_pos len(sample_fields): continue gt sample_fields[gt_pos] # phased基因型用|统一替换成/ gt gt.replace(|, /) alleles gt.split(/) # 只要有一个等位基因不是0且不是.就算携带变异 non_ref [a for a in alleles if a ! 0 and a ! .] if non_ref: sample_counts[sample] 1 return sample_counts if __name__ __main__: result count_snp_per_sample(input.vcf.gz, only_passFalse) for sample, count in sorted(result.items(), keylambda x: x[1], reverseTrue): print(f{sample}\t{count})这段代码默认only_passFalse只统计全部位点因为很多实际VCF的FILTER列很乱PASS、.混着来如果不加判断就把有效位点漏掉了。脚本里已经预留了only_pass开关想只统计PASS位点时改成only_passTrue即可但要看清楚自己文件的FILTER列到底是怎么标记的。4.3 效率与内存gzip逐行读几百万行不虚有朋友会担心Python逐行解析几十万行VCF会不会很慢。以我实测的经验一个包含约50万个位点、88个样本的VCF文件用上面这段纯Python脚本跑完大概需要8到15秒主要时间花在gzip解压和字符串切割上。对于日常交付完全够用。内存方面因为是逐行读取、逐行释放整体内存占用非常低只有样本计数用的字典会随样本数增长通常可以忽略不计。这也是为什么我不建议在这个场景用pandas去read_csv整读VCF——文件一大内存就被吃掉了而且还需要额外处理##注释行收益并不高。如果确实需要极致提速可以考虑两个方向一是用multiprocessing按染色体分块并行解析二是换用下面要讲的cyvcf2这类C扩展库。但90%的场景下上面的纯Python脚本已经够用先保证正确性再谈性能优化这是我一贯的原则。5. 生产环境我更喜欢cyvcf2/pysam为什么5.1 cyvcf2方案代码对比纯Python解析适合理解VCF结构和临时性任务但在生产环境里处理全基因组VCF时我更喜欢用cyvcf2。它是用C语言封装htslib的库解析速度比纯Python快一个数量级而且API设计得比pysam更简洁。from cyvcf2 import VCF def count_snp_with_cyvcf2(vcf_path): vcf VCF(vcf_path) samples vcf.samples counts {s: 0 for s in samples} for variant in vcf: # cyvcf2中FILTER为None表示PASS if variant.FILTER is not None: continue ref variant.REF.upper() alt [a.upper() for a in variant.ALT] # 只保留SNP if len(ref) ! 1 or any(len(a) ! 1 for a in alt): continue # genotype.array()返回每个样本的等位基因索引矩阵-1表示缺失 gts variant.genotype.array() for i, sample in enumerate(samples): a1, a2 gts[i] if (a1 0) or (a2 0): counts[sample] 1 return countsvariant.genotype.array()返回的是一个二维数组形状是(样本数, 2)每个位置的取值是等位基因索引0表示参考等位基因1、2等表示ALT等位基因索引-1表示缺失。判断a1 0 or a2 0的意思是这个样本在该位点上至少携带一个ALT等位基因正好对应GT0/1或1/1这些变异型。为什么这个方案在生产环境更合适关键在于速度。用同一份50万位点、88样本的VCF测试纯Python版本需要10秒左右cyvcf2版本基本在1秒内结束。当你手里的文件从50万位点变成5000万位点时这个差距会从十几秒拉大到十几分钟后者在交互式分析里是难以接受的。5.2 pysam的边界细节pysam同样是htslib的Python绑定但API更底层、更贴近C接口。很多时候pysam是绕不开的因为它支持BAM、CRAM、VCF等多类文件在同一个流程里既能读比对结果又能读变异结果比较方便。import pysam vcf pysam.VariantFile(input.vcf.gz) counts {s.name: 0 for s in vcf.header.samples} for rec in vcf.fetch(): filters rec.filter.keys() if filters and filters ! (PASS,): continue ref rec.ref.upper() alts [a.upper() for a in rec.alts] if rec.alts is not None else [] if len(ref) ! 1 or any(len(a) ! 1 for a in alts): continue for sample in vcf.header.samples: gt rec.samples[sample].get(GT) if gt is not None and any(a is not None and a 0 for a in gt): counts[sample] 1需要注意pysam的几个细节rec.filter.keys()在VCF的FILTER列为.时返回空元组在PASS时返回(PASS,)。所以判断保底逻辑要写成if filters and filters ! (PASS,): continue不能简单地检查PASS in filters。rec.samples[sample].get(GT)返回的是一个元组比如(0, 1)是杂合(0,)表示单倍体位点(None, None)表示缺失。半合子的元组长度是1遍历时如果硬按两个元素去取会越界所以用any()会安全很多。对于染色体X和线粒体等非二倍体样本GT元组长度可能不是2这也是为什么生产脚本里我坚持用any(a is not None and a 0 for a in gt)而不是写死gt[0]和gt[1]。5.3 什么时候用纯Python什么时候上库如果给我一个选择我的分界线是场景推荐方案一次性交付文件不大百万位点以内纯Python脚本最省事无需装额外依赖学习VCF格式、理解字段含义纯Python逐行解析是绝佳教材全基因组VCF几千万行cyvcf2速度差距非常可观需要同时处理BAM和VCFpysam一套API解决两类文件需要按区间、按样本做复杂过滤直接用bcftools命令行更灵活纯Python脚本的另一个不可替代的价值是当你需要自定义统计规则时改起来非常直观。比如有同事问“只统计纯合SNP”或者“只统计支持深度大于10的SNP”加几行判断就行不用去查bcftools的表达式语法。6. 交付统计结果前我总会再查这几个坑6.1 FILTER列是.还是PASS统计口径的差异这个问题坑过不少人。很多VCF在caller输出时并不会给每个位点都标注PASS有的位点FILTER列只是.表示“没有经过任何过滤条件”这不等于通过过滤但在有些统计逻辑里却会被当成通过。反过来有些全外显子组数据分析流程会把所有非PASS位点都过滤掉再交付这时VCF里只会看到PASS位点。我的建议是写脚本时不要默认FILTER一定等于PASS而是先跑一遍描述性统计bcftools query -f %FILTER\n input.vcf.gz | sort | uniq -c看看文件中到底有几种FILTER标记再决定only_pass开不开。如果文件本身是原始VCF且未经过严格过滤我交付统计结果时通常给两个数字全部位点的SNP数和PASS位点的SNP数让下游分析者自己选。6.2 phased基因型、多等位、半合子VCF里的GT字段除了常见的0/0、0/1、1/1外还有几种容易漏处理的场景phased基因型0|1表示两个等位基因分别来自父源和母源在统计有没有变异这件事上|和/没有区别。脚本里统一replace(|, /)即可。多等位位点ALT列有逗号分隔的多个等位基因比如A,G。这种位点的GT可能是1/2表示同时携带两个不同的ALT。判断时只要等位基因索引大于0就算变异不能只判断“是否为0/1”。半合子男性样本的X染色体、线粒体样本可能只有一个等位基因GT元组长度为1。纯Python的split和cyvcf2的any()都能自然处理但pysam里如果写死gt[1]就会越界所以一定要用遍历或any()。还有一个很容易被忽略的点ALT为*或.的位点。*表示缺失等位基因.在ALT列有时代表“没有可报告的变异”这些位点如果恰好被放进来单碱基长度判断通常能过滤掉一部分但保险起见遇到奇怪的字符还是人工看一眼。6.3 和已有报告对不上的排查思路如果你统计出的数字和课题组之前拿到的报告对不上不要把归因停留在“代码写错了”这一步。按以下顺序排查90%的问题能定位先确认统计口径对方要的是SNP还是所有变异含INDEL要的是携带变异的样本数还是位点出现次数这两个口径差出来能到3倍以上。确认是否过滤了PASS原始VCF里每个位点是否都通过质检对方报告里有没有“QC passing”字样确认参考等位基因判断VCF里REF是参考等位基因ALT是变异等位基因有些工具输出时把二者调转统计就歪了。用bcftools验证已知小样本挑两三个样本用bcftools view -s sampleA input.vcf.gz单独提取再用bcftools stats -s sampleA对比小数点后都核对清楚。我经历过最离谱的一次是跑完全部样本发现数字和团队历史报告差了整整一倍。查了很久才发现对方历史流程里用的是“只统计纯合SNP”而我们统计的是“所有携带变异的SNP”。纯合变异的等位基因两个都是ALT恰好把杂合全部漏掉。这个教训让我从此在写统计脚本前一定会盯着对方问清楚统计口径甚至把口径的具体定义写进交付文档里。最后再分享一个实际操作中的小技巧不管用哪种方案统计结果一定要用tab分隔输出且文件名和列名都带上口径标记比如sample_snp_count_all.tsv和sample_snp_count_pass.tsv。刚开始觉得多此一举但当你一个月后回来看结果文件、或者把文件发给协作者时你会感谢当初这个节约沟通成本的决定。