zzh
/
BWA-FastAlign
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

结果正确了，gen-sam过程应该还有性能提升空间

This commit is contained in:
zzh 2026-01-17 01:00:43 +08:00
parent 67cbb5b8e9
commit e335e79543
3 changed files with 56 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Makefile
View File

@ -1,5 +1,5 @@
CC= gcc
CFLAGS= -g -Wall -Wno-unused-function -mavx2 -mavx512bw #-O3
CFLAGS= -g -Wall -Wno-unused-function -mavx2 -mavx512bw -O3
WRAP_MALLOC=-DUSE_MALLOC_WRAPPERS
SHOW_PERF= -DSHOW_PERF

72
paired_sam.c
View File

@ -363,11 +363,11 @@ static void worker_calc_matesw_avx512_u8(void* data, int idx, int tid) {
static void worker_calc_matesw_avx512_i16(void* data, int i, int tid) {}
// 检测添加新的align
static int check_add_align(kswr_avx_t aln, int min_seed_len, int is_rev, int64_t l_pac, mem_alnreg_t* a, int l_ms, uint8_t* ms, mem_alnreg_v* ma,
int64_t rb) {
static int check_add_align(const mem_opt_t* opt, kswr_avx_t aln, int is_rev, int64_t l_pac, mem_alnreg_t* a, int l_ms, uint8_t* ms,
mem_alnreg_v* ma, int64_t rb) {
int res = 0;
if (aln.score >= min_seed_len && aln.qb >= 0) { // something goes wrong if aln.qb < 0
// int i, tmp;
if (aln.score >= opt->min_seed_len && aln.qb >= 0) { // something goes wrong if aln.qb < 0
int i, tmp;
mem_alnreg_t b;
memset(&b, 0, sizeof(mem_alnreg_t));
b.rid = a->rid;
@ -383,14 +383,15 @@ static int check_add_align(kswr_avx_t aln, int min_seed_len, int is_rev, int64_t
// printf("*** %d, [%lld,%lld], %d:%d, (%lld,%lld), (%lld,%lld) == (%lld,%lld)\n", aln.score, rb, re, is_rev,
// is_larger, a->rb, a->re, ma->a[0].rb, ma->a[0].re, b.rb, b.re);
kv_push(mem_alnreg_t, *ma, b); // make room for a new element
// // move b s.t. ma is sorted
// for (i = 0; i < ma->n - 1; ++i) // find the insertion point
// if (ma->a[i].score < b.score)
// break;
// tmp = i;
// for (i = ma->n - 1; i > tmp; --i) ma->a[i] = ma->a[i - 1];
// ma->a[i] = b;
// move b s.t. ma is sorted
for (i = 0; i < ma->n - 1; ++i) // find the insertion point
if (ma->a[i].score < b.score)
break;
tmp = i;
for (i = ma->n - 1; i > tmp; --i) ma->a[i] = ma->a[i - 1];
ma->a[i] = b;
res = 1;
ma->n = mem_sort_dedup_patch(opt, 0, 0, 0, ma->n, ma->a);
}
return res;
}
@ -572,33 +573,48 @@ static void workder_gen_sam(void* data, int idx, int tid) {
endIdx = w->n_reads;
int id = 0, i = 0, j = 0, k = 0;
mem_alnreg_v bb[2];
kv_init(bb[0]); kv_init(bb[1]);
for (id = startIdx; id < endIdx; id += 2) {
// 初始化变量
bseq1_t* s = &w->seqs[id];
mem_alnreg_v* a = &w->regs[id];
seq_sam_t* ss = &w->sams[id];
int cmp = strcmp("ERR194147.42331121", s[0].name) == 0;
#if 1
int cmp = strcmp("ERR194147.1629456", s[0].name) == 0;
// int cmp = (id + w->n_processed) == 202924;
if (cmp == 1) {
fprintf(stderr, "%s\n", s[0].name);
fflush(stderr);
//exit(0);
}
#endif
int n[2] = {0};
// , nn[2] = {0};
//
_clear_vec(bb[0]);
_clear_vec(bb[1]);
for (i = 0; i < 2; ++i)
for (j = 0; j < s[!i].msw.n; ++j) {
msw_task_t* t = s[!i].msw.a[j];
kv_push(mem_alnreg_t, bb[i], a[i].a[t->aj]);
}
int n[2] = {0}, nn[2] = {0};
for (i = 0; i < 2; ++i) {
int si = !i;
int origin_n = a[si].n;
// 这里应该先给task排序，因为u8和i16是分开排序的，需要合在一起
for (k = 0; k < s[si].msw.n; ++k) {
msw_task_t* t = s[si].msw.a[k];
mem_alnreg_t* b = &a[i].a[t->aj];
// mem_alnreg_t* b = &a[i].a[t->aj];
mem_alnreg_t* b = &bb[i].a[k];
mem_alnreg_v* ma = &a[si];
uint8_t skip = 0;
if (n[si]) { // 添加了新的aln，需要检测skip
for (j = origin_n; j < ma->n; ++j) { // check which orinentation has been found
for (j = 0; j < ma->n; ++j) { // check which orinentation has been found
int64_t dist;
int r;
r = mem_infer_dir(bns->l_pac, b->rb, ma->a[j].rb, &dist);
@ -606,18 +622,25 @@ static void workder_gen_sam(void* data, int idx, int tid) {
skip |= 1 << r;
}
}
//if ((t->skip | skip) != 15) { // 检查新添加的aln是不是把当前的任务skip掉了
kswr_avx_t aln = t->aln; // fprintf(gf[0], "id-%ld score-%d te-%d qe-%d score2-%d te2-%d tb-%d qb-%d\n", s[si].id + w->n_processed, aln.score, aln.te, aln.qe, aln.score2, aln.te2, aln.tb, aln.qb);
if (cmp) {
fprintf(stderr, "id-%ld score-%d te-%d qe-%d score2-%d te2-%d tb-%d qb-%d\n", s[si].id + w->n_processed, aln.score, aln.te, aln.qe, aln.score2, aln.te2, aln.tb, aln.qb);
fflush(stderr);
}
nn[si] += 1;
n[si] += check_add_align(t->aln, opt->min_seed_len, t->is_rev, bns->l_pac, b, s[si].l_seq, (uint8_t*)s[si].seq, ma, t->rb);
// if (skip == 15)
if ((t->skip | skip) == 15)
continue;
// if ((t->skip | skip) != 15) { // 检查新添加的aln是不是把当前的任务skip掉了
#if 1
kswr_avx_t aln = t->aln; // fprintf(gf[0], "id-%ld score-%d te-%d qe-%d score2-%d te2-%d tb-%d qb-%d\n", s[si].id +
// w->n_processed, aln.score, aln.te, aln.qe, aln.score2, aln.te2, aln.tb, aln.qb);
if (cmp) {
fprintf(stderr, "id-%ld score-%d te-%d qe-%d score2-%d te2-%d tb-%d qb-%d\n", s[si].id + w->n_processed, aln.score, aln.te,
aln.qe, aln.score2, aln.te2, aln.tb, aln.qb); fflush(stderr);
}
#endif
// nn[si] += 1;
n[si] += check_add_align(opt, t->aln, t->is_rev, bns->l_pac, b, s[si].l_seq, (uint8_t*)s[si].seq, ma, t->rb);
//}
}
}
// 处理完2个pair read之后再排序，因为上面操作有插入，排序后之前记录的索引就无效了，所以上面不能排序，要两个pair read处理完之后再排序
#if 0
for (i = 0; i < 2; ++i) {
int tmp;
if (nn[i] > 0) {
@ -645,6 +668,7 @@ static void workder_gen_sam(void* data, int idx, int tid) {
// a[i].n = mem_sort_dedup_patch(opt, 0, 0, 0, a[i].n, a[i].a);
}
}
#endif
////////////////////////////////////////////////////////////
// 这里需要传入全局id

14
run.sh
View File

@ -1,11 +1,11 @@
thread=1
thread=64
make clean; make -j 32
#n1=~/data/dataset/D2/n1.fq.gz
#n2=~/data/dataset/D2/n2.fq.gz
n1=~/data/dataset/D2/d1.fq
n2=~/data/dataset/D2/d2.fq
n1=~/data/dataset/D2/n1.fq.gz
n2=~/data/dataset/D2/n2.fq.gz
#n1=~/data/dataset/D2/d1.fq
#n2=~/data/dataset/D2/d2.fq
#n1=/file-data/zzh/dataset/D1/45mr1.fq.gz
#n2=/file-data/zzh/dataset/D1/45mr2.fq.gz
@ -14,8 +14,8 @@ reference=~/data/reference/fmt/human_g1k_v37_decoy.fasta
#out=/dev/null
#out=./oldsam-D2-out.sam
#out=./D2-1.sam
out=./d.sam
out=./D2-1.sam
#out=./d.sam
prog=./fastalign
#prog=/home/zzh/fastbwa