GATK4 germline variant calling 分析流程

通过对基因组高通量测序数据的分析可以从中获取Genome Varient的信息。这个过程被称为Call variant，现在已经有很多成熟的软件和流程来完成这个任务。本笔记主要记录了使用GATK4软件来分析 Variant的方法。

软件安装

为了方便，我直接使用conda的环境来安装Call variant所需要的软件。主要的软件版本入下： * GATK4: 4.1.8，主要的分析软件 * Picard: 2.22.8，SAM/BAM处理软件 * BWA: 0.7.17，Mapping 软件

bioconda的channel中包含这些生信分析的软件，bioconda的添加参考介绍。

因为GATK4使用Java 1.8，为了防止和其他软件冲突，最好还是在一个独立的环境中安装。

$ conda create -n gatk gatk4 picard bwa samtools bcftools

这样，我们就创建了一个gatk的环境，其中安装了我们需要的GATK和Picard等软件。

Call variant 流程

使用GATK4来进行Variant calling主要参考了官方 Germline Variant Calling 流程。这个分析流程是从已经处理好的BAM文件开始的。生成可以使用的BAM文件则参考了官方 Data processing 流程。

测序数据比对到基因组

所得到的原始测序数据是FASTQ格式的，经过预先质量控制（Quality Control）后就可以比对（Mapping）到基因组。GATK4软件推荐的比对软件是由Heng Li开发的BWA。这里用到的参考基因组是由 illumina整理好的hg38/GRCh38 iGenome 人类参考基因组。

REF_BWA=./reference/genome.fa
THREAD=4
FQ1=./data/sample.1.fq
FQ2=./data/sample.2.fq

bwa mem -t ${THREAD} ${REF_BWA} ${FQ1} ${FQ2} | \
    samtools view -@ $((THREAD-1)) -b - > ./data/sample.raw.bam

bwa mem 可以接受双端测序（paired-end），默认会把比对结果输出到标准输出（stdout）。
bwa mem 还有一个 -R 参数，可以加上值 '@RG\tID:foo\tSM:bar' 指定读取分组，如果设置错误GATK处理会有问题。这个RG可以在比对的时候添加，也可以用Picard进行添加，见下面。
使用 samtools view 可以把bwa mem的SAM格式的标准输出转为二进制的BAM格式。
$((THREAD-1))是bash变量的数学计算，这里的实际结果是4 - 1=3。 samtools view 的 -@ 参数接受的是增加的线程数，所以这里要减一。

用Picard修改BAM文件的Read Group

如果在序列比对的时候没有在 bwa mem 中设置Read Group，可以用Picard的 AddOrReplaceReadGroups命令来完成。Read Group主要是表示这些序列是同一次测序出来的结果，在后面的一些矫正是在一次测序内部进行的，所以需要设置一下。如果BAM文件中的序列只有单一来源的话，就只要设置这些Read Group完全一样就可以。下面是通过 Picard 来修改和添加 Read Group。

picard AddOrReplaceReadGroups \
    -I ./data/sample.raw.bam \
    -O ./data/sample.AddOrReplaceReadGroups.bam \
    --RGID sample1 \
    --RGLB lib1 \
    --RGPL illumina \
    --RGPU unit1 \
    --RGSM sample1 \
    --VALIDATION_STRINGENCY LENIENT \
    -Xms1g -Xmx10g -XX:ParallelGCThreads=4

--VALIDATION_STRINGENCY LENIENT这里是为了避免Picard因为检查Header 不不符合规范而报错终止。
-Xms1g 是 Java Virtual Machine的参数，设置最小内存为1G。
-Xmx10g 也是 Java Virtual Machine的参数，设置最大内存为10G。
-XX:ParallelGCThreads=4 还是 Java Virtual Machine的参数，设置并行线程为4。

标注重复序列

在做基因组测序的时候，序列会因为PCR的缘故而扩增太多，产生的许多重复序列（Duplicates）会影响到后面计算比例，所以需要区分出那些序列是由于扩增引起的。这里用Picard的 MarkDuplicates 命令来标注重复序列。这些序列不会从BAM中被剔除，而只会在 FLAG处被标注为重复序列。

picard MarkDuplicates \
    -I ./data/sample.AddOrReplaceReadGroups.bam \
    -O ./data/sample.MarkDuplicates.bam \
    -M ./data/sample.MarkDuplicates.matrics \
    -Xms1g -Xmx10g -XX:ParallelGCThreads=4

对BAM文件进行排序

在标注完重复序列后，需要对BAM文件按照染色体坐标进行排序。

picard SortSam \
    -I ./data/sample.MarkDuplicates.bam \
    -O ./data/sample.SortSam.bam \
    -SO "coordinate" \
    --CREATE_INDEX false \
    --CREATE_MD5_FILE false \
    -Xms1g -Xmx10g -XX:ParallelGCThreads=4

设置其他SAM标签

为了方便后面计算，可以设置SAM文件的几个标签：NM、MD和UQ。NM是编辑距离，即序列可以经过多少编辑变为参考序列。MD是表示错配的字符串。UQ是以Phred 形式表示的匹配到参考序列位置的概率。

REF_FA=./Reference/genome.fa

picard SetNmMdAndUqTags \
    -I ./data/sample.SortSam.bam \
    -O ./data/sample.tagged.bam \
    --CREATE_INDEX true \
    --CREATE_MD5_FILE true \
    --REFERENCE_SEQUENCE ${REF_FA} \
    -Xms1g -Xmx10g -XX:ParallelGCThreads=4

矫正碱基质量得分（Base Quality Score）

碱基质量是对每个碱基在测序时候估计测序错误概率得出的分值。而测序仪在评估碱基质量的时候会由于种种原因出现偏差。因为碱基质量得分在后面的Call Variant过程中会被用做评估碱基的可信程度，所以需要在预处理数据中进行矫正，避免测序仪的系统偏差。这里利用GATK4 的 BaseRecalibrator 和 ApplyBQSR 命令来进行矫正。其中 Baserecalibrator 通过已知的一些位点在测序数据中的质量的情况来建立机器学习模型，然后 ApplyBQSR 把模型应用到所有位点进行矫正。

export REF_FA=./Reference/genome.fa
export REF_dbSNP=./Reference/Homo_sapiens_assembly38.dbsnp138.vcf.gz
export REF_Millsindel=./Reference/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
export REF_knownindel=./Reference/Homo_sapiens_assembly38.known_indels.vcf.gz

gatk BaseRecalibrator \
    -R ${REF_FA} \
    -I ./data/sample.tagged.bam \
    -O ./data/sample.bqsr.table \
    --use-original-qualities \
    --known-sites ${REF_dbSNP} \
    --known-sites ${REF_Millsindel} \
    --known-sites ${REF_knownindel}

gatk ApplyBQSR \
    -R ${REF_FA} \
    -I ./sample.tagged.bam \
    -O ./data/sample.bqsr.bam \
    -bqsr ./data/sample.bqsr.table \
    --static-quantized-quals 10 \
    --static-quantized-quals 20 \
    --static-quantized-quals 30 \
    --add-output-sam-program-record \
    --create-output-bam-md5 \
    --use-original-qualities

其中所用到的已经知道的VCF文件需要使用 bgzip 来压缩，并利用 tabix 来建立索引。

Call Variant

有了经过矫正的BAM文件后,可以通过GATK的HaplotypeCaller命令来获取基本的VCF。

gatk --java-options "-Xmx10G -XX:ParallelGCThreads=4" \
    HaplotypeCaller \
    -R ${REF_FA} \
    -I ./data/sample.bqsr.bam \
    -O ./data/sample.g.vcf.gz \
    -OVI \
    -contamination 0 \
    -G StandardAnnotation \
    -G StandardHCAnnotation \
    -G AS_StandardAnnotation \
    -GQB 10 -GQB 20 -GQB 30 -GQB 40 -GQB 50 -GQB 60 -GQB 70 -GQB 80 -GQB 90 \
    -ERC GVCF \
    -bamout ./data/sample.HaplotypeCaller.bam

HaplotypeCaller 可以直接输出压缩的 GVCF，并且可以建立 index（-OVI）。
HaplotypeCaller 还可以把Variant附近的序列输出到BAM文件（-bamout）。

对 SNP Variant 得分的矫正

为了更准确地过滤 VCF，可以对 Variant 的打分进行矫正。对 Variant 的矫正也类似对碱基的矫正，有两个步骤，需要 GATK 的 VariantRecalibrator 和 ApplyVQSR 两个命令。其中 VariantRecalibrator 也是从已知的 SNP 中建立模型，而在 ApplyVQSR 中进行矫正。需要注意的是，当前的版本只能每次单独对 SNP 进行矫正或者对 Indel 进行矫正，而不能同时矫正两个。这里主要矫正 SNP。

gatk VariantRecalibrator \
    -R ${REF_FA} \
    -V ./data/sample.g.vcf.gz \
    --resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap_3.3.hg38.sites.vcf.gz \
    --resource:omni,known=false,training=true,truth=false,prior=12.0 1000G_omni2.5.hg38.sites.vcf.gz \
    --resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G_phase1.snps.high_confidence.hg38.vcf.gz \
    --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 Homo_sapiens_assembly38.dbsnp138.vcf.gz \
    -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \
    -mode SNP \
    -O ./data/sample.recal \
    --tranches-file ./data/sample.tranches \
    --rscript-file ./data/sample.plots.R

gatk ApplyVQSR \
    -R ${REF_FA} \
    -V ./data/sample.g.vcf.gz \
    -O ./data/sample.recalibrated.g.vcf.gz \
    -OVI \
    --truth-sensitivity-filter-level 99.0 \
    --tranches-file ./data/sample.tranches \
    --recal-file ./data/sample.recal \
    -mode SNP

多任务并行

在有的时候需要处理的测序样品比较多，想要同时处理多个数据，可以考虑使用xargs来同时运行多个数据。

echo -e "sample1 ./bam/sample1.bam" > files.txt

export DIR_WK=$(pwd)

function MarkDuplicates {
    LAB=$1
    BAM=$2
    OUTDIR=${DIR_WK}/bam_Sorted
    TMPDIR=${DIR_WK}/tmp
    mkdir -p ${OUTDIR}
    mkdir -p ${TMPDIR}
    OUTBAM=${OUTDIR}/${LAB}.Markduplicates.bam
    OUTTXT=${OUTDIR}/${LAB}.MarkDuplicates.matrics
    picard MarkDuplicates \
           -I ${BAM} \
           -O ${OUTBAM} \
           -M ${OUTTXT} \
           -Xms1g -Xmx10g -XX:ParallelGCThreads=4
}

export -f MarkDuplicates

cat files.txt | xargs -P 4 -n 1 -I {} bash -c 'MarkDuplicates $@' _ {}

files.txt中保存着需要处理的数据，数据有两列，第一列是样品名称，第二列是BAM文件的路径。
xargs 的 -P参数指定了并行的任务的数量。
使用 xargs来传参数到自定义的方程需要使用bash -c来运行命令行形式的 bash脚本，后面的_ {}中的_是用来站第0位的参数的，这样函数中的参数 $1对应的就是数据文件中的第一列，$2对应的是数据的第二列。
为了使得函数内和bash -c使用变量，就需要export所用到的变量。