最近看了GATK的網(wǎng)站滞谢，以前經(jīng)常用GATK call somatic mutation串稀， 但是用的版本太老了，另外服務(wù)器最近換了盤狮杨，很多代碼路徑都變了母截，網(wǎng)上關(guān)于GATK4的教程比較少，不想做伸手黨橄教，于是嘗試自己寫代碼試一試清寇。

1 Getting start with GATK

GATK是Genome AnalysisToolkit的縮寫，它搜集了一系列分析高通量測(cè)序數(shù)據(jù)的工具护蝶，主要目的是在于發(fā)現(xiàn)變異华烟，這個(gè)工具可以單獨(dú)的用，但是也可以和其他的workflow合在一起持灰。

contents：

1. 安裝GATK

GATK 工具具有相當(dāng)簡單的軟件需求： 一個(gè)Unix-style OS 以及java 1.8 版本以上盔夜，其他的工具需要額外的R和python的依賴。 可以在這個(gè)網(wǎng)址下載： https://github.com/broadinstitute/gatk/releases堤魁，
最新版本的GATK可以按照如下命令下載：

   wget https://github.com/broadinstitute/gatk/releases/download/4.1.9.0/gatk-4.1.9.0.zip

如果需要其他版本喂链，可以自行在在這個(gè)網(wǎng)址下載。 一旦你下載了解壓這個(gè)文件夾妥泉，里面一般會(huì)有這四個(gè)文件：

gatk
gatk-package-[version]-local.jar
gatk-package-[version]-spark.jar
README.md

這個(gè)時(shí)候你可能會(huì)問椭微，為什么會(huì)有兩個(gè)jar， 就像名字里一樣盲链，gatk-package-[version]-spark.jar 是在Spark cluster上專門跑Spark tools 設(shè)計(jì)的蝇率。
而 gatk-package-[version]-local.jar是為了跑其他的設(shè)計(jì)的。那么這是否意味著你必須每次跑程序時(shí)指定你要跑哪一個(gè)匈仗，大可不必瓢剿，在文件夾中有g(shù)atk這個(gè)參數(shù)，這是一個(gè)可執(zhí)行的腳本將會(huì)根據(jù)你的命令行選擇適當(dāng)?shù)膉ar悠轩，如果你想的話，當(dāng)然可以用一個(gè)specific的jar攻泼，但是用gatk會(huì)簡單一些火架。 它會(huì)根據(jù)你的情況幫你設(shè)置一些你需要的參數(shù)。

2. 是否安裝成功

為了測(cè)試你是否成功的安裝了GATK忙菠，在你的終端跑下面的命令行何鸡，

./gatk -- help 

3.跑GATK和picard 的命令行

 gatk [--java-options "jvm args like -Xmx4G go here"] ToolName [GATK args go here]

比如說：

gatk --java-options "-Xmx8G" HaplotypeCaller -R reference.fasta -I input.bam -O output.vcf

4.GATK Best Practices

GATK的輸入文件需要做一些預(yù)處理。

（1）數(shù)據(jù)預(yù)處理

在所有call 點(diǎn)突變的步驟前必須的一個(gè)階段牛欢，數(shù)據(jù)的預(yù)處理骡男。 他包括處理原始的測(cè)序數(shù)據(jù)（FASTQ或者uBAM格式）去產(chǎn)生可以分析的BAM文件。 這些步驟包含了比對(duì)到reference 基因組上傍睹，以及一些數(shù)據(jù)清理的操作隔盛，矯正技術(shù)偏差犹菱，將數(shù)據(jù)更加合適分析。

輸入格式

希望輸入的read data的格式是unmapped BAM (uBAM) 格式吮炕。 轉(zhuǎn)換工具可以將數(shù)據(jù)從FASTQ轉(zhuǎn)換到uBAM腊脱。

（2）主要步驟

首先，我們把這個(gè)測(cè)序reads reads 比對(duì)到ref 基因組上去產(chǎn)生一個(gè)SAM/BAM格式的文件龙亲，這個(gè)SAM/BAM格式的數(shù)據(jù)根據(jù)坐標(biāo)順序排列陕凹。下一步，我們標(biāo)記了重復(fù)序列為了減輕由于數(shù)據(jù)產(chǎn)出步驟而帶來的誤差（PCR擴(kuò)增）鳄炉。 最后我們矯正了堿基的質(zhì)量得分杜耙，因?yàn)閏all 突變的算法嚴(yán)重依賴每條測(cè)序read中每個(gè)堿基的質(zhì)量得分。

1 比對(duì)

在這里需要用到的工具是BWA

$RG_string="@RG ID:H0164.2  PL:illumina PU:H0164ALXX140820.2    LB:Solexa-272222    PI:0    DT:2014-08-20T00:00:00-0400 SM:NA12878  CN:BI"
$RG_LB="sample_name"
bwa mem -t 16 -M -R "$RG_string"  fq.1.gz fq.2.gz -o $RG_LB.sam"

2 mark duplicates (標(biāo)記重復(fù))

gatk --java-options "-Xmx20G" MarkDuplicatesSpark
-I $sample.bam -O ${sample}_marked.bam \
-M $sample.metrics \
1>log.mark 2>&1

3 FixMateInformation

gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" FixMateInformation \
    -I ${sample}_marked.bam \
    -O ${sample}_marked_fixed.bam \
    -SO coordinate \
    1>${sample}_log.fix 2>&1 

4 GATK堿基質(zhì)量重矯正（BQSR）：

下載所需要的參考基因組

wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/*

此時(shí)會(huì)下載如下文件

1000G_phase1.snps.high_confidence.hg38.vcf
dbsnp_146.hg38.vcf
Mills_and_1000G_gold_standard.indels.hg38.vcf

然后對(duì)堿基進(jìn)行矯正

gatk --java-options "-Xmx20G" BaseRecalibrator \
   -I ${sample}_marked_fixed.bam \
   -R reference.fasta \
   --known-sites 1000G_phase1.snps.high_confidence.hg38.vcf\
   --known-sites dbsnp_146.hg38.vcf\
   --known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf\
   -O recal_data.table1

gatk  --java-options "-Xmx20G" ApplyBQSR \
   -R reference.fasta \
   -I ${sample}_marked_fixed.bam \
   --bqsr-recal-file  recal_data.table1 \
   -O ${sample}_marked_fixed_BQSR.bam

gatk --java-options "-Xmx20G" BaseRecalibrator \
   -I ${sample}_marked_fixed_BQSR.bam \
   -R reference.fasta \
   --known-sites 1000G_phase1.snps.high_confidence.hg38.vcf\
   --known-sites dbsnp_146.hg38.vcf\
   --known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf\
   -O recal_data.table2

gatk --java-options "-Xmx20G" AnalyzeCovariates \
     -before recal_data.table1 \
     -after recal_data.table2 \
     -plots AnalyzeCovariates.pdf

5. 用mutect2 call點(diǎn)突變

目前mutect2 v4.1.0.0是支持多個(gè)樣本的綜合分析的拂盯。 里面有三種模式（1）tumor-normal mode佑女。 （2） tumor-only mode。 （3） mitochondrial mode磕仅。

a 在v4.1珊豹，不再需要用-tumor 這個(gè)參數(shù)去特異的設(shè)置tumor 樣本。 僅僅是需要用-normal 這個(gè)參數(shù)去設(shè)置正常樣本榕订，如果你有正常樣本的話店茶。

b 在4.0.4.0 開始，GATK提倡使用默認(rèn)參數(shù)--af-of-alleles-not-in-resource劫恒， 工具中不同的模型的參數(shù)設(shè)置不同贩幻。 tumor-only calling 設(shè)置默認(rèn)的是5e-8. tumor-normal calling 設(shè)置的是1e-6。 以及mitochondrial mode設(shè)置的是4e-3两嘴。 在以前的版本里丛楚，默認(rèn)這是的是0.001， 這個(gè)是人類的平均異質(zhì)性憔辫。 關(guān)于其他的物種趣些，把--af-of-alleles-not-in-resource 改成1。

在這里贰您，假設(shè)你已經(jīng)有了normal.bam和tumor.bam, 參考文件是ref.fasta, 目標(biāo)區(qū)域文件是intervals.interval.list .
(0) 構(gòu)建normal panel
如果你有多個(gè)normal.bam做對(duì)照坏平，在開始之前，可以利用這些normal bam構(gòu)建normal panel作為第二步的原始候選變異檢測(cè)輸入集锦亦。

gatk Mutect2 -R ref.fasta  -I normal1.bam -max-mnp-distance 0 -O normal1.vcf.gz

gatk GenomicsDBImport \
-R reference.fasta \
-L intervals.interval_list \
--genomicsdb-workspace-path pon_db \
-V normal1.vcf.gz

gatk CreateSomaticPanelOfNormals \
-R reference.fasta \
-V gendb://pon_db \
-O pon.vcf.gz

（I）Tumor with matched normal
如果有一個(gè)匹配的normal舶替， mutect2是被設(shè)計(jì)僅僅是用來call somatic mutation的。這個(gè)工具可以自動(dòng)的跳過那些已知的存在于germline中的突變杠园。 這樣可以避免一些計(jì)算資源的浪費(fèi)顾瞪。 如果一個(gè)突變處于germline的邊緣。mutect2會(huì)把這個(gè)變量放到callset數(shù)據(jù)集中，被FilterMutectCalls用來做進(jìn)一步的過濾陈醒，或者review惕橙。

 gatk --java-options "-Xmx20G" Mutect2 \
     -R reference.fa \
     -I tumor.bam \
     -I normal.bam \
     -normal normal_sample_name \
     --germline-resource af-only-gnomad.vcf.gz \
     --panel-of-normals pon.vcf.gz \
     --f1r2-tar-gz f1r2.tar.gz \
     -O somatic.vcf.gz

這個(gè)命令產(chǎn)生的輸出文件還需要用FiterMutectCalls過濾。
同時(shí)在mutect2的v4.1還支持來自于同一個(gè)人的多個(gè)腫瘤和正常樣本孵延。 但是必須要設(shè)置這個(gè)-I和-normal參數(shù)吕漂。

gatk --java-options "-Xmx20G" Mutect2 \
     -R reference.fa \
     -I tumor1.bam \
     -I tumor2.bam \
     -I normal1.bam \
     -I normal2.bam \
     -normal normal1_sample_name \
     -normal normal2_sample_name \
     --germline-resource af-only-gnomad.vcf.gz \
     --panel-of-normals pon.vcf.gz \
     -O somatic.vcf.gz

估計(jì)配對(duì)樣本的交叉污染估計(jì)

gatk Pileup \
-R reference.fasta \
-L intervals.interval_list \
-I normal.bam \
-O normal-pileups.table

gatk Pileup \   
-R reference.fasta \   
-L intervals.interval_list \   
-I tumor.bam \   
-O tumor-pileups.table

gatk CalculateContamination \   
-I tumor-pileups.table \   
-matched normal-pileups.table \   
-O contamination.table

測(cè)序偏好矯正以及過濾

gatk LearnReadOrientationModel \
-I f1r2.tar.gz \
-O read-orientation-model.tar.gz

gatk FilterMutectCalls \
-R reference.fasta \
-V somatic.vcf.gz \
--contamination-table contamination.table \
--ob-priors read-orientation-model.tar.gz \
-O filtered.vcf.gz\

最后得到的filtered.vcf.gz就是過濾好的結(jié)果啦，vep注釋起來尘应，發(fā)現(xiàn)GATK4.1也提供了一個(gè)注釋的工具Funcotator惶凝，有興趣也可以嘗試一下。
（ii）Tumor-only mode
這個(gè)模型是針對(duì)一個(gè)樣本犬钢。

 gatk --java-options "-Xmx20G"  Mutect2 \
   -R reference.fa \
   -I sample.bam \
   -O single_sample.vcf.gz

(iii) Mitochondrial mode

 gatk Mutect2 \
  -R reference.fa \
  -L chrM \
  --mitochondria-mode \
  -I mitochondria.bam \
  -O mitochondria.vcf.gz