BWA源码阅读笔记（三）索引文件pac如何实现信息压缩

在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;

}