GATK变异检测

变异检测

GATK 变异检测

GATK是Genome Analysis Toolkit的缩写，是用来处理高通量测序数据的一套软件。最初，GATK被设计用来分析人类基因组和外显子，主要用来寻找SNP和indel。后开，GATK的功能越来越丰富，增加了short variant calling、计算copy number（CNV）和结构变异（SV）等新功能。同时，GATK也越来越广泛地应用于其他物种的数据分析中。现在，GATK已经成为了基因组和RNA-seq分析过程中，寻找变异的行业标准。

数据类型：Illumina数据

软件版本：Gatk 4.1.8.1、fastp 0.20.0、bwa 0.7.17、samtools 1.9

测试数据：

$ref = ~/database/hg19_BWA /hg19.fa
$tumor1.fq = ~/GATK_passway/Illumina测序文件/ 202011_R1.fq
$tumor2.fq = ~/GATK_passway/Illumina测序文件/ 202011_R2.fq
$normal1.fq =  ~/GATK_passway/Illumina测序文件/ 2020NC_R1.fq
$normal2.fq = ~/GATK_passway/Illumina测序文件/ 2020NC_R2.fq
$bed = ~/GATK_passway/Illumina测序文件/ Illumina.bed

GATK的使用流程：

Gatk的使用流程共3个步骤：原始数据的处理，变异检测，过滤质控

第一部份：原始数据的处理

1. 对原始数据进行质控，去除掉质量较低的reads

Fastp –i ${tumor1.fq} –o ${tumor1.clean.fq} –I ${tumor2.fq} –O ${ tumor2.clean.fq}

2. 建立ref的dict索引，便于对ref文件的查询

后续比对，建立panel normal，bed文件处理等

Samtools dict ${ref} -o ${ref.dict}

文件格式：软件的版本，染色体号，染色体长度，MD5，来源

3. 创建参考基因组的索引，用于bwa的比对

Bwa index ${ref}

建立完成后文件夹下有 hg19.fa.amb、hg19.fa.ann、hg19.fa.bwt 、hg19.fa.pac、hg19.fa.sa，5个文件

4. 将质控后的文件与数据库文件做比对

Bwa men –R ‘@RG\tID:id1\tPL:${type}\tSM:${sm}’ ${ref} ${tumor1.clean.fq1} ${tumor2.clean.fq2} | samtools view –bS > ${tumor.bam}

@RG：reads group，ID：reads的id号， PL:下机数据的平台，例如：Illumina，SM：样本名
这个步骤直接将比对后的sam文件转化为了bam

5. 比对后bam文件排序，统计比对情况

将比对后的文件进行排序

Samtools sort ${out.bam} –o ${out.sorted.bam}

建立索引便于查询和搜索

Samtools index ${out.sorted.bam} 

统计比对后的信息，查看比对情况

Samtools flagstat ${out.sorted.bam}

6. 标记重复序列

Gatk MarkDuplicates –I ${.sorted.bam} –O ${.sorted.marked.bam}-M { .txt}

将pcr过程中的多拷贝或高密度的flowcell使得测序的通量显著提升产生的重复，进行标记

7. 将bed文件转化为GATK的所需的格式，同时进行窗口划分

Gatk PreprocessIntervals –L ${bed} –R${ref} --bin-length 0 --interval-merging-rule OVERLAPPING_ONLY –O ${list}	

–bin-length：为你创建的窗口的大小，interval-merging-rule ：OVERLAPPING_ONLY 有重叠区域时合并窗口

第二部分：变异检测

1. somatic call SNP 和INDEL的流程图

call somatic 变异的流程图
对于成对样本的somatic变异检测，将normal的样本创建为panel of normal ，通过GATK的Mutect2的算法进行tumor样本的变异进行筛选

mutect2 的算法基本原理，根据Normal panel，划分活跃区域，对区域内的序列重新组装，建立相似性矩阵，确定基因型，从而获得变异信息

构建normal panel

对normal样本进行call变异处理

Gatk Mutect2 –R ${ref} –I ${.sorted.marked.bam } –O ${normal.vcf}

本次未引入其他数据库，如有需要自行下载

获取tumor样本中的原始变异信息

获取原始变异时，按照变异的频率进行筛选，过滤掉低频率的变异信息

成对样本tumor and normal的变异检测

Gatk Mutect2 -R ${ref} -I ${tumor.MarkDup.sorted.bam} -I ${normal.MarkDup.sorted.bam} -tumor ${sm1} -normal ${sm2} -pon ${normal.vcf} --minimum-allele-fraction ${int} -L ${list} -O ${out.vcf}

将假阳性的结果进行标记

Gatk FilterMutectCalls -V ${out.vcf} -R ${ref} -O ${out.Filtered.vcf}

grep "^#|PASS" ${out.Filtered.vcf} > ${output.vcf}

tumor-only 的变异检测

Gatk Mutect2 -R ${ref} -I ${tumor.MarkDup.sorted.bam} -tumor ${sm} --minimum-allele-fraction ${int} -pon ${normal.vcf} -O ${out.vcf}

somatic变异检测中对于panel of normal 的构建会帮你去除掉许多的变异,变异结果更为的准确

2. Germline SNP和INDEL的变异检测

1.	寻找活跃区域，就是和参考基因组不同部分较多的区域
2.	通过对该区域进行局部重组装，确定单倍型（haplotypes）
3.	在给定的read数据下，计算单倍型的可能性。
4.	分配样本的基因型

对于Germline的变异检测

Gatk HaplotypeCaller -R ${ref} -I ${.sorted.marked.bam} -O ${out.vcf}

3. 拷贝数变异检测

统计每个窗口中ref的GC含量百分比

Gatk AnnotateIntervals -L ${list} -R ${hg19.fa} --interval-merging-rule OVERLAPPING_ONLY -O ${out.tsv}

contig所在的位置，起止位置，GC复制比率

获取样本的read counts

统计bam文件在每个窗口的reads数，成对样本请处理normal 和tumor

Gatk CollectReadCounts -L ${out.list} -R ${hg19.fa} -I ${normal.MarkDup.sorted.bam} --format HDF5 --interval-merging-rule OVERLAPPING_ONLY -O ${normal.HDF5}

创建正常样本的CNV

Gatk CreateReadCountPanelOfNormals -I ${normal.MarkDup.sorted.bam} --minimum-interval-median-percentile 5.0 -O ${normal.HDF5}

降噪DenoiseReadCounts

进行标准化和降噪处理，去除掉创建normal时引入的噪音

Gatk DenoiseReadCounts -I ${tumor.HDF5} --count-panel-of-normals ${normal.HDF5} --standardized-copy-ratios ${normal.standard.tsv} --denoised-copy-ratios ${normal.Denoised.tsv}

计算单倍体的拷贝数，用于下面的拷贝数变异统计

Gatk CollectAllelicCounts -L ${out.list} -R ${hg19.fa} -I ${normal.MarkDup.sorted.bam} -O ${normal.AllelicCounts.tsv}

ModelSegments

将结果信息存放在file文件夹下，文件的名字以name为字段的一部分

Gatk ModelSegments --denoised-copy-ratios ${normal.Denoised.tsv} --allelic-counts ${tumor.AllelicCounts.tsv}  --normal-allelic-counts ${normal.AllelicCounts.tsv} -O ${file} --output-prefix ${name}

获取拷贝数变异结果

Gatk CallCopyRatioSegments -I ${name.cr.seg} -O ${output.cr.called.seg}

结果中的name.igv.seg文件可以使用igv进行查看

【1】contig所在的染色体 【2】【3】【4】【5】平均拷贝的比率，大于 0.14 就是扩增，小于 -0.15 就是缺失，其他的为正常【6】检测结果，+、-、0，分别代表这增加，减少，正常

第三部分：过滤质控

1. 硬过滤标准

对文件进行过滤分析，按照深度对vcf进行过滤，剔除低深度的变异

Gatk FilterVcf --MIN_DP ${int} -I ${output.vcf} -O ${output.filter.vcf}

低于该深度的会在info中进行相应的展示

过滤掉低深度的变异信息

grep -v "LowDP" ${output.filter.vcf} > ${result.vcf}

按照变异类型将文件中的SNP、INDEL进行筛选 ，将文件拆分为两种类型的变异,拆分过程可能会丢失数据，请在原数据中进行查找

Gatk SelectVariants -select-type ${type} -V ${result.vcf} –O ${result.type.vcf}

不同的变异检测hard filter标准不同

Gatk VariantFiltration –V ${result.type.vcf} --filter-expression “QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0” --filter-name ${filter} –O ${result.type.filter.vcf}

过滤后的vcf文件重新整理为新的vcf

Gatk MergeVcfs -I ${result.snp.filter.vcf} -I ${result.indel.filter.vcf} -O ${result.last.vcf}

2. 软过滤标准

该过滤过程是对结果文件通过机器学习的方式，需多种变异数据库的变异检测结果

构建重新校准模型，为筛选目的对变体质量进行评分

 Gatk VariantRecalibrator -R ${ref} -V ${output.vcf} --resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap_3.3.hg38.sites.vcf.gz -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -mode ${type} -O ${output.recal} --tranches-file ${output.tranches} 

训练集名字，这个名字是可以随便改动的
known：该数据是否作为已知变异数据集，用于对变异数据的标注；
training：该数据是否作为模型训练的数据集，用于训练VQSR模型；
truth：该数据是否作为验证模型训练的真集数据，这个数据同时还是VQSR训练bad model时自动进行参数选择的重要数据；
prior：该数据集在VQSR模型训练中的权重，或者叫Prior likelihood（这里转化为Phred-scale，比如20代表的是0.99）

Gatk ApplyVQSR -R ${ref} -V ${output.vcf} --truth-sensitivity-filter-level 99.0 --tranches-file output.tranches --recal-file ${output.recal} -mode SNP -O ${output.vcf.gz}

附录

由下机数据的fastq格式文件到比对后的sam，转换格式之后的bam文件，call变异之后的VCF文件，各种格式文件的内容及所便是的含义如下

fastq文件的格式

第一行：必须以“@”开头，后面跟着唯一的序列ID标识符，然后跟着可选的序列描述内容，标识符与描述内容用空格分开。
第二行：序列字符。
第三行：必须以“+”开头，后面跟着可选的ID标识符和可选的描述内容，如果“+”后面有内容，该内容必须与第一行“@”后的内容相同。
第四行：碱基质量字符，每个字符对应第二行相应位置碱基或氨基酸的质量，该字符可以按一定规则转换为碱基质量得分，碱基质量得分可以反映该碱基的错误率。这行的字符数与第二行中的字符数必须相同。

sam文件的格式

https://blog.csdn.net/genome_denovo/article/details/78712972

bam文件的格式

共计11个字符段
【1】reads的ID标识符 【2】标记数字 【3】比对到的染色体 【4】 比对到的染色体的位置 【5】比对的质量值 【6】比对结果的CIGAR字符串 M：匹配，I：插入，D：删除 【7】 "="表示正确匹配到序列上、"X"表示错误匹配到序列上 【8】mate在序列上的名称【9】mate在序列上的位置 【10】reads的序列 11.reads序列的碱基质量 【11】AS:i 匹配的得分、XS:i 第二好的匹配的得分、S:i mate 序列匹配的得分、XN:i 在参考序列上模糊碱基的个数、XM:i 错配的个数、XO:i gap open的个数、XG:i gap 延伸的个数、NM:i 经过编辑的序列、YF:i 说明为什么这个序列被过滤的字符串、YT:Z、MD:Z? 代表序列和参考序列错配的字符串

vcf文件的格式

【1】染色体位置 【2】变异所在位置 【3】variant的ID 【4】变异在ref上的信息 【5】突变后的情况，出现染色体号的为染色体的倒置【6】突变后的质量值 【7】按照突变质量的过滤情况 【8】详情信息 【9】格式 【10】样本名