企业官网建设流程全解析

做全基因组测序（WGS）的群体项目时，总有一类问题会被反复问到：这个多样品VCF文件里，每个样本到底有多少个SNP？问的人可能是刚接手数据的同学，也可能是等报告的合作方。问题是同一个，但直接甩一个数字出去很容易翻车——VCF里同一个位点在不同样本中可能是0/1、1/1、0/0，也可能直接是./.；ALT列还可能同时出现A、T、C三个等位基因；不同预处理流程出来的VCF，统计口径还不一样。这篇文章基于我最近跑完的一批WGS联合分析数据，把多样品VCF中每个样本SNP数量的统计方法完整过一遍，内容包括VCF核心字段讲解、命令行工具的适用边界、带详细注释的Python脚本、大规模文件优化办法，以及我实际踩过的几个坑，适合正在做群体WGS、外显子测序或准备出QC报告的同学参考。

1. 多样品VCF统计SNP数的应用场景与需求拆解

1.1 为什么每个样本的SNP数如此关键

先说个多数项目里真实遇到的场景。你从GATK HaplotypeCaller跑完联合genotyping，拿到了一个包含几十甚至上百个样本的VCF，第一件事往往不是直接丢给下游注释，而是先看看这批样本整体做得怎么样。这时候每个样本的SNP数就是最直观的体检指标。

样本间SNP数量差异过大，通常意味着某个样本测序深度不够、文库构建出了问题、样本发生了污染，或者样本标签写错了。比如一个样本如果只有其他样本一半的SNP数，先别急着怀疑生物学差异，大概率是这个样本的测序量崩了或者比对率异常。反过来，某个样本SNP数异常多，可能是混入了其他物种或样本间交叉污染。之前我处理过一个48样本的WGS项目，其中一个样本SNP数比其他样本高了快一倍，后来检查发现是存在跨样本污染，最终在数据质控阶段就把这个问题拦住了，没有污染下游的所有分析。

SNP数还不只是质控指标。在做家系分析时，父母和子代的SNP数通常符合孟德尔遗传规律；在做群体结构分析、构建进化树之前，确认每个样本的变异数量在一个合理区间，是避免后续分析被单个离群样本带偏的关键一步；即使是给合作方出交付报告，一张样本名加SNP数的小表格，也比一堆复杂的统计图更能说明“这批数据能用来干什么”。

所以别小看这个“数数”工作，它被高频问到，背后是整个WGS分析流水线质量的门面，也是最重要的一道情报哨兵。

1.2 从VCF结构理解“行存位点、列存样本”的信息路径

要统计每个样本的SNP数，绕不开对VCF格式的理解。VCF的全称是Variant Call Format，设计思路非常简单：文件里的每一行描述一个变异位点，每一列描述一个维度信息。前8列描述这个位点本身，第9列描述样本字段的格式，第10列开始，每一个样本占一列。

这就是统计SNP数的“信息路径”：你要遍历所有数据行，先判断这个位点是不是SNP，然后依次进入每个样本列，解析该样本在这个位点上的基因型（GT），如果基因型非缺失且携带非参考等位基因，就给这个样本记一个变异SNP位点。

不过VCF有两层复杂性。第一层是同一个位点在不同样本里，基因型可能完全不同：位点是真实存在的变异，但样本A可能是0/1杂合，样本B可能是1/1纯合，样本C可能是0/0纯合参考，样本D可能直接是./.缺失；第二层是VCF里所谓的“位点”形态并不统一，有的位点只有双等位基因，有的位点存在多等位基因，有的位点混着SNP和indel，如果没有预先把统计口径定清楚，同一个VCF在不同人手里能数出完全不同的结果。

所以接下来的内容，我会先把VCF格式里统计需要的核心字段讲透，再给具体方法。先看懂字段，再谈统计，这是最不容易走弯路的方式。

2. 统计前必须吃透的VCF核心字段

2.1 元信息行、表头行和数据行的区分

拿一个常见的VCF打开，前几十行通常是以两个#号开头的元信息行，记录的是VCF版本、FILTER定义、INFO字段含义、样本处理流程等。这些行不参与变异统计，解析时直接跳过即可。

关键是从以单个#号开头的表头行开始：

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001 NA00002 NA00003

这一行定义了后面每一个数据列的含义。注意，从第10列开始才是样本名。实际统计前，用下面命令先确认一下样本名，永远是个好习惯：

# 查看VCF中的样本清单，确认没有重复名、没有未知样本 bcftools query -l input.vcf.gz

数据行则是从第1列染色体名称开始的具体位点记录。我从一个真实的WGS项目里截一个典型示例：

chr1 12009 . C T 100 PASS AC=2;AN=8 GT:AD:DP:GQ:PL 0/0:4,0:4:12:0,36,470 0/1:3,10:13:99:130,0,320 1/1:0,8:8:24:290,24,0

这个位点有3个样本，第一个样本是0/0纯合参考，第二个是0/1杂合变异，第三个是1/1纯合变异。统计的时候，第一个样本在这个位点不应该被计入“变异SNP数”，后两个应该计入。

2.2 FORMAT与样本列：GT才是统计的主战场

第9列FORMAT定义了每个样本列里各个字段的顺序和含义。最常见的形式是：

GT:AD:DP:GQ:PL

每个样本列里，冒号分隔的字段数和FORMAT中定义的一一对应。以样本2为例：

0/1:3,10:13:99:130,0,320

对应的就是GT=0/1，AD=3,10（参考等位基因支持数3、变异等位基因支持数10），DP=13（总深度），GQ=99（基因型质量），PL=130,0,320（三种基因型的Phred打分）。统计SNP数时，绝大部分情况下只需要解析GT字段。

GT字段的取值和含义需要记牢：

这里有个容易混的点：基因型缺失（./.）和纯合参考（0/0）在统计时必须分开。很多项目在汇报“SNP数”时，习惯说“该样本有70万个SNP高质量位点”，说的是带变异的位点数，不包括0/0；也有项目会用“非缺失位点数”来评估基因分型完整性，这时候0/0要被算进去。两种口径都有意义，但一定要在报告中写清楚，否则别人看到数字会对不上。

2.3 多样品VCF中同一位点如何对应多个样本

多样品VCF是多个样本单独call变异后，通过GATK GenotypeGVCFs或bcftools merge等步骤合并出来的。合并后，同一个基因组位置的所有等位基因会被汇总到一行里。

举个例子，样本A在chr1:1000检测到C>T，样本B在同位置检测到C>A，合并后这行就是：

chr1 1000 . C T,A ... GT:AD 0/1:5,3,0 0/2:4,0,6

ALT列同时有T和A两个等位基因，样本A的GT是0/1（参考+第一个ALT的T），样本B的GT是0/2（参考+第二个ALT的A）。这种多等位基因位点如果不提前定好规则，统计时极易出错：是按整个位点给两个样本各记1个SNP，还是按等位基因拆分后给样本B记一个A变异、给样本A记一个T变异？两种做法的总数不同，跟其它工具输出的数字也会对不上。

我自己的建议是：在做“每个样本的SNP数”这种简单指标时，统一按“位点”计数，即一个位点只要该样本存在任意一个非参考等位基因，就计1个变异SNP。如果后续要更细致的等位基因频率分析，再单独拆成等位基因计数的口径，不要混在一起。

3. 先用现成工具快速出数：grep、vcftools、bcftools的实测对比

3.1 纯文本处理方案能做到什么程度

最直接的方法是用grep把注释行去掉，然后awk处理。很多人第一版脚本是这样：

# 粗略统计位点总数，但不区分样本，也不解析GT grep -v "^#" input.vcf.gz 2>/dev/null | wc -l # 注意：上面命令对.gz压缩文件无效，需要先解压或使用zgrep zgrep -v "^#" input.vcf | awk '{n++; if (length($4)==1 && length($5)==1) snp++} END {print "total_var=" n, "snp_site=" snp}'

这个方案只能统计这个文件里“总共有多少个变异位点”“其中多少是SNP”，完全做不到按样本拆分。它适合快速瞄一眼文件的量级，不适合回答“每个样本有多少SNP”这个问题。

另外要提醒一个反直觉的点：VCF是文本文件，但它通常以bgzip压缩存储，普通grep直接读是乱码。虽然grep本身能识别部分文本，但正规写法是使用zgrep或zcat，或者先bcftools view再处理。我见过不止一个同事因为没处理.gz压缩，统计结果总是莫名其妙弹出二进制错误。

3.2 vcftools能算但不算直接

vcftools是立了多年的老工具，常用在VCF过滤和格式转换上。跟样本SNP数统计相关的是下面这组命令：

# 计算每个样本的缺失率 vcftools --gzvcf input.vcf.gz --missing-indv --out wgs_qc # 计算每个样本的等位基因频率 vcftools --gzvcf input.vcf.gz --freq2 --out wgs_qc

wgs_qc.imiss里有一列F_MISS，表示每个样本的基因型缺失比例。如果一个VCF固定有N个SNP位点，理论上可以用“N×(1-F_MISS)”算出每个样本的非缺失位点数。但这个数字里包含了0/0纯合参考，并不是大家通常说的“样本的SNP变异数”。

更重要的问题是，--freq2给的是位点级别的等位基因频率，不是样本级别的计数，没法直接把每个样本的SNP变异数拆出来。所以vcftools在实践中更适合做位点过滤或缺失率检查，用来回答“每个样本有多少个SNP”就比较别扭。如果你的需求刚好是算缺失率，它是个利器；但如果你已经明确要“每个样本的变异SNP数”，直接往下看bcftools或Python。

3.3 bcftools query加awk的推荐组合

bcftools是目前处理VCF最顺手的工具，它的query子命令可以把大矩阵文件裁成任意想要的格式。先是过滤出SNP，再按“样本+GT”逐行展开，配合awk聚合，就能得到每个样本的统计表：

# 保留双等位SNP位点，输出样本名和GT，每行一个样本一个基因型 bcftools view -v snps -m2 -M2 input.vcf.gz -Ou | \ bcftools query -f '[%SAMPLE\t%GT\n]' - | \ gawk ' { cnt[$1]++; # 该样本的总基因型记录数 split($2, a, /[/|]/); # 按/或|拆分GT if (a[1] == "." || a[2] == ".") { miss[$1]++; next } if (a[1] != a[2]) { het[$1]++; next } if (a[1] != "0") { homalt[$1]++; next } homref[$1]++; } END { print "sample\thomref\thet\thomalt\tmiss\tgenotyped"; for (s in cnt) print s "\t" homref[s] "\t" het[s] "\t" homalt[s] "\t" miss[s] "\t" cnt[s]-miss[s]; }' | column -t

这条命令里，-v snps只保留SNP，-m2 -M2保留双等位基因位点，避免多等位基因位点的计数歧义。gawk脚本里，het是杂合变异数，homalt是纯合变异数，两者之和就是该样本在这个VCF里的SNP变异数；homref是0/0位点数，miss是缺失数。

实测下来，中等规模（几十个样本、每个染色体一个VCF）的文件，这条管道命令几秒就能跑完。它的优点是快、不写中间文件，缺点是定制性差一点，如果你想在统计时加入复杂的过滤条件，比如“同时要求DP>=10且GQ>=20”，用awk写起来就比较繁琐，这时候就该用自写Python脚本。

4. Python精确统计方案：完整代码与逐段注释

4.1 什么情况下必须自己写脚本

bcftools方案虽快，但它有几个场景不好覆盖：一是要按VCF的FILTER、QUAL、DP、GQ组合过滤；二是要同时输出纯合、杂合、缺失、参考位点数，并形成固定格式的交付CSV；三是处理多等位基因位点时，想保留更精细的计数规则；四是团队里其他人不太会用命令行，希望有个一键运行的脚本。

这些需求在真实的WGS交付项目里几乎都会遇到。尤其当你要把统计结果直接合并到交付报告里，或者需要跟样本测序深度表、比对率表做关联分析时，一个能控制每列含义的Python脚本比每次手写awk要可靠得多。

4.2 完整脚本：SNP判定、GT解析、缺失处理与CSV输出

下面是我在项目里实际使用的脚本精简版。它的逻辑是：读入VCF（支持.gz压缩），跳过元信息行，判断位点是否为SNP，解析每个样本的GT，按“位点”口径统计每个样本的纯合参考、杂合变异、纯合变异、缺失数，并输出CSV。

#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ 多样品VCF文件中每个样品的SNP数统计 用法: python count_sample_snp.py input.vcf.gz python count_sample_snp.py input.vcf.gz --pass-only 输出: sample_snp_counts.csv """ import gzip import sys from collections import defaultdict def open_maybe_gzip(path): """自动识别.gz文件，统一按文本模式打开。""" if path.endswith('.gz'): return gzip.open(path, 'rt', encoding='utf-8') return open(path, 'r', encoding='utf-8') def parse_gt(gt_field): """ 从样本列中解析GT字段，返回两个等位基因编号。 如果不是合法GT，或等位基因缺失，返回None。 支持 0/1、0|1、1/2、./. 等常见写法。 """ gt = gt_field.split(':')[0] # 只取GT，忽略AD/DP/GQ if gt in ('.', './.', '.|.'): return None if '/' in gt: a1, a2 = gt.split('/', 1) elif '|' in gt: a1, a2 = gt.split('|', 1) else: return None if a1 == '.' or a2 == '.': return None return a1, a2 def is_snp(ref, alt_list): """判断位点是否为SNP：REF长度必须为1，且所有ALT长度也必须为1。""" if len(ref) != 1: return False for alt in alt_list: if alt in ('.', '*'): return False if len(alt) != 1: return False return True def main(vcf_path, pass_only=False): samples = [] counts = defaultdict(lambda: { 'het': 0, 'hom_alt': 0, 'hom_ref': 0, 'missing': 0, 'total_alt': 0 }) with open_maybe_gzip(vcf_path) as f: for line in f: if line.startswith('##'): continue if line.startswith('#CHROM'): cols = line.strip().split('\t') samples = cols[9:] # 第10列开始才是样本名 for s in samples: counts[s] # 预初始化，避免后续KeyError continue parts = line.strip().split('\t') if len(parts) < 10: continue # 可选：只统计FILTER列为PASS或.的位点 if pass_only: if parts[6] not in ('PASS', '.'): continue ref = parts[3] alt_list = parts[4].split(',') if not is_snp(ref, alt_list): continue # 定位GT在FORMAT中的位置 fmt_fields = parts[8].split(':') gt_idx = fmt_fields.index('GT') if 'GT' in fmt_fields else 0 for i, sample in enumerate(samples): sample_fields = parts[9 + i].split(':') allele_pair = parse_gt(sample_fields[gt_idx]) if allele_pair is None: counts[sample]['missing'] += 1 continue a1, a2 = allele_pair if a1 != a2: counts[sample]['het'] += 1 counts[sample]['total_alt'] += 1 else: if a1 == '0': counts[sample]['hom_ref'] += 1 else: counts[sample]['hom_alt'] += 1 counts[sample]['total_alt'] += 1 # 输出CSV out_name = 'sample_snp_counts.csv' with open(out_name, 'w', encoding='utf-8') as out: header = ['sample', 'total_alt', 'het', 'hom_alt', 'hom_ref', 'missing', 'genotyped', 'missing_rate'] out.write(','.join(header) + '\n') for s in samples: c = counts[s] genotyped = c['het'] + c['hom_alt'] + c['hom_ref'] total = genotyped + c['missing'] miss_rate = c['missing'] / total if total > 0 else 1.0 row = [s, str(c['total_alt']), str(c['het']), str(c['hom_alt']), str(c['hom_ref']), str(c['missing']), str(genotyped), f'{miss_rate:.4f}'] out.write(','.join(row) + '\n') print(f'结果已写入 {out_name}') if __name__ == '__main__': if len(sys.argv) < 2: print('Usage: python count_sample_snp.py <input.vcf.gz> [--pass-only]') sys.exit(1) main(sys.argv[1], '--pass-only' in sys.argv[2:])

脚本里几个关键点值得单独说明。

第一，parse_gt函数只关心GT字段，不碰AD和DP。很多初学者容易在这里出问题：样本列里的字段顺序虽然大部分是GT:AD:DP:GQ:PL，但不同流程产出的VCF，FORMAT字段顺序可能有差异，比如有的流程是GT:GQ:DP，有的会把PL放在前面。只取冒号分割后的第一个字段，并确认它是GT，而不是靠固定的split(':')[0]硬编码，能避免很多兼容性问题。我给脚本里做了双重保障——先找到GT在FORMAT中的索引，再去样本列取对应位置。

第二，is_snp函数里对ALT为.和*的处理。VCF规范允许ALT为.表示“该位点没有替代等位基因”，*表示“在参考中有间隙的缺失等位基因”。这两种情况严格来说都不是标准SNP，统计时要排除。否则你会出现莫名其妙的“SNP数偏大”问题，而且很难定位。

第三，输出的CSV里有两列数字要区分清楚：total_alt是“该样本携带非参考等位基因的SNP位点数”，也就是大家通常说的“这个样本有多少个SNP”；genotyped是“该样本在所有SNP位点中非缺失的基因型数量”，它包含了0/0位点，用来评估样本的基因分型完整性。

4.3 脚本输出字段解释与口径定义

脚本跑完后，生成的CSV大概是这样的结构：

这个表里，total_alt是核心指标。如果项目汇报里写“样本NA00002有325万个SNP”，说的就是这个数。missing_rate用来判断样本质量，一般WGS项目里缺失率在1%以下算正常，如果一个样本缺失率超过5%，要么是测序深度不足，要么是样本处理有问题，需要重点排查。

关于口径，建议每个人在交付结果时都写一句说明，比如：“统计的是VCF中所有SNP位点，每个样本携带至少一个非参考等位基因的位点数，不保留多等位位点展开，缺失基因型不计入。”这句话会省掉后续大量沟通成本。

5. 大规模样本与超大VCF的优化方案

5.1 用cyvcf2代替纯Python逐行解析

跑小型VCF时，上面的纯Python脚本完全够用。但到了真正的WGS项目，动辄几百GB的VCF，纯Python逐行解析的速度会让你怀疑人生。实测经验是，纯Python脚本处理一个几GB的VCF可能要几十分钟，而用cyvcf2会快5到10倍以上。

cyvcf2是一个基于Cython封装的VCF解析库，处理速度极快。同样是数SNP，它的基本用法是这样：

from cyvcf2 import VCF vcf = VCF('input.vcf.gz') samples = vcf.samples counts = {s: {'het': 0, 'hom_alt': 0, 'hom_ref': 0, 'missing': 0, 'total_alt': 0} for s in samples} for variant in vcf: if not variant.is_snp: # 直接排除indel和多等位复杂位点 continue gt_types = variant.gt_types # 0=HOM_REF, 1=HET, 2=HOM_ALT, 3=UNKNOWN for i, s in enumerate(samples): if gt_types[i] == 3: counts[s]['missing'] += 1 elif gt_types[i] == 0: counts[s]['hom_ref'] += 1 else: counts[s]['total_alt'] += 1 if gt_types[i] == 1: counts[s]['het'] += 1 else: counts[s]['hom_alt'] += 1 # 之后输出CSV的逻辑与纯Python版一致

variant.is_snp是cyvcf2内置的快速判断属性，variant.gt_types直接把所有样本的基因型类型打包成数组，省去了你手动解析GT字符串的开销。如果你的环境能装这个依赖，大规模数据处理首选它。

5.2 按染色体拆分后用多进程并行统计

即使有了cyvcf2，如果只有一个CPU在跑，几百GB的文件依然慢。并行思路并不复杂：VCF本身就是按染色体排布的，把文件按染色体拆开，每个进程统计一条染色体，再把结果合并起来。

mkdir -p chr_split for chr in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 X Y MT; do bcftools view -r ${chr} input.vcf.gz -Ou -o chr_split/chr${chr}.vcf & done wait

然后对chr_split下的每个文件分别跑统计脚本，最后在Python里按样本名汇总。注意三点：

• 不要一次开几十个进程，普通机械磁盘或网络文件系统的IO会成为瓶颈，反而更慢。我实测最好效果是4到8个并行进程。
• 拆分前确认VCF是按坐标排好序的。大多数联合calling流程产出的VCF都排过序，但如果你手动合并过多个来源的VCF，必须先用bcftools sort排序，否则-r按区域抽取会漏位点。
• 拆分后的临时文件记得用bcftools index建索引，这样后续按区域查询会更快。

如果不想拆文件，也可以考虑用tabix按区域切片。但tabix只在VCF是bgzip压缩且建过索引的前提下才能用，普通文本VCF不适用，这点要记牢。

5.3 统计结果如何与测序深度合并成QC报告

统计完每个样本的SNP数，只是第一步。真正有价值的用法是跟其他QC指标合并起来看。我通常会把三张表join在一起：一张是样本SNP数统计表，一张是测序深度表（可以用samtools depth或mosdepth算出平均深度），一张是对比样本间亲缘关系的IBD矩阵。

一张典型的QC汇总表长这样：

看表格的时候，如果一个样本的total_alt_snp明显偏低，同时mean_depth也低，那就是测序深度的问题；如果SNP数正常但缺失率高，可能是call变异时的过滤参数太严格，或者样本纯度不够；如果SNP数偏高，同时duplication率也高，要警惕样本污染。

这个合并动作很简单，用pandas几行就能搞定：

import pandas as pd snp_df = pd.read_csv('sample_snp_counts.csv') depth_df = pd.read_csv('depth_summary.csv') qc = snp_df.merge(depth_df, on='sample', how='left') qc.to_csv('all_sample_qc.csv', index=False)

合并之后的表，就是可以直接交给下游分析或合作方的最终QC报告。

6. 我在实际跑数中踩过的坑与排查链路

6.1 多等位基因位点重复计数的问题

有一次我给一个100样本的外显子数据集统计SNP数，发现样本之间的SNP数差异比预期大，有些样本甚至比别人多出几万个位点。排查后定位到根因：这些样本的基因组区域刚好落在一些已知的多等位基因热点区域，比如HLA区域，一个位点经常出现三个甚至四个ALT等位基因。

当时我的统计脚本没有对多等位基因做明确处理，导致同一个位点在部分样本里被重复计了多次。修正后的规则很简单：按“位点”计数，不按“等位基因”计数。一个样本在一个多等位位点上，只要存在任意一个非参考基因型，就只计1个SNP。同时我建议，如果项目不需要多等位基因信息，直接用bcftools view -m2 -M2或bcftools norm -m-预处理VCF，把多等位位点拆成双等位位点，再从源头避免口径混乱。

注意：bcftools view -m2 -M2和bcftools norm -m-的语义不同。前者是过滤掉多于两个等位基因的位点，后者是把一个多等位位点拆成多行双等位位点。拆开后，“位点数”会变多，这是正常的，但不是每个人都知道。如果你要用拆分后的文件汇报SNP数，一定要注明“按展开的等位基因计数”。

6.2 重复样本列污染统计结果

另一个比较隐蔽的坑来自样本命名。我接过一个VCF，用bcftools query -l一看，48个样本名里有2个重复的样本，只是其中一个末尾多了_R后缀。查下来发现是该样本在两个批次里都做了测序，联合分析时没有去重，导致统计结果里出现两个高度相似的样本名，让交付报告显得很混乱。

排查链路是这样的：

1. 用bcftools query -l列出所有样本名；
2. 用sort | uniq -d找出重复名；
3. 确认重复样本是否确实是同一个人/同一个DNA文库重测序，而不是样本名编码错误；
4. 如果是重复样本，决定是去掉其中一列，还是重命名后保留两个列分别统计。

处理重复列可以用bcftools reheader --samples sample_names.txt，一行一个样本名，重新定义样本名，把这个隐患在源头解决。

6.3 一件事先说清：VCF里的“SNP文件”和射频仿真SnP文件不是一回事

最后说个容易让新手懵的搜索问题。有时候你去网上搜“snp文件”，会看到大量跟射频仿真相关的结果，尤其是HFSS这类电磁仿真软件导出的.s1p、.s2p等文件。那边叫“SnP文件”，全称是S-parameter Network Parameter，用于描述射频器件端口网络中的散射参数，是一种完全不同的技术领域概念，文件格式长这样：

# MHz S RI R 50 1000 0.8 -0.2 0.3 0.1 2000 0.7 -0.3 0.2 0.2

而生物信息学里的SNP，全称是Single Nucleotide Polymorphism，指单核苷酸多态性，记录在VCF文件中。两者只是英文缩写撞了车，内容没有任何关系。如果你在WGS分析中拿到一个.s1p或.s2p文件，先确认它是不是Touchstone格式——如果是，那就不是变异位点文件，别拿去跑注释或统计变异数，方向完全不对。

如果是从VCF里导出SNP位点，通常做法是bcftools view -v snps过滤之后，再用bcftools query提取位点表：

bcftools view -v snps input.vcf.gz -Ou | \ bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\n' - > snp_sites.tsv

这才是标准的“从VCF导出SNP列表”的姿势。

回到统计这件事上，我最后想分享一个习惯：每次跑这类统计前，我会把“输入VCF的版本、是否经过VQSR、是否经过norm拆分、多等位位点怎么处理”写在脚本的注释里，或者在结果表里附一个备注列。因为一个月或半年后你再看这个结果时，文档比记忆靠谱得多。希望这篇文章能把你在多样品VCF里统计每个样本SNP数时遇到的问题都覆盖到，少走几趟弯路。