提取VCF文件中样本清单的方法

今天分享的笔记是对VCF文件中的样本ID进行提取和替换的方法，主要是用linux系统下的bcftools和Python实现。

VCF（Variant Call Format）文件通常用于存储基因组变异数据。要从VCF文件中获取样本ID清单，可以使用以下方法：

bcftools

bcftools 是一个常用的VCF文件处理工具，可以使用以下命令快速获取样本ID清单，该命令能快速提取出所有样本ID：

bcftools query -l <your_file.vcf>

Python

可以使用pysam库来读取VCF文件并提取样本ID，以下是个简单示例：

import pysam

# 打开VCF文件
vcf_file = pysam.VariantFile("your_file.vcf")

# 获取样本ID清单
sample_ids = list(vcf_file.header.samples)

# 打印样本ID
for sample_id in sample_ids:
    print(sample_id)

手动查看

打开VCF文件并查看文件头部（以##开头的行），样本ID通常位于以#CHROM开头的行之后。

##fileformat=VCFv4.2
##...
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample1 Sample2 Sample3

如果发现样本的名称不符合要求，如何进行修改替换呢？

方法一：bcftools

可以使用bcftools reheader命令来更改VCF文件中的样本ID。首先，需要创建一个包含新的样本ID的文件，然后使用该文件进行替换。

创建包含新样本ID的文件

假设文件名为new_samples.txt：

Sample1    NewSample1
Sample2    NewSample2
Sample3    NewSample3

bcftools reheader替换

bcftools reheader -s new_samples.txt -o new_file.vcf original_file.vcf

方法二：Python

可以使用pysam库编写Python脚本，读取VCF文件，替换样本ID并保存为新文件。首先安装一下pysam包：

pip install pysam

以下是替换样本ID的脚本：

import pysam

# 原始和新的样本ID映射
sample_map = {
    "Sample1": "NewSample1",
    "Sample2": "NewSample2",
    "Sample3": "NewSample3"
}

# 打开原始VCF文件
vcf_in = pysam.VariantFile("original_file.vcf", "r")

# 创建新的VCF文件
vcf_out = pysam.VariantFile("new_file.vcf", "w", header=vcf_in.header)

# 替换样本ID
for old_sample, new_sample in sample_map.items():
    vcf_in.header.samples.add(new_sample)
    vcf_in.header.samples.remove(old_sample)

# 写入新的VCF文件
for record in vcf_in:
    vcf_out.write(record)

# 关闭文件
vcf_in.close()
vcf_out.close()

方法三：文本处理神器sed

可以使用文本处理工具直接编辑VCF文件头部，但这种方法更适合处理小文件或简单替换。

sed 's/Sample1/NewSample1/; s/Sample2/NewSample2/; s/Sample3/NewSample3/' original_file.vcf > new_file.vcf

此方法直接在文件头部进行字符串替换，适用于简单场景。对于复杂替换或大文件，推荐使用bcftools或Python脚本。根据具体需求和环境，选择合适的方法进行批量修改替换VCF文件中的样本ID。

本文由 mdnice 多平台发布