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

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

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

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

备注:

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

  • 发表于:
  • 原文链接:https://kuaibao.qq.com/s/20180621G1LRJ300?refer=cp_1026
  • 腾讯「云+社区」是腾讯内容开放平台帐号(企鹅号)传播渠道之一,根据《腾讯内容开放平台服务协议》转载发布内容。

扫码关注云+社区

领取腾讯云代金券