在BWA源码阅读笔记（二）索引文件amb/ann/pac文件是什么?中，我们了解了amb是记录模糊碱基的位置，ann是记录基因组染色体的信息，而pac则是存放压缩后碱基的信息。本篇主要解读bntseq.c中的两行代码, 解读bwa如何实现基因组信息压缩，也就是pac文件的是如何产生的。
#define _set_pac(pac, l, c) ((pac)[(l)>>2] |= (c)<<((~(l)&3)<<1))
#define _get_pac(pac, l) ((pac)[(l)>>2]>>((~(l)&3)<<1)&3)
这两行代码的运算部分全部采用位运算，保证了计算效率
-
>>右移 -
<<左移 -
~: 按位取反 -
&: 按位与 -
|: 按位或
我们先以编码#define _set_pac(pac, l, c) ((pac)[(l)>>2] |= (c)<<((~(l)&3)<<1))为例讲解编码过程
首先对位置按位取反，例如00000000会变成11111111, 00000001会变成11111110
接着x&3相当于求模，例如 7&3 = 7 % 4=1, 那么11111111就变成了00000011, 11111110就变成了00000010
经过这两步，0->3, 1->2, 2->1, 1-> 0, 后续都是3,2,1,0一组的循环。紧接着的<<1是在原先的基础上乘2, 得到了6,4,2,0一组的循环。
c的取值是0,1,2,3, 根据位置的不同每个数字就有不同的表现
- 0: 0, 0, 0, 0
- 1: 64, 16, 4, 1
- 2: 128, 32, 8, 2
- 3: 192, 48, 12 ,3
l>>2是将数字除以4, 保证最终大小就是原来的四分之一，|=将当前编码结果和之前的编码结果进行按位或操作。
举个例子, 第0位编码1，第2位编码3，也就是一开始pac[0]=64,接着pac[0] = 64 | 32 = 96.
接下来讲解#define _get_pac(pac, l) ((pac)[(l)>>2]>>((~(l)&3)<<1)&3)
(pac)[(l)>>2]的操作就是获取压缩后的位置上的编码。通过((~(l)&3)<<1)计算出该编码在二进制的位置。之前是左移操作，此时通过右移，就将对应的二进制位移动到开头前两个位置。通过&3提取这两个位置，就是最终的编码。
根据我的理解，编码和解码的按位取反并没有必要，原来是靠前的移动到最后，去掉按位取反，就是将靠前的防在前面，靠后的移动到后面。 为了验证我的想法，我更改bwa代码中的相应为之后，重新编译，然后测试前后两个版本的bwt索引是否一致。当然最终结果也证明我是对的，并且少一步，速度还会比原来的快一点点。
#3G的人类基因组
[before] Real time: 36.836 sec; CPU: 31.836 sec
[after] Real time: 35.493 sec; CPU: 30.678 sec
同样的逻辑，我自己也写一个函数了，运行效率上会比原来的代码慢，但是会更容易理解一些。
#include <stdio.h>
#define _set_pac(pac, l, c) ((pac)[(l)>>2] |= (c)<<((~(l)&3)<<1))
#define _get_pac(pac, l) ((pac)[(l)>>2]>>((~(l)&3)<<1)&3)
typedef unsigned int uint8_t;
uint8_t encode(uint8_t c, int p){
int reminder = (p % 4)*2;
c = c << reminder;
return c;
}
uint8_t decode(uint8_t *arr, int p){
uint8_t c = arr[p >> 2];
p = (p % 4)*2;
return (c >> p)&3;
}
void test1(){
uint8_t orgin[12] = {1,2,3,3,0,1,2,1,0,3,1,2};
uint8_t res[3] = {0};
for (int i = 0; i < 12; i++){
printf("%u ",orgin[i]);
}
printf("\n");
for (int i = 0; i < 12; i++){
res[i>>2] |= encode(orgin[i],i);
}
for (int i = 0; i < 3; i++){
printf("%u ", res[i]);
}
printf("\n");
for (int i = 0; i < 12; i++){
printf("%u ", decode(res, i));
}
printf("\n");
}
void test2(){
uint8_t orgin[12] = {1,2,3,3,0,1,2,1,0,3,1,2};
uint8_t res[3] = {0};
for (int i = 0; i < 12; i++){
printf("%u ",orgin[i]);
}
printf("\n");
for (int i = 0; i < 12; i++){
_set_pac(res, i, orgin[i]);
}
for (int i = 0; i < 3; i++){
printf("%u ", res[i]);
}
printf("\n");
for (int i = 0; i < 12; i++){
printf("%u ", _get_pac(res, i));
}
printf("\n");
}
int main(int argc, char *argv[]){
test1();
printf("--------------------\n");
test2();
return 0;
}
网友评论