WGS多样品VCF文件：如何准确统计每个样本的SNP数量

做WGS项目时手头有一份多样品联合基因分型得到的VCF文件想把每个样本各自检出的SNP位点数统计出来这个需求看起来简单实际操作起来坑真不少。先说个容易混淆的点平时搜“SNP相关操作”时会碰上射频仿真圈子里HFSS导出SNP文件那个SNP是Touchstone/sNp格式n代表端口数跟基因组里的单核苷酸多态性SNP完全是两码事。这篇文章只围绕WGS数据中的VCF文件讲清楚多样品VCF里每个样品SNP数统计的方法细节。你可能会遇到直接用vcftools --counts2出来的却是所有样本等位基因频率用脚本硬撸VCF结果几十G文件把内存撑爆统计口径没对齐不同人跑出来的数字差一倍。这些坑我都会展开聊并提供可直接运行的代码和注释。适合做群体重测序、WGS家系分析、以及刚接触VCF处理的读者参考。1. 需求拆解多样品VCF中统计每个样品的SNP数到底在统计什么1.1 这个任务从哪来WGS项目里的真实场景先想想什么时候会需要这种统计。最常见的是样本QC一批样本刚完成全基因组测序和比对变异检测需要快速排查哪些样本测序深度不够、比对率低导致检出的SNP数量明显偏少或者怀疑样本污染、样本混淆时也可以通过SNP数量做初步判断。比如有个项目里20个样本客户只想知道每个样本的somatic SNP数量本质上就是把VCF按样本拆开数一遍。还有一种场景是后续分析的前置步骤构建系统发育树、群体结构分析、亲缘关系评估时需要知道每个样本有多少可用于分析的SNP才能决定过滤阈值和缺失率标准。这类统计在WGS项目中几乎绕不开但很多朋友第一次做时是懵的因为VCF文件动辄几十GB样本列多到屏幕塞不下根本不知道该用哪个命令行工具更不知道该按什么规则数。所以第一步不是写脚本而是把任务本身拆清楚。1.2 定义统计口径是位点数不是基因型个数要统计之前先把“SNP数”的定义固定下来。默认情况下这里说的是“每个样品中基因型不是纯合参考、且不是缺失的SNP位点数目”。换句话说某个样本在某SNP位点的基因型若是0/0这个位点不计数是0/1、1/1或1/2这种非参考基因型就计数。容易误会的点在于这是按位点计数不是按等位基因计数。比如某个样本在一个位点是1/1纯合突变如果按等位基因数统计就变成2个SNP按位点统计只有1个。绝大多数下游分析需要的都是位点数因为一个位点无论杂合还是纯合都只代表一个变异事件。另一种口径是只统计“可判定的”位点排除缺失基因型./.的位点。还有更严格的口径要求位点在所有样本里缺失率低于某个阈值或者最小等位基因频率MAF大于0.05。这些口径不同出来的数字相差很大。所以拿到任务时先问清需求别急着跑命令。我建议在动手前先列一份简单的口径清单SNP还是包含InDel是否多等位位点要拆是否过滤缺失率是否过滤MAF这五个问题想清楚了后面写脚本才不会翻车。2. 工具选型解析批量统计SNP数哪条路最靠谱2.1 命令行三件套 vs 写脚本 vs R语言统计方法无非三条路纯命令行工具组合、Python脚本、R语言。纯命令行适合一次性处理、文件规模中等、只输出一张表的场景优点是方便缺点是灵活性差。写脚本适合需要反复调整统计口径、要做复杂过滤、还要把结果跟样本元数据整合的场景优点是完全可控缺点是代码容易写错。R语言里的VariantAnnotation包适合统计完紧接着做PCA或曼哈顿图的情况但R读大VCF本身就是性能瓶颈大规模WGS数据我一般不推荐。如果你数据量很小几十个样本外显子级别R完全没问题但如果是几百个全基因组样本还是老老实实用命令行加Python。2.2 几个关键工具的适用场景对比我日常最常用的组合是bcftools加awk外加一个Python库cyvcf2。bcftools负责VCF的读取、过滤、样本操作awk做行式统计cyvcf2负责在Python里快速解析VCF特别适合要额外算深度、质量等指标的时候。vcftools里的--missing-indv、--freq等也能间接得到一些信息但直接统计每个样本SNP数并不高效因为它通常按位点输出需要自己再聚合。工具组合优点缺点适用场景bcftools awk轻量、快、可管道化多等位和缺失处理要自己写逻辑日常快速统计样本数多也不怕vcftools参数简单生成多种QC信息按样本计数不方便输出字段有限过滤前查看缺失率、频率Python cyvcf2灵活、快、可编程可加复杂逻辑需要装Python包写代码有门槛需要同时输出多个统计维度R VariantAnnotation统计和可视化一体大文件读取慢内存占用高小规模数据、下游作图的场景工具选型没有唯一标准关键看你的数据量多大、后续要不要联动别的表格以及你更熟哪套生态。如果只是临时帮同事看一眼直接用bcftools加awk就好如果这是周期性任务建议一步到位写Python脚本。3. 核心细节解析VCF里藏着哪些“坑”影响了SNP计数3.1 VCF格式的GT字段怎么读VCF的主体是每行一个变异位点前8列是CHROM、POS、ID、REF、ALT、QUAL、FILTER、INFO后面每一列是一个样本。样本列里最重要的子字段是GT也就是样本在该位点的基因型。GT用斜杠或竖线分隔等位基因数字对应REF和ALT列表中的等位基因索引0表示REF1表示第一个ALT2表示第二个ALT。举个例子某行VCF长这样1 1000 rs123 A G . PASS . GT:DP:GQ 0/1:30:99 1/1:28:99 0/0:35:99这里是三个样本GT分别是0/1、1/1、0/0。如果REFAALTG那么0/0就是A/A纯合参考0/1是A/G杂合1/1是G/G纯合突变。如果ALT不止一个比如REFAALTG,C那么1/2表示G/C这时也是一个变异位点但计数时通常仍然只算1个位点除非你特别关心双等位还是多等位。此外还有./.或.|.这种缺失基因型通常是因为测序深度太低、GQ不达标被判定为“无基因型”。统计时你可以选择排除缺失也可以把它当成“未检出”。这个决策直接决定数量必须在脚本里明确处理。我见过一个朋友用Excel打开VCF手工筛选1/1和0/1结果漏掉了一半样本因为Excel把几百列看花眼了。总之GT字段是计数的核心理解它才能真正理解脚本。3.2 SNP过滤和质控参数别乱设很多人拿到VCF直接就开始数结果把大量质量很差的位点也算进去了最后数字虚高。标准的做法是先做一轮基础过滤保留SNP去掉InDel过滤未通过的FILTER列按需要过滤QUAL或DP。bcftools里一条命令bcftools view -v snps -f PASS input.vcf.gz -Oz -o pass.snp.vcf.gz这里-v snps表示只保留SNP-f PASS表示保留FILTER为PASS的位点。但注意有的VCF里FILTER列是.表示没有被标记为异常这时-f PASS会把它过滤掉需要改成-f PASS|.或者用bcftools filter重新标记。接着可以用bcftools view -i QUAL30 INFO/DP10这类表达式做更细的过滤。关于“SNP”定义bcftools的-v snps是指REF和ALT都是单个碱基的变异。如果某个位点既有SNP又有InDel比如ALT是AT这种复杂位点要慎用一般直接丢掉因为很难确定它是不是一个干净的SNP事件。我在实际项目中会连续跑两条命令一条去掉InDel一条去掉多等位这样统计结果最干净bcftools view -v snps input.vcf.gz -Oz -o tmp.snp.vcf.gz bcftools view -m2 -M2 tmp.snp.vcf.gz -Oz -o biallelic.snp.vcf.gz tabix -p vcf biallelic.snp.vcf.gz这里-m2 -M2表示保留等位基因计数从2到2的位点也就是只保留双等位位点。注意这会丢掉多等位SNP如果项目需要保留多等位就不要加这一步。3.3 多等位位点和缺失基因型的处理多样品VCF中经常出现多等位位点比如ALTG,C有些样本是1/1有些是1/2有些是2/2。这类位点如果保留你很难用简单脚本判断它是否“变异”。比较稳妥的办法是先用bcftools view -m2 -M2把多等位位点变成双等位或者用bcftools norm拆分为多个双等位位点。注意-m2 -M2会丢掉多等位位点可能损失一部分信息如果项目不关心多等位这样处理最简单统计结果也更稳定。缺失基因型方面如果你做群体数据我建议用bcftools进行缺失率过滤比如保留缺失率低于20%的位点bcftools view -i F_MISSING0.2 input.vcf.gz -Oz -o miss0.2.vcf.gz这样每个样本的SNP数统计出来相对公平不会因为个别样本在公共位点上缺失太多导致数量异常。当然如果目的是评估样本本身的真实变异检出量就不要过滤缺失率因为缺失少本身就说明样本质检好。这里有个小原则过滤越严格样本间的可比性越强但也越可能洗掉真实生物学差异。具体阈值要根据测序深度和物种多样性来定别照搬别人参数。4. 实操过程与核心环节实现从VCF到每个样本的SNP数4.1 准备阶段检查文件、建索引、提取样本名拿到VCF文件第一件事不是数数而是先确认文件完整性和压缩格式。如果你的文件名是.vcf.gz记得看看旁边有没有.tbi索引没有的话先建索引tabix -p vcf input.vcf.gz如果文件名是.vcf未压缩可以直接用但大文件建议转成bgzip压缩后再操作bgzip -c input.vcf input.vcf.gz tabix -p vcf input.vcf.gz然后提取样本名列表这一步后面循环或核对时都会用到bcftools query -l input.vcf.gz samples.txt这个文件每行一个样本名顺序和VCF中样本列一致。我自己习惯先wc -l samples.txt确认样本数量跟实验设计对得上再看一眼前几行防止样本名有奇怪的符号。样本名里如果带着/、空格、中文后面awk按列拆分时不一定会错但容易在文件命名和shell循环里出问题最好提前用sed把特殊字符替换掉。还有个小习惯执行任何脚本前先head -5 samples.txt和bcftools query -l | head看看别等跑完发现样本名对不上。4.2 方法一bcftools awk 快速统计推荐日常使用假设你已经过滤好了一个双等位SNP的VCFbiallelic.snp.vcf.gz。要统计每个样本有多少个“非参考且非缺失”的位点最直接的方式是用bcftools query把样本基因型按位点展开再用awk按样本聚合。示例脚本如下# 样本列表 bcftools query -l biallelic.snp.vcf.gz samples.txt # 统计%GT会输出每个样本的GT字符串以\t分隔 # [%GT\t] 是每样本输出GT和一个tab bcftools query -f %CHROM\t%POS\t[\t%GT]\n \ biallelic.snp.vcf.gz \ | awk BEGIN{FSOFS\t} NRFNR{samples[NR]$0; next} { # 第1列CHROM、第2列POS不计入样本计数 for(i3; iNF; i){ gt$i # 只有非缺失且不等于0/0或0|0的才计数 if(gt!./. gt!.|. gt!0/0 gt!0|0){ cnt[i-2] } } } END{ for(j1; jlength(samples); j){ print samples[j], cnt[j]0 } } samples.txt - sample_snp_counts.txt这段脚本里有个关键点bcftools query在输出样本GT时如果某些样本在该位点没有信息会输出缺失的./.awk里过滤掉0/0和缺失剩下0/1、1/1、1/2这些都算作一个变异位点。如果VCF已经被过滤成双等位SNP基本只有0/1和1/1两种情况。对于超大规模文件我更常用下面这个更省流的写法直接让bcftools每行输出样本名和GT再用awk按样本累加避免输出宽表bcftools view -H -v snps biallelic.snp.vcf.gz | \ awk { # 样本名通过bcftools query -l 提前读入samples.txt for(i10; iNF; i){ split($i, a, :); gta[1]; # 统计非参考非缺失 if(gt!./. gt!.|. gt!0/0 gt!0|0){ cnt[i-9]; } } } END{ while((getline name samples.txt) 0){ idx; print name, cnt[idx]0; } } - sample_snp_counts.txt这个脚本里假设样本从第10列开始并且GT是每个样本列的第一个子字段。如果你的VCF样本列之前还有固定列数通常8列信息列所以第10列开始没问题。不过更推荐用bcftools query awk的第一种因为bcftools query能够只输出GT不受INFO列多少的影响。另一种完全不需要手动写awk统计的方法是用bcftools stats的--samples参数bcftools stats --samples sample_list.txt biallelic.snp.vcf.gz stats.txt这个命令会输出Samples部分包含每个样本的记录数但要注意它统计的是该样本在所有位点中“至少有一条非缺失等位基因”的数量跟你自己脚本的定义可能略有差异所以用来做数量级参考可以精确对齐口径时还是用自己脚本保险。4.3 方法二Python cyvcf2 灵活定制适合复杂场景如果除了SNP数你还要同时统计转换/颠换比、杂合/纯合比、每个样本的平均深度或者需要根据基因型值做更多判断建议用Python。这里推荐cyvcf2它底层是htslib读取速度接近C比纯Python快几十倍。安装很简单pip install cyvcf2下面这个示例会读入VCF按位点逐行扫描对每个样本判断基因型并输出每个样本的SNP数、杂合数、纯合数# -*- coding: utf-8 -*- from cyvcf2 import VCF # 打开VCF vcf VCF(biallelic.snp.vcf.gz) samples vcf.samples # 统计容器总SNP数、杂合数、纯合数 snp_total {s: 0 for s in samples} het_total {s: 0 for s in samples} hom_total {s: 0 for s in samples} for variant in vcf: # 如果之前没有按SNP过滤这里可以加类型判断 # 只统计SNPREF和ALT长度都为1 if not variant.is_snp: continue # variant.gt_types返回样本基因型类型数组 # 0HOM_REF, 1HET, 2UNKNOWN, 3HOM_ALT gt_types variant.gt_types for idx, sample_name in enumerate(samples): gt gt_types[idx] if gt 2: # UNKNOWN即./. continue if gt 0: # 0/0纯合参考 continue snp_total[sample_name] 1 # 只要是变异就加1 if gt 1: het_total[sample_name] 1 # 0/1 elif gt 3: hom_total[sample_name] 1 # 1/1 # 输出结果表 with open(sample_snp_counts_cyvcf2.txt, w) as f: f.write(sample\tsnp_total\thet_total\thom_total\n) for s in samples: total snp_total[s] het het_total[s] hom hom_total[s] f.write(f{s}\t{total}\t{het}\t{hom}\n) print(完成结果见 sample_snp_counts_cyvcf2.txt)代码注释已经写得比较细了。这里特别要提醒的是gt_types里的UNKNOWN和HOM_REF都要跳过否则会让计数虚高。另外如果VCF中某个位点是多等位基因且样本基因型是1/2gt_types会怎么标记不同版本情况不同最保险的做法还是先用bcftools把多等位位点拆开或过滤掉。如果你不放心也可以直接判断variant.genotypes[idx]返回的等位基因索引列表自己写逻辑。它返回的格式类似[1, 0, True]表示一个等位基因索引、第二个等位基因索引、是否phased。这样你可以精确区分0/1和1/2。后续需要按染色体分段统计也很容易在这个脚本基础上加一个variant.CHROM判断这就很灵活了。4.4 结果核对拿bcftools stats当裁判不管用什么方法统计完我都建议做一次结果核对。最省事的核对方式是先用bcftools stats看整体信息比较总位点数然后跟自己脚本的结果做交叉验证。我自己常用的核对流程是先跑一下完整版本的bcftools stats看总的SNP位点数SNPs number。bcftools stats biallelic.snp.vcf.gz full.stats.txt把我们脚本得到的所有样本SNP数加起来理论上会大于总位点数因为一个位点可能多个样本都有变异但不会大于“总位点数乘样本数”。如果加起来明显不对多半是GT解析出了问题。用R快速读取结果表画一个箱线图看离群样本这样一眼就能发现哪个样本SNP数量异常少。d - read.table(sample_snp_counts.txt, headerFALSE, col.namesc(sample, snp_count)) boxplot(d$snp_count, mainSNP count per sample)我遇到过一个样本统计出来只有其他人的三分之一后来发现它测序深度只有5X确实是一开始的样子就不好。结果核对这一步花不了几分钟但能挡住很多低级错误。特别是当你把结果发给别人的时候客户第一个问题大概率是“这数字对吗”提前自洽很重要。5. 常见问题与排查技巧实录那些让你统计结果对不上的事5.1 我的结果比人家少了一半是怎么回事这个最常见的原因是没有把0/1和1/1都算进去或者只统计了1/1。有些刚入门的同学数GT时只找“1/1”的位点把杂合全丢了结果自然少一半以上。第二种常见原因是位点经过Hardy-Weinberg过滤或MAF过滤后数量本来就减少了这不是bug是过滤阈值不同。第三种原因是对“SNP”的定义不一致有人把多等位SNP拆开成多个位点计数有人按原始VCF行计数数字当然对不上。所以每次比较结果前必须把“过滤命令统计口径”一起对齐别只对数字。如果对方用的是vcftools跑出来的某个指标你最好先确认对方有没有过滤InDel、有没有过滤缺失率。这里我一般会做一个快速验证直接用bcftools统计总SNP位点数再用我们脚本所有样本数量求和如果求和在“总位点数”和“总位点数×样本数”之间基本说明逻辑没问题。5.2 内存爆了shell卡死怎么办如果你用纯Python写了个循环把每个样本逐个用vcftools --indv统计遇到几百个样本的大型VCF可能会非常慢如果直接用pandas读全量VCF几十GB的文件能直接把内存榨干。建议办法文件压缩成bgzip配合bcftools流式读取一次只处理一行awk脚本里只维护样本个数的计数器不要保留整行数据Python用cyvcf2逐行迭代内存占用始终很低。另外如果只是临时脑洞想看看结果可以先用bcftools view -r chr1把范围缩小到一条染色体跑通了再放全基因组这样排错效率高很多。我自己的习惯是先取一条染色体做一小段比如chr1:1-1000000确认输出格式和计数逻辑没问题再全量跑。不然全量跑半小时才发现样本列数不对很浪费时间。5.3 样本名出现乱码或重复怎么处理VCF的样本名有时会带着路径前缀或平台给的编号出现重复的样本名后果非常严重因为awk统计输出时如果两个样本同名最后的输出会少一行。处理办法是在统计前先检查bcftools query -l input.vcf.gz | sort | uniq -d如果看到重复样本名就要返回上游检查样本表或者用bcftools reheader重新命名样本。还有一点小细节样本名里可能有Windows换行符\r在Linux下awk匹配时可能会出错最好先执行sed -i s/\r$// samples.txt这一类问题不容易发现但会导致最终结果文件里某个样本数量变成0排查起来很抓狂。我遇到过一位同事他的样本名带了路径前辍统计出来的结果跟另一个样本完全一样一开始我还以为脚本有bug后来发现是样本名重复改用samples.txt重新映射后问题才解决。5.4 补充一个独家小技巧按染色体区间并行统计如果项目数据量巨大又不想等脚本跑一两个小时可以按染色体拆开来并行算。思路是先利用bcftools query获取染色体列表bcftools query -f %CHROM\n input.vcf.gz | sort -u chroms.txt然后用xargs循环每个染色体分别执行统计脚本最后把结果合并cat chroms.txt | xargs -P 8 -I {} bash -c bcftools view -r {} biallelic.snp.vcf.gz \ | awk ... tmp.{}.txt 并行时要注意CPU和内存别一次开几十个任务把服务器搞死。合并且累加结果时可以用awk按样本名汇总所有tmp文件。这个方法我在200个WGS样本上实测过8线程并行大概能快五六倍值得一试。如果你用的是cyvcf2也可以在Python里用multiprocessing按contig并行思路类似只是需要把每个contig的结果保存下来再合并。5.5 搜索资料时注意区分不同领域的“SNP”最后提醒一句如果你在搜索“SNP统计”时看到射频仿真里的Touchstone文件.snp尤其是HFSS导出SNP文件相关内容那是S参数文件和本文讨论的基因组单核苷酸多态性SNP完全不是一回事。那边说的“snp”是“SnP”n是端口数文件里保存的是不同频率下的S参数而VCF里的SNP是“Single Nucleotide Polymorphism”指的是基因组上单个碱基的变异。这两个领域共用“SNP”这个缩写容易在搜索时造成困扰但处理方法和工具完全不同别把两边内容混在一起看。回到WGS这边统计多样品VCF中每个样本的SNP数核心就是先把统计口径定清楚再选择合适的工具最后用交叉验证确保结果合理。我个人在实际操作中最深的体会是这类任务往往不复杂但因为VCF数据量大、样本列多出错的概率非常高所以每一步都要留下可追溯的命令、样本列表和版本信息这样出了问题才能快速定位。如果你按上面的方法跑完一遍再回头看看最初那份VCF会发现其实就是“定义清楚格式理解流式处理”这十二个字的事。