SNP Calling（一）

        最近和Соня共同忙于准种进化分析，从上一篇的推文的数学模型中可以看出DNA病毒在复制过程中也可以产生类似于准种的基因组群体。但是因为DNA病毒突变速度慢，基因保守。所以在复制过程中无法像RNA病毒那样产生足够量的突变的基因组群体。尽管如此，我们必须承认，在DNA病毒复制过程中突变的基因组也会出现的。我们可以借鉴准种的概念。

        要分析复制过程中产生的突变，分析SNP和SNV是一个很好的方法，因为SNV分析的部分我有保密内容，不能分享，我只会分享出我自己写的SNP分析的脚本。在这个过程中，Соня给了我很大的提示，正因为她的帮助，我才能顺利完成这个脚本！

        对SNP的正确分析要建立在参考基因组的正确组装和reads的正确mapping上。这里说明一下这句话的含义，首先将亲本病毒侵染细胞或者细菌，侵染后由于病毒在细胞或者细菌内复制过程中自身发生突变，细胞或者细菌内的免疫系统又给予复制的病毒自然选择。那些发生了能够抵抗细胞或者细菌免疫的突变的病毒才有可能继续侵染其他细胞或细菌。那么要分析发生了什么突变，就可以用亲本病毒基因组为参考基因组，将测序后突变的病毒测序的reads mapping回亲本病毒基因组，分析变异情况。

          而实际的分析中可能我们要选择作为参考基因组的病毒并没有完整的基因组，需要de novo组装，这个组装过程中需要对测序数据进行质控，过滤，纠错。质控的目的是在于检测测序质量，过滤是将质量低的碱基或者reads cut掉。纠错是通过一定的算法将reads上一些可能是测序错误的位点校正。质控通常采用FastQC软件。过滤可以自己使用Shell命令过滤，也可以使用Trimmomatic过滤，如果使用Trimmomatic需要保证你的测序数据中的接头序列能在Trimmomatic的adapters文件中找到。如果没有，则需要我们公众号前面的教程通过R语言爬虫（一）的方法找接头序列（我将本之前的那个教程放入原文链接中）。对于纠错，SPAdes中的–only-error-correction模式是非常好的纠错模式，它能直接调用bwa比对软件进行纠错。

        最后，在进行SNP Calling前需要将变异后的病毒reads mapping到参考基因组上，将mapping后的数据排序，使用samtools就能对比对到参考基因组上同一区域的reads进行变异分析统计了，对于那些数量少的能mapping到参考基因组上的reads，这些reads在某些位点上存在变异且这个变异频率在1%以上，就能被统计为SNP。如下图所示

图中灰色部分的彩色的小竖线就是可能为变异的位点，灰色为reads，最下面的彩色条形部分为参考基因组。如果这些变异位点在出现频率上满足条件，就可以被统计为SNP。下面我分享出我做SNP Calling的脚本。（脚本中用到的软件需要自行安装）

运行完以上脚本，可以得到SNP的变异统计数据，用Excel打开可以查看SNP信息

例如第二列代表SNP在参考基因组上位置，第三列是参考基因组该位点碱基，第四列是该位点突变为的碱基，后面的列代表这里SNP的统计学信息。例如，第五列表示Phred格式(Phred_scaled)的质量值，可以理解为所call出来的变异位点的质量值。表示在该位点存在variant的可能性；该值越高，则variant的可能性越大，计算方法：Phred值Q = -10 * lg (1-p) ，p为variant存在的概率。因为SNP是依据一定的统计学模型计算出来的，不一定是真实的SNP，所以FILTER列表示使用模型对对变异位点进行过滤，过滤完了之后，在FILTER一栏都会留下过滤记录，如果是通过了过滤标准，那么这些通过标准的好的变异位点的FILTER一栏就会注释一个PASS，如果没有通过过滤，就会在FILTER这一栏提示除了PASS的其他信息。如果这一栏是一个“.”的话，就说明没有进行过任何过滤。由于过滤SNP通常会取决于你的研究目的，所以我在这里并没有过滤SNP。

阅读原文