序列比对碱基数统计的简单算法

评估两条核酸序列相似性的参数主要有两个：比对的覆盖度和比对序列的相似性。今天小编就和大家来聊聊如何高效准确地计算比对的覆盖度。

一般流程是：首先计算一条序列比对上参考基因组的总碱基数，然后除以该序列的总长，得到的数值即为此次比对的覆盖度。小编先后写过不同的脚本统计比对总碱基数，在此和大家分享一种相对简单的算法。对于一些没有编程基础的朋友可以从中了解一些编程的知识。对于常做生信分析的朋友，希望内附的源码能帮助您做覆盖度的分析统计。

A和B两序列比对，通过blast得到下列表格。现要统计A序列比对上B序列的总碱基数。

首先想到的做法是用A比对的终止位置(q.end)减去A比对的起始位置(q.start)，得到比对长度值，再把每行的比对长度进行累加。即如下公式：

但是，问题就出在这，在初步进行blast的结果中往往会有区域之间的重叠，例如上述列表红色标记的区域之间有重叠部分，蓝色标记的区域之间也有重叠部分。如果直接累加比对长度，会造成重叠区域重复统计，使计算值比实际值偏大。因此统计总比对长度，首先要得到没有重叠的比对区间，再累加各区间长度值。

我们把列表的q.start和q.end 这两列数据写成这种格式：

[[20508,22209],[51106,78043],.....]

按区间的第一位进行从小到大排序：

[[3,9991],[10560,28514],.....]

对排序后的列表，按如下算法进行处理：

主要算法

相邻的2个区间(a和b)，如果b区间起始值大于a区间结束值，则说明a和b没有重叠。

如果b区间起始值小于a区间结束值，则说明a区间和b区间有重叠，需要把a和b融合成一个新的区间，此时如果b区间结束值小于a区间结束值，说明b区间包含在a区间内部，新区间即为a。否则新区间为a的起始值到b的结束值。

以一个实际案例进行说明：原始的列表格式如下：

最终的得到的数值（15）就是去冗余之后，两条序列比对的总碱基数。这个思路同样可以应用的我们的实际统计中去，听完小编的讲解是不是瞬间感觉so easy！

最后给大家附上小编写好的源码（Python）。附. 源码（Python）：

运算结果为：

76662

备注：

上述源码省去了二维列表生成的代码。