GATK(全称The Genome Analysis Toolkit)是Broad Institute开发的用于二代重测序数据分析的一款软件，里面包含了很多有用的工具。

需要的文件

参考基因组：reference.fasta

测序文件：bam 格式比对结果文件

所需软件： samtools 、gatk 、 Picard

---

一、创建gatk所需的索引文件

GATK需要reference序列是经过index的，而且需要两个index文件，一个是后缀名为.fai的，另外一个是后缀名称为.dict的，缺少这些文件，或者两个文件中的内容不一致都可能导致程序报错。

samtools faidx reference.fasta     #生成后缀名为.fai

gatk CreateSequenceDictionary-R reference.fasta  -O reference.dict   #生成后缀名为.dict

二、 增加GATk 要求的read group的格式

GATK要求输入的bam文件包含Read groups，如果没有就会报错。

Read group是@RG开始，包括以下几个部分：

ID= Read group identifier

每一个Read group独有的ID；

PU= Platform Unit

PL= Platform/technology used to produce the read

测序使用的平台： ILLUMINA, SOLID, LS454, HELICOS and PACBIO。

LB= DNA preparation library identifier

对一个read group的reads进行重复序列标记时，需要使用LB来区分reads来自那条lane;有时候，同一个库可能在不同的lane上完成测序;为了加以区分，同一个或不同库只要是在不同的lane产生的reads都要单独给一个ID。

SM= Sample

reads属于的样品名，自由设定

可以在BWA比对时增加read group:

bwa mem -R '@RG\tID:group\tLB:library\tPL:illumina\tPU:unit1\tSM:676R' ~/ref/reference.fasta read1.fq read2.fq > bulk.sam

或者使用Picard增加：

Picard: 它是目前最著名的组学研究中心-Broad研究所开发的一款强大的NGS数据处理工具，功能方面和Samtools有些重叠，但更多的是互补，它是由java编写的，我们直接下载最新的.jar包就行了。

下载链接：wget https://github.com/broadinstitute/picard/releases/download/2.25.5/picard.jar

安装    

java -jar picard.jar

执行命令：

java -jar ~/biosoft/picard.jar AddOrReplaceReadGroups I=bulk.bam O=bulk.add.bam RGID=4 RGLB=library1 RGPL=illumina RGPU=unit1 SORT_ORDER=coordinate RGSM=M05

三、去除PCR重复序列

gatk MarkDuplicates -I bulk.add.bam -O bulk.marked.bam -M bu.metrics 1>log.mark 2>&1

四、VCF输出

#先生成gvcf格式文件     gvcf可记录所有位点的变异情况

gatk HaplotypeCaller -R  ~/ref/reference.fasta   -I bulk.marked.bam -O output.g.vcf.gz -ERC GVCF

#然后在从gvcf提取变异情况

gatk GenotypeGVCFs -R  ~/ref/reference.fasta   -V output.g.vcf.gz -O output.vcf.gz

 或者直接生成vcf文件

gatk HaplotypeCaller -R ~/ref/reference.fasta -I bulk.marked.bam -O   out.vcf

五、执行筛选

gatk VariantFiltration \

    -V  out.vcf \

    -filter "QD < 2.0" --filter-name "QD2" \

    -filter "QUAL < 30.0" --filter-name "QUAL30" \

    -filter "FS > 60.0" --filter-name "FS60" \

    -filter "MQ < 40.0" --filter-name "MQ40" \

    -filter "MQRankSum < -12.5" --filter-name "MQRankSum-12.5" \

    -filter "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum-8" \

    -O  out_prefix.vcf

六、合并VCF文件

第一种方法：

gatk -T CombineVariants -V file1.vcf.gz -V file2.vcf.gz -o merge.vcf.gz -R ref.fa

第二种方法：

bcftools merge file1.vcf.gz  file2.vcf.g  -o merge_bcftools.vcf