Zenn
🛠️

vcfファイルのCHROM列の染色体番号を0埋めしてSortする

2024/04/01に公開
#
bioinformatics
#
vcf
#
genome
tech

はじめに

vcfファイルの染色体名を0埋め(ゼロパディング:zero padding)していないと後続の解析で面倒なことが頻発する。

ここでの0埋めしていない状態というのは、chr1やscaffold1といった状態のことを指し、これをchr01やscaffold01という用に数字の前に0入れて桁数を合わせることで、文字列ソートする際、01→10→20といった期待する並び替えができるようになるようにする事を目的とする。

CHROM列のリネームとSort

リネームには、bcftoolsのannotate --rename-chrsとshellスクリプトを、SortにはbcftoolsのSortを使用する。

Step1 : bcftoolsのインストール

condaでインストールする場合は以下のように、仮想環境を作成してインストールする。既に使用している仮想環境がある場合にはそこにインストールする。

mamba install -c bioconda -c conda-forge bcftools

最新バージョンをソースコードからインストールする場合は以下の通り。gitコマンドを使用するので、なければインストールしておく。

https://samtools.github.io/bcftools/howtos/install.html

# htslibのリポジトリのクローン
git clone --recurse-submodules https://github.com/samtools/htslib.git
# bcftoolsのリポジトリのクローン
git clone https://github.com/samtools/bcftools.git

# 移動とコンパイル
cd bcftools
make -j 32

# 確認
./bcftools

# 今回は使わないがプラグインのパスの追加 (永続的に反映する場合は、~/.bashrcに記載する)
export BCFTOOLS_PLUGINS=$PWD/plugins

Step2 : CHROM列の情報の抽出

bcftools queryでCHROM列の情報を抽出する。vcfファイルの名前をinput.vcfと仮定すると、以下のコマンドでCHROM列の情報が抽出できる。

すべての行の情報が表示されてしまうので、uniq関数で重複を削除する。

input_vcf="input.vcf"

bcftools query -f "%CHROM\n" "$input_vcf" | uniq

Step3 : 受け取ったCHROM列の情報を0埋めして、bcftools annotate --rename-chrsに渡せる形式のファイルに変換

先程のbcftoolsの出力を使って0埋めしつつ、1列目には抽出元の情報、2列目には0埋めしたCHROM列に割り当てたいChromosome名を記載した対応表を作成する。bcftools queryの出力はパイプで渡す。

# 入力用vcf
input_vcf="input.vcf"
# 変換ルールを保存する一時ファイル
rename_rules="rename_rules.txt"

bcftools query -f "%CHROM\n" "$input_vcf" | uniq | while IFS= read -r line; do
  # プレフィックス（chrまたはscaffold）と数字を分離
  prefix=$(echo "$line" | grep -oE '^(chr|scaffold)')
  number=$(echo "$line" | grep -oE '[0-9]+$')

  # 数字を2桁のゼロ埋め形式に変換
  formatted_number=$(printf "%02d" "$number")

  # 変換ルールをファイルに書き込む
  echo "${line} ${prefix}${formatted_number}" >> "$rename_rules"
done

これで、rename_rules.txtファイルに対応表が出力される(下記はsortがうまくできていないvcfから生成した例)。

cat rename_rules.txt | grep -i "chr" | head -n 15
# chr10 chr10
# chr11 chr11
# chr12 chr12
# chr13 chr13
# chr14 chr14
# chr15 chr15
# chr16 chr16
# chr17 chr17
# chr18 chr18
# chr19 chr19
# chr1 chr01
# chr20 chr20
# chr21 chr21
# chr22 chr22
# chr23 chr23

cat rename_rules.txt | grep -i "scaffold" | head -n 15
# scaffold10 scaffold10
# scaffold11 scaffold11
# scaffold12 scaffold12
# scaffold13 scaffold13
# scaffold14 scaffold14
# scaffold15 scaffold15
# scaffold16 scaffold16
# scaffold17 scaffold17
# scaffold18 scaffold18
# scaffold19 scaffold19
# scaffold1 scaffold01
# scaffold20 scaffold20
# scaffold21 scaffold21
# scaffold22 scaffold22
# scaffold23 scaffold23

Step4 : 対応表を使ってvcfファイルのCHROM列をrenameしてsort

--rename-chrsオプションの後ろに、対応表、rename対象のvcfと続ける。また、sortコマンドをパイプで繋ぎ、CHROM列の情報で昇順の並び替えを行う。出力ファイルは-oオプションの後ろにoutput用のvcfファイルを記載する。

input_vcf="input.vcf"
# 出力ファイル名
output_vcf="renamed_${input_vcf}"

# 染色体名を置換
bcftools annotate --rename-chrs "$rename_rules" "$input_vcf" | bcftools sort -Oz -o "$output_vcf"

確認

output_vcf="renamed_merged.vcf.gz"
# 変換後のCHROM列を表示
bcftools query -f "%CHROM\n" ${output_vcf} | uniq

これで0埋めされたChromosome名を使ってsortされたrenamed_${input_vcf}が得られた。

後は作成したvcfファイルを使用するとストレス少なく後続の統計解析などが実施できる。

ただ、はじめからreference.fastaの配列情報は0埋めしておくとこういう事態は発生しないので次回からは注意する。

自動的に判断するような条件分岐を記載してもいいかもしれない。

Information

Discussion