企业官网建设流程全解析

上周隔壁课题组找我帮忙，他们的变异检测结果已经拿到了VCF文件，但问题很实际：每个样本到底有多少个SNP？听起来就是一行统计的事，可真上手去做，不少人会在这里卡上一整天。用Excel打开VCF，几十万行直接卡死；用bcftools stats跑完，发现它默认按位点整体统计，不按样本拆分；自己写Python脚本，又搞不清楚FORMAT里GT、AD、DP这些字段到底哪个该用来判断“这个样本在这个位点有没有变异”。这篇我就把平时处理这类需求的完整流程写出来，从bcftools安装到纯Python解析，再到生产环境用的cyvcf2方案，附带一堆实测踩坑记录，适合刚拿到VCF文件、需要快速给课题组交统计结果的生信同学参考。

1. 拿到VCF别急着写代码：先弄清楚要统计什么、口径怎么定

1.1 VCF这张表的结构：前8列是位点，后面才是样本

VCF（Variant 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 DP=14 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这一列用冒号分隔了多个字段，其中GT（Genotype）一定排在最前面，是判断样本有没有变异的根本依据。很多新手一上来就盯着INFO里的DP或AD类字段看，那是测序深度和等位基因支持数，不等于基因型判定结果。

1.2 三种常见统计口径，先跟提需求的人对齐

“每个样本的SNP统计”这句话有歧义，我在实际工作中至少遇到三种不同理解，如果不先对齐，后面代码白写：

大多数课题组问“每个样本多少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 bcftools=1.19 conda activate bio

如果觉得conda装软件时解析依赖太慢，建议换成mamba：

conda install -n base -c conda-forge mamba mamba create -n bio -c conda-forge -c bioconda bcftools=1.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.13，CentOS则可能停留在1.9。老版本不是不能用，而是部分新功能缺失，比如bcftools plugin相关命令和新的-Ou压缩输出行为。如果你只是做简单统计，系统包够用；如果你后续要跑比较复杂的过滤、合并或分割功能，建议还是上conda或源码编译，避免因为版本差异浪费时间。

2.3 源码编译路线：依赖处理与经典报错

在没有root权限的服务器上，源码编译是唯一可行方案。bcftools的源码包可以从GitHub官方仓库下载，推荐下载release版的tar.bz2：

wget 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

编译过程最大的坑是系统缺少基础依赖库，我列几个常见报错及对应解决方式：

在没有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 | head

bcftools 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 sampleA=0/1 sampleB=0/0 sampleC=1/1 chr1 100078 C T sampleA=0/0 sampleB=./. sampleC=0/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 (i=1; i<=NF; i++) { split($i, a, "=") gt = a[2] gsub(/\|/, "/", gt) n = split(gt, g, "/") has_var = 0 for (j=1; j<=n; 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 (k=1; k<=n_alt; k++) if (length(a[k]) != 1) is_snp = 0 if (!is_snp) next for (i=3; i<=NF; i++) { split($i, b, "=") gt = b[2] gsub(/\|/, "/", gt) n = split(gt, g, "/") has_var = 0 for (j=1; j<=n; 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' 'NR<=10 {print}'

挑前10个位点，人工数一下每个样本的变异型数，再和3.2节的统计结果比对。这个小动作看似笨，但能把“命令写错了但看着像那么回事”的风险降到最低。我在带新人时反复强调：管道的每一环都可能有隐性问题，不要信任一条没验证过的统计命令的输出。

4. Python逐行解析VCF：不依赖第三方库的统计脚本

4.1 脚本设计：三步走

虽然bcftools已经能解决统计问题，但很多场景下大家还是希望有一个不依赖外部工具的Python脚本，原因不外乎几点：一是最终要集成到已有的Python分析流程中，二是需要对统计结果做更多自定义加工，三是团队环境里不一定有权限装bcftools。既然标题是以Python实战为主，这一步才是重头戏。

纯Python解析VCF的思路可以拆成三步：

1. 逐行读取文件，跳过##开头的元信息；
2. 遇到#CHROM行时解析样本名列表；
3. 对每个非表头行，先判断是否为SNP位点，再解析每个样本的GT字段并计数。

4.2 完整代码及各段逻辑注释

我用gzip模块读取压缩的VCF，这样无需手动解压，也不会占用大量磁盘空间。判断SNP位点的方法是看REF和ALT是否为单碱基：REF必须长度1；ALT可能有多等位基因，用逗号分隔，所有ALT都必须长度1。GT的解析则要注意phased的竖线符号|，统一替换成斜杠/再判断。

import gzip from collections import defaultdict def count_snp_per_sample(vcf_path, only_pass=True): 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() # 判断SNP：REF必须是单碱基，所有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_pass=False) for sample, count in sorted(result.items(), key=lambda x: x[1], reverse=True): print(f"{sample}\t{count}")

这段代码默认only_pass=False只统计全部位点，因为很多实际VCF的FILTER列很乱，PASS、.混着来，如果不加判断就把有效位点漏掉了。脚本里已经预留了only_pass开关，想只统计PASS位点时改成only_pass=True即可，但要看清楚自己文件的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 counts

variant.genotype.array()返回的是一个二维数组，形状是(样本数, 2)，每个位置的取值是等位基因索引，0表示参考等位基因，1、2等表示ALT等位基因索引，-1表示缺失。判断a1 > 0 or a2 > 0的意思是这个样本在该位点上至少携带一个ALT等位基因，正好对应GT=0/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脚本的另一个不可替代的价值是：当你需要自定义统计规则时，改起来非常直观。比如有同事问“只统计纯合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%的问题能定位：

1. 先确认统计口径：对方要的是SNP还是所有变异（含INDEL）？要的是携带变异的样本数还是位点出现次数？这两个口径差出来能到3倍以上。
2. 确认是否过滤了PASS：原始VCF里每个位点是否都通过质检？对方报告里有没有“QC passing”字样？
3. 确认参考等位基因判断：VCF里REF是参考等位基因，ALT是变异等位基因，有些工具输出时把二者调转，统计就歪了。
4. 用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。刚开始觉得多此一举，但当你一个月后回来看结果文件、或者把文件发给协作者时，你会感谢当初这个节约沟通成本的决定。