实操 | 合并VCF文件的几种方法及注意事项
背 景
在基因组分析领域的很多不同场景中,需要合并VCF文件。
VCF文件。简单来说,就是记录样本基因型的文件。但多数VCF文件不只记录了基因型,也包含有关该基因型的来源的细节。
其它文件。VCF文件的上游是BAM文件,主要记录Reads与参考基因组的比对信息;更上游的,就是FASTQ测序数据,以及物种的参考基因组。
不同类型的VCF文件。VCF文件有单样本的、多样本的;也有普通VCF文件 (只记录变异,未测到的、野生型的都不记录),及GVCF文件 (野生型的、变异的都记录,未测到的不记录)。
一个典型的GVCF文件
因此,本篇文章的使用者,需要首先了解GVCF文件与普通VCF文件的不同。因为二者对应的生信处理方法也非常不同。但具体有哪些不同,这里不再继续讲,可自行在网络资源上按照相应的关键词检索。
合并VCF文件需要注意的问题。VCF文件有时是多个样本、多个个体、或多个病例的合并;有时是不同染色体区域的VCF合并。上述每个场景都涉及不同的软件、程序,甚至算法,需要非常小心、谨慎地操作。
合并:不同样本的GVCF文件
GATK的CombineGVCFs + GenotypeGVCFs
上面2个程序是一套组合,不可拆分,不可单独使用。
那为什么GATK开发者将二者分开呢,推测有两个原因:① 二者分别有些特定的参数;② 第1个程序非常耗时,第2个程序相对较快、但算法复杂。这个问题对使用者也无关紧要。
# 合并队列中每个样本的每一个变异 (GVCF文件)
gatk CombineGVCFs \-R $ref \$(for i in `tail -n +2 metadata.txt | cut -f 1 `; do echo "--variant ${i}.hg38.raw.g.vcf " ;done) \-O cohort.g.vcf.gz
# 获取具体基因型,完成变异Calling
gatk GenotypeGVCFs \-R $ref \--dbsnp ${dbSNP} \--variant cohort.g.vcf.gz \-O Genotype.cohort.dbSNP.g.vcf
当样本很多、数据量也大时,CombineGVCF程序很消耗内存,并且一旦中断(文件不全)就得重新来。其"-L"参数 (如"-L chr1:1-10000") 也不推荐使用 (否则GenotypeGVCFs步骤可能报错),但"-L chr1"没问题。
解决方法:① 限制内存:用"--java-option Xmx20g"等;② 分染色体,用"-L chrX"等。③ 组合使用①和②。
需要了解的是,GenomicsDB可以替代:CombineGVCFs + GenotypeGVCFs,将多个样本GVCF处理生成一起的工作空间。两种方案各有各的优缺点。
根据GATK官网的描述,GenomicsDB更适用于几百个样本以上的情形。
合并:不同染色体区域的VCF文件
cat chr.list.25
# chr1
# chr2
# chr3
# ...
# chr22
# chrX
# chrY
# chrM
gatk MergeVcfs \$(for i in `cat chr.list.25 | cut -f 1 `; do echo "-I Genotype.cohort.${i}.dbSNP.g.vcf " ;done) \-O Genotype.cohort.dbSNP.g.vcf
MergeVcfs与CombineGVCFs不同。前者用于单纯地合并:样本相同、位点独立的VCF文件。如:同一个(或一组)样本的不同染色体的结果。
不像CombineGVCFs,MergeVcfs不做"gVCF block"的计算。此外MergeVcfs会检测两个VCF文件里的样本名是否相互"match"。
如果只查看MergeVcfs程序的介绍,根本看不出来它的用法的特点 (例如:对GVCF文件的合并可能无效),那么必然容易踩坑:
MergeVcfs (Picard) - Combines multiple variant files into a single variant file.
事实上,MergeVcfs及其等价的程序 (GatherVcfs) 不可用于合并不同样本的GVCF文件。但用来合并不同基因组区域的文件非常方便。
此场景除了MergeVcfs、GatherVcfs外的其它程序:
vcf-concat sample1.chr1.vcf sample.chr2.vcf ... > sample1.chrAll.vcf
bcftools concat sample1.chr1.vcf sample.chr2.vcf ... -o sample1.chrAll.vcf
其中,通过conda安装的vcftools,可能不带vcf-concat等程序。从这一点看,bcftools更方便。
1个经验是:既然有GATK的MergeVcfs可用,那就尽量不用vcftools或bcftools,否则可能踩到另一个坑:不同程序对VCF文件的索引格式的要求不同、VCF的"FORMAT"列等也可能改变。
合并:不同样本的普通VCF文件
普通VCF文件只记录变异,即:① 无0/0基因型 (测序测到了、但未变异,即"Wild type") ;② 无"./."基因型 (即"缺失",测序未测到,即"No call") 。
对不同样本的⽂件合并,共有位点会合并、统计;非共有位点若在某1个样本中无变异,则会⾃动记为缺失 ("./.") 。
1个典型的普通VCF文件 (只查看了第9、10列)
vcftools和bcftools在使用之前一般都需要对VCF文件:压缩、索引 (略)。
vcf-merge # 略 (对于连软件安装都麻烦的程序)
bcftools sample1.vcf sample2.vcf ... -o sample.all.vcf
vcf-merge重新计算了AC、AN等指标的值
合并分型质量 ("QUAL"列) 时,vcf-merge取了平均取值,bcftools取了最⼤值, (下图的)gatk CombineVariants (不是CombineGVCFs,也不是MergeVcfs/GatherVcfs)取了最⼩值 (gatk4)
图片来源:https://wenku.baidu.com/view/a0ecad5602f69e3143323968011ca300a7c3f643.html
gatk CombineVariants (GATK4已无此程序)
# 压缩、索引单个VCF文件
ls sorted.*.vcf | while read id;dobcftools view $id -Oz -o $id.gzbcftools index $id.gz
done
# 合并
bcftools merge --threads 8 -m id -O z sorted.*.vcf.gz \-o bcftools.merged.103samps.vcf.gz# -m, --merge (关于多等位基因)Allow multiallelic records for <snps|indels|both|all|none|id>, see man page for details [both]# -O, --output-type <b|u|z|v> 'b' compressed BCF; 'u' uncompressed BCF; 'z' compressed VCF; 'v' uncompressed VCF [v]# gVCF参数(可能不适用于gVCF文件):# -g, --gvcf <-|ref.fa> merge gVCF blocks, INFO/END tag is expected. Implies -i QS:sum,MinDP:min,I16:sum,IDV:max,IMF:max
# 测试某个位点# zcat bcftools.merged.103samps.vcf.gz | grep '1 12626 '# 有返回结果。但无DP等信息
# 索引合并后的文件
bcftools index bcftools.merged.103samps.vcf.gz &
bcftools merge虽然有"--gvcf"参数,但根据之前的测试,可能不适用于对gVCF文件的合并。
总 结
总之,① 合并VCF文件要区分其文件类型,如:是否为gVCF文件,是否为基因组的不同区域,其内部的样本名称等;② 考虑到整个流程的兼容性和流畅性,建议当GATK有相应的工具时,优先使用之;③ 其它场景可依次考虑:bcftools、vcftools。
往期精品(点击图片直达文字对应教程)
机器学习
后台回复“生信宝典福利第一波”或点击阅读原文获取教程合集
实操 | 合并VCF文件的几种方法及注意事项相关推荐
- 合并BIN文件的两种方法
合并BIN文件的两种方法 在单片机的开发过程中,经常需要将两个单独的BIN文件合并成一个文件,方便烧写和生产.下面结合STM32的IAP Bootloader Code和Application Cod ...
- 合并HEX文件的一种方法
<合并BIN文件的两种方法>介绍了如何合并BIN文件,在这个基础之上再配合hex2bin.exe和bin2hex.exe这两个小工具就可以很方便的将两个HEX文件合并成一个了.当然,最终目 ...
- 使用ffmpeg合并视频文件的三种方法
ffmpeg合并视频的方法有三种.国内大多数仅介绍了其中之一.于是觉得有必要翻译一下.其实在ffmpeg的 FAQ文档中有比较详细的说明. 使用concat协议进行视频文件的合并 这种方式的适用场景是 ...
- Java合并PDF文件的几种方法
最近需要做一个把多个pdf报告合并成一个以方便预览的需求,下面总结一下自己用的方法和遇到的一些问题, 第一种方法: 此方法引用了itextpdf.jar包: private static void m ...
- android选择多个文件_一分钟合并多个Excel、PDF文件,3种方法任你选择,好用到没朋友...
一分钟合并多个Excel.PDF文件,3种方法任你选择,好用到没朋友 前情提要: Excel.PDF多个文件怎样合并成一个文件?需求场景: PDF文件合并 当一份完成的PDF资料分为很多份的时候,我们 ...
- 多个PDF文件如何合并成一个?两种方法轻松get
在日常学习生活中,如果你需要将多个文档整合为一个完整的文件,比如说多篇文章.多张图片.多个表格等等,这时候就需要将这些文档合并成一个PDF文件.如何将多个PDF文件如何合并成一个?两种方法轻松帮你搞定 ...
- android写入文件方法,Android 追加写入文件的三种方法
一.使用FileOutputStream 使用FileOutputStream,在构造FileOutputStream时,把第二个参数设为true public static void method1 ...
- 用python下载文件的若干种方法汇总
压缩文件可以直接放到下载器里面下载的 you-get 连接 下载任意文件 重点 用python下载文件的若干种方法汇总 写文章 用python下载文件的若干种方法汇总 zhangqibot发表于Met ...
- Htaccess文件是什么以及Windows下自由创建.htaccess文件的N种方法
.htaccess是什么 概述来说,htaccess文件是Apache服务器中的一个配置文件,它负责相关目录下的网页配置. 通过htaccess文件,可以帮我们实现:网页301重定向.自定义404错误 ...
最新文章
- Guice系列之用户指南(五)
- DataGrid中添加背景
- glut编译问题 (程序无法运行)
- AWS 之于 K8s,如同 Windows 之于 Linux!
- erlang的简单模拟半包的产生
- opencv基础:罗德里格斯旋转公式(Rodrigues' rotation formula)推导 rodrigues()函数原理
- unity 1 学习 物体旋转和通过脚本调用单击事件函数,find函数找物体的方法
- android无法实例化服务器,android – 无法实例化类型PagerAdapter
- java消息头_java中怎么进行头消息校验
- Hybrid App实现原理
- Nordic nRF52840 入门学习
- printf()输出格式大全(附 - 示例代码)
- VMware安装win7操作系统
- kaggle 电商数据分析
- ADO.NET 概述
- Modern PHP
- 中国乡镇企业会计杂志中国乡镇企业会计杂志社中国乡镇企业会计编辑部2022年第12期目录
- 按字寻址和按字节寻址
- hdu2639(01背包变形-第k大背包)
- 微信的原创保护机制到底是如何实现的?