GATK 4.0 全外显子call variant

测试数据:KPGP的WES测序数据,下载地址ftp://ftp.kobic.re.kr/pub/KPGP/2017_release_candidate/WES/,分别下载了KPGP-00265KPGP-00266KPGP-00267KPGP-00270KPGP-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捕获区域最适合的文件,但是凑合用用还是可以的

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转载请注明出处