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

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