mem_chain GPU并行化重写
mem_chain 主要分两步，第一步：收集最大精准匹配，第二步：生成种子链，这里考虑收集最大精准匹配实现，假如使用一个线程块协助处理一个seq。
收集最大精准匹配BWA源码实现
static void mem_collect_intv(const mem_opt_t *opt, const bwt_t *bwt, int len, const uint8_t *seq, smem_aux_t *a)
{
int i, k, x = 0, old_n;
int start_width = 1;
int split_len = (int)(opt->min_seed_len * opt->split_factor + .499);
a->mem.n = 0;
// first pass: find all SMEMs
while (x < len) {
if (seq[x] < 4) {
x = bwt_smem1(bwt, len, seq, x, start_width, &a->mem1, a->tmpv);
for (i = 0; i < a->mem1.n; ++i) {
bwtintv_t *p = &a->mem1.a[i];
int slen = (uint32_t)p->info - (p->info>>32); // seed length
if (slen >= opt->min_seed_len)
kv_push(bwtintv_t, a->mem, *p);
}
} else ++x;
}
// second pass: find MEMs inside a long SMEM
old_n = a->mem.n;
for (k = 0; k < old_n; ++k) {
bwtintv_t *p = &a->mem.a[k];
int start = p->info>>32, end = (int32_t)p->info;
if (end - start < split_len || p->x[2] > opt->split_width) continue;
bwt_smem1(bwt, len, seq, (start + end)>>1, p->x[2]+1, &a->mem1, a->tmpv);
for (i = 0; i < a->mem1.n; ++i)
if ((uint32_t)a->mem1.a[i].info - (a->mem1.a[i].info>>32) >= opt->min_seed_len)
kv_push(bwtintv_t, a->mem, a->mem1.a[i]);
}
// third pass: LAST-like
if (opt->max_mem_intv > 0) {
x = 0;
while (x < len) {
if (seq[x] < 4) {
if (1) {
bwtintv_t m;
x = bwt_seed_strategy1(bwt, len, seq, x, opt->min_seed_len, opt->max_mem_intv, &m);
if (m.x[2] > 0) kv_push(bwtintv_t, a->mem, m);
} else { // for now, we never come to this block which is slower
x = bwt_smem1a(bwt, len, seq, x, start_width, opt->max_mem_intv, &a->mem1, a->tmpv);
for (i = 0; i < a->mem1.n; ++i)
kv_push(bwtintv_t, a->mem, a->mem1.a[i]);
}
} else ++x;
}
}
// sort
ks_introsort(mem_intv, a->mem.n, a->mem.a);
}
第一步find all SMEMs看上去并行困难，每次循环内依靠上一次的x，实际上可并行。
bwt_smem1内部调用bwt_smem1a，其实现：给定 BWT 索引、一个查询序列 q 和一个起始位置 x，找出以 q[x] 为起点的 SMEM 匹配（既往前尽量扩展，又往后尽量扩展），并把结果保存到 mem 向量里。
结合bwt_smem1a实现，find all SMEMs实际上是0到len区间的每个正常碱基位置进行左右扩展，bwt_smem1a只是将X位置后的向前扩展合并。因此find all SMEMs实际上是可并行的，从0到len每个碱基位置向后扩展到极限，如果扩展到同一个位置，保留位置更大的（位置小的被包含），再往前扩展一格，去除区间未变化的位置（未变化说明被包含）。
// NOTE: $max_intv is not currently used in BWA-MEM
int bwt_smem1a(const bwt_t *bwt, int len, const uint8_t *q, int x, int min_intv, uint64_t max_intv, bwtintv_v *mem, bwtintv_v *tmpvec[2])
{
bwt_set_intv(bwt, q[x], ik); // the initial interval of a single base
ik.info = x + 1;
for (i = x + 1, curr->n = 0; i < len; ++i) { // forward search
if (ik.x[2] < max_intv) { // an interval small enough
kv_push(bwtintv_t, *curr, ik);
break;
} else if (q[i] < 4) { // an A/C/G/T base
c = 3 - q[i]; // complement of q[i]
bwt_extend(bwt, &ik, ok, 0);
if (ok[c].x[2] != ik.x[2]) { // change of the interval size
kv_push(bwtintv_t, *curr, ik);
if (ok[c].x[2] < min_intv) break; // the interval size is too small to be extended further
}
ik = ok[c]; ik.info = i + 1;
} else { // an ambiguous base
kv_push(bwtintv_t, *curr, ik);
break; // always terminate extension at an ambiguous base; in this case, i<len always stands
}
}
if (i == len) kv_push(bwtintv_t, *curr, ik); // push the last interval if we reach the end
swap = curr; curr = prev; prev = swap;
for (i = x - 1; i >= -1; --i) { // backward search for MEMs
c = i < 0? -1 : q[i] < 4? q[i] : -1; // c==-1 if i<0 or q[i] is an ambiguous base
for (j = 0, curr->n = 0; j < prev->n; ++j) {
bwtintv_t *p = &prev->a[j];
if (c >= 0 && ik.x[2] >= max_intv) bwt_extend(bwt, p, ok, 1);
if (c < 0 || ik.x[2] < max_intv || ok[c].x[2] < min_intv) { // keep the hit if reaching the beginning or an ambiguous base or the intv is small enough
if (curr->n == 0) { // test curr->n>0 to make sure there are no longer matches
if (mem->n == 0 || i + 1 < mem->a[mem->n-1].info>>32) { // skip contained matches
ik = *p; ik.info |= (uint64_t)(i + 1)<<32;
kv_push(bwtintv_t, *mem, ik);
}
} // otherwise the match is contained in another longer match
} else if (curr->n == 0 || ok[c].x[2] != curr->a[curr->n-1].x[2]) {
ok[c].info = p->info;
kv_push(bwtintv_t, *curr, ok[c]);
}
}
if (curr->n == 0) break;
swap = curr; curr = prev; prev = swap;
}
}
第二步 find MEMs inside a long SMEM 大致同步骤1，多了查找需要额外拆减的SMEM的步骤。
第3步 third pass 难以并行,其内部调用bwt_seed_strategy1单向执行
int bwt_seed_strategy1(const bwt_t *bwt, int len, const uint8_t *q, int x, int min_len, int max_intv, bwtintv_t *mem)
{
int i, c;
bwtintv_t ik, ok[4];
memset(mem, 0, sizeof(bwtintv_t));
if (q[x] > 3) return x + 1;
bwt_set_intv(bwt, q[x], ik); // the initial interval of a single base
for (i = x + 1; i < len; ++i) { // forward search
if (q[i] < 4) { // an A/C/G/T base
c = 3 - q[i]; // complement of q[i]
bwt_extend(bwt, &ik, ok, 0);
if (ok[c].x[2] < max_intv && i - x >= min_len) {
*mem = ok[c];
mem->info = (uint64_t)x<<32 | (i + 1);
return i + 1;
}
ik = ok[c];
} else return i + 1;
}
return len;
}
关于bwt_extend(bwt, p, ok, 1)，在查找EME的3个步骤中都会在线程内部频繁调用
其内部主要调用bwt_occ4,该函数主要是在BWT上查找扩展,BWT存储于全局显存，而且也无法加载到共享内存中，通过一个线程块共同协助处理一个seq不会有访问存储上的提升，在生成最大精准匹配的过程中bwt_extend是最为频繁的，即使能够将seq加载到共享内存，由于BWT无法加载到共享内存（甚至无法做到访问内存合并，出现大量随机全局内存访问），其整体上的提升将不会太明显，因此在查找最大精准匹配的过程使用一个线程管理一个seq会更合适，并行CODE会更容易编写，考虑合并内存访问会更好，只有在特别的地方才考虑一个线程块处理一个SEQ。
//bwt_occ4 BWA源码实现
#define __occ_aux4(bwt, b) \
((bwt)->cnt_table[(b)&0xff] + (bwt)->cnt_table[(b)>>8&0xff] \
+ (bwt)->cnt_table[(b)>>16&0xff] + (bwt)->cnt_table[(b)>>24])
void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4])
{
bwtint_t x;
uint32_t *p, tmp, *end;
if (k == (bwtint_t)(-1)) {
memset(cnt, 0, 4 * sizeof(bwtint_t));
return;
}
k -= (k >= bwt->primary); // because $ is not in bwt
p = bwt_occ_intv(bwt, k);
memcpy(cnt, p, 4 * sizeof(bwtint_t));
p += sizeof(bwtint_t); // sizeof(bwtint_t) = 4*(sizeof(bwtint_t)/sizeof(uint32_t))
end = p + ((k>>4) - ((k&~OCC_INTV_MASK)>>4)); // this is the end point of the following loop
for (x = 0; p < end; ++p) x += __occ_aux4(bwt, *p);
tmp = *p & ~((1U<<((~k&15)<<1)) - 1);
x += __occ_aux4(bwt, tmp) - (~k&15);
cnt[0] += x&0xff; cnt[1] += x>>8&0xff; cnt[2] += x>>16&0xff; cnt[3] += x>>24;
}
浙公网安备 33010602011771号