GATK软件操作生物信息DNA-seq学习

gatk4使用总结

2020-09-19  本文已影响0人  xiaoji_hb

昨天看了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

#如有侵权,请告知删除#
#如有错误,欢迎指正#

上一篇下一篇

猜你喜欢

热点阅读