GATK4.1 call SNP


 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

 

  

歡迎交流

  

GATK4.0全基因組數據分析實戰

GATK - Read groups


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM