昨天看了gatk的官網(wǎng)，從2018年發(fā)布正式版的4.0.0開始，到現(xiàn)在已經(jīng)更新到4.1.8幽七，在速度和準(zhǔn)確度上都有了大幅的提升。gatk4除了整合picard軟件之外溅呢，在使用上與gatk3基本相同澡屡，只不過是在命令運行、功能劃分及運行速度上進行了調(diào)整咐旧。gatk4軟件的安裝及使用驶鹉，網(wǎng)上有很多帖子，可以自行查閱铣墨。今天主要給大家介紹gatk4與gatk3明顯不同的地方梁厉，這些使得gatk4的應(yīng)用更加方便。

1. 輸入文件的變化

變異數(shù)據(jù)是從比對結(jié)果中獲得的踏兜，gatk之前就已經(jīng)支持cram的輸入词顾，但是數(shù)據(jù)去重部分的picard軟件并不支持，運行時會報錯碱妆。至于cram文件肉盹，也是一種比對結(jié)果文件，只不過文件體積更小疹尾，即對于存儲受限的分析人員里說是一個值得考慮的方向上忍，具體的信息可查看公眾號文章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之后對其進行了優(yōu)化(猜測繁成，具體是不是整合前picrad進行了更新不清楚吓笙，有興趣的可以去測測看)，現(xiàn)在是支持cram的輸入和去重的巾腕。

需要注意的是面睛，軟件運行結(jié)果是bam文件，而不是cram文件尊搬。

2. gatk4變異檢測模塊調(diào)整

使用gatk4進行變異檢測時叁鉴，主要包括以下5個模塊

一般的變異使用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主要應(yīng)用在群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包的運行模式，命令行相對簡單鸟召。當(dāng)然胆绊，對于非常熟悉gatk3運行并且不喜歡gatk4的分析人員來說，也可以直接使用安裝包下面的jar包來運行程序欧募，使用方法與GATK3基本相同压状。

3. 群call策略

群call的最終目的是獲得所有樣本的分型結(jié)果，最終會使用GenotypeGVCFs模塊轉(zhuǎn)化為vcf，但是由于輸入文件個數(shù)的限制种冬，需要先使用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中新出現(xiàn)的娱两，GenomicsDBImport莺匠，該方法先對gvcf文件進行處理，輸出一個類似于數(shù)據(jù)庫的目錄十兢，后面可以直接使用該數(shù)據(jù)目錄作為輸入文件進行變異檢測趣竣。

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參數(shù)指定數(shù)據(jù)庫的輸出路徑，該路徑在運行之前不能被創(chuàng)建旱物，否則會報錯遥缕。合并完成之后可以進行位點變異的檢測。

gatk GenotypeGVCFs \
    -R data/ref/ref.fasta \
    -V gendb://my_database \
    -G StandardAnnotation -newQual \
    -O test_output.vcf

群call最大的優(yōu)勢在于我們可以添加樣本后重新分析宵呛，獲取群體的變異位點单匣。如果使用GenomicsDBImport進行分析，若要添加新樣本的變異數(shù)據(jù)烤蜕，只要將新樣本的gvcf信息添加到已有的數(shù)據(jù)庫中即可封孙。

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

數(shù)據(jù)庫中含有的文件是二進制的文件，若想要多的數(shù)據(jù)庫中存在的染色體信息讽营，或者想要提取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和準(zhǔn)確度的優(yōu)勢糖荒，特別是大樣本量的差異分析杉辙，gatk的優(yōu)勢還是很明顯的。

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

#如有侵權(quán)捶朵，請告知刪除#
#如有錯誤蜘矢，歡迎指正#