title: bioinfo-GATK_Call_SNP index_img: /img/bioinformatics.png banner_img: /img/bioinformatics.png tags:
- Bioinformatics categories:
- Bioinformatics
- Pipelines_RNA-seq date: 2020-05-02 16:34:14
参考：
- https://blog.csdn.net/ccArtermices/article/details/101023390
- https://mp.weixin.qq.com/s?__biz=MzAxOTUxOTM0Nw==&mid=2649798425&idx=1&sn=ae355ed362848578e5c853413f23dfd7&chksm=83c1d505b4b65c13124c9acd210356c4364ec9f5498bbd16fa4475be29811213abb64ea9720f&scene=21#wechat_redirect
- https://mp.weixin.qq.com/s?__biz=MzAxOTUxOTM0Nw==&mid=2649798455&idx=1&sn=67a7407980a57ce138948eb46992b603&chksm=83c1d52bb4b65c3dde31df94e9686654bf616166c7311b531213ebf0010f67a32ce827e677b1&scene=21#wechat_redirect
- SNP是指在基因组上单个核苷酸的变异,包括置换、颠换、缺失和插入。
- SNP在基因组中分布相当广泛，近来的研究表明在很多物种基因组中每300bp就出现一次。
- 大量存在的SNP位点，使人们有机会发现与各种疾病，包括肿瘤相关的基因组突变。
- 既然SNP那么广泛存在，获得SNP就变成一项重要的任务。在经历了样本收集、测序、质控和mapping后，我们输出了bam格式的数据。之后，我们就要尝试利用GATK call SNP了。
安装GATK4、JAVA、Python、R等环境
配置java环境
# java -version # 如果是1.8就跳过这一步骤，不过比这个版本低的话就继续阅读本步骤的内容。 conda install -c cyclus java-jdk=8.45
GATK要求安装python 3.6或以上，我这里直接用的conda python=3.7的环境，不再赘述了
配置R，最新版本即可
conda install -c pkgs/r r=3.6
# 我这里用的3.6版本
- 安装GATK4
由于conda中的GATK没有4的版本
这次选择安装编译好的GATK4直接解压运行，从github的仓库，https://github.com/broadinstitute/gatk/releases/ 选择最新版本下载。
我这里下载的是4.1.7.0版本，unzip之后运行gatk即可，或者选择带local字样的jar文件，这样子我们使用java -jar的时候方便指定内存和调用多核，gatk可以加到环境变量中方便使用。
配置picard https://github.com/broadinstitute/picard/tree/2.22.4
git clone https://github.com/broadinstitute/picard.git cd picard/ ./gradlew shadowJar
其实最新的GATK4已经将picard整合进去了，可以很方便的调用，具体可以查看GATK4的官方文档
GATK4 call SNP流程
使用snakemake流程
# _*_ coding: UTF-8 _*_
########################################################################
# ZHAO Huanan
# 2020-05-02
# 293T RNASeq analysis pipeline
########################################################################
# run on abyss
# /home/zhaohuanan/zhaohn_HD/3.project/2.repeat_David_Liu_2020_NBT/2.GATK_CALL_SNP
# --------------------------------------------------------------->>>>>>>
# pipeline
# --------------------------------------------------------------->>>>>>>
# 1. STAR alingment √
# 2. picard change RG √
# 3. picard mark duplicates
# 4. SplitNCigarReads??
# 5. call SNP
# 6. hard filter of SNP
# --------------------------------------------------------------->>>>>>>
# software
# --------------------------------------------------------------->>>>>>>
PICARD = "/gpfs/user/zhaohuanan/1.apps/picard/picard.jar"
GATK4 = "/home/zhaohuanan/zhaohn_HD/1.apps/GATK4/gatk-4.1.7.0/gatk-package-4.1.7.0-local.jar"
HG38_FA = "/home/zhaohuanan/zhaohn_HD/2.database/bwa_hg38/hg38_only_chromosome.fa"
JAVA = "/home/zhaohuanan/zhaohn_HD/1.apps/anaconda3/envs/py37/bin/java"
# --------------------------------------------------------------->>>>>>>
# vars
# --------------------------------------------------------------->>>>>>>
SAMPLES = [
"batch1_293T-RNASeq-KJY-All-EMX1_rep1",
"batch1_293T-RNASeq-KJY-All-EMX1_rep2",
"batch1_293T-RNASeq-KJY-APODel_rep1",
"batch1_293T-RNASeq-KJY-APODel_rep2",
"batch2_293T-RNASeq-BE4-All_rep1",
"batch2_293T-RNASeq-Mock_rep1",
"batch2_293T-RNASeq-Vector_rep1"
]
##############################
# test
##############################
# SAMPLES = [
# 'test'
# ]
##############################
##############################
# all
##############################
rule all:
input:
expand("bam/{sample}_Aligned_sort_MarkDup.bam",sample=SAMPLES),
expand("bam/{sample}_Aligned_sort_MarkDup.bam.bai",sample=SAMPLES),
expand("vcf/{sample}_gatk4_HaplotypeCaller_raw_variants.vcf",sample=SAMPLES),
expand("vcf/{sample}_gatk4_HaplotypeCaller_filtered_variants_merge.vcf",sample=SAMPLES)
##############################
# cut adapt
##############################
##############################
# mapping
##############################
##############################
# add RG tag
##############################
# rule add_RG_tag:
# input:
# "0.sort_bam/{sample}.bam"
# output:
# "0.sort_bam/{sample}_Aligned.out.fix_RG.bam"
# params:
# tag = "'@RG\\tID:{sample}\\tSM:{sample}\\tPL:ILLUMINA'"
# shell:
# """
# srun -T 24 -c 24 \
# samtools addreplacerg -@ 24 \
# -r {params.tag} \
# -O BAM -o {output} \
# --reference {HG38_FA} \
# {input}
# """
##############################
# sort bam
##############################
# rule BAM_sort_by_position:
# input:
# "0.sort_bam/{sample}_Aligned.out.fix_RG.bam"
# output:
# "0.sort_bam/{sample}_Aligned_sort.bam"
# shell:
# """
# srun -T 24 -c 24 \
# samtools sort -@ 24 \
# -O BAM -o {output} \
# -T {output}.temp \
# -m 4G \
# {input}
# """
##############################
# mark duplicate
##############################
rule BAM_mark_duplicate:
input:
"0.sort_bam/{sample}_Aligned_sort.bam"
output:
"bam/{sample}_Aligned_sort_MarkDup.bam",
"bam/{sample}_Aligned_sort_MarkDup.matrix"
log:
"bam/{sample}_Aligned_sort_MarkDup.log"
shell:
"""
srun -T 24 -c 24 \
{JAVA} -Xms90g -Xmx90g -XX:ParallelGCThreads=24 \
-jar {PICARD} \
MarkDuplicates \
I={input} \
O={output[0]} \
M={output[1]} \
VALIDATION_STRINGENCY=SILENT \
ASO=queryname \
CREATE_MD5_FILE=true 2>{log}
"""
##############################
# samtools bam index -> bai
##############################
rule BAM_index:
input:
"bam/{sample}_Aligned_sort_MarkDup.bam"
output:
"bam/{sample}_Aligned_sort_MarkDup.bam.bai"
shell:
"""
srun -T 24 -c 24 \
samtools index -@ 24 \
{input} \
{output}
"""
##############################
# HaplotypeCaller -> g.vcf
##############################
rule GATK_HaplotypeCaller:
input:
"bam/{sample}_Aligned_sort_MarkDup.bam"
output:
"vcf/{sample}_gatk4_HaplotypeCaller_raw_variants.g.vcf"
log:
"vcf/{sample}_gatk4_HaplotypeCaller_raw_variants.g.vcf.log"
shell:
"""
srun -T 24 -c 24 \
{JAVA} -Xms90g -Xmx90g -XX:ParallelGCThreads=24 \
-jar {GATK4} \
HaplotypeCaller \
-R {HG38_FA} \
--emit-ref-confidence GVCF \
-I {input} \
-O {output} 2>{log}
"""
##############################
# GenotypeGVCFs -> raw variants
##############################
rule GATK_GenotypeGVCFs:
input:
"vcf/{sample}_gatk4_HaplotypeCaller_raw_variants.g.vcf"
output:
"vcf/{sample}_gatk4_HaplotypeCaller_raw_variants.vcf"
log:
"vcf/{sample}_gatk4_HaplotypeCaller_raw_variants.log"
shell:
"""
srun -T 24 -c 24 \
{JAVA} -Xms90g -Xmx90g -XX:ParallelGCThreads=24 \
-jar {GATK4} \
GenotypeGVCFs \
-R {HG38_FA} \
-V {input} \
-O {output} 2>{log}
"""
##############################
# SelectVariants SNP
##############################
rule GATK_SelectVariants_SNP:
input:
"vcf/{sample}_gatk4_HaplotypeCaller_raw_variants.vcf"
output:
"vcf/{sample}_gatk4_HaplotypeCaller_select_variants_SNP.vcf"
log:
"vcf/{sample}_gatk4_HaplotypeCaller_select_variants_SNP.log"
shell:
"""
srun -T 24 -c 24 \
{JAVA} -Xms90g -Xmx90g -XX:ParallelGCThreads=24 \
-jar {GATK4} \
SelectVariants \
-select-type SNP \
-V {input} \
-O {output} 2>{log}
"""
##############################
# VariantFiltration SNP
##############################
rule GATK_VariantFiltration_SNP:
input:
"vcf/{sample}_gatk4_HaplotypeCaller_select_variants_SNP.vcf"
output:
"vcf/{sample}_gatk4_HaplotypeCaller_filtered_variants_SNP.vcf"
log:
"vcf/{sample}_gatk4_HaplotypeCaller_filtered_variants_SNP.log"
shell:
"""
srun -T 24 -c 24 \
{JAVA} -Xms90g -Xmx90g -XX:ParallelGCThreads=24 \
-jar {GATK4} \
VariantFiltration \
--filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
--filter-name "PASS" \
-V {input} \
-O {output} 2>{log}
"""
##############################
# SelectVariants INDEL
##############################
rule GATK_SelectVariants_INDEL:
input:
"vcf/{sample}_gatk4_HaplotypeCaller_raw_variants.vcf"
output:
"vcf/{sample}_gatk4_HaplotypeCaller_select_variants_INDEL.vcf"
log:
"vcf/{sample}_gatk4_HaplotypeCaller_select_variants_INDEL.log"
shell:
"""
srun -T 24 -c 24 \
{JAVA} -Xms90g -Xmx90g -XX:ParallelGCThreads=24 \
-jar {GATK4} \
SelectVariants \
-select-type INDEL \
-V {input} \
-O {output} 2>{log}
"""
##############################
# VariantFiltration INDEL
##############################
rule GATK_VariantFiltration_INDEL:
input:
"vcf/{sample}_gatk4_HaplotypeCaller_select_variants_INDEL.vcf"
output:
"vcf/{sample}_gatk4_HaplotypeCaller_filtered_variants_INDEL.vcf"
log:
"vcf/{sample}_gatk4_HaplotypeCaller_filtered_variants_INDEL.log"
shell:
"""
srun -T 24 -c 24 \
{JAVA} -Xms90g -Xmx90g -XX:ParallelGCThreads=24 \
-jar {GATK4} \
VariantFiltration \
--filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
--filter-name "PASS" \
-V {input} \
-O {output} 2>{log}
"""
##############################
# MergeVcfs -> vcf
##############################
rule GATK_MergeVcfs:
input:
"vcf/{sample}_gatk4_HaplotypeCaller_filtered_variants_SNP.vcf",
"vcf/{sample}_gatk4_HaplotypeCaller_filtered_variants_INDEL.vcf"
output:
"vcf/{sample}_gatk4_HaplotypeCaller_filtered_variants_merge.vcf"
log:
"vcf/{sample}_gatk4_HaplotypeCaller_filtered_variants_merge.log"
shell:
"""
srun -T 24 -c 24 \
{JAVA} -Xms90g -Xmx90g -XX:ParallelGCThreads=24 \
-jar {GATK4} \
MergeVcfs \
-I {input[0]} \
-I {input[1]} \
-O {output} 2>{log}
"""
上述流程的细节
目前GATK4.0的正式版本已经发布，它的使用方式与之前相比有着一些差异（变得更加简单，功能也更加丰富了），增加了结构性变异检测和很多Spark、Cloud-Only的功能，并集成了MuTect2和picard的所有功能（以及其他很多有用的工具），这为我们减少了许多额外的工具，更加有利于流程的构建和维护，4.0之后的GATK是一个新的篇章，大家最好是掌握这一个版本！
### 质控
质控是必须做的，我们需要完整认识原始的测序数据质量到底如何，该步骤不能省略。（WGS系列第三节）
### 比对
首先是比对。所谓比对就是把测序数据定位到参考基因组上，确定每一个read在基因组中的位置。这里，我们依然用目前使用最广的BWA来完成这个工作。在正式比对之前，需要先为参考序列构建BWA比对所需的FM-index（比对索引）。
使用bwa完成比对，用samtools完成BAM格式转换、排序并标记PCR重复序列
```shell # 参考，这个我用snakemake
# 简单解释一下这个shell所做的事情：首先利用bwa mem比对模块将E.coli K12质控后的测序数据定位到其参考基因组上（我们这里设置了4个线程来完成比对，根据电脑性能可以适当调大），同时通过管道（'|' 操作符）将比对数据流引到samtools转换为BAM格式（SAM的二进制压缩格式），然后重定向('>'操作符)输出到文件中保存下来。
#1 比对
time /Tools/common/bin/bwa mem -t 4 -R '@RG\tID:foo\tPL:illumina\tSM:E.coli_K12' /Project/201802_wgs_practice/input/E.coli/fasta/E.coli_K12_MG1655.fa /Project/201802_wgs_practice/input/E.coli/fastq/SRR1770413_1.fastq.gz /Project/201802_wgs_practice/input/E.coli/fastq/SRR1770413_2.fastq.gz | /Tools/common/bin/samtools view -Sb - > /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam && echo " bwa mapping done "
# 接着用samtools对原始的比对结果按照参考序列位置从小到大进行排序（同样是4个线程），只有这个步骤完成之后才可以继续往下。
#2 排序 time /Tools/common/bin/samtools sort -@ 4 -m 4G -O bam -o /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.bam /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam && echo "** BAM sort done"rm -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam
# 然后，我们使用GATK标记出排完序的数据中的PCR重复序列。这个步骤完成后，如无特殊需要，我们就可以直接删除前面那两个BAM文件了（原始比对结果和排序后的结果）——后续几乎不会再用到那两份文件了。关于标记PCR重复序列的操作比较简单，不再细说（如果希望了解更多有关重复序列特征的信息可以回看WGS系列第四节中的内容）。
#3 标记PCR重复 time /Tools/common/bin/gatk/4.0.1.2/gatk MarkDuplicates -I /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.bam -O /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup.bam -M /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup_metrics.txt && echo " markdup done "
#4 删除不必要文件(可选) rm -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam rm -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.bam
# 最后，我们再用samtools为E_coli_K12.sorted.markdup.bam创建索引。我认为不论是否有后续分析，为BAM文件创建索引应该作为一个常规步骤，它可以让我们快速地访问基因组上任意位置的比对情况，这一点非常有助于我们随时了解数据。
#5 创建比对索引文件 time /Tools/common/bin/samtools index /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup.bam && echo " index done "
# 至于每个步骤最前面的time，则是用于记录执行时间的，有助于我们清楚地知道每一个分析过程都花了多少时间，当需要优化流程的时候这个信息会很有用。
```
> -R 设置Read Group信息，虽然我在以前的文章中已经反复强调过它的重要性，但这里还是再说一次，它是read数据的组别标识，并且其中的ID，PL和SM信息在正式的项目中是不能缺少的(如果样本包含多个测序文库的话，LB信息也不要省略)，另外由于考虑到与GATK的兼容关系，PL（测序平台）信息不能随意指定，必须是：ILLUMINA，SLX，SOLEXA，SOLID，454，LS454，COMPLETE，PACBIO，IONTORRENT，CAPILLARY，HELICOS或UNKNOWN这12个中的一个。### 变异检测
接下来是用GATK完成变异检测。但在开始之前之前我们还需要先为E.coli K12的参考序列生成一个.dict文件，这可以通过调用CreateSequenceDictonary模块来完成(这是原来picard的功能)。
$ /Tools/common/bin/gatk/4.0.1.2/gatk CreateSequenceDictionary -R E.coli_K12_MG1655.fa -O E.coli_K12_MG1655.dict && echo "** dict done **"
唯一需要注意的是.dict文件的名字前缀需要和fasta的一样，并跟它在同一个路径下，这样GATK才能够找到。
OK，现在我们就可以进行变异检测了，同样使用GATK 4.0的HaplotypeCaller模块来完成。由于我们只有一个样本，要完成这个工作其实很简单，直接输入比对文件和参考序列就行了，但是考虑到实际的情况，我想告诉大家一个更好的方式（虽然这会多花些时间），就是：先为每个样本生成一个GVCF，然后再用GenotypeGVCFs对这些GVCF进行joint calling，如下 ，我把命令都写在gatk.sh中，并执行。
```shell
#1 生成中间文件gvcf
time /Tools/common/bin/gatk/4.0.1.2/gatk HaplotypeCaller \
-R /Project/201802_wgs_practice/input/E.coli/fasta/E.coli_K12_MG1655.fa \
--emit-ref-confidence GVCF \
-I /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup.bam \
-O /Project/201802_wgs_practice/output/E.coli/E_coli_K12.g.vcf && echo "** gvcf done **"
#2 通过gvcf检测变异
time /Tools/common/bin/gatk/4.0.1.2/gatk GenotypeGVCFs \
-R /Project/201802_wgs_practice/input/E.coli/fasta/E.coli_K12_MG1655.fa \
-V /Project/201802_wgs_practice/output/E.coli/E_coli_K12.g.vcf \
-O /Project/201802_wgs_practice/output/E.coli/E_coli_K12.vcf && echo "** vcf done **"
```
很快我们就获得了E.coli K12这个样本初步的变异结果——E_coli_K12.vcf。之所以非要分成两个步骤，是因为我想借此告诉大家，**变异检测不是一个样本的事情，有越多的同类样本放在一起joint calling结果将会越准确，而如果样本足够多的话，在低测序深度的情况下也同样可以获得完整并且准确的结果**，而这样的分步方式是应对多样本的好方法。
最后，我们用bgzip对这个VCF进行压缩，并用tabix为它构建索引，方便以后的分析。
```shell
#1 压缩
time /Tools/common/bin/bgzip -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.vcf
#2 构建tabix索引
time /Tools/common/bin/tabix -p vcf /Project/201802_wgs_practice/output/E.coli/E_coli_K12.vcf.gz
```
bgzip压缩完成之后，原来的VCF文件会被自动删除。
如果大家仔细看过WGS系列第四节的话，会发现我这里缺少了两个步骤：重比对和BQSR。没有执行BQSR是因为E.coli K12没有那些必须的known变异集（或者有但我没找到），所以无法进行；但没有重比对，则是因为我在GATK 4.0中没发现IndelRealigner这个功能，虽然我们使用GATK HaplotypeCaller或者Mutect2的话确实可以省略这个步骤，但如果是其他软件来进行变异检测那么该步骤依然十分重要，我目前不太清楚为何GATK 4.0没有将这个功能单独分离出来。
后面要谈到的就是变异的质控了。很遗憾我们这个E.coli K12的变异结果并不适合通过VQSR来进行过滤，原因上面也提到了一些，它不像人类的基因组数据，有着一套适合用来训练过滤模型的已知变异集（dbSNP，1000G，Hapmap和omini等）。其实这种情况有时候我们在工作中也会碰到，比如有些捕获测序（Panel测序数据，甚至外显子测序）的数据，由于它的区域较小，获得的变异也不多，导致最终没法满足VQSR进行模型训练时所需的最低变异数要求，那时你也不能通过这个方式协助变异质控。那么碰到这种情况的时候该怎么办？我将这部分的内容放在了下一篇文章中，在那里我们再来讨论这个问题。我也会告诉大家变异质控的基本逻辑，而不是简单罗列一个命令，同时也会再用NA12878这个人的数据来进一步告诉大家如何比较和评估变异结果。小结
至此，这个篇文章的上半部分就到此为止了。除了那些重要的内容之外，在上文中，你会看到我反复提到了创建“索引”这个事情，比如为fasta，为BAM，为VCF。我为什么非要反复强调这个事情不可呢？因为我发现许多初学者并不知道索引的作用，当被问到如何从巨大的比对文件或者变异文件中提取某个信息时，总是要走弯路——努力写程序去提取，既慢又费力，结果还不一定好，甚至有些有一定经验的同学也不知道使用bgzip和tabix的好处，因此我才反复在文章里提及。
GATK4 call SNP的质量控制
在上一篇文章中我已经用例子仔细跟大家分享了WGS从原始数据到变异数据（Fastq->VCF）的具体执行过程。那么，在这一篇文章里，我们就来好好谈谈后续非常重要的一个环节——也是本次实践分析的最后一个部分——变异的质控。
什么是质控？我们不妨先给它下个定义：
质控的含义和目的是指通过一定的标准，最大可能地剔除假阳性的结果，并尽可能地保留最多的正确数据。有了这么一个定义之后，我们也就能够更加清晰地知道接下来该做些什么了。
VQSR
在上次的文章里我已经说到在GATK HaplotypeCaller之后，首选的质控方案是GATK VQSR，它通过机器学习的方法利用多个不同的数据特征训练一个模型（高斯混合模型）对变异数据进行质控，然而不幸的是使用VQSR需要具备以下两个条件：
第一，需要一个精心准备的已知变异集，它将作为训练质控模型的真集。比如，对于我们人来说，就有Hapmap、OMNI，1000G和dbsnp等这些国际性项目的数据，这些可以作为高质量的已知变异集。GATK的bundle主要就是对这四个数据集做了精心的处理和选择，然后把它们作为VQSR时的真集位点。这里我强调一个地方：是真集的『位点』而不是真集的『数据』！还请大家多多注意。因为，VQSR并不是用这些变异集里的数据来训练的，而是用我们自己的变异数据。这个对于刚接WGS的同学来说特别容易搞混，不要因为VQSR中用了那四份变异数据，就以为是用它们的数据来训练模型。
实际上，这些已知变异集的意义是告诉我们群体中哪些位点存在着变异，如果在其他人的数据里能观察到落入这个集合中的变异位点，那么这些被已知集包括的变异就有很大的可能是正确的。也就是说，我们可以从数据中筛选出那些和真集『位点』相同的变异，把它们当作是真实的变异结果。接着，进行VQSR的时候，程序就可以用这个筛选出来的数据作为真集数据来训练，并构造模型了。
关于VQSR的内在原理，前不久在我的知识星球中我做过简单的回答（下图），这里就不展开了，感兴趣的同学看下图的内容基本上也是足够的，虽然不详细，但应该可以帮你建立一个关于VQSR的基本认识。对于希望深入理解算法细节的同学来说，我的建议是直接阅读GATK这一部分的代码。但要注意，类似这样的过滤算法实际上还可以用很多不同的机器学习算法来解决，比如SVM，或者用深度学习来构造这个质控模型也都是OK的。
第二，要求新检测的结果中有足够多的变异，不然VQSR在进行模型训练的时候会因为可用的变异位点数目不足而无法进行。
由于条件1的限制，会导致很多非人的物种在完成变异检测之后没法使用GATK VQSR的方法进行质控。而由于条件2，也常常导致一些小panel甚至外显子测序，由于最后的变异位点不够，也无法使用VQSR。这个时候，我们就不得不选择硬过滤的方式来质控了。
那什么叫做硬过滤呢？所谓硬过滤其实就是通过人为设定一个或者若干个指标阈值（也可以叫数据特征值），然后把所有不满足阈值的变异位点采用一刀切掉的方法。
那么如何执行硬过滤？首先，需要我们确定该用哪些指标来评价变异的好坏。这个非常重要，选择对了事半功倍，选得不合理，过滤的结果有时还不如不过滤的。如果把这个问题放在从前，我们需要做比较多的尝试才能确定一些合适的指标，但现在就方便很多了，可以直接使用GATK VQSR所用的指标——毕竟这些指标都是经过精挑细选的。我想这应该不难理解，既然VQSR就是用这些指标来训练质控模型的，那么它们就可以在一定程度上描述每个变异的质量，我们用这些指标设置对应的阈值来进行硬过滤也将是合理的。VQSR使用的数据指标有6个（这些指标都在VCF文件的INFO域中，如果不是GATK得到的变异，可能会有所不同，但知道它们的含义之后也是可以自己计算的），分别是：
QualByDepth（QD）
FisherStrand (FS)
StrandOddsRatio (SOR)
RMSMappingQuality (MQ)
MappingQualityRankSumTest (MQRankSum)
ReadPosRankSumTest (ReadPosRankSum)
指标有了，那么阈值应该设置为多少？下面我想先给出一个硬过滤的例子，然后再逐个来对其进行分析，以便大家能够更好地理解变异质控的思路。值得注意的是不同的数据，有不同的情况，它的阈值有时是不同的。不过不用担心，当你掌握了如何做的思路之后完全有能力根据具体的情况举一反三。
执行硬过滤
首先是硬过滤的例子，这个过程我都用最新的GATK来完成。GATK 4.0中有一个专门的VariantFiltration模块（继承自GATK 3.x），它可以很方便地帮我们完成这个事情。不过，过滤的时候，需要分SNP和Indel这两个不同的变异类型来进行，它们有些阈值是不同的，需要区别对待。在下面的例子里，我们还是用上一篇文章中最后得到的变异数据（E_coli_K12.vcf.gz）为例子，这是具体的执行命令：
# 使用SelectVariants，选出SNP
time /Tools/common/bin/gatk/4.0.1.2/gatk SelectVariants \
-select-type SNP \
-V ../output/E.coli/E_coli_K12.vcf.gz \
-O ../output/E.coli/E_coli_K12.snp.vcf.gz
# 为SNP作硬过滤
time /Tools/common/bin/gatk/4.0.1.2/gatk VariantFiltration \
-V ../output/E.coli/E_coli_K12.snp.vcf.gz \
--filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
--filter-name "PASS" \
-O ../output/E.coli/E_coli_K12.snp.filter.vcf.gz
# 使用SelectVariants，选出Indel
time /Tools/common/bin/gatk/4.0.1.2/gatk SelectVariants \
-select-type INDEL \
-V ../output/E.coli/E_coli_K12.vcf.gz \
-O ../output/E.coli/E_coli_K12.indel.vcf.gz
# 为Indel作过滤
time /Tools/common/bin/gatk/4.0.1.2/gatk VariantFiltration \
-V ../output/E.coli/E_coli_K12.indel.vcf.gz \
--filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
--filter-name "PASS" \
-O ../output/E.coli/E_coli_K12.indel.filter.vcf.gz
# 重新合并过滤后的SNP和Indel
time /Tools/common/bin/gatk/4.0.1.2/gatk MergeVcfs \
-I ../output/E.coli/E_coli_K12.snp.filter.vcf.gz \
-I ../output/E.coli/E_coli_K12.indel.filter.vcf.gz \
-O ../output/E.coli/E_coli_K12.filter.vcf.gz
# 删除无用中间文件
rm -f ../output/E.coli/E_coli_K12.snp.vcf.gz* ../output/E.coli/E_coli_K12.snp.filter.vcf.gz* ../output/E.coli/E_coli_K12.indel.vcf.gz* ../output/E.coli/E_coli_K12.indel.filter.vcf.gz*与上一篇文章的目录逻辑一样，我们把这些shell命令都写入到bin目录下一个名为“variant_filtration.sh”的文件中，然后运行它。最后，只要符合了上面任意一个阈值的变异都会被留下来。流程的最后，我们需要把分开质控的SNP和Indel结果重新合并在一起，然后再把那些不必要的中间文件删除掉。
在具体的项目中，你如果需要使用硬过滤的策略，这个例子中的参数可以作为参考，特别是对于高深度数据而言。接下来我结合GATK所提供的资料与大家分享如何理解这些指标以及得出这些阈值的思路。
如何理解硬过滤的指标和阈值的计算
考虑到SNP和Indel在判断指标和阈值方面的思路是一致的，因此没必要重复说。所以，下面我只以SNP为例子，告诉大家设定阈值的思路。强调一下，为了更具有通用价值，这些阈值是借用NA12878（来自GIAB）的高深度数据进行计算获得的，所以如果你的数据（或者物种）相对比较特殊（不是哺乳动物），那么就不建议直接套用了，但可以依照类似的思路去寻找新阈值。
QualByDepth（QD）
QD是变异质量值（Quality）除以覆盖深度（Depth）得到的比值。这里的变异质量值就是VCF中QUAL的值——用来衡量变异的可靠程度，这里的覆盖深度是这个位点上所有含有变异碱基的样本的覆盖深度之和，通俗一点说，就是这个值可以通过累加每个含变异的样本（GT为非0/0的样本）的覆盖深度（VCF中每个样本里面的DP）而得到。举个例子：
1 1429249 . C T 1044.77 . . GT:AD:DP:GQ:PL 0/1:48,15:63:99:311,0,1644 0/0:47,0:47:99:392,0,0 1/1:0,76:76:99:3010,228,0这个位点是1:1429249，VCF格式，但我把FILTER和INFO的信息省略了，它的变异质量值QUAL=1044.77。我们可以从中看到一共有三个样本，其中一个是杂合变异（GT=0/1），一个纯合的非变异（GT=0/0），最后一个是纯合的变异（GT=1/1）。每个样本的覆盖深度都在其各自的DP域上，分别是：63，47和76。按照定义，这个位点的QD值就应该等于质量值除以另外两个含有变异的样本的深度之和（排除中间GT=0/0这个不含变异的样本），也就是:
QD = 1044.77 / (63+76) = 7.516QD这个值描述的实际上就是单位深度的变异质量值，也可以理解为是对变异质量值的一个归一化，QD越高一般来说变异的可信度也越高。在质控的时候，相比于QUAL或者DP（深度）来说，QD是一个更加合理的值。因为我们知道，原始的变异质量值实际上与覆盖的read数目是密切相关的，深度越高的位点QUAL一般都是越高的，而任何一个测序数据，都不可避免地会存在局部深度不均的情况，如果直接使用QUAL或者DP都会很容易因为覆盖深度的差异而带来有偏的质控结果。
在上面『执行硬过滤』的例子里面，我们看到认为好的SNP变异，QD的值不能低于2，但问题是为什么是2，而不是3或者其它数值呢？
要回答这个问题，我们可以通过利用NA12878 VQSR质控之后的变异数据和原始的变异数据来进行比较，并把它说明白。
首先，我们可以先把所有变异位点的QD值都提取出来，然后画一个密度分布图（Y轴代表的是对应QD值所占总数的比例，而不是个数），看看QD值的总体分布情况（如下图，来自NA12878的数据）。
从这个图里，我们可以看到QD的范围主要集中在0~40之间。同时，我们可以明显地看到有两个峰值（QD=12和QD=32）。这两个峰所反映的恰恰是杂合变异和纯合变异的QD值所集中的地方。这里大家可以思考一下，哪一个是代表杂合变异的峰，哪一个是代表纯合变异的峰呢？
回答是，第一个峰（QD=12）代表杂合，而第二峰（QD=32）代表纯合，为什么呢？因为对于纯合变异来说，贡献于质量值的read是杂合变异的两倍，同样深度的情况下，QD会更大。对于大多数的高深度测序数据来说，QD的分布和上面的情况差不多，因此这个分布具有一定的代表性。
接着，我们同时画出VQSR之后所有可信变异（FILTER=Pass）和不可信变异的QD分布图，如下，浅绿色代表可信变异的QD分布图，浅红色代表不可信变异的QD分布图。
你可以看到，大多数Fail的变异，都集中在左边的低QD区域，而且红波峰恰好是QD=2的地方，这就是为什么硬过滤时设置QD>2的原因了。
可是在上面的图里，我想你也看到了，有很多Fail的变异它的QD还是很高的，有些甚至高于30，通过这样的硬过滤参数所得到的结果中就会包含这部分本该要过滤掉的坏变异；而同样的，在低QD（<2）区域其实也有一些是好的变异，但是却被过滤掉了。这其实也是硬过滤的一大弊端，它不会像VQSR那样，通过多个不同维度的数据构造合适的高维分类模型，从而能够更准确地区分好和坏的变异，而仅仅是一刀切。
当你理解了上面有关QD的计算和阈值选择的过程之后，要弄懂后面的指标和阈值也就容易了，因为用的也都是同样的思路。
FisherStrand（FS）
FS是一个通过Fisher检验的p-value转换而来的值，它要描述的是测序或者比对时对于只含有变异的read以及只含有参考序列碱基的read是否存在着明显的正负链特异性（Strand bias，或者说是差异性）。这个差异反应了测序过程不够随机，或者是比对算法在基因组的某些区域存在一定的选择偏向。如果测序过程是随机的，比对是没问题的，那么不管read是否含有变异，以及是否来自基因组的正链或者负链，只要是真实的它们就都应该是比较均匀的，也就是说，不会出现链特异的比对结果，FS应该接近于零。
这里多说一句，在VCF的INFO中有时除了FS之外，有时你还会看到SB或者SOR。它们实际上是从不同的层面对链特异的现象进行描述。只不过SB给出的是原始各链比对数目，而FS则是对这些数据做了精确Fisher检验；SOR原理和FS类似，但是用了不同的统计检验方法计算差异程度，它更适合于高覆盖度的数据。
与QD一样，我们先来看一下质控前所有变异的FS总体密度分布图（如下）。很明显与QD相比，FS的范围更加的大，从0到好几百的都有。不过从图中也可以看出，绝大部分的变异还是在100以下的。
下面这一个图则是经过VQSR之后，画出来的FS分布图。跟上面的QD一样，浅绿色代表好变异，浅红色代表坏变异。我们可以看到，大部分好变异的FS都集中在0~10之间，而且坏变异的峰值在60左右的位置上，因此过滤的时候，我们把FS设置为大于60。其实设置这么高的一个阈值是比较激进的（留下很多假变异），但是从图中你也可以看到，不过设置得多低，我们总会保留下很多假的变异，既然如此我们就干脆选择尽可能保留更多好的变异，然后祈祷可以通过『执行硬过滤』里其他的阈值来过滤掉那些无法通过FS过滤的假变异。
StrandOddsRatio（SOR）
关于SOR在上面讲到FS的时候，我就在注释里提及过了。它同样是对链特异（Strand bias）的一种描述，但是从上面我们也可以看到FS在硬过滤的时候并不是非常给力，而且由于很多时候read在外显子区域末端的覆盖存在着一定的链特异（这个区域的现象其实是正常的），往往只有一个方向的read，这个时候该区域中如果有变异位点的话，那么FS通常会给出很差的分值，这时SOR就能够起到比较好的校正作用了。计算SOR所用的统计检验方法也与FS不同，它用的是symmetric odds ratio test，数据是一个2×2的列联表（如下），公式也十分简单，我把公式进行了简单的展开，从中可以清楚地看出，它考虑的其实就是ALT和REF这两个碱基的read覆盖方向的比例是否有偏，如果完全无偏，那么应该等于1。
sor = (ref_fwd/ref_rev) / (alt_fwd/alt_rev) = (ref_fwd * alt_rev) / (ref_rev * alt_fwd)OK，那么同样的，我们先看一下这个值总体的密度分布情况（如下）。总的来说，这个分布明显集中在0~3之间，这也和我们的预期比较一致。不过也有比较明显的长尾现象，这个时候我们也没必要定下太过明确的阈值预期，先看VQSR的分布结果。
下面这个图就是在VQSR之后，区分了好和坏变异之后，SOR的密度分布。很明显，好的变异基本就在1附近。结合这个分布图，我们在上面的例子里把它的阈值定为3基本上也不会过损失好的变异了，虽然单靠这个阈值还是会保留下不少假的变异，但是至少不合理的长尾部分可以被砍掉。
RMSMappingQuality（MQ）
MQ这个值是所有比对至该位点上的read的比对质量值的均方根（先平方、再平均、然后开方，如下公式）。
它和平均值相比更能够准确地描述比对质量值的离散程度。而且有意思的是，如果我们的比对工具用的是bwa mem，那么按照它的算法，对于一个好的变异位点，我们可以预期，它的MQ值将等于60。
下面是所有未过滤的变异位点上MQ的密度分布图。基本上就只在60的地方存在一个很瘦很高的峰。可以说这是目前为止这几个指标中图形最为规则的了，在这个图上，我们甚至就可以直接定出MQ的阈值了，比如所有小于50的就可以过滤掉了。
但是，理性告诉我们还是要看一下VQSR的对比结果（下图）。
你会发现似乎所有好的变异都紧紧集中在60旁边了，其它地方就都是假的变异了，所以MQ的阈值设置为50也是合理的。但是同样要注意到的地方是，60这个范围实际上依然有假的变异位点在那里，我们把这个区域放大来看，如下图，这里你就会发现其实假变异的密度分布图也覆盖到60这个范围了。
考虑到篇幅的问题，接下来MappingQualityRankSumTest（MQRankSum）和ReadPOSRankSumTest（ReadPOSRankSum）的阈值设定原理，我不打算再细说下去，思路和上面的4个是完全一样的。都是通过比较VQSR之后的密度分布图，最后确定了硬过滤的阈值。
但请不要以为这只是适用于GATK得到的变异，实际上，只要我们弄懂了这些指标选择的原因和过滤的思路，那么通过任何其他的变异检测工具也是依旧可以适用的，区别就在于GATK帮我们把这些要用的指标算好了。
同样地，这些指标也不是一成不变的，可以根据实际的情况换成其他，或者我们自己重新计算。
Ti/Tv处于合理的范围
Ti/Tv的值是物种在与自然相互作用和演化过程中在基因组上留下来的一个统计标记，在物种中这个值具有一定的稳定性。因此，一般来说，在完成了以上的质控之后，还会看一下这些变异位点Ti/Tv的值是多少，以此来进一步确定结果的可靠程度。
Ti（Transition）指的是嘌呤转嘌呤，或者嘧啶转嘧啶的变异位点数目，即A<->G或C<->T；
Tv（Transversion）指的则是嘌呤和嘧啶互转的变异位点数目，即A<->C，A<->T，G<->C和G<->T。（如下图）
另外，在哺乳动物基因组上C->T的转换比较多，这是因为基因组上的胞嘧啶C在甲基化的修饰下容易发生C->T的转变。
说了这么多Ti/Tv的比值应该是多少才是正常的呢？如果没有选择压力的存在，Ti/Tv将等于0.5，因为从概率上讲Tv将是Ti的两倍。但现实当然不是这样的，比如对于人来说，全基因组正常的Ti/Tv在2.1左右，而外显子区域是3.0左右，新发的变异（Novel variants）则在1.5左右**。**
最后多说一句，Ti/Tv是一个被动指标，它是对最后质控结果的一个反应，我们是不能够在一开始的时候使用这个值来进行变异过滤的。
小结
虽然本文一直在谈论的是如何做好硬过滤，但不管我们的指标和阈值设置的多么完美，硬过滤的硬伤都是一刀切，它并不会根据多个维度的数据自动设计出更加合理的过滤模式。硬过滤作为一个简单粗暴的方法，不到不得已的时候不推荐使用，即使使用了，也一定要明白每一个过滤指标和阈值都意味着什么。
最后，我也希望通过这一篇文章能够完整地为你呈现一个变异质控的思路。