獲得高質(zhì)量的SNP是進(jìn)行GWAS分析经柴，群體遺傳學(xué)分析的前提捉邢。GATK是比較流行的call snp方法播急。另外還有samtools和bcftools結(jié)合的方法再沧。實(shí)際的項(xiàng)目開(kāi)展根據(jù)需要而進(jìn)行庇配。
本篇的數(shù)據(jù)來(lái)源于2014年的一篇關(guān)于桃的重測(cè)序文章斩跌，直達(dá)連接
https://www.nature.com/articles/ncomms13246
這是文章中給與的call snp方法
方法上是將GATK和samtools+bcftools，這兩種策略結(jié)合在一起進(jìn)行捞慌。
那么耀鸦，總結(jié)一些文章的策略：
1，BWA比對(duì)到參考基因組啸澡，得到bam文件（這個(gè)當(dāng)然需要排序和去除PCR重復(fù)袖订，文章沒(méi)說(shuō)）；
2嗅虏，使用GATK進(jìn)行局部區(qū)域重比對(duì)：重比對(duì)的過(guò)程分為兩步：
第一步洛姑，RealignerTargetCreator ，定位出所有需要進(jìn)行序列重比對(duì)的目標(biāo)區(qū)域皮服；
第二步楞艾，IndelRealigner，對(duì)所有在第一步中找到的目標(biāo)區(qū)域運(yùn)用算法進(jìn)行序列重比對(duì)龄广，最后得到新的重比對(duì)bam文件硫眯。
3, 用2中獲得的重比對(duì)的bam文件進(jìn)行GATK和samtools兩種方法進(jìn)行call SNP
文章中使用的GATK 2.4版本。
2019年的文章的方法
1择同，BWA比對(duì)到參考基因組两入，得到bam文件（這個(gè)當(dāng)然需要排序和去除PCR重復(fù)，文章沒(méi)說(shuō)）奠衔；再用samtools去除quality <20的谆刨。
2塘娶，使用GATK進(jìn)行局部區(qū)域重比對(duì)：重比對(duì)的過(guò)程分為兩步：
第一步，RealignerTargetCreator 痊夭，定位出所有需要進(jìn)行序列重比對(duì)的目標(biāo)區(qū)域刁岸；
第二步，IndelRealigner她我，對(duì)所有在第一步中找到的目標(biāo)區(qū)域運(yùn)用算法進(jìn)行序列重比對(duì)虹曙，最后得到新的重比對(duì)bam文件。
3,使用PrintReads 做SNP檢測(cè)
4番舆，使用HaplotypeCaller做call snp
為得到高質(zhì)量的SNP值酝碳，設(shè)置一下參數(shù)：
-stand_call_conf 30 -stand_emit_conf 40.
5，硬過(guò)濾
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.
有些地方不是太明白疏哗。不過(guò)這個(gè)方法和2014年的文章基本差不多
關(guān)于-stand_call_conf 、-stand_emit_conf禾怠， 這兩個(gè)參數(shù)
-stand_emit_conf：在變異檢測(cè)過(guò)程中返奉，所容許的最小質(zhì)量值。只有大于等于這個(gè)設(shè)定值的變異位點(diǎn)會(huì)被輸出到結(jié)果中吗氏。
-stand_call_conf：在變異檢測(cè)過(guò)程中芽偏，用于區(qū)分低質(zhì)量變異位點(diǎn)和高質(zhì)量變異位點(diǎn)的閾值。只有質(zhì)量值高于這個(gè)閾值的位點(diǎn)才會(huì)被視為高質(zhì)量的弦讽。低于這個(gè)質(zhì)量值的變異位點(diǎn)會(huì)在輸出結(jié)果中標(biāo)注LowQual污尉。在千人基因組計(jì)劃第二階段的變異檢測(cè)時(shí)，利用35x的數(shù)據(jù)進(jìn)行snp calling的時(shí)候往产，當(dāng)設(shè)置成50時(shí)被碗，有大概10%的假陽(yáng)性。
參考;https://www.plob.org/article/7023.html
SNP calling的閾值簡(jiǎn)化為1個(gè)捂齐。
3.7以前使用兩個(gè)：-stand_call_conf 蛮放、-stand_emit_conf缩抡，
3.7中去掉了-stand_emit_conf奠宜，同時(shí)把-stand_call_conf的默認(rèn)值由30將為10。
參考：https://zhuanlan.zhihu.com/p/26262338
我又搜到了一篇GATK的流程：https://gencore.bio.nyu.edu/variant-calling-pipeline/瞻想。我們使用這個(gè)流程重復(fù)一下文章的一個(gè)數(shù)據(jù)压真。
在GATK中有一步是進(jìn)行BQSR。人類是有已知的變異位點(diǎn)信息蘑险，可以直接下載來(lái)用滴肿。但是非人類物種很多是沒(méi)有這些位點(diǎn)信息的。GATK原本就是為了人類數(shù)據(jù)設(shè)計(jì)的佃迄。但是泼差，這套流程針對(duì)的是非人類數(shù)據(jù)贵少。采用的策略是在call snp后，進(jìn)行過(guò)濾堆缘，得到高質(zhì)量的SNP滔灶，然后利用這個(gè)過(guò)濾的SNP進(jìn)行BQSR部分的操作。流程給的代碼非常清楚吼肥，先占個(gè)坑录平，這幾天實(shí)踐一下。下面詳細(xì)操作一下缀皱。
Overview of the pipeline
