GATK变异检测

大家好，欢迎来到IT知识分享网。

变异检测

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测序文件/ _R1.fq $tumor2.fq = ~/GATK_passway/Illumina测序文件/ _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的使用流程：

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

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} 

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} 

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的流程图

构建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} 

获取样本的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} 

第三部分：过滤质控

1. 硬过滤标准

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

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

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} 

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/ 

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】样本名