GATK4.0 和之前的版本相比還是有較大的不同,更加趨於流程化。
軟件安裝
1 wget https://github.com/broadinstitute/gatk/releases/download/4.1.5.0/gatk-4.1.5.0.zip 2 unzip gatk-4.1.5.0.zip
GATK 簡單說明
1 ## 幫助信息 2 gat --help 3 4 ## 列出所有的工具 5 gatk --list 6 7 ## 工具的說明,比如以VariantAnnotator 為例 8 gatk VariantAnnotator --help
GATK分析簡要流程
-
所需數據 : ref.fa
-
- reads1.fq
- reads2.fq
-
-
建立索引
1 bwa index ref.fa 2 samtools faidx ref.fa 3 gatk CreateSequenceDictionary -R ref.fa -O ref.dict 4 5 ## 6 -R Input reference fasta or fasta.gz Required 7 -O 輸出文件
-
比對
1 ## bwa 比對 2 bwa mem -t 4 -R '@RG\tID:id1\tPL:illumina\tSM:test' ref.fa test_1.fq test_2.fq | samtools view -bS - >test.bam 3 4 ##參數 5 -R 設置reads group,gatk必須要的信息,其中ID,PL和SM信息是必須要的 6 7 ## 排序 8 samtools sort -@ 3 -o test.sorted.bam test.bam 9 rm test.bam
GATK 要求read group的格式
ID = Read group identifier
每一個read group 獨有的ID,每一對reads 均有一個獨特的ID,可以自定義命名;
PL = Platform
測序平台;ILLUMINA, SOLID, LS454, HELICOS and PACBIO,不區分大小寫;
SM = sample
reads屬於的樣品名;SM要設定正確,因為GATK產生的VCF文件也使用這個名字;
LB = DNA preparation library identifier
對一個read group的reads進行重復序列標記時,需要使用LB來區分reads來自那條lane;有時候,同一個庫可能在不同的lane上完成測序;為了加以區分,
同一個或不同庫只要是在不同的lane產生的reads都要單獨給一個ID. 一般無特殊說明,成對兒read屬於同一庫,可自定義,比如:library1
若是忘記添加read group信息還以通過 AddOrReplaceReadGroups 添加
1 gatk AddOrReplaceReadGroups -I .bam -O .add.bam -LB library1 -PL illumina -PU pl1 -SM name 2 3 ##參數 4 -I Input file (BAM or SAM or a GA4GH url); 5 -O Output file (BAM or SAM); 6 -LB Read-Group library; 7 -PL Read-Group platform (e.g. ILLUMINA, SOLID); 8 -PU Read-Group platform unit (eg. run barcode); 9 -SM Read-Group sample name
-
標記重復序列
2 gatk MarkDuplicates -I test.sorted.bam -O test.sorted.markdup.bam -M test.sorted.markdup_metrics.txt
3 ##參數 4 -I 排序后的一個或者多個bam或者sam文件 5 -M 輸出重復矩陣 6 -O 輸出文件 7 8 ## 建立索引 9 samtools index test.sorted.markup.bam
-
檢測變異
1 ##兩種方法 2 3 ##(1)多樣本一起call,此次只有一個樣本,若有多個樣本,則繼續用 -I 參數添加即可 4 gatk --java-options -Xmx4G HaplotypeCaller -I test.sorted.markup.bam -O test.gvcf1 -R ref.fa 5 6 ## (2)單個樣本call,然后在合並 7 ## 生成中間文件gvcf 8 gatk --java-options -Xmx4G HaplotypeCaller -I test.sorted.markup.bam -O test.gvcf -R ref.fa --emit-ref-confidence GVCF 9 10 ##通過gvcf檢測變異, -V 添加上步得到的gvcf 11 gatk GenotypeGVCFs -R ref.fa -V test.gvcf -O test.vcf 13 14 ##參數 15 -I BAM/SAM/CRAM file 16 -O 輸出文件 17 -R 參考基因組 18 --java-options: 若設置java則需要添加 19 -Xmx4G:內存為4G,防止內存太大 20 -V A VCF file containing variants
-
提取SNP,INDEL
1 ## 提取SNP 2 gatk SelectVariants -V test.vcf -O test.snp.vcf --select-type-to-include SNP 3 4 ## 提取INDEL 5 gatk SelectVariants -V test.vcf -O test.indel.vcf --select-type-to-include INDEL 6 7 ##參數 8 -O 輸出vcf文件 9 -V 輸入vcf文件 10 --select-type-to-include 選擇提取的變異類型{NO_VARIATION, SNP, MNP, INDEL, 11 SYMBOLIC, MIXED}
-
對vcf文件進行過濾
1 gatk VariantFiltration -O test.snp.fil.vcf.temp -V test.snp.vcf --filter-expression 'QUAL < 30.0 || QD < 2.0 || FS > 60.0 || SOR > 4.0' \ 2 --filter-name lowQualFilter --cluster-window-size 10 --cluster-size 3 --missing-values-evaluate-as-failing 3 4 ## 參數 5 -O 輸出filt.vcf文件 6 -V 輸入vcf文件 7 --filter-expression 過濾條件, VCF INFO 信息 8 --cluster-window-size 以10個鹼基為一個窗口 9 --cluster-size 10個鹼基為窗口,若存在3以上個則過濾 10 --filter-name 被過濾掉的SNP不會刪除,而是給一個標簽, 比如 Filter 11 --missing-values-evaluate-as-failing 當篩選標准比較多的時候,可能有一些位點沒有篩選條件當中的一條或幾條,例如下面的這個表達式;QUAL < 30.0 || QD < 2.0 || FS > 60.0 || MQ < 40.0 || HaplotypeScore > 13.0 並不一定所有位點都有這些信息,這種情況下GATK運行的時候會報很多WARNING信息,用這個參數可以把這些缺少某些FLAG的位點也給標記成沒有通過篩選的。
-
篩選PASS的SNP,INDEL
1 ## 根據FILTER那列信息進行篩選 2 grep PASS test.snp.fil.vcf.temp > test.snp.fil.vcf
歡迎交流