将核碱基表示形式从 ASCII 转换为 UCSC .2 位

Question

明确的 DNA 序列仅由核碱基腺嘌呤 (A)、胞嘧啶 (C)、鸟嘌呤 (G)、胸腺嘧啶 (T) 组成。对于人类消费，基数可以用相应的char大写或小写字母表示：A, C, G, T, 或a, c, g, t。然而，当需要存储长序列时，这种表示方式效率低下。由于只需要存储四个符号，因此可以为每个符号分配一个2位代码。UCSC指定的常用.2bit格式正是这样做的，使用以下编码：T = 、C = 、A = 、G = 。0b000b010b100b11

下面的 C 代码显示了为清晰起见而编写的参考实现。转换表示为序列的基因组序列的各种开源软件通常使用由每个序列char索引的 256 条目查找表。char这也与 的内部表示隔离char。然而，即使访问的是片上高速缓存，存储器访问也非常昂贵，并且通用表查找很难进行 SIMDize。因此，如果可以通过简单的整数算术来完成转换，则将是有利的。鉴于 ASCII 是占主导地位的char编码，我们可以对此进行限制。

将 ASCII 字符给出的核碱基转换为其.2bit表示形式的有效计算方法是什么？

/* Convert nucleobases A, C, G, T represented as either uppercase or lowercase
   ASCII characters into UCSC .2bit-presentation. Passing any values other than
   those corresponding to 'A', 'C', 'G', 'T', 'a', 'c', 'g', 't' results in an
   indeterminate response.
*/
unsigned int ascii_to_2bit (unsigned int a)
{
    const unsigned int UCSC_2BIT_A = 2;
    const unsigned int UCSC_2BIT_C = 1;
    const unsigned int UCSC_2BIT_G = 3;
    const unsigned int UCSC_2BIT_T = 0;

    switch (a) {
    case 'A':
    case 'a':
        return UCSC_2BIT_A;
        break;
    case 'C':
    case 'c':
        return UCSC_2BIT_C;
        break;
    case 'G':
    case 'g':
        return UCSC_2BIT_G;
        break;
    case 'T':
    case 't':
    default:
        return UCSC_2BIT_T;
        break;
    }
}

Answer 1

如果专心地盯着核碱基 ASCII 字符的二进制代码，就会清楚地看到，位 1 和位 2 提供了唯一的两位代码：A= 0b01000001-> 0b00、C= 0b01000011-> 0b01、G= 0b01000111-> 0b11、T= 0b01010100-> 0b10。与小写 ASCII 字符类似，仅在第 5 位上有所不同。不幸的是，这个简单的映射与 -encoding 不太匹配.2bit，因为 A 和 T 的代码被交换了。解决此问题的一种方法是使用存储在变量中的简单四条目排列表，可能在优化后分配给寄存器（“寄存器内查找表”）：

unsigned int ascii_to_2bit_perm (unsigned int a)
{
    unsigned int perm = (2 << 0) | (1 << 2) | (0 << 4) | (3 << 6);
    return (perm >> (a & 6)) & 3;
}

另一种方法通过观察 和 的位 1 为 0 ，而 和 的位 1 为A1 ，纠正通过简单的位操作即时提取位 1 和 2 生成的“原始”代码。因此，我们可以通过将输入的第 1 位与初始代码的第 1 位进行异或来交换 A 和 T 的编码：TCG

unsigned int ascii_to_2bit_twiddle (unsigned int a)
{
    return ((a >> 1) & 3) ^ (~a & 2);
}

该版本对于具有快速位域提取指令的处理器和没有桶形移位器的低端处理器来说是有利的，因为只需要右移一位位置。在乱序处理器上，这种方法比排列表提供更多的指令级并行性。由于在所有字节通道中使用相同的移位计数，因此适应 SIMD 实现似乎也更容易。

在我专心研究相关 ASCII 字符的二进制编码之前，我研究过使用简单的数学计算。对小乘数和除数进行简单的暴力搜索得到：

unsigned int ascii_to_2bit_math (unsigned int a)
{
    return ((18 * (a & 31)) % 41) & 3;
}

该乘法器18对于没有快速整数乘法器的处理器很友好。现代编译器可以有效地处理带有编译时常量除数的模计算，并且不需要除法。即便如此，我注意到即使是最好的可用编译器也难以利用非常有限的输入和输出范围，因此我必须手动处理它以简化代码：

unsigned int ascii_to_2bit_math (unsigned int a)
{
    unsigned int t = 18 * (a & 31);
    return (t - ((t * 25) >> 10)) & 3;
}

即使以这种形式并假设单周期乘法的可用性，这似乎通常与之前的两种方法没有竞争力，因为它产生更多的指令和更长的依赖链。然而，在 64 位平台上，整个计算可以替换为 64 位、32 项查找表，如果该 64 位表可以有效地放入寄存器中，则可以提供有竞争力的性能（x86 就是这种情况） 64，它作为立即数加载。

unsigned int ascii_to_2bit_tab (unsigned int a)
{
    uint64_t tab = ((0ULL <<  0) | (2ULL <<  2) | (0ULL <<  4) | (1ULL <<  6) |
                    (3ULL <<  8) | (0ULL << 10) | (2ULL << 12) | (3ULL << 14) |
                    (1ULL << 16) | (3ULL << 18) | (0ULL << 20) | (2ULL << 22) |
                    (3ULL << 24) | (1ULL << 26) | (2ULL << 28) | (0ULL << 30) |
                    (1ULL << 32) | (3ULL << 34) | (1ULL << 36) | (2ULL << 38) |
                    (0ULL << 40) | (1ULL << 42) | (3ULL << 44) | (0ULL << 46) |
                    (2ULL << 48) | (0ULL << 50) | (1ULL << 52) | (3ULL << 54) |
                    (0ULL << 56) | (2ULL << 58) | (3ULL << 60) | (1ULL << 62));
    return (tab >> (2 * (a & 31))) & 3;
}

我附上我的测试框架以供参考：

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>

#define ORIGINAL_MATH (1)

unsigned int ascii_to_2bit_perm (unsigned int a)
{
    unsigned int perm = (2 << 0) | (1 << 2) | (0 << 4) | (3 << 6);
    return (perm >> (a & 6)) & 3;
}

unsigned int ascii_to_2bit_twiddle (unsigned int a)
{
    return ((a >> 1) & 3) ^ (~a & 2);
}

unsigned int ascii_to_2bit_math (unsigned int a)
{
#if ORIGINAL_MATH
    return ((18 * (a & 31)) % 41) & 3;
#else // ORIGINAL_MATH
    unsigned int t = 18 * (a & 31);
    return (t - ((t * 25) >> 10)) & 3;
#endif // ORIGINAL_MATH
}

unsigned int ascii_to_2bit_tab (unsigned int a)
{
    uint64_t tab = ((0ULL <<  0) | (2ULL <<  2) | (0ULL <<  4) | (1ULL <<  6) |
                    (3ULL <<  8) | (0ULL << 10) | (2ULL << 12) | (3ULL << 14) |
                    (1ULL << 16) | (3ULL << 18) | (0ULL << 20) | (2ULL << 22) |
                    (3ULL << 24) | (1ULL << 26) | (2ULL << 28) | (0ULL << 30) |
                    (1ULL << 32) | (3ULL << 34) | (1ULL << 36) | (2ULL << 38) |
                    (0ULL << 40) | (1ULL << 42) | (3ULL << 44) | (0ULL << 46) |
                    (2ULL << 48) | (0ULL << 50) | (1ULL << 52) | (3ULL << 54) |
                    (0ULL << 56) | (2ULL << 58) | (3ULL << 60) | (1ULL << 62));
    return (tab >> (2 * (a & 31))) & 3;
}

/* Convert nucleobases A, C, G, T represented as either uppercase or lowercase
   ASCII characters into UCSC .2bit-presentation. Passing any values other than
   those corresponding to 'A', 'C', 'G', 'T', 'a', 'c', 'g', 't' results in an
   inderminate response.
*/
unsigned int ascii_to_2bit (unsigned int a)
{
    const unsigned int UCSC_2BIT_A = 2;
    const unsigned int UCSC_2BIT_C = 1;
    const unsigned int UCSC_2BIT_G = 3;
    const unsigned int UCSC_2BIT_T = 0;

    switch (a) {
    case 'A':
    case 'a':
        return UCSC_2BIT_A;
        break;
    case 'C':
    case 'c':
        return UCSC_2BIT_C;
        break;
    case 'G':
    case 'g':
        return UCSC_2BIT_G;
        break;
    case 'T':
    case 't':
    default:
        return UCSC_2BIT_T;
        break;
    }
}

int main (void)
{
    char nucleobase[8] = {'A', 'C', 'G', 'T', 'a', 'c', 'g', 't'};
    printf ("Testing permutation variant:\n");
    for (unsigned int i = 0; i < sizeof nucleobase; i++) {
        unsigned int ref = ascii_to_2bit (nucleobase[i]);
        unsigned int res = ascii_to_2bit_perm (nucleobase[i]);
        printf ("i=%2u %c res=%u ref=%u %c\n", 
                i, nucleobase[i], res, ref, (res == ref) ? 'p' : 'F');
    }
    printf ("Testing bit-twiddling variant:\n");
    for (unsigned int i = 0; i < sizeof nucleobase; i++) {
        unsigned int ref = ascii_to_2bit (nucleobase[i]);
        unsigned int res = ascii_to_2bit_twiddle (nucleobase[i]);
        printf ("i=%2u %c res=%u ref=%u %c\n", 
                i, nucleobase[i], res, ref, (res == ref) ? 'p' : 'F');
    }
    printf ("Testing math-based variant:\n");
    for (unsigned int i = 0; i < sizeof nucleobase; i++) {
        unsigned int ref = ascii_to_2bit (nucleobase[i]);
        unsigned int res = ascii_to_2bit_math (nucleobase[i]);
        printf ("i=%2u %c res=%u ref=%u %c\n", 
                i, nucleobase[i], res, ref, (res == ref) ? 'p' : 'F');
    }
    printf ("Testing table-based variant:\n");
    for (unsigned int i = 0; i < sizeof nucleobase; i++) {
        unsigned int ref = ascii_to_2bit (nucleobase[i]);
        unsigned int res = ascii_to_2bit_tab (nucleobase[i]);
        printf ("i=%2u %c res=%u ref=%u %c\n", 
                i, nucleobase[i], res, ref, (res == ref) ? 'p' : 'F');
    }
    return EXIT_SUCCESS;
}

Answer 2

一种选择如下：

unsigned int ascii_to_2bit (unsigned int a)
{
    return ((0xad - a) & 6) >> 1;
}

这样做的优点是它只需要 8 位，不会溢出，并且不包含非常量移位，因此即使没有特定的 SIMD 指令，例如以 64 位输入 8 个字符，它也可以立即用于并行化。单词

unsigned int ascii_to_2bit_8bytes (uint64_t a)
{
    return ((0xadadadadadadadad - a) & 0x0606060606060606) >> 1;
}

将两个输出位保留在每个字节的底部。