526互联

生信： 一起读官方文档 GATK2.1版本 篇

发布时间 2023-08-29 17:53:53作者: 仗剑天涯横刀笑

一起读官方文档 GATK2.1版本 篇

参考文章:

GATK使用：https://www.plob.org/article/7070.html

GATK介绍

GATK做什么的？

它主要用于从sequencing 数据中进行variant calling，包括SNP、INDEL。比如现在风行的exome sequencing找variant，一般通过BWA+GATK的pipeline进行数据分析。

BWA流程上一篇文章已经讲完了，这一篇主要讲一下GATK2版本的使用。

GATK参数 -- GATK2.1版本

这里只讲述了在BWA + GATK流程中的GATK2.1版本应用

STEP01: 局部重新对齐

分为两个小步骤1.RealignerTargetCreator 和2. IndelRealigner

  1. GATK RealignerTargetCreator

RealignerTargetCreator 用于找到那些可能需要进行局部重新对齐的区域。

-Xmx256g 分配了最多 256GB 的内存给 Java 虚拟机

-T RealignerTargetCreator 使用 RealignerTargetCreator 工具。

-R 指定参考基因组

-I 指定输入的 BAM 文件

-o输出文件，其中包含了需要进行局部重新对齐的目标区域。这通常是一个 .intervals 文件

# 实例
java -Xmx256g -jar GenomeAnalysisTK-2.1/GenomeAnalysisTK.jar -T RealignerTargetCreator -R $ref -I bam/$file.bam -o bam/$file.realn.intervals
  1. GATK IndelRealigner

IndelRealigner根据上面的RealignerTargetCreator对这些局部区域进行重新对齐

-T IndelRealigner使用 IndelRealigner 工具。

-targetIntervals指定由上一步 RealignerTargetCreator 生成的 .intervals 文件路径，主要表示哪些区域需要进行重新对齐

java -Xmx256g -jar GenomeAnalysisTK-2.1/GenomeAnalysisTK.jar -T IndelRealigner -R $ref -targetIntervals bam/$file.realn.intervals -I bam/$file.bam -o bam/$file.realn.bam

注意RealignerTargetCreator和IndelRealigner是一块使用的，这在GATK4版本中已经可以使用HaplotypeCaller

STEP02: picard进行MarkDuplicates

MarkDuplicates 工具对GATK的IndelRealigner后的$file.realn.bam进行 标记和处理重复的序列读取。

-Xmx4g: 为 Java 虚拟机分配最多 4GB 的内存

-XX:ConcGCThreads=2: 设置并发垃圾回收器线程的数量为 2

-XX:ParallelGCThreads=12: 设置用于并行垃圾回收的线程数量为 12

MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000: 设置用于存储 read ends 信息的最大文件句柄数

TMP_DIR=TMP: 设置临时目录。

INPUT=bam/$file.realn.bam: 输入 BAM 文件

OUTPUT=bam/$file.realn.marked.bam: 包含标记了重复读取信息的BAM文件

METRICS_FILE=bam/$file.realn.marked.txt: 输出的文本文件，其中包含有关重复标记的各种度量和统计信息。

VALIDATION_STRINGENCY=LENIENT: 验证严格性的设置。LENIENT 表示在遇到不完全符合标准的行时发出警告，而不是终止。

java -Xmx4g -XX:ConcGCThreads=2 -XX:ParallelGCThreads=12 -jar picard.jar MarkDuplicates MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 \
TMP_DIR=TMP \
INPUT=bam/$file.realn.bam \
OUTPUT=bam/$file.realn.marked.bam \
METRICS_FILE=bam/$file.realn.marked.txt \
VALIDATION_STRINGENCY=LENIENT 

# samltools对bam文件建立索引后续使用
samtools index bam/$file.realn.marked.bam

STEP03: GATK UnifiedGenotyper进行基因变异检测

UnifiedGenotyper同时检测单核苷酸变异（SNPs）和短插入/缺失（indels）。最新版本GATK4已经使用更先进的工具HaplotypeCaller，之后会出一期专门的文章阅读GATK4教程。下面是对该GATK2.1中的UnifiedGenotyper参数讲解。

-Xmx256g: 为 Java 虚拟机分配最多 256GB 的内存。

-T UnifiedGenotyper: 指定使用 UnifiedGenotyper 工具。

--num_threads 30: 指定使用 30 个线程进行计算。

-R $ref: 指定参考基因组的路径（这里由 shell 变量 $ref 保存）。

-I bam/$file.realn.marked.bam -I bam/WT.realn.marked.bam: 指定输入的 BAM 文件。这里有两个输入文件：一个是样本（由 $file 表示），另一个是野生型（WT）。

-o VCF/$file-WT.GATK2.1.vcf: 指定输出的 VCF 文件的路径和名称。

-stand_call_conf 30.0: 设置变异检测的置信度阈值为 30。

-glm BOTH: 指定同时检测 SNPs 和 indels。

-rf BadCigar: 使用BadCigar读取过滤器来过滤具有不合规或损坏的 CIGAR 字符串的读取。

java -Xmx256g -jar GenomeAnalysisTK-2.1/GenomeAnalysisTK.jar -T UnifiedGenotyper \
--num_threads 30 \
-R $ref \
-I  bam/$file.realn.marked.bam -I bam/WT.realn.marked.bam \
-o VCF/$file\-WT.GATK2.1.vcf \
-stand_call_conf 30.0 \
-glm BOTH \
-rf BadCigar

使用bcftools查看一下生成的vcf文件 : bcftools view W356-WT.GATK2.1.vcf.gz | less -SN

image-20230829170008725