生信（九）生信代码中的位操作

本文介绍了生信代码中用到的一些位操作。

我们知道，0和1构成的二进制充斥着计算机语言的世界。一般来说，我们对二进制可以操作的最小单位就是一个bit（位）了，一个bit要么是0，要么是1。在编写代码的过程中，如果我们能了解一点位操作，有时可以简化代码、提高效率。

这一点对于生信的编程同样适用。

应用一：列举k-mer

比如，在《算法（三）列举所有k-mer的组合》一文中，笔者曾经分享过一段代码，意在解决NGS数据分析中时常会碰到的列举k-mer的问题：

“如何打印出特定长度的全部 k-mer 呢？”

当然可以用递归的方法（见前文，此处略过），但是下面的代码更简洁：

#include <stdio.h>

int PrintKMer(const int k) {  // for k <= 31
    int i;
    unsigned long long x, y;

    for (x = 0; x < 1ULL<<(2 * k); ++x) {
        for (i = 0, y = x; i < k; ++i, y >>= 2)
            putchar("ATGC"[y & 3]);
        putchar('\n');
    }
    return 0;
}

int main(void) {
    PrintKMer(2);
    return 0;
}

这个方法巧妙地运用位操作解决了问题，读者可以仔细品读，该段代码源自lh3在biostars网站上的分享。效果如下：

应用二：寻找最接近的2的幂

在NGS领域著名的kseq.h这个头文件中，我们可以看到lh3另一段运用位操作的代码：

#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))

这段代码想要实现的功能是“返回不小于x并且最接近x的2的幂”。当然这是有应用限制的，只能对小于等于$2^{30}$（有符号整数）或者小于等于$2^{31}$（无符号整数）的整数起作用。

代码中加上--(x)的效果是当x是2的幂时，会返回原值（比如，当x=8时，会返回8）。如果去掉--(x)这一小句，那么当x是2的幂时，会返回下一个2的幂（比如，当x=8时，会返回16）。

我们运行测试代码：

#include <stdio.h>
#include <stdint.h>
#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))

int main(void) {
    int32_t a[] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1073741823, 1073741824, 1073741825};
    int n = 14;
    int i, k, j;
    for (i = 0; i < n; i++) {
        k = j = a[i];
        printf("%d --> %d\n", k, kroundup32(j));
    }
    return 0;
}

效果是：

更多关于位操作的技巧

从上面两个应用来看，位运算的确可以应用于生信领域的代码中。那么为什么要用位操作呢？一般有两个原因：一是很多时候运用位操作可以简化代码；二是运用位操作通常可以提高代码运行效率（比起乘法和除法操作来说）。

如果你想了解更多位操作的技巧，可以参考Bit Twiddling Hacks这个网站，其实上文“寻找最接近的2的幂”的代码也出现在了该网站的小节中。

除此以外，里面还有很多经过验证的实用的位操作。比如：

如何判断一个数是不是2的幂

可以这样做：

unsigned int v; // we want to see if v is a power of 2
bool f;         // the result goes here 
f = (v & (v - 1)) == 0;

//Note that 0 is incorrectly considered a power of 2 here. To remedy this, use:
f = v && !(v & (v - 1));

有时间的时候去看看，说不定就可以获取一些启发。