作为一个plink的忠实依赖者不得不说有时候得换换思维补充学习一下处理vcf文件的软件常用方法。一、bcftools简介bcftools是SAMtools项目的一部分专门用于处理VCFVariant Call Format和BCFBinary Call Format文件。作为生物信息学分析中的核心工具它提供了查看、过滤、统计、格式转换、合并与比较等多种功能是基因组变异数据分析的必备工具。1. MAF过滤使用-i MAF阈值表达式2. 缺失率过滤使用-i F_MISSING阈值表达式3. 复合过滤支持多个条件组合如MAF0.01 F_MISSING0.14. 灵活扩展可结合质量值、深度、HWE等其它过滤条件5. 效率优势在VCF层面直接过滤避免格式转换开销与PLINK相比bcftools在VCF文件处理上更加轻量高效特别适合在流程前期进行快速质控能够有效减少下游分析的数据量。二、安装方法1. Conda安装推荐conda install -c bioconda bcftools2. 源码编译安装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 make sudo make install3. 验证安装bcftools --version三、常用命令详解1. 查看与基本信息提取查看头信息# 查看完整的头信息 bcftools view -h file.vcf.gz # 查看特定元信息 bcftools view -h file.vcf.gz | grep -E ^##(INFO|FORMAT|contig|FILTER)提取样本名称列表# 提取所有样本名到文件 bcftools query -l file.vcf.gz samples.txt # 直接显示样本列表 bcftools query -l file.vcf.gz提取样本# 保留单个样本 bcftools view -s SAMPLE_NAME input.vcf -o output.vcf # 保留多个样本用逗号分隔无空格 bcftools view -s SAMPLE1,SAMPLE2,SAMPLE3 input.vcf -o output.vcf # 从文件中读取样本列表每行一个样本名 bcftools view -S samples_to_keep.txt input.vcf -o output.vcf查看基因型信息# 提取基本变异信息 bcftools query -f %CHROM\t%POS\t%REF\t%ALT[\t%GT]\n file.vcf.gz # 提取深度信息 bcftools query -f %CHROM\t%POS\t%REF\t%ALT[\t%DP]\n file.vcf.gz2. 过滤与筛选常用质控表达式详解过滤条件 表达式 说明MAF过滤 MAF0.01 保留MAF1%的位点缺失率过滤 F_MISSING0.1 保留缺失率10%的位点质量值过滤 QUAL30 保留质量值30的位点深度过滤 INFO/DP10 保留深度10的位点哈迪-温伯格平衡 INFO/HWE1e-6 保留HWE p-value1e-6的位点杂合度过滤 INFO/Het0.2 INFO/Het0.8 保留杂合度在20%-80%的位点与PLINK区别质控项目 bcftools命令 PLINK等效参数MAF过滤 -i MAF0.01 --maf 0.01位点缺失率 -i F_MISSING0.1 --geno 0.1样本缺失率 需结合--missing-indv --mind 0.1哈迪-温伯格 -i HWE1e-6 --hwe 1e-6最小等位基因数 -i MAC3 --mac 3基于质量值过滤# 保留QUAL≥30的位点 bcftools filter -i QUAL30 file.vcf.gz -o filtered.vcf.gz # 保留PASS的位点 bcftools view -f PASS file.vcf.gz -o pass.vcf.gz基于样本基因型过滤# 保留至少有一个样本为杂合子的位点 bcftools view -g het file.vcf.gz -o het.vcf.gz # 保留所有样本均为纯合子的位点 bcftools view -g hom file.vcf.gz -o hom.vcf.gz基于区域过滤# 提取特定染色体区域 bcftools view -r chr1:1000000-2000000 file.vcf.gz -o region.vcf.gz # 使用BED文件提取区域 bcftools view -R regions.bed file.vcf.gz -o bed_regions.vcf.gz基于MAF质控# 过滤MAF 0.01的位点 bcftools view -i MAF0.01 -Oz -o filtered_maf.vcf.gz input.vcf.gz # 过滤MAF 0.05的位点 bcftools view -i MAF0.05 -Oz -o filtered_maf_0.05.vcf.gz input.vcf.gz ############### 方法2先计算再质控 ####### # 计算MAF并添加到INFO字段 bcftools fill-tags input.vcf.gz -Oz -o with_maf.vcf.gz -- -t MAF # 然后过滤 bcftools view -i MAF0.01 with_maf.vcf.gz -Oz -o filtered.vcf.gz基于缺失率质控方法1基于样本缺失率 # 过滤缺失率 10% 的位点 bcftools view -i F_MISSING0.1 -Oz -o filtered_missing.vcf.gz input.vcf.gz # 更严格的过滤缺失率 5% bcftools view -i F_MISSING0.05 -Oz -o filtered_missing_0.05.vcf.gz input.vcf.gz 方法2基于基因型缺失计数 # 过滤缺失基因型数 5 的位点 bcftools view -i N_MISSING5 -Oz -o filtered_missing_count.vcf.gz input.vcf.gz综合质控MAF 缺失率 其他指标复合条件过滤 # 同时过滤MAF、缺失率和质量值 bcftools view -i MAF0.01 F_MISSING0.1 QUAL30 \ -Oz -o qc_filtered.vcf.gz \ input.vcf.gz # 更全面的质控条件 bcftools view \ -i MAF0.01 F_MISSING0.1 QUAL30 INFO/DP10 \ -Oz -o strict_qc.vcf.gz \ input.vcf.gz 完整质控流程 #!/bin/bash # 文件名vcf_qc_pipeline.sh # 功能VCF文件完整质控流程 INPUT_VCFinput.vcf.gz OUTPUT_VCFqc_filtered.vcf.gz # 步骤1计算统计信息可选 bcftools stats $INPUT_VCF input.stats.txt # 步骤2添加MAF注解如果不存在 bcftools fill-tags $INPUT_VCF -Oz -o with_tags.vcf.gz -- -t MAF,F_MISSING # 步骤3应用质控过滤器 bcftools view \ -i MAF0.01 F_MISSING0.1 QUAL20 \ -Oz -o temp_filtered.vcf.gz \ with_tags.vcf.gz # 步骤4移除单态位点可选 bcftools view -c 1 temp_filtered.vcf.gz -Oz -o $OUTPUT_VCF # 步骤5创建索引 bcftools index $OUTPUT_VCF # 步骤6验证结果 bcftools stats $OUTPUT_VCF qc_filtered.stats.txt echo 质控完成 echo 输入文件变异数: $(bcftools view -H $INPUT_VCF | wc -l) echo 输出文件变异数: $(bcftools view -H $OUTPUT_VCF | wc -l) # 清理临时文件 rm with_tags.vcf.gz temp_filtered.vcf.gz3. 统计与质量评估基本统计# 生成基本统计报告 bcftools stats file.vcf.gz stats.txt # 查看统计摘要 bcftools stats file.vcf.gz | grep ^SN | cut -f 3-等位基因频率统计# 计算等位基因频率 bcftools fill-tags file.vcf.gz -o with_af.vcf.gz -- -t AF # 提取AF信息 bcftools query -f %CHROM\t%POS\t%REF\t%ALT\t%AF\n with_af.vcf.gz4. 格式转换与处理VCF格式转换# VCF转BCF格式 bcftools convert -O b file.vcf.gz -o file.bcf # BCF转VCF格式 bcftools convert -O v file.bcf -o file.vcf.gz # 压缩VCF文件 bcftools view file.vcf -Oz -o file.vcf.gz提取特定字段# 提取INFO字段 bcftools query -f %CHROM\t%POS\t%INFO/DP\t%INFO/AF\n file.vcf.gz # 提取FORMAT字段 bcftools query -f %CHROM\t%POS[\t%GT:%DP:%GQ]\n file.vcf.gz5. 合并与比较合并多个VCF文件注意这个在重测序数据中将特别重要因为文件太大所需内存大大部分情况下需要按照染色体拆分# 合并相同样本的不同染色体文件 bcftools concat chr1.vcf.gz chr2.vcf.gz chr3.vcf.gz -o merged.vcf.gz # 合并不同样本的VCF文件 bcftools merge sample1.vcf.gz sample2.vcf.gz sample3.vcf.gz -o all_samples.vcf.gz比较两个VCF文件# 找出两个文件的差异 bcftools isec -p output_dir file1.vcf.gz file2.vcf.gz # 输出共有位点 bcftools isec -p output_dir -n2 file1.vcf.gz file2.vcf.gz提取指定染色体标记/拆分VCF# 命令1提取指定染色体的VCF数据 bcftools view -r ${CHR} -Oz -o chr${CHR}.vcf.gz input.vcf.gz # 命令2为提取的文件创建索引 bcftools index chr${CHR}.vcf.gz bcftools view 查看/提取VCF/BCF文件 bcftools的核心子命令用于查看、筛选、转换VCF/BCF文件 -r ${CHR} 指定区域/染色体 从输入文件中提取指定染色体或基因组区域的数据 ${CHR}是变量如1、chr1、2:1000-5000等 -Oz 输出压缩格式 组合参数 -O指定输出格式 z表示bgzip压缩的VCF格式.vcf.gz -o chr${CHR}.vcf.gz 输出文件 指定输出文件名通常以染色体编号命名 input.vcf.gz 输入文件 输入的VCF文件必须是bgzip压缩且已建立索引6. 注释与功能预测添加注释信息# 使用VEP注释需先安装VEP bcftools annotate -a annotations.vcf.gz -c INFO file.vcf.gz -o annotated.vcf.gz # 添加自定义注释 bcftools annotate -x INFO/OLD_TAG -a new_info.vcf.gz -c INFO/NEW_TAG file.vcf.gz四、实战示例示例1从原始VCF到过滤后的分析文件# 1. 查看样本列表 bcftools query -l raw.vcf.gz samples.txt # 2. 过滤低质量位点 bcftools filter -i QUAL20 DP10 raw.vcf.gz -o filtered.vcf.gz # 3. 提取特定区域 bcftools view -r chr1:1000000-5000000 filtered.vcf.gz -o target_region.vcf.gz # 4. 统计过滤效果 bcftools stats raw.vcf.gz raw_stats.txt bcftools stats target_region.vcf.gz filtered_stats.txt示例2批量处理多个样本#!/bin/bash # 批量处理脚本 for vcf in *.vcf.gz do # 提取样本名 sample$(bcftools query -l $vcf) # 过滤并重命名 bcftools filter -i QUAL30 FORMAT/DP10 $vcf -o ${sample}_filtered.vcf.gz # 统计 bcftools stats ${sample}_filtered.vcf.gz ${sample}_stats.txt done示例3创建分析报告#!/bin/bash # 生成分析报告 echo VCF文件分析报告 report.txt echo 生成时间: $(date) report.txt echo report.txt # 基本信息 echo 1. 文件基本信息 report.txt echo 文件: $1 report.txt echo 样本数: $(bcftools query -l $1 | wc -l) report.txt echo 位点数: $(bcftools view $1 | grep -v ^# | wc -l) report.txt echo report.txt # 质量分布 echo 2. 质量值分布 report.txt bcftools query -f %QUAL\n $1 | awk BEGIN {print QUAL范围\t计数} {if($110) a; else if($120) b; else if($130) c; else d} END {print 0-9\t a \n10-19\t b \n20-29\t c \n30\t d} report.txt五、实用技巧与注意事项1. 性能优化# 使用多线程加速 bcftools view --threads 8 large.vcf.gz -o output.vcf.gz # 设置缓存大小 bcftools view --max-mem 4G huge.vcf.gz -o output.vcf.gz2. 常见问题解决内存不足使用--max-mem参数限制内存使用文件过大先按染色体分割处理再合并结果特别是重测序数据几百G那种编码问题确保VCF文件使用UTF-8编码3. 最佳实践建议始终使用压缩格式.vcf.gz比.vcf节省大量空间建立处理流水线将常用命令封装为脚本记录处理步骤保存所有命令到日志文件验证结果使用bcftools stats比较处理前后数据六、总结bcftools作为VCF文件处理的瑞士军刀其功能强大且灵活。掌握这些常用命令可以显著提高基因组变异数据分析的效率。建议初学者从查看和过滤命令开始逐步学习更高级的功能。在实际工作中结合具体需求灵活组合这些命令可以构建出高效的数据处理流程。学习资源推荐官方文档bcftools(1)GitHub仓库GitHub - samtools/bcftools: This is the official development repository for BCFtools. See installation instructions and other documentation here http://samtools.github.io/bcftools/howtos/install.html · GitHubBiostars论坛https://www.biostars.org/上的bcftools相关讨论希望这篇指南能帮助你在生物信息学分析中更高效地使用bcftools如果有任何问题或建议欢迎在评论区留言讨论。