从reads到SNP-GATK流程（待续）

获得高质量的SNP是进行GWAS分析，群体遗传学分析的前提。GATK是比较流行的call snp方法。另外还有samtools和bcftools结合的方法。实际的项目开展根据需要而进行。

本篇的数据来源于2014年的一篇关于桃的重测序文章，直达连接
https://www.nature.com/articles/ncomms13246

这是文章中给与的call snp方法

image.png

方法上是将GATK和samtools+bcftools，这两种策略结合在一起进行。

那么，总结一些文章的策略：
1，BWA比对到参考基因组，得到bam文件（这个当然需要排序和去除PCR重复，文章没说）；
2，使用GATK进行局部区域重比对：重比对的过程分为两步：
第一步，RealignerTargetCreator ，定位出所有需要进行序列重比对的目标区域；
第二步，IndelRealigner，对所有在第一步中找到的目标区域运用算法进行序列重比对，最后得到新的重比对bam文件。
3, 用2中获得的重比对的bam文件进行GATK和samtools两种方法进行call SNP
文章中使用的GATK 2.4版本。

2019年的文章的方法

image.png

1，BWA比对到参考基因组，得到bam文件（这个当然需要排序和去除PCR重复，文章没说）；再用samtools去除quality <20的。
2，使用GATK进行局部区域重比对：重比对的过程分为两步：
第一步，RealignerTargetCreator ，定位出所有需要进行序列重比对的目标区域；
第二步，IndelRealigner，对所有在第一步中找到的目标区域运用算法进行序列重比对，最后得到新的重比对bam文件。
3,使用PrintReads 做SNP检测
4，使用HaplotypeCaller做call snp
为得到高质量的SNP值，设置一下参数：
-stand_call_conf 30 -stand_emit_conf 40.
5，硬过滤
QUAL < 40, QD < 2.0, FS > 60.0, MQ < 40.0, MQRankSum < − 12.5, ReadPosRankSum < − 8.0
6，To further reduce false positives, the SNP number per 10 bp was limited to three using the following settings: --clusterWindowSize 10, --clusterSize 3.
有些地方不是太明白。不过这个方法和2014年的文章基本差不多

关于-stand_call_conf 、-stand_emit_conf， 这两个参数
-stand_emit_conf：在变异检测过程中，所容许的最小质量值。只有大于等于这个设定值的变异位点会被输出到结果中。
-stand_call_conf：在变异检测过程中，用于区分低质量变异位点和高质量变异位点的阈值。只有质量值高于这个阈值的位点才会被视为高质量的。低于这个质量值的变异位点会在输出结果中标注LowQual。在千人基因组计划第二阶段的变异检测时，利用35x的数据进行snp calling的时候，当设置成50时，有大概10%的假阳性。
参考;https://www.plob.org/article/7023.html

SNP calling的阈值简化为1个。
3.7以前使用两个：-stand_call_conf 、-stand_emit_conf，
3.7中去掉了-stand_emit_conf，同时把-stand_call_conf的默认值由30将为10。
参考：https://zhuanlan.zhihu.com/p/26262338

我又搜到了一篇GATK的流程：https://gencore.bio.nyu.edu/variant-calling-pipeline/。我们使用这个流程重复一下文章的一个数据。

在GATK中有一步是进行BQSR。人类是有已知的变异位点信息，可以直接下载来用。但是非人类物种很多是没有这些位点信息的。GATK原本就是为了人类数据设计的。但是，这套流程针对的是非人类数据。采用的策略是在call snp后，进行过滤，得到高质量的SNP，然后利用这个过滤的SNP进行BQSR部分的操作。流程给的代码非常清楚，先占个坑，这几天实践一下。下面详细操作一下。

Overview of the pipeline

image.png