GATK4测试学习
GATK代表GenomeAnalysisToolkit。 它是用于分析高通量测序数据的命令行工具的集合，主要侧重于突变体发现。
- 这些工具可以单独使用，也可以链接成完整的工作流程。
- 从4.0版开始，GATK包含Picard工具包。
Contents
- Quick start for the impatient
- Get GATK
- Install it
- Test that it works
- Run GATK and Picard commands
快速入门
- 在Linux或MacOSX上运行； 不支持MS Windows
- 确保使用的是Java 8 / JDK 1.8
- 下载GATK软件包，或获取Docker映像（不知道的东西太多了！之后学习下Docker映像）
- 由于一些原因有两个jar文件，通过gatk包装器脚本调用GATK，而不是直接调用jar
- 基本语法是gatk [--java-options“ -Xmx4G”] ToolName [GATK args]
- 使用Terra工作空间来测试管道并了解每种工具的功能(并不懂这句)
获取GATK
下载并解压缩软件包（名为gatk- [version]）后，在结果目录中找到文件：
gatk
gatk-completion.sh
gatkcondaenv.yml
GATKConfig.EXAMPLE.properties
gatkdoc
gatk-package-4.1.3.0-local.jar
gatk-package-4.1.3.0-spark.jar
gatkPythonPackageArchive.zip
README.md
scripts
- gatk-package- [version] -spark.jar是在Spark集群上运行Spark工具的jar
- gatk-package- [version] -local.jar是用于其他所有项目的jar（ 包括在本地（即在常规服务器或群集上）运行Spark工具
- 每次执行不需要指定运行哪一个，gatk文件是调用的可执行包装脚本，它将根据命令行的其余部分选择适当的jar。 如果需要，仍然可以调用特定的jar，但是使用gatk会更容易，并且它还将帮助设置一些必须手动指定的参数。
安装GATK
不需要安装，添加到环境变量就可以用啦！确保在路径中包含最后一个“/”！！！
> PATH="/path/to/gatk-package/:$PATH"
#测试它是否有效
> gatk --help
运行GATK和Picard命令
在“工具索引”部分列出并详细描述了可用的工具以及可用的选项。 调用任何GATK或Picard工具的基本语法如下：
> gatk [--java-options "jvm args like -Xmx4G go here"] ToolName [GATK args go here]
# for example, a simple GATK command would look like:
> gatk --java-options "-Xmx8G" HaplotypeCaller -R reference.fasta -I input.bam -O output.vcf
在此处找到有关GATK命令行语法的更多信息。
Picard工具的语法
- 在GATK中使用时，所有Picard工具都使用与GATK相同的语法。
- I = input.bam 转换为执行-I input.bam。
- 因此，例如，一个简单的Picard命令如下所示：
> gatk ValidateSamFile -I input.bam -MODE SUMMARY
😀 😃 😄 😁 😆 😅 😂 🤣 😊 😇 🙂 分割线 🙃 😉 😌 😍 😘 😗 😋 😛 😝 😜
Base Quality Score Recalibration (BQSR)
BQSR代表“基本质量得分重新校准”。 简而言之，这是一个数据预处理步骤，它在估计每个碱基检出的准确性时检测测序机产生的系统性错误。
注意，此基本重新校准过程（BQSR）不应与变量重新校准（VQSR）混淆，后者是一种复杂的过滤技术。
Contents
- Overview
- Base recalibration procedure details
- Important factors for successful recalibration
- Examples of pre- and post-recalibration metrics
- Recalibration report
概述
基本质量得分是测序机发出的每个碱基的错误估计值； 他们表示机器对每次调用正确的基准的信心如何。 例如，假设机器读取了一个A核苷酸，并分配了质量分数Q20（以Phred表示），这意味着它有99％的把握可以正确识别碱基。 这看似很高，但这确实意味着我们可以预期在100个案例中有1个是错误的。 因此，如果我们有数十亿个碱基调用（我们在30倍的基因组中获得约900亿个碱基），那么按照该速率，机器将在9亿个碱基中进行错误的调用-这是很多错误的碱基。
为什么这有关系？ 因为我们的简短变体调用算法在很大程度上依赖于在每个读取序列中分配给各个碱基调用的质量得分。 这是因为质量得分告诉我们，我们可以信任该特定观察多少信息，以告知我们该碱基所在位置的生物学真相。 如果我们有一个质量得分较低的碱基检出，则意味着我们不确定我们是否正确地读取了A，也可能是其他原因。 因此，我们不会像其他具有更高质量的基本调用那样信任它。 换句话说，我们使用该分数来权衡我们对特定位点存在的变异等位基因的支持或反对的证据。
基本的重新校准过程涉及两个关键步骤：首先，BaseRecalibrator工具基于输入数据和一组已知变量构建协变量模型，从而生成重新校准文件。 然后ApplyBQSR工具会根据模型调整数据中的基本质量得分，从而生成一个新的BAM文件。 已知变体用于掩盖真实（预期）变异位点的碱基，以避免将真实变异算作错误。 在被屏蔽的站点之外，每个不匹配项都计为一个错误。 其余大部分是会计。有一个可选的但强烈建议的步骤，其中涉及建立第二个模型并生成绘图前后的图表，以可视化重新校准过程的效果。 这对于质量控制很有用。
基本重新校准程序详细信息
BaseRecalibrator建立模型
要构建重新校准模型，第一个工具将检查输入BAM文件中的所有读数，并以表格形式列出有关以下碱基特征的数据：
- 读取所属的读取组
- 机器报告的质量得分
- 机器周期产生此基数（第N个周期=从读取开始起第N个基数）
- 当前碱基+先前碱基（二核苷酸）
根据已知的变异资源(通常是 dbSNP) ，我们计算每个 bin 中的碱基数，以及这些碱基与参考碱基不匹配的频率，但不包括已知在人群中发生变异的位点。这些信息以 GATKReport 格式输出到一个重新校准的文件中。
ApplyBQSR调整分数
第二个工具再次执行所有读取，使用重新校准文件根据每个基地落入的箱子调整每个基地的分数。 因此，新的质量分数实际上是：
- 报告质量分数和经验质量之间的全局差值之和。
- 加上质量箱特定班次。
- 加上环x等效和二核苷酸x等效。
经过重新校正后，阅读质量分数比以前更接近他们的经验分数。 这意味着它们可以以统计上可靠的方式用于下游处理，例如变量调用。 此外，通过考虑循环和序列上下文的质量变化，我们可以在读数中识别真正高质量的碱基，经常找到碱基的子集，即使在最初没有碱基被标记为Q30的情况下也是如此。
重新校准成功的重要因素
Read groups
重新校准系统支持读取组，这意味着它使用@RG标记按读取组对数据进行分区。这允许它对每个读组执行重新校准，这反映读属于哪个库，以及它在flowcell上的哪个lane上排序。我们知道，系统偏差可能发生在一条车道上，但不会发生在另一条车道上，或者一个库可能发生在另一个库上，因此能够在每个序列数据单元内重新校准，可以使建模过程更加准确。作为一个必然结果，这意味着可以对具有多个读组的BAM文件运行BQSR。但是，请注意，内存需求随文件中的读组数量呈线性增长，因此具有许多读组的文件可能需要大量RAM来存储所有协变量数据。
Amount of data
重新校准质量的一个关键决定因素是观察到的碱基的数量和每个容器中的不匹配。这个过程在少量对齐读取时不能很好地工作。我们通常预计每个读组有超过1亿个碱基;根据经验，数字越大越好。
No excuses
应该总是对排序数据执行重新校准。在人类数据中，考虑到拥有的详尽的变异数据库，几乎所有剩余的不匹配项（甚至在癌症中）都将是错误，所以为你的数据确定一个准确的错误模型超级容易，这对下游分析是至关重要的。
你真的需要跳过BQSR的主要例子是当你有太少的数据（一些小的基因panel有这个问题），或者你正在处理一个非常奇怪的有机体，它显示出疯狂的变异量。
重新校准报告
重新校准报告包含以下5个表格：
- Arguments Table -- a table with all the arguments and its values
该表包含用于运行此数据集的BQSR的所有参数。 - Quantization Table
GATK量化程序使用一种统计方法来确定最佳合并系统，该系统将通过合并特定数据集中存在的不同质量而导致的误差最小化。 运行BQSR时，将生成一个包含每种基本质量的基本计数的表，并生成一个“默认”量化表。 如果要量化质量得分，此表是GATK中任何其他工具的必需参数。 - ReadGroup Table
此表包含每个阅读组不匹配插入和删除的经验质量分数。 - Quality Score Table
该表包含每个阅读组的经验质量得分和原始质量得分，用于错配的插入和缺失。 - Covariates Table
该表具有数据集中使用的每个协变量的经验性质。 默认协变量是循环和上下文。 在当前实现中，上下文具有固定大小（默认为6）。 每个上下文和每个周期在此表上都有一个按阅读组和原始质量得分分层的条目。
分割线分割线分割线分割线分割线分割线分割线分割线分割线
GATK-4.1.3.0全基因组和全外显子组分析实战
[参考地址]https://cloud.tencent.com/developer/news/180158
注意点：
- 取消了RealignerTargetCreator、IndelRealigner，HaplotypeCaller继承了这部分功能。
- GATK4使用了新的设计模式，做了很多功能的整合，picard已经完全整合。
1.原始数据质控
2.数据预处理
序列比对
bwa mem -t 4 -R '@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name' /path/to/human.fasta read_1.fq.gz read_2.fq.gz | samtools view -S -b - > sample_name.bam
- 2.fq是非强制性的，只有是双末端测序的数据时才需要添加。
- -t：线程数，我们在这里使用4个线程；
- -R：Read Group的字符串信息，以@RG开头，用来将比对的read进行分组。不同的组之间测序过程被认为是相互独立的，这个信息对于后续对比对数据进行错误率分析和Mark duplicate时非常重要；
- BWA MEM比对模块是有一定适用范围的：它是专门为长read比对设计的，目的是为了解决，第三代测序技术这种能够产生长达几十kb甚至几Mbp的read情况。一般只有当read长度≥70bp的时候，才推荐使用，如果比这个要小，建议使用BWA ALN模块。
排序
gatk SortSam -I data.bam -O NGS171171.sort.bam -SO coordinate
- 也有用samtools进行排序的，但是我直接用了gatk里的。
去除重复序列（或者标记重复序列）
gatk MarkDuplicates -I NGS171171.sort.bam -O NGS171171.sort.dedup.bam -M NGS171171_metrics.txt
# 创建索引文件
gatk BuildBamIndex -I NGS171171.sort.dedup.bam
# 同个样本的所有比对结果合并成一个大的BAM文件
# samtools merge <out.bam> <in1.bam> [<in2.bam> ... <inN.bam>]
- Picard整合在gatk4中，简化成直接输入Picard工具名MarkDuplicates就行！
重新校正碱基质量值（BQSR）
# 下载参考文件
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/dbsnp_138.hg19.vcf.gz &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/dbsnp_138.hg19.vcf.idx.gz &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.idx.gz &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/ucsc.hg19.fasta.gz &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/ucsc.hg19.fasta.fai.gz &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/ucsc.hg19.dict.gz &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.idx.gz &
- 第一步，BaseRecalibrator，这里计算出了所有需要进行重校正的read和特征值，然后把这些信息输出为一份校准表文件。（sample_name.recal_data.table）
gatk BaseRecalibrator -I NGS171171.sort.dedup.bam -R ucsc.hg19.fasta --known-sites 1000G_phase1.snps.high_confidence.hg19.sites.vcf --known-sites Mills_and_1000G_gold_standard.indels.hg19.sites.vcf -O NGS171171.sort.dedup.recal_data.table --disable-sequence-dictionary-validation
- 第二步，ApplyBQSR，这一步利用第一步得到的校准表文件（wes.recal_data.table）重新调整原来BAM文件中的碱基质量值，并使用这个新的质量值重新输出一份新的BAM文件。
- 注意，这一步gatk4与3是不一样的，原先为PrintReads。
gatk ApplyBQSR -R ucsc.hg19.fasta -I NGS171171.sort.dedup.bam -bqsr NGS171171.sort.dedup.recal_data.table -O NGS171171.sort.dedup.BQSR.bam --disable-sequence-dictionary-validation
- 变异检测前确定bam文件是否符合GATK要求。如果显示 no error，则可以用HaplotypeCaller call SNP/Indel。
gatk ValidateSamFile -I NGS171171.sort.dedup.BQSR.bam
- 建立索引。
samtools index NGS171171.sort.dedup.BQSR.bam NGS171171.sort.dedup.BQSR.bai
samtools mpileup -B -q 30 -Q 30 -f ucsc.hg19.fasta NGS171171.sort.dedup.BQSR.bam > NGS171171.sort.dedup.BQSR.bam.mpileup
- 统计去除PCR重复之后的mapping率和重复率、coverage。
> samtools flagstat NGS171171.sort.dedup.BQSR.bam > NGS171171.sort.dedup.BQSR.bam.DataStat
> bedtools2/bin/genomeCoverageBed -i NGS171171.sort.dedup.BQSR.bam -bga -max 100 > NGS171171.sort.dedup.BQSR.bam.coverage
> perl $Perlpath/statWGScoverage.pl -i NGS171171.sort.dedup.BQSR.bam.coverage -o NGS171171.sort.dedup.BQSR.bam.coverage.stat
HaplotypeCaller突变检测
有两种做法，区别在于是否生成中间文件gVCF
#（1）直接进行HaplotypeCaller，这适合于单样本，只执行一次HaplotypeCaller。
> nohup gatk HaplotypeCaller -R ucsc.hg19.fasta -I NGS171171.sort.dedup.BQSR.bam -D dbsnp_138.hg19.vcf -O NGS171171.g.vcf --disable-sequence-dictionary-validation &
# 如果有多个样本那么可以继续用-I参数输入-I sample1.bam [-I sample2.bam ...]
# 为了提高效率我们可以按照染色体一条条来独立执行这个步骤，最后再把结果合并起来就好了。但其实完整的算也挺快哈哈～
# (2）每个样本先各自生成gVCF，然后再进行群体joint-genotype。
gVCF全称是genome VCF，是每个样本用于变异检测的中间文件，格式类似于VCF，它把joint-genotype过程中所需的所有信息都记录在这里面，文件无论是大小还是数据量都远远小于原来的BAM文件。这样一旦新增加样本也不需要再重新去读取所有人的BAM文件了，只需为新样本生成一份gVCF，然后重新执行这个joint-genotype就行了。我用的第一种方法～
变异质控和过滤
质控的含义和目的是指通过一定的标准，最大可能地剔除假阳性的结果，并尽可能地保留最多的正确数据。
第一种方法
GATK VQSR，它通过机器学习的方法利用多个不同的数据特征训练一个模型（高斯混合模型）对变异数据进行质控，使用VQSR需要具备以下两个条件：
- 第一，需要一个精心准备的已知变异集，它将作为训练质控模型的真集。比如，Hapmap、OMNI，1000G和dbsnp等这些国际性项目的数据，这些可以作为高质量的已知变异集。
- 第二，要求新检测的结果中有足够多的变异，不然VQSR在进行模型训练的时候会因为可用的变异位点数目不足而无法进行。适合全基因组分析。
- 此方法要求新检测的结果中有足够多的变异，不然VQSR在进行模型训练的时候会因为可用的变异位点数目不足而无法进行。可能很多非人的物种在完成变异检测之后没法使用GATK VQSR的方法进行质控，一些小panel、外显子测序，由于最后的变异位点不够，也无法使用VQSR。全基因组分析或多个样本的全外显子组分析适合用此方法。
第二种方法 - 通过过滤指标过滤。
突变注释
MuTect
v4.1开始，不再需要使用-tumor指定肿瘤样本名称。 如果包括普通样品，则只需用-normal指定普通样品名称。
gatk Mutect2 -R ucsc.hg19.fasta -I NGS171171.sort.dedup.BQSR.bam -O NGS171171.vcf.gz --disable-sequence-dictionary-validation
