如何快速简化vcf信息?

发布时间 2023-07-26 11:41:34作者: 米源MY

需求描述

vcf是标准的基因型格式文件,其中包含的信息可多可少。主要在于INFO可无限扩展特征,以及每个样本的FORMAT信息,会大大增加vcf文件的大小。一般来说,GATK等软件得到的基因型都会有这些信息,初始变异我们最好保留它们,因为这是过滤位点/样本的依据。但是当我们确定了最终用于后续下游分析的位点时,除了每个样本的基因型,其他信息我们往往不会再关注。这时,仍然在一个巨大文件中保留这些冗余信息,显然是没有必要的。

如何快速简化vcf,只保留位点及基因型等基本信息?

可能存在错误的做法

vcf和plink文件之间进行相互转化:vcf——>plink——>vcf,这样一番操作,就将vcf简化了。

plink --vcf snp.pass.vcf.gz --make-bed --out snp --allow-extra-chr --double-id
plink --bfile snp --recode vcf-iid --out snp --allow-extra-chr

但这样操作存在一些潜在的问题:

一是染色体名称和位点的变化。chr/Chr等前缀会直接去掉(导致后续要做注释时,需要对注释gff/gtf做一些额外的操作),而且非数字染色体上位点可能直接被删掉,--allow-extra-chr可能也起不到作用;

二是样本名称的潜在变化。plink默认用下划线对样本名进行分隔,分隔的两个字段分别作为ped文件中的family id和sample id, 如果vcf中的样本名含有多个下划线,无法正确进行划分,软件会报错,此时可以修改--id-delim参数,该参数设定了分隔符,默认是下划线,可以设置成其他字符,以达到正确区分的目的。为保险起见,这点可以在转化格式时加个--double-id参数。

三是REF和ALT可能调换。这是最严重的,格式转回来后,有些位点的REF和ALT碱基直接调换了,大部分则是正常的。这个坑我踩过,目前仍然百思不得其解。

以上潜在的问题,可能导致文件格式后的位点基因型与原来不一致,而你可能察觉不到。

更靠谱的做法

我更倾向于推荐用bcftools来简化vcf文件。

bcftools annotate --remove QUAL,FILTER,INFO,^FORMAT/GT snp.pass.vcf.gz -Oz -o snp.simply.vcf.gz

解释:
--remove QUAL,FILTER,INFO,^FORMAT/GT: 删除 QUAL,FILTER,INFO字段信息,以及FORMAT 字段中除 GT(基因型)开头的字段外的所有字段。它实际上只保留了基因型信息,而删除了其他所有信息。

原文件:
image.png

简化后文件:

image.png

当然还有其他的方法,比如自己写脚本,这里不讨论。

Ref:
http://samtools.github.io/bcftools/howtos/annotate.html