2025.06.18【WGS】|基于Python与bcftools的VCF文件SNP频率分布可视化实战（附代码）

1. 从VCF到SNP频率分布为什么需要可视化在群体遗传学研究中VCF文件就像一本记录所有样本遗传变异的百科全书。但当我们面对数百个样本、数百万个SNP位点时如何快速把握数据特征这时候SNP频率分布可视化就像给数据装上显微镜能直观展现群体遗传结构。我处理过的一个真实案例某水稻群体测序项目中30个样本的VCF文件包含280万个SNP。用Excel打开直接卡死用文本编辑器查看全是天书。直到用Python画出频率分布图才在5秒内发现样本间存在明显的频率差异模式——这正是后续GWAS分析的关键线索。频率分布图的价值一眼识别异常样本如污染或混样发现群体分层现象评估测序数据质量验证实验设计合理性2. 工具准备bcftoolsPython黄金组合2.1 bcftoolsVCF处理的瑞士军刀bcftools是处理VCF文件的命令行神器相比纯Python解析它有三大优势内存效率高流式读取大文件不卡顿速度快C语言编写的底层实现功能全内置统计、过滤、转换等功能安装方法conda环境conda install -c bioconda bcftools常用频率统计命令# 计算全局等位基因频率 bcftools query -f %CHROM\t%POS\t%REF\t%ALT\t%AF\n input.vcf # 按群体计算频率需先定义群体样本 bcftools view -S population1.txt input.vcf | bcftools query -f %AF\n2.2 Python生态从数据到图表推荐这套工具链pysam高效读取VCF/BCF文件pandas数据清洗与统计matplotlib/seaborn可视化核心plotly交互式图表可选环境配置conda create -n vcf-viz python3.9 conda install -c bioconda pysam conda install -c conda-forge pandas matplotlib seaborn3. 实战四部曲从数据到图表3.1 数据提取bcftools精准狙击假设我们要研究某人群chr1区域的SNP频率先用bcftools提取关键信息# 提取指定区域的AF信息1KG项目格式示例 bcftools query -r chr1:1000000-2000000 -f %CHROM\t%POS\t%REF\t%ALT\t%AF\n ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz chr1_af.tsv提取结果示例chr1 1000001 T C 0.123 chr1 1000002 A G 0.456 chr1 1000003 C T 0.7893.2 数据清洗pandas妙手回春原始数据常存在缺失值、异常值等问题需要清洗import pandas as pd df pd.read_csv(chr1_af.tsv, sep\t, headerNone, names[CHROM,POS,REF,ALT,AF]) # 处理缺失值 print(f原始数据量: {len(df)}) df df[df[AF] ! .] df[AF] df[AF].astype(float) print(f有效数据量: {len(df)}) # 过滤低质量位点 df df[(df[AF] 0.01) (df[AF] 0.99)] print(fQC后数据量: {len(df)})3.3 频率分布统计直方图的秘密不同分箱策略会揭示不同模式我常用两种方法import numpy as np # 自动分箱适合快速探索 auto_bins np.linspace(0, 1, 21) # 手动分箱突出稀有变异 custom_bins [0, 0.001, 0.01] list(np.linspace(0.05, 0.95, 19)) [0.99, 0.999, 1]统计代码freq_counts pd.cut(df[AF], binsauto_bins).value_counts().sort_index()3.4 可视化进阶让图表会说话基础直方图import matplotlib.pyplot as plt plt.figure(figsize(10,6)) plt.hist(df[AF], binsauto_bins, edgecolorblack) plt.xlabel(Allele Frequency) plt.ylabel(SNP Count) plt.title(Global SNP Frequency Distribution) plt.grid(True, alpha0.3) plt.savefig(basic_af_dist.png, dpi300, bbox_inchestight)高级版本添加群体对比# 假设我们已经加载了不同群体的AF数据 populations { EUR: eur_af, EAS: eas_af, AFR: afr_af } plt.figure(figsize(12,7)) for pop, data in populations.items(): sns.kdeplot(data, labelpop, fillTrue, alpha0.3) plt.xlim(0,1) plt.legend(titlePopulation) plt.xlabel(Allele Frequency) plt.ylabel(Density) plt.title(Allele Frequency Spectrum by Population) plt.savefig(afs_by_pop.png, dpi300)4. 完整代码示例从VCF到出版级图表#!/usr/bin/env python3 VCF SNP频率分布可视化工具 输入VCF文件 输出等位基因频率分布图 import subprocess import pandas as pd import matplotlib.pyplot as plt import seaborn as sns from pathlib import Path def extract_af(vcf_path, regionNone): 使用bcftools提取等位基因频率 cmd fbcftools query -f %CHROM\t%POS\t%REF\t%ALT\t%AF\n {vcf_path} if region: cmd fbcftools query -r {region} -f %CHROM\t%POS\t%REF\t%ALT\t%AF\n {vcf_path} result subprocess.run(cmd, shellTrue, capture_outputTrue, textTrue) if result.returncode ! 0: raise RuntimeError(fbcftools执行失败: {result.stderr}) df pd.read_csv(StringIO(result.stdout), sep\t, names[CHROM,POS,REF,ALT,AF]) return df def plot_afs(df, output_path, titleNone): 绘制频率分布图 plt.figure(figsize(12,7)) # 主分布图 sns.histplot(df[AF], bins50, kdeTrue, edgecolornone, alpha0.7) # 添加统计信息 median_af df[AF].median() plt.axvline(median_af, colorred, linestyle--) plt.text(median_af0.02, plt.ylim()[1]*0.9, fMedian AF: {median_af:.3f}, colorred) # 图表修饰 plt.xlabel(Allele Frequency, fontsize12) plt.ylabel(Density, fontsize12) plt.title(title or Allele Frequency Spectrum, fontsize14) plt.grid(True, alpha0.2) plt.savefig(output_path, dpi300, bbox_inchestight) plt.close() if __name__ __main__: import sys from io import StringIO vcf_file sys.argv[1] output_prefix Path(vcf_file).stem try: print(f处理文件: {vcf_file}) df extract_af(vcf_file) # 数据清洗 df df[df[AF] ! .] df[AF] df[AF].astype(float) df df[(df[AF] 0) (df[AF] 1)] # 输出统计 print(f总SNP数: {len(df)}) print(fAF范围: {df[AF].min():.3f}-{df[AF].max():.3f}) # 绘图 plot_afs(df, f{output_prefix}_afs.png, titlefAllele Frequency Spectrum\n{vcf_file}) except Exception as e: print(f错误发生: {str(e)}) sys.exit(1)使用方式python vcf_afs_plot.py input.vcf.gz5. 常见问题与性能优化5.1 内存爆炸试试这些技巧处理大型VCF时遇到过内存不足的问题总结出这些经验分染色体处理添加-r chr1等区域参数流式处理使用bcftools view -H配合Python生成器Dask替代pandas超过1亿个位点时使用# 流式处理示例 def stream_vcf(vcf_path): cmd fbcftools view -H {vcf_path} process subprocess.Popen(cmd, shellTrue, stdoutsubprocess.PIPE, textTrue) for line in process.stdout: fields line.strip().split(\t) af float(fields[4]) # 假设AF在第5列 if 0 af 1: yield af5.2 图形美化期刊级别的调整要让图表达到发表质量我通常做这些调整plt.style.use(seaborn-whitegrid) # 设置样式 sns.set_palette(husl) # 设置调色板 plt.rcParams[font.family] Arial # 设置字体 # 在plot_afs函数中添加 plt.xticks(fontsize10) plt.yticks(fontsize10) plt.grid(True, whichboth, linestyle--, alpha0.5)5.3 交互式探索plotly实战对于需要深入探索的数据推荐交互式可视化import plotly.express as px fig px.histogram(df, xAF, nbins100, titleInteractive AF Distribution) fig.update_layout( hovermodex unified, xaxis_titleAllele Frequency, yaxis_titleCount ) fig.show()6. 扩展应用从基础到进阶6.1 群体间频率差异分析比较不同群体的频率分布能发现选择信号# 计算群体间AF差异 def compare_populations(pop1, pop2): merged pd.merge(pop1, pop2, on[CHROM,POS], suffixes(_pop1, _pop2)) merged[delta_af] merged[AF_pop1] - merged[AF_pop2] return merged # 绘制曼哈顿图 plt.scatter(merged[POS], -np.log10(merged[delta_af].abs()))6.2 功能注释整合结合ANNOVAR注释结果分析不同功能区域的频率分布# 假设已加载注释数据 synonymous df[df[FUNCTION] synonymous] nonsynonymous df[df[FUNCTION] nonsynonymous] plt.figure() sns.kdeplot(synonymous[AF], labelSynonymous) sns.kdeplot(nonsynonymous[AF], labelNonsynonymous)6.3 时间序列分析对于古代DNA数据可以绘制频率随时间变化# 按年代分组计算平均AF time_bins [0, 1000, 5000, 10000] df[time_period] pd.cut(df[AGE], binstime_bins) sns.boxplot(datadf, xtime_period, yAF)