基于GATK（Genome Analysis Toolkit）进行SNP calling

GATK（Genome Analysis Toolkit）是进行DNA和RNAseq数据变异检测的常用工具，目前已成为变异检测的“金标准”。
本文提供其与其他软件联合使用进行SNP calling的方法。
所需数据：
参考基因组及其索引：ref_genome.fa
PE reads数据：read1.fq.gz、read2.fq.gz

1 比对

首先使用bwa-mem2进行比对，注意ID的设置必须与样品名一致且唯一。

bwa-mem2 \
	mem \
	-M \
	-t 10 \
	-R "@RG\tID:sample1\tSM:sample\tPL:Illumina" \
	ref_genome.fa \
	read1.fq.gz \
	read2.fq.gz \
	> sample1.sam
1
2
3
4
5
6
7
8
9

2 sort

使用gatk的SortSam工具进行reads的sort：

gatk \
	SortSam \
	-I sample1.sam \
	-O sample1.bam \
	-SO coordinate \
	-VALIDATION_STRINGENCY LENIENT \
	-CREATE_INDEX true \
	-TMP_DIR ./sample1tmp.sort
1
2
3
4
5
6
7
8

3 标记重复序列

比对到参考基因组上的下机数据存在PCR冗余的重复reads，需要将重复的reads进行标记：

gatk \
	MarkDuplicates \
	-I sample1.bam \
	-O sample1.dup.bam \
	-M sample1.dup.txt \
	-REMOVE_DUPLICATES true \
	-VALIDATION_STRINGENCY LENIENT \
	-CREATE_INDEX true \
	-TMP_DIR sample1tmp.dup
1
2
3
4
5
6
7
8
9

4 统计比对情况（可选）

samtools可进行比对情况的统计：

samtools \
	flagstat \
	sample1.dup.bam \
	> sample1.dup.bam.stat
1
2
3
4

5 变异检测

该步骤是真正意义上的snp calling:

gatk \
	HaplotypeCaller \
	-R ref_genome.fa \
	-I sample1.dup.bam \
	-O sample1.dup.vcf
1
2
3
4
5

Ending