GATK-3.8(最新稳定版)遗传突变分析流程(SNPs和IND

作者: 生信杂谈 | 来源:发表于2017-12-14 21:40 被阅读553次

GATK现在最新的稳定版已经到了3.8，测试版是4.0。3.8版和之前的版本还是有比较大的不同的，但核心算法与4.0的差异不大，4.0主要整合GATK和picard工具并实现并行运算，所以4.0更趋向于流程化。
但可以看到网上很多的教程并没有对GATK3.8做出相应的更新,比如IndelRealigner工具在HaplotypeCaller和MuTect2中不再建议被使用了（在UnifiedGenotyper和MuTect2中依然被要求使用）；还有UnifiedGenotyper已经全面被HaplotypeCaller取代了，官方已经弃用了。
本期详细介绍GATK3.8版本的详细的遗传突变分析流程（包括代码）

GATK pipeline

一、使用bwa进行序列比对

在对reads质控后使用bwa的mem算法将reads比对到基因组上，在进行这一步的时候需要注意的一点就是需要添加"read group"信息到每个bam文件，并且这一步在生成bam文件后使用bamtools添加也可以的。
bam文件的"read group"信息可以通过以下命令获取：

samtools view -H sample.bam | grep '@RG'

"read group"信息以@RG开头：

@RG ID:H0164.2  PL:illumina PU:H0164ALXX140820.2    LB:Solexa-272222    PI:0    DT:2014-08-20T00:00:00-0400 SM:NA12878  CN:BI

"read group"信息每部分所代表的信息如下：

ID = Read group的ID，必须唯一
PU = Platform Unit，这个可以没有，因为ID已经是唯一的了,PU的优先级要比ID高,并且其表示的也是flowcell+line+barcode.
SM = Sample样本序号
PL = Platform/technology used to produce the read 测序技术或机构:
LB = DNA preparation library identifier 表示这个文件是文库的哪一部分，同一文库在不同的flowcell中测，在这里可以标识出来。

所以对于一个被分成4部分在不同flowcell中测序的样本，其使用bwa进行序列比对并添加"read group"信息的代码如下：

# 获得@RG的ID:
RG_ID=ID1
# SM:用样本序号来表示:
RG_SM="sample_1"
# PL表示测序技术或机构:
RG_PL="Illumina"
for((j=1;j<5;j++))
    do
        # RG_LB在每个分割的样本里是不一样的.先提取文件的library part ID
        LB_ID=$j
        # LB表示这个文件是library的哪一部分.
        RG_LB="LB_$LB_ID"
        # 合并RG string:
        RG_string="@RG\tID:$RG_ID\tSM:$RG_SM\tPL:$RG_PL\tLB:$RG_LB"
        echo $RG_string
        # 使用bwa进行map:
       bwa mem -t 16 -M -R "$RG_string" ./GRCh38/Homo_38 ./Clean/"$i"/CL10002*_L?_B5EHUMwkyEAAA?AAA-"$LB_ID"_1.fq.gz  ./Clean/"$i"/CL10002*_L?_B5EHUMwkyEAAA?AAA-"$LB_ID"_2.fq.gz -o ./Outcome/bwa_out/Clean/"$i"/"$RG_LB.sam"
    done

二、标记duplicates并合并bam文件：

在制备文库的过程中，由于PCR扩增过程中会存在一些偏差，也就是说有的序列会被过量扩增。这样在比对的时候，这些过量扩增出来的完全相同的序列就会比对到基因组的相同位置。而这些过量扩增的reads并不是基因组自身固有序列，不能作为变异检测的证据。因此，要尽量去除这些由PCR扩增所形成的duplicates, 使用picard-tools来给这些序列设置一个Flag以标志它们，方便给后面做variants calling时的GATK所识别。这里定义的重复序列是这样的：如果两条reads具有相同的长度而且比对到了基因组的同一位置，那么就认为这样的reads是由PCR扩增而来，就会被GATK标记。
对于一个样本在多个Flowcell下测序得到的多个bam文件，需要先sort bam文件，然后mark duplicates再进行合并多个bam文件为一个。代码及注释如下：

# 对于单样本4个bam文件：
for((i=1;i<=4;i++))
do
   # 先对每个bam文件排序：
    java  -jar picard.jar SortSam I=./"$i"/cleaned.bam O=./"$i"/cleaned_sorted.bam SORT_ORDER=coordinate
    # 鉴定重复的reads:
    java -jar picard.jar MarkDuplicates I=./"$i"/cleaned_sorted.bam O=./"$i"/marked_duplicates.bam  M=./"$i"/marked_dup_metrics.txt
done

# 合并bam文件：,USE_THREADING参数会多占20%的CPU,但可以加快速度.
java -jar picard.jar MergeSamFiles USE_THREADING=true  I=./1/marked_duplicates.bam I=./2/marked_duplicates.bam I=./3/marked_duplicates.bam I=./4/marked_duplicates.bam O=./marked_duplicates.bam

顺便还可以看看文库的插入片段（Insert size）是多少：

java -jar  picard.jar CollectInsertSizeMetrics I=./marked_duplicates.bam  O=./insert_size_metrics.txt H=./insert_size_histogram.pdf M=0.5

Insert size Histogram
可以看出这个样本插入片段主要分布在100-300bp之间，插入片段长度的峰在150bp左右。

三、GATK碱基质量值重校正(BQSR)：

这一步是对bam文件里reads的碱基质量值进行重新校正，使最后输出的bam文件中reads中碱基的质量值能够更加接近真实的与参考基因组之间错配的概率。还需要准备两类外部文件。
一类是knownSites,已知的SNP位点数据在GATK Resource Bundle中下载，下载地址为：ftp://ftp.broadinstitute.org/bundle/
另一类是-L参数指定的基因组区域，GATK将在这个参数指定的区域进行variants calling，如果是外显子测序，则此区域为全基因组外显子区域，但GATK并不支持bed格式，所以这里需要使用Picard的BedToIntervalList工具将bed格式转为GATK能够识别的IntervalList文件。注意，这个intervalList文件在后面要用到多次，设置UNIQUE=true对于结果无影响，因为后续的步骤中GATK将自动将intervalList转为UNIQUE（注意后缀必须被命名为.interval_list）。
代码如下：

# 对exon bed文件排序:
sortBed -i ./ref_files/GENCODEv24_exons.bed  >./ref_files/GENCODEv24_exons_sorted.bed

# 将exon的bed格式转为picard格式:
java -jar  picard.jar BedToIntervalList \
I=./ref_files/GENCODEv24_exons_sorted.bed \
O=./ref_files/GENCODEv24_exons.interval_list \
SD=./GRCh38/Homo_38.dict \
UNIQUE=true

接下来进行BQSR:

Base Quality Score Recalibration (BQSR) 第一步：

$JAVA_HOME/bin/java -jar ./GenomeAnalysisTK.jar \
    -T BaseRecalibrator \
    -R ref.fasta \
    -I input.Bam \
    -knownSites 1000G_phase1.snps.high_confidence.hg38.vcf \
    -knownSites dbsnp_146.hg38.vcf \
    -knownSites Mills_and_1000G_gold_standard.indels.hg38.vcf \
    -L GENCODEv24_exons.interval_list \
    -o  BQSR_1_out_grp

Base Quality Score Recalibration (BQSR) 第二步：

$JAVA_HOME/bin/java -jar ./GenomeAnalysisTK.jar \
    -T BaseRecalibrator \
    -R ref.fasta \
    -I input.Bam \
    -knownSites 1000G_phase1.snps.high_confidence.hg38.vcf \
    -knownSites dbsnp_146.hg38.vcf \
    -knownSites Mills_and_1000G_gold_standard.indels.hg38.vcf \
    -L  GENCODEv24_exons.interval_list \
    -o  BQSR_2_out_grp\
    -BQSR  BQSR_1_out_grp

下面这一步是可选的，画图并产生.csv中间结果文件：

 $JAVA_HOME/bin/java -jar ./GenomeAnalysisTK.jar \
    -T AnalyzeCovariates \
    -R ref.fasta \
    -before BQSR_1_out_grp \
    -after BQSR_2_out_grp \
    -csv AnalyzeCovariates.csv \
    -plots AnalyzeCovariates.png

最后一步：将recal后的PrintReads写入bam文件:

$JAVA_HOME/bin/java -jar ./GenomeAnalysisTK.jar \
    -T PrintReads \
    -R ref.fasta \
    -I input.Bam \
    -o PrintReads_out.Bam \
    -BQSR BQSR_1_out_grp

四、使用HaplotypeCaller进行varinats calling：

目前寻找SNP和indel的工具有很多，但对于生殖遗传(germline)的突变位点识别，使用最多、准确率最高并且是GATK最推荐使用的是HaplotypeCaller。HaplotypeCaller寻找活跃区域，也就是和参考基因组差异较多的区域，对该区域进行局部重组装，确定单倍型（haplotypes），在给定的read数据下，计算单倍型的可能性，分配样本的基因型（genotype）。对于每一个样本，代码如下：

 $JAVA_HOME/bin/java -jar ./GenomeAnalysisTK.jar \
     -R ref.fasta \
     -T HaplotypeCaller \
     -I PrintReads_out.Bam \
     --emitRefConfidence GVCF \
     --dbsnp 1000G_phase1.snps.high_confidence.hg38.vcf \
     -L GENCODEv24_exons.interval_list \
     -o HaplotypeCaller_out.g.vcf

五、使用GenotypeGVCFs进行joint genotype

上一步通过HaplotypeCaller产生的.gVCF文件需要合并每个样本的突变信息到单一vcf文件来方便进行下一步的过滤分析。

# 假设我有8个样本：
$JAVA_HOME/bin/java -jar ./GenomeAnalysisTK.jar \
  -T GenotypeGVCFs \
  -R ref.fasta \
  --variant ./1/HaplotypeCaller_out.g.vcf \
  --variant ./2/HaplotypeCaller_out.g.vcf \
  --variant ./3/HaplotypeCaller_out.g.vcf \
  --variant ./4/HaplotypeCaller_out.g.vcf \
  --variant ./5/HaplotypeCaller_out.g.vcf \
  --variant ./6/HaplotypeCaller_out.g.vcf \
  --variant ./7/HaplotypeCaller_out.g.vcf \
  --variant ./8/HaplotypeCaller_out.g.vcf \
  -o GenotypeGVCFs_out.vcf

提取SNV:

$JAVA_HOME/bin/java -jar ./GenomeAnalysisTK.jar \
  -T SelectVariants \
  -R ref.fasta \
  --variant GenotypeGVCFs_out.vcf \
  -o GenotypeGVCFs_out_snp.vcf \
  --selectTypeToInclude SNP

提取INDEL:

$JAVA_HOME/bin/java -jar ~/software_source/GATK3.8/GenomeAnalysisTK.jar \
  -T SelectVariants \
  -R ref.fasta \
  --variant GenotypeGVCFs_out.vcf \
  -o GenotypeGVCFs_out_indel.vcf \
  --selectTypeToInclude INDEL

四、突变位点质量值重矫正(VQSR)

这一步是基于"true sites"(如HapMap 3和Omni的sites等)建立混合高斯模型来对突变位点质量值重新矫正，并评估该位点作为一个真实突变的可能性。这一步最少样本量为30，但对于INDEL来说，即使有30个样本，也有可能由于突变位点太少也有可能导致建模失败，解决方法是从千人基因组中添加样本进行VQSR,笔者在对16个样本的INDEL进行VQSR时，一度增加至40个样本才成功建模。代码如下：

对于SNP:

# 第一步：Variant recalibrat for SNP:
$JAVA_HOME/bin/java -jar ./GenomeAnalysisTK.jar \
  -T VariantRecalibrator \
  -R $ref_fasta \
  -input GenotypeGVCFs_out_snp.vcf \
  -resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap_3.3_grch38_pop_stratified_af.vcf \
  -resource:omni,known=false,training=true,truth=true,prior=12.0 1000G_omni2.5.hg38.vcf \
  -resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G_phase1.snps.high_confidence.hg38.vcf\
  -resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp_146.hg38.vcf \
  -an QD \
  -an FS \
  -an SOR \
  -an MQ \
  -an MQRankSum \
  -an ReadPosRankSum \
  -an InbreedingCoeff \
  -mode SNP  \
  -tranche 100.0 \
  -tranche 99.9 \
  -tranche 99.0 \
  -tranche 90.0 \
  -recalFile VQSR_recal_out_snp \
  -tranchesFile VQSR_tranch_out_snp \
  -rscriptFile VQSR_rscript_out_snp 

# 第二步：ApplyRecalibration for SNP:
$JAVA_HOME/bin/java -jar ./GenomeAnalysisTK.jar \
  -T ApplyRecalibration \
  -R ref.fasta \
  -input GenotypeGVCFs_out_snp.vcf \
  -mode SNP \
  --ts_filter_level 99.0 \
  -recalFile VQSR_recal_out_snp \
  -tranchesFile VQSR_tranch_out_snp \
  -o VQSR_recalibrated_out_snp

对于INDEL:

# 第一步：Variant recalibrat for INDEL: 
$JAVA_HOME/bin/java -jar ./GenomeAnalysisTK.jar \
    -T VariantRecalibrator \
    -R $ref_fasta \
    -input GenotypeGVCFs_out_indel.vcf  \
    -resource:mills,known=false,training=true,truth=true,prior=12.0 dbsnp_146.hg38.vcf \
    -resource:dbsnp,known=true,training=false,truth=false,prior=2.0 Mills_and_1000G_gold_standard.indels.hg38.vcf \
    -an QD \
    -an FS \
    -an SOR \
    -an MQRankSum \
    -an ReadPosRankSum \
    -an InbreedingCoeff \
    -mode INDEL \
    -tranche 100.0 \
    -tranche 99.9 \
    -tranche 99.0 \
    -tranche 90.0 \
    --maxGaussians 4 \
    -recalFile  VQSR_recal_out_indel \
    -tranchesFile VQSR_tranch_out_indel \
    -rscriptFile VQSR_rscript_out_indel

# 第二步：ApplyRecalibration for INDEL:
$JAVA_HOME/bin/java  -jar ~/software_source/GATK3.8/GenomeAnalysisTK.jar \
    -T ApplyRecalibration \
    -R ref.fasta \
    -input GenotypeGVCFs_out_indel.vcf  \
    -mode INDEL \
    --ts_filter_level 99.0 \
    -recalFile VQSR_recal_out_indel \
    -tranchesFile VQSR_tranch_out_indel \
    -o  VQSR_recalibrated_out_indel.vcf

五、对VCF文件进行过滤

如果在VQSR过程中加入了千人基因组的样本，需要通过SelectVariants选择自己样本相关的突变，去除与自己样本不想管的突变：

$JAVA_HOME/bin/java -jar ./GenomeAnalysisTK.jar \
   -T SelectVariants \
   -R ref.fasta \
   -V VQSR_recalibrated_out_snp.vcf \
   -o SelectVariants_out_snp.vcf \
   # 筛选条件为正则表达式，由于我的样本都是sample_1、sample_2这样的，所以正则为smple_.+。
   -se "sample_.+" \
   -env \
   -ef

选择指定染色体的突变位点进行分析，如果只想选择1-24和x+y染色体的突变位点进行分析，可以使用vcftools进行过滤：

#首先使用vcftools提取24个染色体的variants,再进行下一步分析:
vcftools --vcf variantFiltered_Out_snp.vcf --recode --recode-INFO-all --chr chr1 --chr chr2 --chr chr3 --chr chr4 --chr chr5 --chr chr6 --chr chr7 --chr chr8 --chr chr9 --chr chr10 --chr chr11 --chr chr12 --chr chr13 --chr chr14 --chr chr15 --chr chr16 --chr chr17 --chr chr18 --chr chr19 --chr chr20  --chr chr21  --chr chr22 --chr chrX --chr chrY --out vcftools_24chrom_Out_snp.vcf

选择FILTER列为PASS并且reads深度大于10的突变位点，还是使用SelectVariants工具(SelectVariants是支持JEXL语法的)：

$JAVA_HOME/bin/java -jar ~/software_source/GATK3.8/GenomeAnalysisTK.jar \
   -T SelectVariants \
   -R ref.fasta \
   -V vcftools_24chrom_Out_snp.recode.vcf" \
   -o SelectVariants_24_out_snp.vcf \
   -select "vc.isNotFiltered() && DP>10 "

六、对结果进行进评估

可以通过VariantEval工具获得每个样本的TiTv Ratio:

# 对SNP结果进行评估：For SNP:
$JAVA_HOME/bin/java -jar ~/software_source/GATK3.8/GenomeAnalysisTK.jar \
   -T VariantEval \
   -R $ref_fasta \
   -o $VariantEval_Out_snp \
   --eval $SelectVariants_24_out_snp

结果如下：

TiTv Ratio
一般来说全基因组范围内，颠换/转换一般是在2左右，而在编码区域，颠换/转换一般稍微高于3左右，较高的Ts/TV一般是由于密码子的第三个碱基上发生了颠换。

有段时间没有更新了，非常抱歉，立个Flag:以后至少每周一篇。

更多原创精彩视频敬请关注生信杂谈:

GATK-3.8(最新稳定版)遗传突变分析流程(SNPs和IND

一、使用bwa进行序列比对

二、标记duplicates并合并bam文件：

三、GATK碱基质量值重校正(BQSR)：

接下来进行BQSR:

四、使用HaplotypeCaller进行varinats calling：

五、使用GenotypeGVCFs进行joint genotype

四、突变位点质量值重矫正(VQSR)

对于SNP:

对于INDEL:

五、对VCF文件进行过滤

六、对结果进行进评估

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

生物信息杂谈

生物信息学与算法

生信相关

生物信息

生信工具

生信程序员

测序数据分析代码

WGS专题