轉換 VCF 檔染色體命名慣例

發表於

VCF 檔中的 ##contig=<ID=*,lenth=*> 記錄了參考基因體所有染色體的名稱和長度等資訊，而 CHROM 欄位則標示變異所在的染色體。目前，染色體名稱有 UCSC （chr + 染色體編號，例如：chr1和 chrX）和 Ensembl （沒有前綴，僅有染色體編碼，例如：1 和 X）兩種慣例。執行分析時， 如果輸入的 VCF 檔和對應的參考資料採用相異命名慣例，往往會導致分析結果失真。

舉例來說，GIAB 提供的 VCF 檔是評估變異分析常用的資料。這些 VCF 檔採用 Ensembl 慣例，如果拿 UCSC 慣例的 VCF 檔和它們比較，程式會以為沒有找到這些變異。因此，執行分析時，常常要調整染色體名稱，確保分析結果正確。

本文簡介使用 linux sed 和 bcftools 轉換 VCF 檔染色體命名慣例的方式。

範例資料：分別採用 UCSC 和 Ensembl 慣例的 VCF 檔

ucsc.vcf，UCSC 慣例會在染色體名稱加上 chr 前綴，粒線體基因則命名為 chrM。

1
2
3
4
5
6
##fileformat=VCFv4.2
##contig=<ID=chr1>
##contig=<ID=chrM>
#CHROM    POS    ID    REF    ALT    QUAL    FILTER    INFO    FORMAT    test
chr1    100    .    T    C    100    PASS    .    .    .
chrM    100    .    T    C    100    PASS    .    .    .

ensembl.vcf，Ensembl 慣例只標記染色體名稱，粒線體基因則命名為 MT。

1
2
3
4
5
6
##fileformat=VCFv4.2
##contig=<ID=1>
##contig=<ID=MT>
#CHROM    POS    ID    REF    ALT    QUAL    FILTER    INFO    FORMAT    test
1    100    .    T    C    100    PASS    .    .    .
MT    100    .    T    C    100    PASS    .    .    .

除了 CHROM，轉換慣例時也要一併異動 ##contig=<ID=* 的格式。許多軟體（例如 GATK）會檢查兩者是否相符，如果 ##contig=<ID=*> 含有 CHROM 不存在的染色體（反之亦然）就會報錯。

使用 sed 轉換染色體命名慣例

我們可以使用 sed 添加或是刪除 chr 來轉換命名慣例。使用時要注意替換字串的順序和模式。舉例來說，從 UCSC 轉換為 Ensembl 慣例時，如果先把 VCF 檔內所有的 chr 移除，那麼調整粒線體基因名稱時要把 M 改成 MT。然而，M 可能配對到太多字符，所以需要新增其它判斷的根據。

UCSC -> Ensembl

1
2
3
4
5
6
cat ucsc.vcf |\
sed '
  s/##contig=<ID=chrM/##contig=MT/;
  s/^chrM/MT/;
  s/^chr//;
  s/##contig=<ID=chr/##contig=<ID=/'

Ensembl -> UCSC

1
2
3
4
5
cat ensembl.vcf |\
sed 's/^\([0-9XYM]\)/chr\1/;
  s/^chrMT/chrM/;
  s/^##contig=<ID=/##contig=<ID=chr/;
  s/^##contig=<ID=chrMT/##contig=<ID=chrM/'

使用 bcftools 轉換染色體命名慣例

bcftools annotate 會依據用戶提供的命名格式辭典，同時更新 ##contig=<ID=*> 和 CHROM 的紀錄。例如：

ucsc_to_ensembl.tsv

1
2
chr1    1
chrM    MT

以及

ensembl_to_ucsc.tsv

1
2
1    chr1
MT    chrM

UCSC -> Ensembl

1
bcftools annotate --rename-chrs ucsc_to_ensembl.tsv ucsc.vcf

Ensembl -> UCSC

1
bcftools annotate --rename-chrs ensembl_to_ucsc.tsv ensembl.vcf

如果你拿到的是壓縮過的 VCF 檔

前述案例都是處理未壓縮的 VCF 檔，但實務上為了節省空間和時間都會事先壓縮 VCF 檔為 vcf.gz，並且用索引檔加速後續的資料處理。

調整 vcf.gz 的染色體命名慣例的步驟為解壓縮、調整格式、壓縮和索引，程式碼範例如下，記得把 <commands> 部分替換成對應的 sed 指令。

1
2
3
4
5
6
# 解壓縮 vcf.gz 以執行資料處理，隨後再壓縮
zcat <input.vcf.gz> |\
<commands> |\
bcftools -o <output.vcf.gz> -O z
# 為更新後的 vcf.gz 建立索引檔
bcftools index -t <output.vcf.gz>

bcftools annotate 可以直接處理並輸出壓縮的 VCF 檔，但最後還是要自行建立索引檔：

1
2
3
4
5
6
# UCSC to Ensembl
bcftools annotate --rename-chrs ucsc_to_ensembl.tsv ucsc.vcf.gz -O z -o ensembl.vcf.gz
bcftools index -t ensembl.vcf.gz # 會產出 ensembl.vcf.gz.tbi
# Ensembl to UCSC
bcftools annotate --rename-chrs ensembl_to_ucsc.tsv ensembl.vcf.gz -O z -o ucsc.vcf.gz
bcftools index -t ucsc.vcf.gz # 會產出 ucsc.vcf.gz.tbi

延伸資訊