转换参考基因组坐标

Posted on Views: Waline: Word count in article: 1.9k Reading time ≈ 3 mins.

我们有时候在处理多个vcf文件时，他们的参考基因组可能各不相同，这时候就需要把vcf 文件中的变异位置从一个参考基因组映射到另一个参考基因组。

1.picard方法

先picard 建立基因组索引.dict文件，再vcf坐标转换

1
java -jar /home/liuxiao/soft/picard/picard.jar CreateSequenceDictionary R= sequences.fa O=sequences.fa.dict
1
2
3
4
5
6
7
#使用 Picard 的 LiftoverVcf工具
java -jar /home/lx/soft/picard.jar LiftoverVcf \
I=/home/lx/wsl/caid323.vcf.gz  \
O=323.vcf \
CHAIN=/home/lx/wsl/oviAri4Tooarv1.over.chain  \
REJECT=unmap_variants.vcf  \
R=/home/lx/wsl/oarchange/oarv1.fasta

相关的.fa文件，.dict文件，.vcf文件, .chain文件均尽量软连接在同一文件夹

1
ln -s $truepath  #软连接

2、liftOver工具。

步骤 1：准备 VCF 文件

确保你的 VCF 文件格式正确，可以用以下命令查看和处理 VCF 文件：

1
vcftools --vcf input.vcf --recode --stdout > prepared_input.vcf

步骤 2：提取 VCF 文件中的变异位置

需要将 VCF 文件转换为 BED 格式，因为 liftOver 工具处理 BED 格式文件。

1
vcf2bed < prepared_input.vcf > input.bed

可以使用 bedops 工具中的 vcf2bed 命令。如果没有安装 bedops，可以使用以下 Python 脚本进行转换：

1
2
3
4
5
6
7
8
9
10
11
import vcf
vcf_reader = vcf.Reader(open('prepared_input.vcf', 'r'))
with open('input.bed', 'w') as bed_file:
    for record in vcf_reader:
        chrom = record.CHROM
        start = record.POS - 1
        end = record.POS
        name = record.ID
        score = '.'
        strand = '.'
        bed_file.write(f"{chrom}\t{start}\t{end}\t{name}\t{score}\t{strand}\n")

步骤 3：运行 liftOver

使用 liftOver 工具将坐标转换到新的参考基因组版本。

1
./liftOver input.bed oldToNew.chain output.bed unlifted.bed
  • oldToNew.chain 是你下载的链文件，它描述了两个基因组版本之间的坐标转换关系。
  • output.bed 是成功转换的坐标。
  • unlifted.bed 是未能转换的坐标。

步骤 4：将 BED 文件转换回 VCF 文件

使用以下脚本将转换后的 BED 文件转换回 VCF 文件：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import vcf

vcf_reader = vcf.Reader(open('prepared_input.vcf', 'r'))
records = list(vcf_reader)

with open('output.bed', 'r') as bed_file, open('output.vcf', 'w') as vcf_file:
    vcf_writer = vcf.Writer(vcf_file, vcf_reader)
    for line in bed_file:
        fields = line.strip().split('\t')
        chrom = fields[0]
        pos = int(fields[2])
        for record in records:
            if record.CHROM == chrom and record.POS == pos:
                record.POS = pos
                vcf_writer.write_record(record)