我是怎么被基因组坑了10000块钱的
说真的，这都得怪基因组...
事情起因很简单，我有个课题呢就是对40个人病人做了WES，拿到原始测序数据，经过质控 > mapping > GATK > 注释 > 筛选潜在致病位点 > 大队列sanger验证一条龙服务后，我发现有个被验证的突变位点在其他300个该病患者中一个携带者都没验证到（麻蛋，300个样10000块钱啊，一个都没验证到，老板知道了不得踹死我...）
这个位点在前期的40个WES的样本中携带率超过90%，队列验证怎么可能一个携带者都没有，这不科学...难道是我的筛选流程有误？
我怀疑是他们验证实验的问题，最容易出问题的就是扩增的不是目标序列了。于是我要来了跑胶图和引物序列，引物特异性很好，引物序列blast的结果也是目标序列，这说明应该不是扩增的问题。呃..难道是sanger实验有误？还好之前WES的血样还有剩下的，拿去扩增并sanger验证了下，结果让我大吃一惊，这个位点的40个样本的sanger验证一个都没验证到，也就是说这个位点WES的结果和sanger的结果完全不一致..要么公司纯粹在忽悠我，要么就是我的snp calling有误...
我不想背锅啊..在和公司拉锯式撕b的过程中，我还是自查了整个分析流程，但一无所获。索性我重头跑了一遍，为了排除所有可能存在的干扰因素，我这次跑了两个pipeline，一个是用hg19跑的，一个是用hg38跑的（之前使用hg19跑的），万万没想到啊，这个位点在hg19中是被call出来了，但在hg38中并没有被call出来，这TM就尴尬了..随便找个样本在IGV里放大这个位点一看。
该位点在某个样本A中hg19的IGV结果如下，表现很好啊，这颜色，这分布，这深度，简直perfect！
但再看下该位点在某个样本A中hg38的IGV结果， WTF，variant-supporting reads不见了..
是的，同样位置，hg19中variant-supporting reads怎么在hg38中不见了，我随便找了个hg19中携带该突变的read，根据read name在hg38的bam文件中看这个read死哪去了。然后就发现这一对reads在hg38中完美的匹配到chrUn（unplaced sequences，是人的序列，但不知道在哪）上去了，一搜hg19的fasta里果然没有这条序列。
所以呢，就是由于基因组版本差异导致本该map到chrUn上的reads错误的map到了某个基因，然后我就根据错误的mapping结果call出了错误的突变。其实在结果出来之前我已经猜到了可能是由于基因组差异导致的，但看文章里大家对于基因组版本的选择并没有太高的要求，hg19和hg38都ok，我也觉得应该没这么巧吧，结果刚好是版本差异序列产生的假阳性位点。至于我为什么选择hg19的基因组版本，因为hg19的注释数据多啊。。
以上就是我被基因组坑了10000块钱的整个过程，由于数据敏感性的问题，马赛克打的有点多，分析流程和很多信息没有列出来。当然现在比较重要的是我怎么和老板解释....
hg38和hg19的差异还不止在于chrUn，hg38还新增了ALT contig，修正了数千的错误序列，GATK也strongly recommend 使用hg38（https://software.broadinstitute.org/gatk/documentation/article?id=11010），我也建议大家都适用hg38吧..至于外部注释信息，可以用liftover转为hg38的再进行注释VCF。
更多原创精彩视频敬请关注生信杂谈:
