将核碱基表示从ASCII转换为UCSC .2bit

ymzxtsji  于 2024-01-06  发布在  其他
关注(0)|答案(3)|浏览(292)

非抗生素DNA序列仅由核碱基腺嘌呤（A）、胞嘧啶（C）、鸟嘌呤（G）、胸腺嘧啶（T）组成。对于人类消费，碱基可以用相应的char表示，在char或char中：A、C、G、T或a、c、g、t。然而，当需要存储长序列时，这种表示是低效的。由于只需要存储四个符号，每个符号都可以分配一个2位代码。2 UCSC规定的常用.2bit正是这样做的，使用以下编码：T = 0b00，C = 0b01，A = 0b10，G = 0b11。
下面的C代码显示了为清晰起见而编写的参考实现。各种转换以char序列表示的基因组序列的开源软件通常使用由序列中的每个char索引的256条目查找表。这也与char的内部表示隔离。然而，存储器访问在能量上是昂贵的，即使访问的是片上缓存，并且通用表查找难以SIMD化。因此，如果转换可以通过简单的整数运算来完成，这将是有利的。考虑到ASCII是主要的char编码，可以限制于此。
什么是有效的计算方法来转换作为ASCII字符的核碱基到它们的.2bit-表示？

  1. /* Convert nucleobases A, C, G, T represented as either uppercase or lowercase
  2.    ASCII characters into UCSC .2bit-presentation. Passing any values other than
  3.    those corresponding to 'A', 'C', 'G', 'T', 'a', 'c', 'g', 't' results in an
  4.    indeterminate response.
  5. */
  6. unsigned int ascii_to_2bit (unsigned int a)
  7. {
  8.     const unsigned int UCSC_2BIT_A = 2;
  9.     const unsigned int UCSC_2BIT_C = 1;
  10.     const unsigned int UCSC_2BIT_G = 3;
  11.     const unsigned int UCSC_2BIT_T = 0;
  13.     switch (a) {
  14.     case 'A':
  15.     case 'a':
  16.         return UCSC_2BIT_A;
  17.         break;
  18.     case 'C':
  19.     case 'c':
  20.         return UCSC_2BIT_C;
  21.         break;
  22.     case 'G':
  23.     case 'g':
  24.         return UCSC_2BIT_G;
  25.         break;
  26.     case 'T':
  27.     case 't':
  28.     default:
  29.         return UCSC_2BIT_T;
  30.         break;
  31.     }
  32. }

字符串

c

3条答案

按热度按时间
kmynzznz

kmynzznz1#

如果我们仔细观察碱基的ASCII字符的二进制代码，就会发现第1位和第2位提供了一个唯一的两位代码：A = 0b01000001 -> 0b00，C = 0b01000011 -> 0b01，G = 0b01000111 -> 0b11，T = 0b01010100-> 0b10。类似于ASCII字符，仅在第5位不同。不幸的是，这种简单的Map并不完全匹配.2bit编码，其中A和T的代码被交换。解决这个问题的一种方法是用一个简单的四项排列表存储在一个变量中，可能在优化后分配给一个寄存器（“寄存器内查找表”）：

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

字符串
另一种方法通过观察A和T的位1为0，但C和G为1。因此，我们可以通过对 inverse 进行XOR来交换A和T的编码。将输入的位1与初始代码的位1进行比较：

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

型
这个版本在具有快速位域提取指令的处理器和没有桶形移位器的低端处理器上很有优势，因为只需要右移一个位位置。在乱序处理器上，这种方法提供了比置换表更高的并行度。它似乎更容易适应SIMD实现，因为所有字节通道都使用相同的移位计数。
在我仔细研究相关ASCII字符的二进制编码之前，我研究了一个简单的数学计算，对较小的乘数和除数进行简单的蛮力搜索，得到：

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

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

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

型
然而，在64位平台上，这整个计算可以用64位、32条目的查找表来代替，如果这个64位表可以有效地放入寄存器中，则该查找表可以给予有竞争力的性能，这是x86-64的情况，其中它被作为立即数加载。

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

型
我添加了我的测试框架以供参考：

  1. #include <stdio.h>
  2. #include <stdint.h>
  3. #include <stdlib.h>
  5. #define ORIGINAL_MATH (1)
  7. unsigned int ascii_to_2bit_perm (unsigned int a)
  8. {
  9.     unsigned int perm = (2 << 0) | (1 << 2) | (0 << 4) | (3 << 6);
  10.     return (perm >> (a & 6)) & 3;
  11. }
  13. unsigned int ascii_to_2bit_twiddle (unsigned int a)
  14. {
  15.     return ((a >> 1) & 3) ^ (~a & 2);
  16. }
  18. unsigned int ascii_to_2bit_math (unsigned int a)
  19. {
  20. #if ORIGINAL_MATH
  21.     return ((18 * (a & 31)) % 41) & 3;
  22. #else // ORIGINAL_MATH
  23.     unsigned int t = 18 * (a & 31);
  24.     return (t - ((t * 25) >> 10)) & 3;
  25. #endif // ORIGINAL_MATH
  26. }
  28. unsigned int ascii_to_2bit_tab (unsigned int a)
  29. {
  30.     uint64_t tab = ((0ULL <<  0) | (2ULL <<  2) | (0ULL <<  4) | (1ULL <<  6) |
  31.                     (3ULL <<  8) | (0ULL << 10) | (2ULL << 12) | (3ULL << 14) |
  32.                     (1ULL << 16) | (3ULL << 18) | (0ULL << 20) | (2ULL << 22) |
  33.                     (3ULL << 24) | (1ULL << 26) | (2ULL << 28) | (0ULL << 30) |
  34.                     (1ULL << 32) | (3ULL << 34) | (1ULL << 36) | (2ULL << 38) |
  35.                     (0ULL << 40) | (1ULL << 42) | (3ULL << 44) | (0ULL << 46) |
  36.                     (2ULL << 48) | (0ULL << 50) | (1ULL << 52) | (3ULL << 54) |
  37.                     (0ULL << 56) | (2ULL << 58) | (3ULL << 60) | (1ULL << 62));
  38.     return (tab >> (2 * (a & 31))) & 3;
  39. }
  41. /* Convert nucleobases A, C, G, T represented as either uppercase or lowercase
  42.    ASCII characters into UCSC .2bit-presentation. Passing any values other than
  43.    those corresponding to 'A', 'C', 'G', 'T', 'a', 'c', 'g', 't' results in an
  44.    inderminate response.
  45. */
  46. unsigned int ascii_to_2bit (unsigned int a)
  47. {
  48.     const unsigned int UCSC_2BIT_A = 2;
  49.     const unsigned int UCSC_2BIT_C = 1;
  50.     const unsigned int UCSC_2BIT_G = 3;
  51.     const unsigned int UCSC_2BIT_T = 0;
  53.     switch (a) {
  54.     case 'A':
  55.     case 'a':
  56.         return UCSC_2BIT_A;
  57.         break;
  58.     case 'C':
  59.     case 'c':
  60.         return UCSC_2BIT_C;
  61.         break;
  62.     case 'G':
  63.     case 'g':
  64.         return UCSC_2BIT_G;
  65.         break;
  66.     case 'T':
  67.     case 't':
  68.     default:
  69.         return UCSC_2BIT_T;
  70.         break;
  71.     }
  72. }
  74. int main (void)
  75. {
  76.     char nucleobase[8] = {'A', 'C', 'G', 'T', 'a', 'c', 'g', 't'};
  77.     printf ("Testing permutation variant:\n");
  78.     for (unsigned int i = 0; i < sizeof nucleobase; i++) {
  79.         unsigned int ref = ascii_to_2bit (nucleobase[i]);
  80.         unsigned int res = ascii_to_2bit_perm (nucleobase[i]);
  81.         printf ("i=%2u %c res=%u ref=%u %c\n", 
  82.                 i, nucleobase[i], res, ref, (res == ref) ? 'p' : 'F');
  83.     }
  84.     printf ("Testing bit-twiddling variant:\n");
  85.     for (unsigned int i = 0; i < sizeof nucleobase; i++) {
  86.         unsigned int ref = ascii_to_2bit (nucleobase[i]);
  87.         unsigned int res = ascii_to_2bit_twiddle (nucleobase[i]);
  88.         printf ("i=%2u %c res=%u ref=%u %c\n", 
  89.                 i, nucleobase[i], res, ref, (res == ref) ? 'p' : 'F');
  90.     }
  91.     printf ("Testing math-based variant:\n");
  92.     for (unsigned int i = 0; i < sizeof nucleobase; i++) {
  93.         unsigned int ref = ascii_to_2bit (nucleobase[i]);
  94.         unsigned int res = ascii_to_2bit_math (nucleobase[i]);
  95.         printf ("i=%2u %c res=%u ref=%u %c\n", 
  96.                 i, nucleobase[i], res, ref, (res == ref) ? 'p' : 'F');
  97.     }
  98.     printf ("Testing table-based variant:\n");
  99.     for (unsigned int i = 0; i < sizeof nucleobase; i++) {
  100.         unsigned int ref = ascii_to_2bit (nucleobase[i]);
  101.         unsigned int res = ascii_to_2bit_tab (nucleobase[i]);
  102.         printf ("i=%2u %c res=%u ref=%u %c\n", 
  103.                 i, nucleobase[i], res, ref, (res == ref) ? 'p' : 'F');
  104.     }
  105.     return EXIT_SUCCESS;
  106. }

型

展开查看全部
赞(0）分享  回复(0）举报  2024-01-06
zbq4xfa0

zbq4xfa02#

一种选择是：

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

字符串
其优点在于，它只需要8位，不会溢出，并且不包含非常数移位，因此即使没有特定的SIMD指令，例如，在64位字中输入8个字符，也可以立即用于并行化

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

型
在每个字节的底部留下两个输出位。

展开查看全部
赞(0）分享  回复(0）举报  2024-01-06
wz8daaqr

wz8daaqr3#

下面的例程将用ASCII字符串的内容填充uint32_t值的数组，该ASCII字符串由您指定的字符填充，并将保存状态，以便能够附加第二个，第三个，使用它的方法通过一个main()例程来说明，该例程从命令行获取字符串参数，并将它们传递给TOTAL元素。例程的参数在它上面的注解中描述。

  1. #include <ctype.h>
  2. #include <stdint.h>
  3. #include <stdio.h>
  4. #include <stdlib.h>
  5. #include <string.h>
  7. #define UCSC_2BIT_T  (0)
  8. #define UCSC_2BIT_C  (1)
  9. #define UCSC_2BIT_A  (2)
  10. #define UCSC_2BIT_G  (3)
  12. #define PAIRS_PER_INTEGER   16
  14. /**
  15.  * Converts a string of nucleobases in ASCII form into an array
  16.  * of 32bit integers in 2bitform.  As it should be possible to
  17.  * continue the array of integers, a status reference is
  18.  * maintained in order to determine determine where in the
  19.  * integer to start putting the two bit sequences.  For this,
  20.  * a (*state) variable is maintained (initialize it to 0) to
  21.  * remember where to start putting bitpairs in the array.
  22.  *
  23.  * @param state_ref reference to the state variable to be maintained
  24.  *               with the position on which to put the base in the
  25.  *               array entry.
  26.  * @param dna    string with the ASCII chain of bases.
  27.  * @param twobit array reference to start filling.
  28.  * @param sz     size of the array ***in array entries***.
  29.  * @return the number of array elements written.  This allows to
  30.  *         use a pointer and advance it at each function call
  31.  *         with the number of entries consumed on each call.
  32.  */
  33. ssize_t
  34. dna2twobit(
  35.     int *state_ref,
  36.     char *dna,
  37.     uint32_t twobit[],
  38.     size_t sz)
  39. {
  40.     /* local copy so we only change the state in case of
  41.      * successful execution */
  42.     int state = 30 - *state_ref;
  43.     if (state == 30) *twobit = 0;
  44.     ssize_t total_nb = 0; /* total number of array elements consumed */
  45.     while (*dna && sz) {
  46.         char base = toupper(*dna++);
  47.         uint32_t tb;
  48.         switch (base) {
  49.         case 'A': tb = UCSC_2BIT_A; break;
  50.         case 'C': tb = UCSC_2BIT_C; break;
  51.         case 'T': tb = UCSC_2BIT_T; break;
  52.         case 'G': tb = UCSC_2BIT_G; break;
  53.         default: return -1;
  54.         }
  55.         *twobit |= tb << state;
  56.         state -= 2;
  57.         if (state < 0) {
  58.             --sz; ++twobit;
  59.             state += 32;
  60.             *twobit = 0;
  61.             total_nb++;
  62.         }
  63.     }
  64.     *state_ref = 30 - state;
  65.     return total_nb;
  66. }

字符串
这个函数可以单独链接到任何程序中，但我提供了一个main()代码来说明函数的使用。您可以在命令行参数中使用ASCII编码的实际链来调用程序。程序将一个接一个地对它们进行编码，以演示多重转换（16个碱基适合一个32位整数，正如您在问题中发布的格式定义所述，因此，如果没有编码16的倍数，则状态反映了最后一个中覆盖了多少位。
代码继续执行下面的main函数：

  1. #define TOTAL 16
  3. int main(int argc, char **argv)
  4. {
  5.     int state = 0, i;
  6.     uint32_t twobit_array[TOTAL], *p = twobit_array;
  7.     size_t twobit_cap = TOTAL;
  9.     for (i = 1; i < argc; ++i) {
  10.         printf("Adding the following chain: [%s]\n", argv[i]);
  11.         ssize_t n = dna2twobit(
  12.                         &state,
  13.                         argv[i],
  14.                         p,
  15.                         twobit_cap);
  16.         if (n < 0) {
  17.             fprintf(stderr,
  18.                     "argument invalid: %s\n"
  19.                     "continuing to next\n",
  20.                     argv[i]);
  21.             continue;
  22.         }
  23.         twobit_cap -= n;
  24.         p += n;
  25.     }
  26.     if (!state) --p;
  27.     uint32_t *q = twobit_array;
  28.     size_t off = 0;
  29.     for (int j = 0; q <= p; j++) {
  30.         char *sep = "";
  31.         printf("%09zd: ", off);
  32.         for (int k = 0; k < 4 && q <= p; k++) {
  33.             printf("%s%08x", sep, *q);
  34.             sep = "-";
  35.             q++;
  36.             off += 16;
  37.         }
  38.         printf("\n");
  39.     }
  40.     return 0;
  41. }

型

展开查看全部
赞(0）分享  回复(0）举报  2024-01-06
我来回答

相关问题

    查看更多

    热门标签

    更多
    JavaquerypythonNode开发语言requestUtil数据库Table后端算法LoggerMessageElementParser

    最新问答

    更多
    • 蜀ICP备13028337号-1 大数据知识库 https://www.saoniuhuo.com © All rights reserved
    • 本站内容来源互联网,如果侵犯您的权益请联系我们删除, 联系方式：448109455@qq.com