企业官网建设流程全解析

做生信分析的人，几乎都绕不开VCF文件。无论你是做全基因组测序、外显子组测序，还是基因panel，最终拿到的高通量变异结果基本都是VCF格式。最近好几个朋友都在问我同一个问题：手里有一批样本的VCF，怎么快速知道每个样本到底各有多少个SNP？这个需求听起来简单，但真要统计得干净、准确，还要能批量产出结果，还是有不少细节坑的。这篇文章就把我自己实测下来最顺的一套流程完整分享出来：先用bcftools做快速统计打底，再写一个Python脚本把结果整理成干净表格，顺带把bcftools的安装和常见坑都讲清楚，保证你照着做也能5分钟跑出结果。

这套方案适合的人群很广：刚入门的生信学生、做临床样本分析的技术员、需要批量整理变异数据的科研人员都适用。不需要你很懂Python，也不需要你把VCF格式背下来，只要跟着步骤走就行。我会把背后的原理、为什么要这么处理、哪个环节容易出错，都尽量说透，这样你以后遇到类似场景也能自己举一反三。

1. VCF与SNP统计：先搞清楚我们在处理什么

1.1 VCF文件里的核心信息长什么样

VCF的全称是Variant Call Format，是存放基因变异信息的标准文本格式。之所以叫“标准”，是因为它能同时容纳多个样本的变异信息，并且每个位点的等位基因、基因型、质量值、过滤标签都按固定列顺序排列。VCF文件大致分两部分：以##开头的元信息行，和以#CHROM开头的表头行。表头之后每行代表一个变异位点，列顺序是固定的：CHROM、POS、ID、REF、ALT、QUAL、FILTER、INFO、FORMAT，后面跟着的就是每个样本的基因型数据。

我刚开始接触VCF时也觉得这格式又长又枯燥，但只要你理解了它“一列一个样本、每格一个基因型”的结构，后面做统计就顺理成章了。真正重要的是FORMAT列和后面每个样本列之间的对应关系。FORMAT里通常写着GT:AD:DP:GQ:PL，GT是基因型，AD是各等位基因的深度，DP是总深度，GQ是基因型质量，PL是三种基因型的Phred似然值。我们要做SNP统计，主要盯住GT就够了，其他字段只有在过滤低质量位点时才会用到。

1.2 什么是“每个样本的SNP统计”

SNP是单核苷酸多态性，也就是单个碱基位置发生变异的现象。从VCF文件的角度看，只要某个位点的REF和ALT分别是单个碱基，比如A>G、C>T，这个位点就是SNP。相比之下，插入缺失、结构变异都不属于SNP的范畴，统计时要格外注意区分。

“每个样本的SNP统计”这句话再拆开，意思就是：在这个VCF文件包含的所有样本中，逐个计算每个样本在多少SNP位点上携带了非参考等位基因。这里有几种常见计数口径：可以只数非参考纯合位点个数，可以只数杂合位点个数，也可以两者加起来作为该样本的SNP总数。实际项目里，我通常会把总数、杂合数、纯合突变数、缺失数一起统计出来，看起来只是一个数字，但多拆几列会给后续分析省下大量时间。

不同环境和不同研究目的下，这个统计结果的含义也会有差异。比如医学外显子项目里，样本的SNP总数异常偏高，可能提示样本污染或测序深度不均；育种研究里，不同品系间的SNP数量差别大，通常是亲缘关系或基因组多态性差异的直观体现。所以不要小看这一步基础统计，很多时候它是判断数据质量的第一道关口。

1.3 哪些场景最需要这个统计

最常见的场景是质控。拿到一批新样本的VCF之后，我会先按样本统计SNP数量，如果某个样本明显偏离群体平均水平，比如是其他样本的三倍，那就要警惕是不是样本建库或分析流程出了岔子。另一个高频场景是样本分组比较，比如突变体与野生型、用药组与对照组，想快速看组间是否有整体差异，这时候一张“样本名+SNP数”的表格就是最直接的依据。

这个统计还有一个很实用的场景，就是筛选代表样本。有些项目不需要分析全部样本，只想挑几个遗传信息最丰富的个体做重测序或单细胞验证，这时候每个样本的SNP总数就能帮你快速排序，把多态性最高的几个样本挑出来。我自己做群体遗传分析时，也经常先用这个统计给样本集做个“体检”，确认样本编号、分组信息是否和数据本身匹配。

2. bcftools：最好用的VCF处理瑞士军刀

2.1 为什么选bcftools而不是自己写脚本硬解

我最早做VCF统计的时候，确实写过一段纯Python遍历VCF的脚本，后来数据量一上来就发现不行。几百个样本、上百万个位点的VCF，纯Python逐行解析一夜都不一定能跑完。bcftools是samtools团队维护的VCF/BCF处理工具，底层用C实现，性能远高于我们自己写的Python解析逻辑，而且生态成熟，几乎所有的生信流程都会用到它。

bcftools最大的优势在于它把所有常用的VCF操作都封装成了子命令：view按区域或样本过滤、query灵活提取字段、stats直接生成统计报告、isect处理交集，还有index、norm等配套工具。你做SNP统计时，其实只需要其中两三个子命令就能完成大部分工作。更关键的是，bcftools对VCF.gz格式支持得非常好，不需要手动解压，这一点在处理几个GB级别的文件时简直救命。

2.2 三种安装方式，实测哪种最省心

bcftools的安装方式很多，我按实际体验排个序：conda最省心，系统包管理器次之，源码编译适合特殊需求。

如果你已经装了conda或mamba，只需要一句话：

conda install -c bioconda -c conda-forge bcftools

这个方式会帮你自动处理依赖和版本兼容，尤其适合不太想折腾环境变量的朋友。装完直接bcftools --version验证，如果提示找不到命令，检查一下你当前的conda环境是否激活。

如果用Linux系统自带的包管理器，Ubuntu和Debian系列可以试：

sudo apt update sudo apt install bcftools

CentOS系列则是：

sudo yum install bcftools

这类方式的好处是系统集成度高，缺点是版本可能偏旧。旧版bcftools在部分过滤语法和stats输出格式上和8.x版本有差异，如果你发现自己的命令和网上教程对不上，先看版本。

最后一种是从源码编译。到GitHub的samtools/bcftools仓库下载源码，按INSTALL文档配置，适合需要特定编译选项、或者要在没有管理员权限的服务器上安装的情况。源码编译对初学者不太友好，可能踩到依赖库缺失的坑，但如果你的服务器环境特殊，这也是唯一能走通的路。任何方式装完之后，都建议立刻验证：

bcftools --version

看到类似bcftools 1.19这样的输出，说明安装成功。如果提示bcftools: command not found，先确认有没有把conda环境或安装目录的bin路径加进PATH。

2.3 和统计样本SNP最相关的几个子命令

bcftools query是我最常用的子命令之一。它像一个字段提取器，能按你需要的方式从VCF里抽出指定的列，然后转成自由格式的文本。做每个样本的SNP统计时，query -l可以列出所有样本名称，query -f可以自定义输出格式，比如同时输出样本名、基因型、染色体和位置。

bcftools stats是另一个核心命令，它会对整个VCF文件生成一份非常详细的统计报告，里面包含了变异类型分布、转换/颠换比、每个样本的SNP计数，甚至能按质量值和深度分布给出统计。我们待会要用的“SNP counts by sample”就在这里。bcftools view则负责过滤，比如只留下PASS位点、只留下SNP类型、只留下某几个样本。

这三个子命令配合使用，几乎能覆盖所有日常统计需求。把这个组合理解清楚，后面做起分析来非常顺手。

3. 五分钟起步：用bcftools完成每个样本的SNP计数

3.1 先想清楚统计口径，再做命令

很多人在这一步翻车就是因为没想清楚要统计什么。同样是“每个样本的SNP数量”，可以指样本携带的非参考等位基因位点数，也可以指样本中所有变异位点里SNP类型的总数，还可以要求只统计通过了质量过滤的位点。所以我建议你在跑命令之前先明确三点：只看SNP还是包含Indel？只统计PASS位点还是全部位点？杂合和纯合是分开统计还是加总？这三点直接决定了你的命令参数怎么加。

以最常用的口径为例：我只统计通过FILTER标签为PASS的SNP位点，并且把每个样本的杂合、纯合突变和总数分别列出来。这个需求用bcftools一行就能摸清样本列表，再用stats直接出样本统计。

先列出VCF文件里有哪些样本：

bcftools query -l your_file.vcf.gz

看到返回的样本名列表后，执行：

bcftools stats your_file.vcf.gz > stats.txt

然后从stats.txt里找到下面这段：

# SNP counts by sample: # [1] id [2] number of SNPs [3] number of transitions (ts) [4] number of transversions (tv) [5] number of ts/tv sample_A 1532 1020 512 1.99 sample_B 1478 998 480 2.08

stats默认按VCF里每个样本分别统计，输出的第二列就是每个样本的SNP数量。这个结果里的SNP数量是统计了所有位点的，如果你想先过滤到PASS，可以先用view做一次过滤再管道传给stats：

bcftools view -i 'FILTER="PASS"' your_file.vcf.gz | bcftools stats - > stats_pass.vcf.gz

注意管道传给bcftools stats -时，从stdin读取的VCF必须是未压缩文本格式，否则会报错。

3.2 让bcftools只统计某一组样本

有些场景不需要所有样本，只想看某个子集。假设样本名单存在samples_to_check.txt文件里，一行一个样本名，可以用--samples参数指定：

bcftools stats --samples samples_to_check.txt your_file.vcf.gz > sub_stats.txt

这个操作只对给定样本做统计，运行速度也快很多。你还可以在stats文档里看到--samples-file的用法，实际上就是从文件读取样本子集。配合bcftools view -s同样能实现“先筛样本后统计”，但stats的--samples更直接，不会改动原始文件，也不会因为样本名顺序不同产生额外文件。

如果你手头不是VCF.gz而是普通的VCF，bcftools也能读，只是大文件不压缩读起来会很慢。建议所有大VCF都先压缩并建索引：

bgzip your_file.vcf bcftools index -t your_file.vcf.gz

这一步几乎不会出错，但很多人容易漏掉。没有索引的文件在很多bcftools子命令里会直接报错，比如“Failed to open index”。

3.3 理解stats输出，别把Indel混进SNP里

bcftools stats的输出很长，很多人一看就眼晕。这里我最想强调的一点是：SNP counts by sample这一段统计的是SNP，不是所有变异。bcftools内部会区分SNP与Indel，所以直接看这一段就能拿到干净的SNP计数，不需要你再另外过滤。

如果你觉得自己用Python写会更可控，也可以先用bcftools query把每个样本的GT字段全部导出：

bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE=%GT]\n' your_file.vcf.gz > per_sample_gt.txt

然后让Python去读这个精简文件。这样就把“解析VCF”这个重活交给bcftools，Python只做纯文本统计，速度和稳定性都会好很多。这种方法尤其适合那种需要做定制化统计、但我又想避开复杂VCF解析逻辑的时候。

4. Python批量处理：把统计结果变成干净表格

4.1 为什么已经有了bcftools，还要写Python

bcftools stats确实能输出每个样本的SNP计数，但它输出的是文本报告，不是直接可用的数据分析表格。你想做后续的样本分组比较、画图、筛选异常样本，都得先把结果整理成CSV或DataFrame。这时候Python就有优势了：整合多个统计维度、按组求均值、生成可视化图表都非常方便。

另一个用Python的原因是自定义统计逻辑。比如你想同时统计每个样本的杂合SNP数、纯合SNP数、缺失率，还想去掉某些低质量位点，bcftools stats给的信息不够细，但写Python脚本能完全按你的规则来。我这里分享一个我自己常用的脚本，它不依赖pysam，只读文本VCF，结构简单，容易读懂和修改。

import gzip from collections import defaultdict def parse_vcf_snp_stats(vcf_path, pass_only=True): # 记录每个样本的统计结果 stats = defaultdict(lambda: {"total_snp": 0, "het": 0, "hom_alt": 0, "missing": 0}) opener = gzip.open if vcf_path.endswith(".gz") else open with opener(vcf_path, "rt") as fin: sample_names = [] for line in fin: if line.startswith("##"): continue if line.startswith("#CHROM"): header = line.strip().split("\t") sample_names = header[9:] continue fields = line.strip().split("\t") chrom, pos, ref, alt, filt = fields[0], fields[1], fields[3], fields[4], fields[6] # 只统计SNP：REF和所有ALT都必须是单碱基 if len(ref) != 1: continue if any(len(a) != 1 for a in alt.split(",")): continue if pass_only and filt != "PASS": continue # 从第10列开始是每个样本的数据 sample_data = fields[9:] for idx, sample in enumerate(sample_names): gt_field = sample_data[idx].split(":")[0] gt_format = gt_field.replace("|", "/") alleles = gt_format.split("/") if "." in alleles: stats[sample]["missing"] += 1 continue if len(alleles) == 2: a1, a2 = int(alleles[0]), int(alleles[1]) if a1 == 0 and a2 == 0: # 参考纯合，不计入SNP总数 continue elif a1 != a2: stats[sample]["het"] += 1 stats[sample]["total_snp"] += 1 else: stats[sample]["hom_alt"] += 1 stats[sample]["total_snp"] += 1 return sample_names, stats if __name__ == "__main__": vcf_file = "your_file.vcf.gz" samples, result = parse_vcf_snp_stats(vcf_file) print("sample\ttotal_snp\thet\thom_alt\tmissing") for s in samples: r = result[s] print(f"{s}\t{r['total_snp']}\t{r['het']}\t{r['hom_alt']}\t{r['missing']}")

这个脚本里有一个细节值得注意：GT字段里可能用|分隔等位基因，也可能用/分隔，前者通常表示已经分型（phased），后者表示未分型。统计时我都会先统一成/再处理，避免漏掉。另外，面对多等位基因位点比如A>C,G，脚本的SNP判断逻辑会依据ALT是否都是单碱基，但真正的基因型计数仍依赖于样本的GT是0/1还是0/2等情况，目前我把它统一视为非参考等位基因存在，所以只要不是0/0就计入总数，这在绝大多数场景是够用的。

4.2 脚本运行起来，输出长这样

假设你的VCF里有sample_A和sample_B两个样本，运行上面的脚本后，终端会输出：

sample total_snp het hom_alt missing sample_A 1532 1008 524 12 sample_B 1478 976 502 35

这样一张表格比bcftools stats输出更直观，可以让你一眼看出不同样本的杂合与纯合突变比例。如果你还想继续做样本间比较，可以把输出重定向到CSV文件，或者在Python里直接生成pandas DataFrame。修改一下print部分，或者把samples和result直接喂给pd.DataFrame都非常容易。

这里额外说一个统计口径的坑：有的项目里，样本在某个位点的GT是0/0，但位点本身在群体里是一个已知SNP，这个“0/0”是不是要算进该样本的SNP里？我的习惯是不算，因为这个统计看的是“该样本携带的非参考等位基因”，0/0表示该样本在这个位点没有变异。但如果你的项目关心的是“该样本在多少已知SNP位点有基因型数据”，那统计逻辑就完全不同了。所以写脚本之前，先明确统计口径，这一步比任何代码优化都重要。

4.3 用Python给统计结果画图，直观看异常样本

拿到每个样本的SNP统计结果后，另外一个很有价值的操作是画箱线图或柱状图，快速发现离群样本。我这里用matplotlib画一个最基础的分组柱状图：

import matplotlib.pyplot as plt samples = ["sample_A", "sample_B", "sample_C", "sample_D"] total_snps = [1532, 1478, 2015, 1490] hets = [1008, 976, 1350, 999] hom_alts = [524, 502, 665, 491] x = range(len(samples)) plt.figure(figsize=(10, 6)) plt.bar(x, hets, label="het", color="#4C72B0") plt.bar(x, hom_alts, bottom=hets, label="hom_alt", color="#DD8452") plt.xticks(x, samples, rotation=45) plt.ylabel("SNP count") plt.title("SNP counts per sample") plt.legend() plt.tight_layout() plt.savefig("snp_counts_per_sample.png", dpi=150)

从图上一眼就能看出sample_C的SNP总数明显偏高，这时候我就知道要去查这个样本是否有样本污染、测序覆盖度是否异常，或者是否来自遗传背景差异较大的个体。画图看起来只是锦上添花，但在我实际项目里，它经常能第一时间暴露问题。

4.4 大规模VCF文件下的性能优化思路

如果你处理的VCF文件动辄几百个样本、上千万个位点，纯Python逐行解析会变得很慢。我实测过一个300样本、约800万位点的VCF，用上面的纯Python脚本大概要跑七八分钟，勉强能接受。如果还想再提速，有两条路可以走。

第一条是先用bcftools把文件里的样本数和位点数缩减，只保留需要分析的样本，还有通过质量过滤的位点，再交给Python处理。比如：

bcftools view -S samples_selected.txt -i 'FILTER="PASS"' your_file.vcf.gz > filtered.vcf.gz

然后再对filtered.vcf.gz跑Python脚本，速度能快好几倍。第二条路是使用pysam库，它是一个Python和HTSlib的绑定库，能直接读取BCF/VCF压缩文件，速度接近C工具。如果你的分析逻辑很复杂，需要频繁查询位点，pysam是更好的选择。

import pysam vcf_in = pysam.VariantFile("your_file.vcf.gz") sample_names = list(vcf_in.header.samples) stats = {s: {"total": 0, "het": 0, "hom_alt": 0} for s in sample_names} for rec in vcf_in: if len(rec.ref) != 1: continue if any(len(a) != 1 for a in rec.alts or []): continue if rec.filter.keys() and "PASS" not in rec.filter.keys(): continue for s in sample_names: gt = rec.samples[s].get("GT") if gt is None: continue if -1 in gt: continue a1, a2 = gt if a1 == 0 and a2 == 0: continue if a1 != a2: stats[s]["het"] += 1 stats[s]["total"] += 1 else: stats[s]["hom_alt"] += 1 stats[s]["total"] += 1

这种方式熟悉之后非常顺手，尤其适合后续要做更多变异注释、信息提取的项目。我的建议是：一次性任务用纯Python更安全，重复性的、在流程里要跑很多次的任务，直接把pysam版本封装成函数。

5. 实操中高频遇到的坑与排查技巧

5.1 样本ID和文件编码问题

第一个容易坑人的点是样本名里的隐藏字符。有时候VCF里的样本ID是从外部工具传过来的，带着不可见字符或换行符，用bcftools query -l看不出来，但落到脚本里就成了奇怪的样本名，导致最后的统计表出现一串莫名其妙的多余行。处理方式很简单，拿到样本列表后，先检查长度和字符，如果发现异常就做strip和清洗。

另外，VCF文件本身的编码也值得留意。如果是从Windows环境传过来的文件，可能存在换行符问题，最好在Linux下用dos2unix整理一下再分析。压缩的VCF.gz文件则建议规范压缩名和索引名，不要一会儿叫xxx.vcf.gz一会儿又叫xxx.vcf.gz.tbi，路径写错是再低级但也最常见的错误。

5.2 FILTER和INFO筛选陷阱

很多刚接触bcftools的人会默认stats输出的就是“最终可靠”的SNP数量，其实不然。如果你的VCF里包含了大量低质量位点，而且FILTER列不是PASS，这些位点默认都会被算进去。所以我自己的习惯是，做大范围质控时先过滤：

bcftools view -i 'FILTER="PASS"' input.vcf.gz | bcftools stats - > stats_pass.txt

在Python脚本里，我也在开头就加入了pass_only参数，默认只统计PASS位点。另一个容易被忽视的筛选项是INFO字段里的variant类型。有些VCF会写明INFO/ANN或INFO/CSQ等注释信息，但如果你只关心SNP，还是直接用REF和ALT长度来判断更可靠，因为有些注释工具的过滤标注并不严格等同于变异类型。

5.3 多等位基因位点与GT格式的细节

VCF里的位点不一定是双等位基因，比如0/1和0/2出现在同一个位点，说明这个位置有多条ALT序列。统计样本SNP的时候，如果样本基因型是0/2，那它确实携带非参考等位基因，应当计入该样本的SNP总数。但如果你的目标是“样本有多少个SNP位点发生了变异”，那么无论该位点有多少条ALT，一个样本最多计一次。

我的脚本里就采用了“一个位点一个样本只计一次”的逻辑。这个逻辑的取舍要根据你的项目定义来，我强烈建议在脚本注释里写明这一条，不然过两个星期回来看脚本，自己都会犹豫当时是怎么数数的。

5.4 bcftools提示打不开文件或索引缺失

在bcftools操作中，最经典也最让人着急的报错是“Failed to open index”。出现这个报错，几乎都是因为你的VCF是压缩格式，但没有生成对应的索引文件，或者索引文件是旧版本。解决办法是重新生成索引：

bcftools index -t your_file.vcf.gz

如果你拿到的是一个未压缩的VCF，但文件非常大，也建议先压成bgzip格式再建索引，否则后续很多操作都会变慢。还有一种情况是，你用管道把bcftools view的输出传给bcftools stats，但忘了在stats后面加-，bcftools会认为你要打开一个名为空字符串的文件，直接报错。这个坑我踩过不只一次，命令行里的小短线经常被忽略。

我把这部分常用问题整理成表，方便你排查：

5.5 关于bcftools stats结果的一个小心得

我再分享一个实际项目里非常有用的小技巧：bcftools stats报告里的Ts/Tv（转换/颠换比）也能辅助判断数据质量。正常人类全基因组SNP的Ts/Tv比值大约在2.0左右，如果某个样本的Ts/Tv明显低于1.5，那可能是测序错误率高、或者样本混合出了问题。做每个样本的SNP统计时顺手看一下这个值，比只看数量多了一个质量的判断维度。

6. 把统计结果用起来：扩展方向与进阶思路

6.1 合并多个VCF文件的样本统计结果

有些项目按染色体或区域分析，每个文件只包含部分位点。你想得到全基因组尺度的每个样本SNP总数，最顺滑的做法是先合并VCF，再做统计：

bcftools concat chr1.vcf.gz chr2.vcf.gz chr3.vcf.gz -Oz -o all_chr.vcf.gz bcftools index -t all_chr.vcf.gz bcftools stats all_chr.vcf.gz > all_stats.txt

注意concat要求所有输入文件的样本集合完全一致，位置信息不能重叠，否则会报错。合并后的VCF再做Python脚本统计，结果就是每个样本全基因组范围的SNP完整计数。如果不想合并，也可以在多个文件上分别跑统计，然后写个小脚本把结果按样本累加，但这种分散处理方式容易漏掉某些样本，远不如先合并再统计省心。

6.2 在Python里做样本分组比较

拿到每个样本的SNP统计表格之后，最常见的下一步是分组比较。比如你想比较对照组和处理组的SNP数量是否存在显著差异，直接用pandas和scipy就能做：

import pandas as pd from scipy import stats df = pd.read_csv("snp_counts.csv") group1 = df[df["group"] == "control"]["total_snp"] group2 = df[df["group"] == "treatment"]["total_snp"] t_stat, p_value = stats.ttest_ind(group1, group2) print(f"t={t_stat:.4f}, p={p_value:.4f}")

这个逻辑对很多场景都适用，但要注意先检查数据是否符合正态分布、方差是否齐性。如果不满足，就用Mann-Whitney U检验。这一步看似简单，但会直接影响结论的可靠性。

6.3 将统计结果输出为报告

如果你要给团队成员或合作方交付结果，纯文本或CSV可能不够直观。我一般会把统计结果做成一个包含表格和柱状图的PDF报告。Python里用reportlab或直接生成Markdown再转PDF都可以。重点是把样本ID、SNP总数、杂合数、纯合数、缺失数、Ts/Tv以及分组信息都列全，这样拿到报告的人不需要自己再翻原始文件。

6.4 后续还能做哪些分析

每个样本的SNP统计是整个变异分析链条的起点，接下来通常能接很多方向。如果你有表型数据，可以做GWAS；如果你想研究样本间的关系，可以基于这些SNP计算遗传距离；如果你想做进化分析，可以筛选高多态性SNP位点用于构建系统发育树。这篇文章讲的统计方法虽然基础，但它是后面一切分析的基石，把这个环节做得干净准确，后面每一步都会顺畅很多。

我自己在实操中对这个流程最深的感受是：不要让工具替代你的思考。bcftools帮你快速出数，Python帮你整理和画图，但真正的统计口径、过滤标准和后续分析方向，都需要你根据项目需求来决定。遇到问题时，先把“我要统计什么”想清楚，再去找命令和代码，往往比瞎试工具更高效。

最后再分享一个我个人的小习惯：每次做完SNP统计，我都会把样本列表、原始VCF版本、bcftools版本、Python脚本、输出表格这5样东西放在同一个目录下，并把过滤条件写在脚本注释里。这样即使过了很久再回来复现结果，也不会因为忘了口径而抓狂。