转载

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-00273 5组数据

用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捕获区域最适合的文件，但是凑合用用还是可以的

GTAK对于是否需要在一些工具中用-L参数给出了一定的说明 When should I use -L to pass in a list of intervals? 以及 When should I restrict my analysis to specific intervals?

现在可以在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文件

本文出自于 http://www.bioinfo-scrounger.com 转载请注明出处

原文 http://www.bioinfo-scrounger.com/archives/642

正文到此结束