GATK4+mutect2 call somatic mutation

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

1 Getting start with GATK

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

contents：

1. 安裝GATK

GATK 工具具有相當簡單的軟件需求： 一個Unix-style OS 以及java 1.8 版本以上，其他的工具需要額外的R和python的依賴。 可以在這個網(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

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

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

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

2. 是否安裝成功

為了測試你是否成功的安裝了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的輸入文件需要做一些預處理。

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

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

輸入格式

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

（2）主要步驟

首先，我們把這個測序reads reads 比對到ref 基因組上去產(chǎn)生一個SAM/BAM格式的文件，這個SAM/BAM格式的數(shù)據(jù)根據(jù)坐標順序排列。下一步，我們標記了重復序列為了減輕由于數(shù)據(jù)產(chǎn)出步驟而帶來的誤差（PCR擴增）。 最后我們矯正了堿基的質量得分，因為call 突變的算法嚴重依賴每條測序read中每個堿基的質量得分。

1 比對

在這里需要用到的工具是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 (標記重復)
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堿基質量重矯正（BQSR）：

下載所需要的參考基因組

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

此時會下載如下文件

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

然后對堿基進行矯正

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點突變

目前mutect2 v4.1.0.0是支持多個樣本的綜合分析的。 里面有三種模式（1）tumor-normal mode。 （2） tumor-only mode。 （3） mitochondrial mode。

a 在v4.1，不再需要用-tumor 這個參數(shù)去特異的設置tumor 樣本。 僅僅是需要用-normal 這個參數(shù)去設置正常樣本，如果你有正常樣本的話。
b 在4.0.4.0 開始，GATK提倡使用默認參數(shù)--af-of-alleles-not-in-resource， 工具中不同的模型的參數(shù)設置不同。 tumor-only calling 設置默認的是5e-8. tumor-normal calling 設置的是1e-6。 以及mitochondrial mode設置的是4e-3。 在以前的版本里，默認這是的是0.001， 這個是人類的平均異質性。 關于其他的物種，把--af-of-alleles-not-in-resource 改成1。

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

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
如果有一個匹配的normal， mutect2是被設計僅僅是用來call somatic mutation的。這個工具可以自動的跳過那些已知的存在于germline中的突變。 這樣可以避免一些計算資源的浪費。 如果一個突變處于germline的邊緣。mutect2會把這個變量放到callset數(shù)據(jù)集中，被FilterMutectCalls用來做進一步的過濾，或者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

這個命令產(chǎn)生的輸出文件還需要用FiterMutectCalls過濾。
同時在mutect2的v4.1還支持來自于同一個人的多個腫瘤和正常樣本。 但是必須要設置這個-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

估計配對樣本的交叉污染估計

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

測序偏好矯正以及過濾

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就是過濾好的結果啦，vep注釋起來，發(fā)現(xiàn)GATK4.1也提供了一個注釋的工具Funcotator，有興趣也可以嘗試一下。
（ii）Tumor-only mode
這個模型是針對一個樣本。

 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
最后編輯于 ：
?著作權歸作者所有,轉載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

相關閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容