生物的DNA都是由四种碱基构成的，这四种碱基A,T,G和C的排列顺序包含了物种的所有遗传信息。这些碱基的任意一个排列可以叫做一个k-mer，其中k指的是这个排列的长度，例如 AATT就是一个4-mer。在分析海量的NGS数据时，有时会需要统计序列里的k-mer分布情况来估计测序数据的质量和覆盖度，检测基因组的重复度等等。举个例子，给定一个15个碱基长的序列，它一共会有多少种碱基排列方式呢？这个计算不难，因为每个位置都会有可能是四种碱基里的任意一种，那么一共就会有4的15次方，也就是十亿种不同的碱基排列组合方式。那么有没有一种方法可以短时间内对NGS数据里十亿多个不同的15-mer的出现次数进行统计呢?

哈希算法(hashing)是一种非常快速高效的字符串索引和查找算法，它有到目前为止最好的的查找复杂度O(1)，理论上在一次查找之后就可以发现要找的东西是不是在数据库里。举个例子，如果一位生物老师想要在有两百位同学的名单上寻找一个人的名字，他应该会从名单的第一个开始看起，直到他发现这个名字为止，这样很正常但是会浪费时间。但如果他是一个比较geek的算法老师，他就有可能事先把这个名字通过某种变换变成一个数字N，然后把这个名字记在这个花名册的第N位，第二次需要寻找这个名字的时候就可以直接对这个名字进行相同的变换来得到它在花名册中的位置N了。虽然对两百个人这样做好像没有必要，但是如果名单上有一亿人的话，这样做就显现出它的优势来了。哈希算法的思路就是这样，详细的哈希算法可以去看看这篇文章

大自然用四种碱基来组成所有物种的遗传信息，这似乎和计算机里用0和1来进行所有运算有着异曲同工之妙，实际上它们确实是相通的。既然构成DNA的碱基只有四种，我们完全可以把碱基组成的序列看作是一个四进制数。假设在这个四进制数里A为0, T为1, C为2, G为3. 那么TTAAG就可以写成11003。如果我们把所有的4-mer从AAAA，AAAT，AAAC，AAAG，AATA等等一直列举到GGGG，其实就等同于把一个数从0数到了4的4次方罢了。当然也可以把每一个碱基转化为一个二进制数，那么ATGC就可以分别变成 00, 01, 10和11. 刚才的TTAAG就变成01`01`00`00`11了。这时你会发现：既然四进制数和二进制数是一一对应的（当然表示大小的数字本身不管在几进制里都会是一一对应的），那么对于一个确定的k，列举出的所有k-mer也会与二进制数一一对应，并且根据这个二进制数还能推导出原来的碱基序列。所以，我们的哈希函数到这里也已经写好了，那个与k-mer一一对应的二进制数就正好可以作为这个碱基序列在哈希表中的索引。

在一般情况下，哈希表一直有两个令人头疼的问题，第一个是表的空间常常不能有效利用，二是字符串的哈希值有时会发生碰撞，导致两个不同的字符串会被存储到哈系表中的同一位置。但是，在用上面的哈希函数来转换碱基序列后，这两个问题都不存在。这个哈希函数不仅转换起来相当简便，甚至还能毫不费力地推回到原来的序列，这样我们就不需要在内存里存储序列本身，这对节省内存空间和运算量带来了很大的余地，这就是为什么我们要选择用二进制数来代表K-mer了。实际上所有的数字在计算机里都是以二进制的形式来储存的，这样的好处是，不仅通过简单的位运算就可以把K-mer换成相对应的二进制数，而且在后来将K-mer在一条序列上移动的时候，只需要进行两次位运算就可以得到一个新的K-mer的索引，方便并且直观。

有了哈希函数，下一步就是做哈希表了。这个表实际上是一个数列。拿15-mer来举例，我们有4的15次方也就是十亿条不同的K-mer，这就需要有一个十亿位的数组来对应。虽然我们不需要存储序列本身，（因为数组的索引本身就代表着一个不同的K-mer）但是数组的每一位还是需要占用四个字节来统计它对应K-mer的出现次数，那么这个十亿位的数组总共就要占用4G的内存。这样的内存占用在普通电脑里就已经有点吃力了，当然这也是哈希表用空间换时间的一个缺点。

所以大概的思路就是，在一开始我们可以把这个数组的每一位都赋为0, 代表到目前为止还没有见到过对应的K-mer，然后对于每一个在序列里遇到的15-mer，我们把它逐位转换为二进制数，然后直接在数列中将以这个二进制数为索引对应的值加1, 表示我们见到过了一次这个15-mer，等遍历完所有的序列以后，这个数列也就被填满了。

在遍历一条序列的时候，假设序列长度为N，那么这条序列就含有N-K+1个K-mer。如果我们逐个截取长度为15的小片段，然后再将之转化为二进制数的话，对于每条序列，我们就需要进行(N-K+1)*15次转换，实际上，在这个过程中发生了很多重复运算。因为两个相邻15-mer之间有14个碱基是重叠的，也就是说，前一个15-mer所对应的30位二进制数中的后28位是与后一个K-mer所对应的二进制数的前28位是完全一样的，我们没有必要整个重新计算下一个K-mer的索引，只需要读进下一个碱基，然后把这个碱基对应的两位二进制数加在这个二进制数的末尾，然后顶掉最头上的两位数字就OK了。两次位运算就可以得到一个新的K-mer的存储位置，这要比老老实实转换15次方便得多，这样的好处是，整体运算量并不会随着K的增大而增大一丁点儿。结果对于长度为N的序列，不管K是多少，只需要进行N次字母到二进制数的转换就能完成这条序列中所有的K-mer的统计，就象下面这样：

可是在实际应用中，我们一般会将一条K-mer和它的反向互补序列作为同一条序列来看待，如何才能办到这一点呢？如果等全部K-mer的频率都统计出来之后再将所有反向互补序列出现的次数分别加起来的话，计算量是相当可观的。其实我们可以在哈希函数上做一些改进，让两个反向互补的序列都能哈希到同一个哈希值，这样就等于在计算的过程中顺带地解决了这个问题。有一个办法是，将每个K-mer的最终哈希值定为它和它反向互补序列的哈希值的和，这样虽然行得通，但不是最好的办法。因为两数相加的哈希值必然会超过原来的最大值，那就需要更大的数列和内存去装它，这个办法对于大的K-mer一点也不合用。一个更好的办法是：在计算一个K-mer的哈希值时，也同时对它的反向互补序列做哈希，然后永远取两数中小的那个哈希值作为它们在数组里的共同索引。这个方法的好处在于，它不仅不会要求更大的空间，用得巧的话反而能够正好减小一半的内存使用。但是这需要巧妙地设置ATGC四个字母所对应的二进制数：假设A用11, T用01的话， 那么以A开头的序列的索引和它反向互补序列的索引一定是在索引最大值二分之一的两旁。也就是说，在索引小于中点的序列的反向互补序列的索引一定大于这个中点，这只是在AT和GC分别取11,10和00,01中任意一个的情况下才成立，如果AT分别取11和10的话可就不行了。

用上面这个方法，刚才统计15-mer时原本需要的4G内存就可以减小到2G，但是在遍历每条序列的时候， 需要同时计算K-mer和它的反向互补序列的索引然后比较大小。对每条序列来说，就意味着多做N次转换运算和N-K+1次比大小运算，但是和减小一半的内存比起来，这样做还是值得的，这是这个算法的python代码。

当我们得到所有的K-mer和它反向互补序列的出现次数后，下一步就是报告出现频率最高的K-mers了。这一步貌似容易，直接对这个数列排序不就可以了吗？但实际上这样是行不通的。暂且不考虑直接排序的时间复杂度，等到把这个数列排好序之后你会发现，一个很重要的信息丢失了。因为排序的过程中需要换位，但是位置一换，原来的索引就乱了，这样就不能通过原来的索引确定它所对应的K-mer，排序的意义就不是很大了。实际上在现实应用中，我们感兴趣的常常只是出现频率最高的那几个K-mer，在这种情况下，我们可以考虑建立一个小数列来遍历这个数组，随时更新记录下这个数组里的最大值和它所对应的索引就OK了。这时，堆一定是一个再好不过的选择。堆是一个设计相当精巧的二叉树，在这个二叉树中，父节点的值总是大于等于或小于等于它两个子节点所对应的值。所以一个数列里的最大值或者最小值总是保存在这个树的根部。这样做的好处是，我们在任何时候都只需要一次操作就可以提取出这个堆里的最大值或最小值。堆有两种，父节点总是大于等于它的子节点的叫做最大堆，反之就是最小堆，在最大堆里，提取在根部的最大值非常方便，但是寻找最小值需要花很长时间，最小堆和它恰好相反，寻找最小值很方便，但是最大值就不是那么明显了。在这里，我们只需要保存一个固定大小的最小堆（最小值在根部的堆），然后遍历整个数列，把每次遇到的新值和这个堆里的最小值作对比，然后来确定是否替换掉最小值，对这个堆重新排序就好了。

在这个过程中，我们只需要对这个最小堆进行两个操作，找最小值和替换最小值然后重新排列，而对堆来说，这两个操作所需要的时间分别是O(1)和O(logN), 所以用堆来遍历这个数列找到TOP m的时间复杂度就只是Nlog(m), 和直接对数组进行排序的O(Nlog(N))的复杂度要进步不少。所以说堆应用在这里可以说最适合不过了。这个堆的python代码可以在这里找到。

从整体来看，整个运算时间无疑是比较快的。遍历序列的时间复杂度是O(N),而找出Top m的时间复杂度是O(Nlog(m)),加起来就是O(Nlog(m))，因为m通常很小，一般是10左右，所以说几乎可以把m看作是常数。但是从空间上来看就不是那么乐观了，这也是这个算法的缺点，一个15-mer的运算都要消耗至少2G的内存，这就是为什么一些k-mer计数软件比如说Sylamer只支持到15mer的原因了。如果对哈希函数做一些小的修饰，就会达到更节省空间目的，但是这也就会相应地牺牲一些查找速度，这就需要具体问题具体分析了。