java queue GATK_GATK 4.0 全外显子call variant

测试数据：KPGP的WES测序数据，下载地址ftp://ftp.kobic.re.kr/pub/KPGP/2017_release_candidate/WES/，分别下载了KPGP-00265，KPGP-00266，KPGP-00267，KPGP-00270和KPGP-002735组数据 用fastp过滤下raw data数据，以KPGP-00265为例：

~/biosoft/fastp/fastp -i KPGP-00266_L005_R1.fq.gz -I KPGP-00266_L005_R2.fq.gz -o KPGP-00266_L005_R1.clean.fq -O KPGP-00266_L005_R2.clean.fq

然后根据之前WGS的流程GATK 4.0 WGS germline call variant，WES从比对到BQSR其实都是差不多的，可以写到一个简单的shell脚本中，如：

#!/bin/bash

fq1=$1

fq2=$2

sample=$3

index=/home/anlan/reference/index/bwa/gatk_hg38/gatk_hg38

genome=/home/anlan/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta

GATK=/home/anlan/biosoft/GATK4.0/gatk-4.0.5.1/gatk

dbsnp=/home/anlan/annotation/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz

dbsnp1000G=/home/anlan/annotation/GATK/resources/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz

dbindel100G=/home/anlan/annotation/GATK/resources/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz

## BWA Alignment

bwa mem -t 4 -M -R "@RG\tID:$sample\tPL:illumina\tLB:WES\tSM:$sample" $index $fq1 $fq2 1>$sample.sam 2>/dev/null

echo bwa `date`

## Sam to bam and sort

$GATK --java-options "-Xmx12G -Djava.io.tmpdir=/home/anlan/tmp" SortSam \

-I $sample.sam -O $sample.sorted.bam \

-SO coordinate --CREATE_INDEX true

echo gatk-SortSam `date`

## MarkDuplicate

$GATK --java-options "-Xmx12G -Djava.io.tmpdir=/home/anlan/tmp" MarkDuplicates \

-I $sample.sorted.bam -O $sample.sorted.marked.bam \

-M $sample.metrics

echo gatk-MarkDuplicates `date`

## BQSR

$GATK --java-options "-Xmx12G -Djava.io.tmpdir=/home/anlan/tmp" BaseRecalibrator \

-R $genome -I $sample.sorted.marked.bam -O recal_data.table \

--known-sites $dbsnp --known-sites $dbsnp1000G --known-sites $dbindel100G

$GATK --java-options "-Xmx12G -Djava.io.tmpdir=/home/anlan/tmp" ApplyBQSR \

-R $genome -I $sample.sorted.marked.bam --bqsr-recal-file recal_data.table \

-O $sample.sorted.marked.BQSR.bam

echo gatk-BQSR `date`

接下来则是call variant步骤，WES可以像WGS一样先call个g.vcf文件，然后再到VCF文件，如:

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk --java-options "-Xmx12G -Djava.io.tmpdir=/home/anlan/tmp" HaplotypeCaller \

-R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta \

-I KPGP-00265.sorted.marked.BQSR.bam \

--dbsnp ~/annotation/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz \

-ERC GVCF -O KPGP-00265.g.vcf

但从GATK 4.0版本起，GenotypeGVCFs只支持a single single-sample GVCF，a single multi-sample GVCF created by CombineGVCFs 以及a GenomicsDB workspace created by GenomicsDBImport

所以我们需要在用GenotypeGVCFs前需要将多个样本的g.vcf文件用CombineGVCFs方式或者GenomicsDBImport方式合并成一个文件，前者(比较传统)是一个总的g.vcf文件，后者是一个GenomicsDB(XX.db)

CombineGVCFs方法：

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk CombineGVCFs \

-R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta $(for i in $(ls *.vcf);do echo "--variant $i";done) \

-O combined.g.vcf

用GTAK的GenotypeGVCFs工具对g.vcf文件进行 joint genotyping，如：

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk --java-options "-Xmx12G -Djava.io.tmpdir=/home/anlan/tmp" \

GenotypeGVCFs \

-R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta \

-V combined.g.vcf \

-O combined.vcf

GenomicsDBImport方法(这里需要注意的是其必须要输入一个区间One or more genomic intervals，所以可以选择分染色体进行)：

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk GenomicsDBImport $(for i in $(ls *.vcf);do echo "-V $i";done) \

-L chr1 \

--genomicsdb-workspace-path my_database_chr1

接着也是跟上面一样，用GenotypeGVCFs工具来joint genotyping，只是参数变成-V gendb://my_database_chr1

但WES在建库过程中有个捕获的过程，这里简单的说就是只测了人类外显子的区域，那么在call variant过程中我们不必要对全基因组范围进行计算，只需要人为定义一个区域(比如捕获区域的bed文件，或者外显子数据库的bed文件也行)

一般厂家的人全外显子版本的目标捕获区域是整合了多个数据库的，如安捷伦V7版本的全外是基于最新的RefSeq、GENCODE、CCDS 和 UCSC Known Genes 数据库

但我的是测试数据，也不知道其当初的捕获区域是哪个，所以简单的以CCDS数据库为例，首先需要先下载个最新的CCDS文件 (CCDS.20180614.txt) ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human

然后提取其外显子的起始/终止碱基位置，然后以150bp测序话，可以再加减150bp区间，由于bed格式是以0-based for the start coordinates，所以还要再减去1，最终做成一个bed格式文件，bed文件需要有3列信息chr，start位置和stop位置

其实GATK -L参数可支持的文件不仅仅是bed文件，还有Picard-style.interval_list，GATK-style.list，具体可以参照Intervals and interval lists

说个巨坑遭遇：

起初我的将CCDS文件处理成bed文件格式作为GATK HaplotypeCaller工具-L参数输入文件，但是GATK却意外提示了，没有丝毫ERROR等报错信息，网上搜了好久也没有查到有用的信息，但至少确定是bed文件的问题，但是什么问题却不晓得；最终我决定再仔细看了GATK官方文档，才查到上面网站所说的3中输入格式；因此我抱着试试的想法用了gatk-style作为输入文件，结果GATK报错了！！！终于等到了报错的提示；看了下，原来是我的intervals文件的坐标超出了参考基因组的序列长度，我将一些超过区间长度的坐标从intervals文件删除后就正常运行了！！！

这里放一个GATK-style的intervals文件exon.list，虽然不是WES捕获区域最适合的文件，但是凑合用用还是可以的

现在可以在HaplotypeCaller生成g.vcf文件的时候就指定好区间intervals(这样可以只检测出捕获区域以内的突变，并减少运行时间)，如：

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk --java-options "-Xmx8G -Djava.io.tmpdir=/home/anlan/tmp" HaplotypeCaller \

-R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta \

-L exon.list \

-I KPGP-00265.sorted.marked.BQSR.bam \

-ERC GVCF -O KPGP-00265.g.vcf

最后再重复上述步骤合并g.vcf文件，然后joint genotyping生成vcf文件