Bcftools:从VCF文件中提取样本ID,利用Python脚本批量替换样本编号

今天分享的笔记是对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。

LICENSED UNDER CC BY-NC-SA 4.0 素材来源于互联网公开资料,如有侵权请联系后台删除
Comment