GATK4.1 call SNP

阅读原文时间：2023年07月12日阅读：1

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

1 ## 帮助信息
2 gat --help
3
4 ## 列出所有的工具
5 gatk --list
6
7 ## 工具的说明，比如以VariantAnnotator 为例
8 gatk VariantAnnotator --help

所需数据： 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