gatk4使用总结
昨天看了gatk的官网,从2018年发布正式版的4.0.0开始,到现在已经更新到4.1.8,在速度和准确度上都有了大幅的提升。gatk4除了整合picard软件之外,在使用上与gatk3基本相同,只不过是在命令运行、功能划分及运行速度上进行了调整。gatk4软件的安装及使用,网上有很多帖子,可以自行查阅。今天主要给大家介绍gatk4与gatk3明显不同的地方,这些使得gatk4的应用更加方便。
1. 输入文件的变化
变异数据是从比对结果中获得的,gatk之前就已经支持cram的输入,但是数据去重部分的picard软件并不支持,运行时会报错。至于cram文件,也是一种比对结果文件,只不过文件体积更小,即对于存储受限的分析人员里说是一个值得考虑的方向,具体的信息可查看公众号文章cram输出及使用。
picard原装软件使用cram运行时报错信息。
Exception in thread "main" java.lang.IllegalStateException: A valid CRAM reference was not supplied and one cannot be acquired via the property settings reference_fasta or use_cram_ref_download
at htsjdk.samtools.cram.ref.ReferenceSource.getDefaultCRAMReferenceSource(ReferenceSource.java:105)
at htsjdk.samtools.SamReaderFactory$SamReaderFactoryImpl.open(SamReaderFactory.java:406)
at picard.sam.markduplicates.util.AbstractMarkDuplicatesCommandLineProgram.openInputs(AbstractMarkDuplicatesCommandLineProgram.java:220)
at picard.sam.markduplicates.MarkDuplicates.buildSortedReadEndLists(MarkDuplicates.java:535)
at picard.sam.markduplicates.MarkDuplicates.doWork(MarkDuplicates.java:264)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:295)
at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:103)
at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:113)
不能给定基因组,就意味着不能使用cram文件。而gatk4在整合了picard之后对其进行了优化(猜测,具体是不是整合前picrad进行了更新不清楚,有兴趣的可以去测测看),现在是支持cram的输入和去重的。
需要注意的是,软件运行结果是bam文件,而不是cram文件。
2. gatk4变异检测模块调整
使用gatk4进行变异检测时,主要包括以下5个模块
Name | Summary |
---|---|
CombineGVCFs | Merges one or more HaplotypeCaller GVCF files into a single GVCF with appropriate annotations |
GenomicsDBImport | Import VCFs to GenomicsDB |
GenotypeGVCFs | Perform joint genotyping on one or more samples pre-called with HaplotypeCaller |
HaplotypeCaller | Call germline SNPs and indels via local re-assembly of haplotypes |
Mutect2 | Call somatic SNVs and indels via local assembly of haplotypes |
一般的变异使用HaplotypeCaller进行检测,输出格式可以是vcf,也可以是gvcf
# GVCF输出
gatk --java-options "-Xmx4g" HaplotypeCaller \
-R Homo_sapiens_assembly38.fasta \
-I input.bam \
-O output.g.vcf.gz \
-ERC GVCF
# VCF输出
gatk --java-options "-Xmx4g" HaplotypeCaller \
-R Homo_sapiens_assembly38.fasta \
-I input.bam \
-O output.vcf.gz \
-bamout bamout.bam
vcf适合单样本变异检测,GVCF主要应用在群call变异检测上。对于gvcf文件的处理,可以使用GenotypeGVCFs模块输出变异位点的信息。
gatk --java-options "-Xmx4g" GenotypeGVCFs \
-R Homo_sapiens_assembly38.fasta \
-V input.g.vcf.gz \
-O output.vcf.gz
该方法在使用上有限制,只能输入一个gvcf文件,对于多个样本群call的方式,需要使用另外两个软件对gvcf文件进行处理。
另外还有一点,gatk4在命令上有很大的改进,弃用了直接使用java包的运行模式,命令行相对简单。当然,对于非常熟悉gatk3运行并且不喜欢gatk4的分析人员来说,也可以直接使用安装包下面的jar包来运行程序,使用方法与GATK3基本相同。
3. 群call策略
群call的最终目的是获得所有样本的分型结果,最终会使用GenotypeGVCFs模块转化为vcf,但是由于输入文件个数的限制,需要先使用GenomicsDBImport和GenotypeGVCFs2个模块处理。
CombineGVCFs模块是用来合并多个gvcf文件。
gatk CombineGVCFs \
-R reference.fasta \
--variant sample1.g.vcf.gz \
--variant sample2.g.vcf.gz \
-O cohort.g.vcf.gz
另外一种可以合并gvcf方式是gatk4中新出现的,GenomicsDBImport,该方法先对gvcf文件进行处理,输出一个类似于数据库的目录,后面可以直接使用该数据目录作为输入文件进行变异检测。
gatk --java-options "-Xmx4g -Xms4g" GenomicsDBImport \
-V data/gvcfs/mother.g.vcf.gz \
-V data/gvcfs/father.g.vcf.gz \
-V data/gvcfs/son.g.vcf.gz \
--genomicsdb-workspace-path my_database \
--tmp-dir=/path/to/large/tmp \
-L 20
--genomicsdb-workspace-path参数指定数据库的输出路径,该路径在运行之前不能被创建,否则会报错。合并完成之后可以进行位点变异的检测。
gatk GenotypeGVCFs \
-R data/ref/ref.fasta \
-V gendb://my_database \
-G StandardAnnotation -newQual \
-O test_output.vcf
群call最大的优势在于我们可以添加样本后重新分析,获取群体的变异位点。如果使用GenomicsDBImport进行分析,若要添加新样本的变异数据,只要将新样本的gvcf信息添加到已有的数据库中即可。
gatk GenomicsDBImport \
-V data/gvcfs/mother.g.vcf \
-V data/gvcfs/father.g.vcf \
-V data/gvcfs/son.g.vcf \
--genomicsdb-update-workspace-path existing_database
数据库中含有的文件是二进制的文件,若想要多的数据库中存在的染色体信息,或者想要提取gvcf文件,可以使用如下命令
# 提取gvcf文件
gatk SelectVariants \
-R data/ref/ref.fasta \
-V gendb://my_database \
-O combined.g.vcf
# 提取interval信息
gatk GenomicsDBImport \
--genomicsdb-update-workspace-path existing_database
--output-interval-list-to-file /path/to/output/file
总体来说,gatk4更新之后还是很好用的,相比于gatk3,只要选择合适的分析流程和策略,时间上也会大幅缩短。不过相比于samtools这些软件,gatk分析时间还是最大的问题,不过如果考虑到群call和准确度的优势,特别是大样本量的差异分析,gatk的优势还是很明显的。
4. 参考资料
[1] https://gatk.broadinstitute.org/hc/en-us/articles/360036712151-HaplotypeCaller
[2] https://gatk.broadinstitute.org/hc/en-us/articles/360042914991-GenotypeGVCFs
[3] https://gatk.broadinstitute.org/hc/en-us/articles/360036713211-CombineGVCFs
[4] https://gatk.broadinstitute.org/hc/en-us/articles/360036712071-GenomicsDBImport
[5] https://gatk.broadinstitute.org/hc/en-us/articles/360035891051-genomicsdb
#如有侵权,请告知删除#
#如有错误,欢迎指正#