一定要使用包含Decoy的参考基因组啊！

缘起

今天有人咨询我，PRSS1基因NM_002769.5转录本的 c.86A>T，p.Asn29Ile 突变（rs111033566，位于chr7，hg19坐标为chr7:142458451，hg38坐标为chr7: 142750600）是否应该报告？这个基因组在ClinVar中被认为Pathogenic，据文献报道跟Hereditary pancreatitis等遗传疾病相关。但是，该突变接连出现在多个健康人的检测样品中，致病突变似乎不应该有如此高的人群频率！

于是进入dbSNP，看看该突变在人群中的频率，一查发现，该突变在不同队列中，仿佛分为两个“流派”！

首先，在诸如gnomAD-Genomes、1000Genomes等项目中，该突变以A>C为主，且人群突变频率极低：

然而在Exac中，该位点又变成了A>T为主，且人群突变频率接近50% ！

基因组中添加Decoy

于是，我随机打开先前分析过的几个全外数据，发现这个位点上几乎都是A>T突变：

而且，带有突变的Reads，在上游也存在T>C以及下游的几处点突变（若覆盖到）。

随机点开一个带有突变的Reads，发现该Reads比对分数不高，且比对到其他坐标！

偶然间，发现另外一家公司提供的bam文件中，该位置没有发现这类Reads，于是查看bam文件头上的操作信息，发现他们使用的参考基因组为hs37d5。

相比我自己使用的GATK b37，hs37d5只是多出了两条病毒序列，以及李恒大神提供的Decoy序列。莫非Decoy序列会帮我们把比对结果弄干净了？

于是赶紧下载hs37d5的序列，构建bwa的index，重新比对后终于发现：

带有Decoy的hs37d5比对结果中，该突变的Reads消失了（IGV截图上半截）！

Decoy咋来的？

在hs37d5下载的ftp上，附有李恒大神的讲解幻灯片：

这群Decoy序列在1000genome项目时代已经被发现！此外还添加了部分重复序列（422kb）和其他无法比对到人参考基因组上的序列。最终去完冗余，这对序列一共34.5Mb。

添加了Decoy，在李恒大神的幻灯片中，效果更是立竿见影：

看来Exac的频率也不那么可信。。。

联想到，在没有Decoy的比对结果中，chr7:142458371 T>C（rs367779270）总伴随着带有rs111033566的reads出现，顺便也查查rs367779270在dbSNP中的突变记录：

人群频率还是一般泾渭分明！

Exac项目的Nature论文中提及的方法，是将Reads比对到gencode v19的 hg19 参考基因组版本上的。有没有另外加Decoy我不好揣测。但现在看，Exac中给出的部分突变人群频率，或许要存疑了！