修改vcf文件中的染色体名

问题描述

在生物信息分析过程中，常常需要修改 vcf 染色体名称，即删除染色体名中的chr字符或添加chr字符到vcf文件中。

这里使用shell脚本和bcftools两种方式实现。大致比较了一下，shell脚本在速度上稍微快一些。

1. 通过 shell 脚本实现

下面的shell脚本要求输入vcf为解压缩（因为我的数据是没有压缩的…），如果是压缩文档，自行修改一下代码即可。

通过 awk 添加 chr 字符

awk '{ 
        if($0 !~ /^#/) 
            print "chr"$0;
        else if(match($0,/(##contig=<ID=)(.*)/,m))
            print m[1]"chr"m[2];
        else print $0 
      }' no_chr.vcf > with_chr.vcf

通过 sed 删除 chr 字符

sed 's/##contig=<ID=chr/##contig=<ID=/g' with_chr.vcf | sed 's/^chr//g'  > no_chr.vcf

代码封装

因为需要经常用到上面的两段代码，所以将它们封装到一个名为change_vcf_name.sh的脚本中。

#! /user/bin/bash

echo "usage: change_vcf_name.sh a|r input_vcf_file output_vcf_file"
echo "a (add) : 添加chr字符到vcf文件中"
echo "r (remove): 删除vcf文件中的chr字符"

if [ $1 = "a" ]
then
awk '{
        if($0 !~ /^#/)
            print "chr"$0;
        else if(match($0,/(##contig=<ID=)(.*)/,m))
            print m[1]"chr"m[2];
        else print $0
      }' $2 > $3
elif [ $1 = "r" ]
then
sed 's/##contig=<ID=chr/##contig=<ID=/g' $2 | sed 's/^chr//g'  > $3
else
echo "输入的第一个参数不合适，请输入a或r。"
fi

代码调用：

# 删除vcf中的chr字符
bash change_vcf_name.sh r input.vcf output.vcf

# 加入chr字符到vcf中
bash change_vcf_name.sh a input.vcf output.vcf

2. 通过 bcftools 程序实现

需先准备一个bcftools的配置文件 （这里命名为chr_name_change.txt，内容如下，由tab分隔成两列。

bcftools会根据这个配置文件修改vcf的名称（将左边的名称替换为右边，例如这里会将chr1替换为1、chr2替换为2，…）

chr1 1
chr2 2
chr3 3
chr4 4
chr5 5
chr6 6
chr7 7
chr8 8
chr9 9
chr10 10
chr11 11
chr12 12
chr13 13
chr14 14
chr15 15
chr16 16
chr17 17
chr18 18
chr19 19
chr20 20
chr21 21
chr22 22
chrX X
chrY Y
chrMT M

程序运行代码

• 调用bcftools的annotate子程序
• 调用annotate子程序的rename-chrs模块
• 调用bgzip进行压缩

bcftools annotate --rename-chrs chr_name_change.txt cohart.filter.vcf.gz | bgzip -c > cohart.filter.nochr.vcf.gz