From 463f7da13805a91cd7d760acf7111db4ce4e4d00 Mon Sep 17 00:00:00 2001 From: zzh Date: Wed, 7 Feb 2024 22:08:51 +0800 Subject: [PATCH] =?UTF-8?q?=E5=B0=86smem1=E5=87=BD=E6=95=B0=E7=94=A8fmt?= =?UTF-8?q?=E7=BB=93=E6=9E=84=E5=AE=9E=E7=8E=B0=E4=BA=86=EF=BC=8C=E7=BB=93?= =?UTF-8?q?=E6=9E=9C=E5=9F=BA=E6=9C=AC=E6=AD=A3=E7=A1=AE?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .gitignore | 1 + .vscode/launch.json | 8 +- .vscode/settings.json | 11 +- bwa.c | 59 +++- bwa.h | 2 +- bwamem.c | 45 ++- bwamem.h | 3 +- bwt.c | 42 ++- fastmap.c | 21 +- fmt_idx.c | 739 +++++++++++++++++++++++++++++++++++++++++- fmt_idx.h | 116 +++++-- run.sh | 19 +- utils.h | 5 + 13 files changed, 1006 insertions(+), 65 deletions(-) diff --git a/.gitignore b/.gitignore index 0dc1011..035c401 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ *.[oa] +*.txt bwa test test64 diff --git a/.vscode/launch.json b/.vscode/launch.json index 54cb5f4..39dbc60 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -13,13 +13,13 @@ "args": [ "mem", "-t", - "12", + "1", "-M", "-R", "'@RG\\tID:normal\\tSM:normal\\tPL:illumina\\tLB:normal\\tPG:bwa'", - "/home/zzh/data/reference/human_g1k_v37_decoy.fasta", - "/home/zzh/data/fastq/n_s1.fq", - "/home/zzh/data/fastq/n_s2.fq", + "~/reference/human_g1k_v37_decoy.fasta", + "~/fastq/sn_r1.fq", + "~/fastq/sn_r2.fq", "-o", "/dev/null" ], diff --git a/.vscode/settings.json b/.vscode/settings.json index 1728226..77e1001 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -1,5 +1,14 @@ { "files.associations": { - "random": "c" + "random": "c", + "bwt.h": "c", + "bwa.h": "c", + "*.tcc": "c", + "functional": "c", + "string_view": "c", + "istream": "c", + "limits": "c", + "bit": "c", + "numeric": "c" } } \ No newline at end of file diff --git a/bwa.c b/bwa.c index 347e8a9..ed5f771 100644 --- a/bwa.c +++ b/bwa.c @@ -28,6 +28,7 @@ #include #include #include +#include #include "bntseq.h" #include "bwa.h" #include "ksw.h" @@ -248,9 +249,10 @@ char *bwa_idx_infer_prefix(const char *hint) int l_hint; FILE *fp; l_hint = strlen(hint); - prefix = malloc(l_hint + 3 + 4 + 1); + prefix = malloc(l_hint + 3 + 4 + 1 + 10); strcpy(prefix, hint); - strcpy(prefix + l_hint, ".64.bwt"); + //strcpy(prefix + l_hint, ".64.bwt"); + strcpy(prefix + l_hint, ".ne.bwt"); if ((fp = fopen(prefix, "rb")) != 0) { fclose(fp); prefix[l_hint + 3] = 0; @@ -277,15 +279,59 @@ bwt_t *bwa_idx_load_bwt(const char *hint) if (bwa_verbose >= 1) fprintf(stderr, "[E::%s] fail to locate the index files\n", __func__); return 0; } + fprintf(stderr, "zzh-1\n"); tmp = calloc(strlen(prefix) + 5, 1); strcat(strcpy(tmp, prefix), ".bwt"); // FM-index bwt = bwt_restore_bwt(tmp); + fprintf(stderr, "zzh-1\n"); strcat(strcpy(tmp, prefix), ".33.4.sa"); // partial suffix array (SA) bwt_restore_sa(tmp, bwt); + fprintf(stderr, "zzh-after-sa\n"); free(tmp); free(prefix); return bwt; } +FMTIndex *bwa_idx_load_fmt(const char *hint) +{ + char *fmt_idx_fn, *kmer_idx_fn, *sa_fn; + FMTIndex *fmt; + char suffix[32]; + int l_hint = strlen(hint); + fmt_idx_fn = malloc(l_hint + 32); + kmer_idx_fn = malloc(l_hint + 32); + sa_fn = malloc(l_hint + 32); + sprintf(suffix, ".256.%d.fmt", FMT_MID_INTERVAL); + strcpy(fmt_idx_fn, hint); + strcpy(fmt_idx_fn + l_hint, suffix); + sprintf(suffix, ".%d.kmer", KMER_LEN); + strcpy(kmer_idx_fn, hint); + strcpy(kmer_idx_fn + l_hint, suffix); + + if (access(fmt_idx_fn, F_OK) != 0 || access(kmer_idx_fn, F_OK) != 0) + { + if (bwa_verbose >= 1) + fprintf(stderr, "[E::%s] fail to locate the index files\n", __func__); + return 0; + } + fprintf(stderr, "zzh-fmt-1\n"); + fmt = fmt_restore_fmt(fmt_idx_fn); + fprintf(stderr, "%s\n", kmer_idx_fn); + fmt->kmer_entry = fmt_restore_kmer_idx(kmer_idx_fn); + + fprintf(stderr, "zzh-fmt-2\n"); + + strcpy(sa_fn, hint); + sprintf(suffix, ".33.%d.sa", SA_INTV); + strcpy(sa_fn + l_hint, suffix); // partial suffix array (SA) +// fmt_restore_sa(sa_fn, fmt); + + free(fmt_idx_fn); + free(kmer_idx_fn); + free(sa_fn); + + return fmt; +} + bwaidx_t *bwa_idx_load_from_disk(const char *hint, int which) { bwaidx_t *idx; @@ -297,7 +343,14 @@ bwaidx_t *bwa_idx_load_from_disk(const char *hint, int which) } idx = calloc(1, sizeof(bwaidx_t)); if (which & BWA_IDX_BWT) idx->bwt = bwa_idx_load_bwt(hint); - if (which & BWA_IDX_BNS) { + if (which & BWA_IDX_BWT) { + idx->fmt = bwa_idx_load_fmt(hint); + idx->fmt->sa = idx->bwt->sa; + idx->fmt->n_sa = idx->bwt->n_sa; + idx->fmt->sa_intv = idx->bwt->sa_intv; + } + if (which & BWA_IDX_BNS) + { int i, c; idx->bns = bns_restore(prefix); for (i = c = 0; i < idx->bns->n_seqs; ++i) diff --git a/bwa.h b/bwa.h index 90e18d2..aa155b6 100644 --- a/bwa.h +++ b/bwa.h @@ -48,7 +48,7 @@ typedef struct { bwt_t *bwt; // FM-index - bwtd_t *bwtd;// FMT-index + FMTIndex *fmt;// FMT-index bntseq_t *bns; // information on the reference sequences uint8_t *pac; // the actual 2-bit encoded reference sequences with 'N' converted to a random base diff --git a/bwamem.c b/bwamem.c index 80e801b..448bc4d 100644 --- a/bwamem.c +++ b/bwamem.c @@ -41,6 +41,7 @@ #include "kvec.h" #include "ksort.h" #include "utils.h" +#include "fmt_idx.h" #ifdef USE_MALLOC_WRAPPERS # include "malloc_wrap.h" @@ -137,32 +138,49 @@ static void smem_aux_destroy(smem_aux_t *a) free(a); } -static void mem_collect_intv(const mem_opt_t *opt, const bwt_t *bwt, int len, const uint8_t *seq, smem_aux_t *a) +static void mem_collect_intv(const mem_opt_t *opt, const bwt_t *bwt, const FMTIndex *fmt, 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); + int max_seed_len = 0; a->mem.n = 0; + // char *BASE = "ACGT"; + // for (i = 0; i < len; ++i) + // if (seq[i] < 4) + // fprintf(stderr, "%c", BASE[seq[i]]); + // else + // fprintf(stderr, "N"); + // fprintf(stderr, "\n"); + // first pass: find all SMEMs + //fprintf(fp1, "seq: %ld\n", dn++); + // dn ++; while (x < len) { if (seq[x] < 4) { #ifdef SHOW_PERF int64_t tmp_time = realtime_msec(); #endif - x = bwt_smem1(bwt, len, seq, x, start_width, &a->mem1, a->tmpv); + //x = bwt_smem1(bwt, len, seq, x, start_width, &a->mem1, a->tmpv); + x = fmt_smem(fmt, len, seq, x, start_width, &a->mem1, a->tmpv); #ifdef SHOW_PERF tmp_time = realtime_msec() - tmp_time; __sync_fetch_and_add(&time_bwt_smem1a, tmp_time); #endif 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 + //fprintf(fp1, "%ld %ld %ld %ld %d\n", p->x[0], p->x[1], p->x[2], p->info >> 32, (uint32_t)p->info); + int slen = (uint32_t)p->info - (p->info >> 32); // seed length + max_seed_len = fmax(max_seed_len, slen); if (slen >= opt->min_seed_len) kv_push(bwtintv_t, a->mem, *p); } } else ++x; } // second pass: find MEMs inside a long SMEM + //if (max_seed_len == len) + goto collect_intv_end; + old_n = a->mem.n; for (k = 0; k < old_n; ++k) { bwtintv_t *p = &a->mem.a[k]; @@ -204,6 +222,8 @@ static void mem_collect_intv(const mem_opt_t *opt, const bwt_t *bwt, int len, co } else ++x; } } + +collect_intv_end: // sort ks_introsort(mem_intv, a->mem.n, a->mem.a); } @@ -295,7 +315,7 @@ void mem_print_chain(const bntseq_t *bns, mem_chain_v *chn) } } -mem_chain_v mem_chain(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t *bns, int len, const uint8_t *seq, void *buf) +mem_chain_v mem_chain(const mem_opt_t *opt, const bwt_t *bwt, const FMTIndex *fmt, const bntseq_t *bns, int len, const uint8_t *seq, void *buf) { int i, b, e, l_rep; int64_t l_pac = bns->l_pac; @@ -308,7 +328,7 @@ mem_chain_v mem_chain(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t *bn tree = kb_init(chn, KB_DEFAULT_SIZE); aux = buf? (smem_aux_t*)buf : smem_aux_init(); - mem_collect_intv(opt, bwt, len, seq, aux); + mem_collect_intv(opt, bwt, fmt, len, seq, aux); for (i = 0, b = e = l_rep = 0; i < aux->mem.n; ++i) { // compute frac_rep bwtintv_t *p = &aux->mem.a[i]; int sb = (p->info>>32), se = (uint32_t)p->info; @@ -1120,7 +1140,7 @@ void mem_reg2sam(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, } } -mem_alnreg_v mem_align1_core(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t *bns, const uint8_t *pac, int l_seq, char *seq, void *buf) +mem_alnreg_v mem_align1_core(const mem_opt_t *opt, const bwt_t *bwt, const FMTIndex *fmt, const bntseq_t *bns, const uint8_t *pac, int l_seq, char *seq, void *buf) { int i; mem_chain_v chn; @@ -1129,7 +1149,7 @@ mem_alnreg_v mem_align1_core(const mem_opt_t *opt, const bwt_t *bwt, const bntse for (i = 0; i < l_seq; ++i) // convert to 2-bit encoding if we have not done so seq[i] = seq[i] < 4? seq[i] : nst_nt4_table[(int)seq[i]]; - chn = mem_chain(opt, bwt, bns, l_seq, (uint8_t*)seq, buf); + chn = mem_chain(opt, bwt, fmt, bns, l_seq, (uint8_t*)seq, buf); chn.n = mem_chain_flt(opt, chn.n, chn.a); mem_flt_chained_seeds(opt, bns, pac, l_seq, (uint8_t*)seq, chn.n, chn.a); if (bwa_verbose >= 4) mem_print_chain(bns, &chn); @@ -1233,6 +1253,7 @@ mem_aln_t mem_reg2aln(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t * typedef struct { const mem_opt_t *opt; const bwt_t *bwt; + const FMTIndex *fmt; const bntseq_t *bns; const uint8_t *pac; const mem_pestat_t *pes; @@ -1247,12 +1268,12 @@ static void worker1(void *data, int i, int tid) worker_t *w = (worker_t*)data; if (!(w->opt->flag&MEM_F_PE)) { if (bwa_verbose >= 4) printf("=====> Processing read '%s' <=====\n", w->seqs[i].name); - w->regs[i] = mem_align1_core(w->opt, w->bwt, w->bns, w->pac, w->seqs[i].l_seq, w->seqs[i].seq, w->aux[tid]); + w->regs[i] = mem_align1_core(w->opt, w->bwt, w->fmt, w->bns, w->pac, w->seqs[i].l_seq, w->seqs[i].seq, w->aux[tid]); } else { if (bwa_verbose >= 4) printf("=====> Processing read '%s'/1 <=====\n", w->seqs[i<<1|0].name); - w->regs[i<<1|0] = mem_align1_core(w->opt, w->bwt, w->bns, w->pac, w->seqs[i<<1|0].l_seq, w->seqs[i<<1|0].seq, w->aux[tid]); + w->regs[i<<1|0] = mem_align1_core(w->opt, w->bwt, w->fmt, w->bns, w->pac, w->seqs[i<<1|0].l_seq, w->seqs[i<<1|0].seq, w->aux[tid]); if (bwa_verbose >= 4) printf("=====> Processing read '%s'/2 <=====\n", w->seqs[i<<1|1].name); - w->regs[i<<1|1] = mem_align1_core(w->opt, w->bwt, w->bns, w->pac, w->seqs[i<<1|1].l_seq, w->seqs[i<<1|1].seq, w->aux[tid]); + w->regs[i<<1|1] = mem_align1_core(w->opt, w->bwt, w->fmt, w->bns, w->pac, w->seqs[i<<1|1].l_seq, w->seqs[i<<1|1].seq, w->aux[tid]); } } @@ -1274,7 +1295,7 @@ static void worker2(void *data, int i, int tid) } } -void mem_process_seqs(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t *bns, const uint8_t *pac, int64_t n_processed, int n, bseq1_t *seqs, const mem_pestat_t *pes0) +void mem_process_seqs(const mem_opt_t *opt, const bwt_t *bwt, const FMTIndex *fmt, const bntseq_t *bns, const uint8_t *pac, int64_t n_processed, int n, bseq1_t *seqs, const mem_pestat_t *pes0) { extern void kt_for(int n_threads, void (*func)(void*,int,int), void *data, int n); worker_t w; @@ -1287,7 +1308,7 @@ void mem_process_seqs(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t *bn w.regs = malloc(n * sizeof(mem_alnreg_v)); w.opt = opt; w.bwt = bwt; w.bns = bns; w.pac = pac; w.seqs = seqs; w.n_processed = n_processed; - w.pes = &pes[0]; + w.pes = &pes[0]; w.fmt = fmt; w.aux = malloc(opt->n_threads * sizeof(smem_aux_t)); for (i = 0; i < opt->n_threads; ++i) w.aux[i] = smem_aux_init(); diff --git a/bwamem.h b/bwamem.h index 0a0e3bb..6f92c5e 100644 --- a/bwamem.h +++ b/bwamem.h @@ -30,6 +30,7 @@ #include "bwt.h" #include "bntseq.h" #include "bwa.h" +#include "fmt_idx.h" #define MEM_MAPQ_COEF 30.0 #define MEM_MAPQ_MAX 60 @@ -158,7 +159,7 @@ extern "C" { * @param pes0 insert-size info; if NULL, infer from data; if not NULL, it should be an array with 4 elements, * corresponding to each FF, FR, RF and RR orientation. See mem_pestat() for more info. */ - void mem_process_seqs(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t *bns, const uint8_t *pac, int64_t n_processed, int n, bseq1_t *seqs, const mem_pestat_t *pes0); + void mem_process_seqs(const mem_opt_t *opt, const bwt_t *bwt, const FMTIndex *fmt, const bntseq_t *bns, const uint8_t *pac, int64_t n_processed, int n, bseq1_t *seqs, const mem_pestat_t *pes0); /** * Find the aligned regions for one query sequence diff --git a/bwt.c b/bwt.c index aa41749..be958b0 100644 --- a/bwt.c +++ b/bwt.c @@ -314,6 +314,11 @@ int bwt_match_exact_alt(const bwt_t *bwt, int len, const ubyte_t *str, bwtint_t void bwt_extend(const bwt_t *bwt, const bwtintv_t *ik, bwtintv_t ok[4], int is_back) { +#ifdef SHOW_PERF +#if 1 + int64_t tmp_time = realtime_msec(); +#endif +#endif bwtint_t tk[4], tl[4]; int i; bwt_2occ4(bwt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], tk, tl); @@ -325,6 +330,12 @@ void bwt_extend(const bwt_t *bwt, const bwtintv_t *ik, bwtintv_t ok[4], int is_b ok[2].x[is_back] = ok[3].x[is_back] + ok[3].x[2]; ok[1].x[is_back] = ok[2].x[is_back] + ok[2].x[2]; ok[0].x[is_back] = ok[1].x[is_back] + ok[1].x[2]; +#ifdef SHOW_PERF +#if 1 + tmp_time = realtime_msec() - tmp_time; + __sync_fetch_and_add(&time_bwt_extend, tmp_time); +#endif +#endif } static void bwt_reverse_intvs(bwtintv_v *p) @@ -355,12 +366,30 @@ int bwt_smem1a(const bwt_t *bwt, int len, const uint8_t *q, int x, int min_intv, 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 + // int flag = 0; + // int small_intv_qlen = 0; + // int small_intv_val = 0; + // int intv_19 = 0; + for (i = x + 1, curr->n = 0; i < len; ++i) + { // forward search + //if (ik.x[2] <= 2) + //{ + // if (!flag) { + // small_intv_qlen = i - x; + // small_intv_val = ik.x[2]; + // flag = 1; + // } + //} + //if (i-x == 19) + //{ + // intv_19 = ik.x[2]; + //} 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); @@ -372,6 +401,12 @@ int bwt_smem1a(const bwt_t *bwt, int len, const uint8_t *q, int x, int min_intv, break; // always terminate extension at an ambiguous base; in this case, ia[0].info; // this will be the returned value @@ -462,7 +497,8 @@ void bwt_dump_sa(const char *fn, const bwt_t *bwt) static bwtint_t fread_fix(FILE *fp, bwtint_t size, void *a) { // Mac/Darwin has a bug when reading data longer than 2GB. This function fixes this issue by reading data in small chunks - const int bufsize = 0x1000000; // 16M block + //const int bufsize = 0x1000000; // 16M block + const int bufsize = 0x4000000; //64M block bwtint_t offset = 0; while (size) { int x = bufsize < size? bufsize : size; @@ -488,7 +524,7 @@ void bwt_restore_sa(const char *fn, bwt_t *bwt) bwt->n_sa = (bwt->seq_len + bwt->sa_intv) / bwt->sa_intv; bwt->sa = (uint8_t*)malloc(SA_BYTES(bwt->n_sa)); - + fprintf(stderr, "zzh-read-sa %ld\n", SA_BYTES(bwt->n_sa)); fread_fix(fp, SA_BYTES(bwt->n_sa), bwt->sa); err_fclose(fp); } diff --git a/fastmap.c b/fastmap.c index 1bf0b61..d31340b 100644 --- a/fastmap.c +++ b/fastmap.c @@ -48,10 +48,15 @@ int64_t time_ksw_extend2 = 0, time_ksw_global2 = 0, time_ksw_align2 = 0, time_bwt_smem1a = 0, + time_bwt_extend = 0, time_bwt_occ4 = 0, time_bwt_sa = 0, time_bwt_sa_read = 0; +int64_t dn = 0; + +FILE *fp1; + #endif extern unsigned char nst_nt4_table[256]; @@ -110,18 +115,19 @@ static void *process(void *shared, int step, void *_data) fprintf(stderr, "[M::%s] %d single-end sequences; %d paired-end sequences\n", __func__, n_sep[0], n_sep[1]); if (n_sep[0]) { tmp_opt.flag &= ~MEM_F_PE; - mem_process_seqs(&tmp_opt, idx->bwt, idx->bns, idx->pac, aux->n_processed, n_sep[0], sep[0], 0); + mem_process_seqs(&tmp_opt, idx->bwt, idx->fmt, idx->bns, idx->pac, aux->n_processed, n_sep[0], sep[0], 0); for (i = 0; i < n_sep[0]; ++i) data->seqs[sep[0][i].id].sam = sep[0][i].sam; } if (n_sep[1]) { tmp_opt.flag |= MEM_F_PE; - mem_process_seqs(&tmp_opt, idx->bwt, idx->bns, idx->pac, aux->n_processed + n_sep[0], n_sep[1], sep[1], aux->pes0); + mem_process_seqs(&tmp_opt, idx->bwt, idx->fmt, idx->bns, idx->pac, aux->n_processed + n_sep[0], n_sep[1], sep[1], aux->pes0); for (i = 0; i < n_sep[1]; ++i) data->seqs[sep[1][i].id].sam = sep[1][i].sam; } free(sep[0]); free(sep[1]); - } else mem_process_seqs(opt, idx->bwt, idx->bns, idx->pac, aux->n_processed, data->n_seqs, data->seqs, aux->pes0); + } else + mem_process_seqs(opt, idx->bwt, idx->fmt, idx->bns, idx->pac, aux->n_processed, data->n_seqs, data->seqs, aux->pes0); aux->n_processed += data->n_seqs; return data; } else if (step == 2) { @@ -154,6 +160,8 @@ static void update_a(mem_opt_t *opt, const mem_opt_t *opt0) int main_mem(int argc, char *argv[]) { + fp1 = fopen("./test_out.txt", "w"); + mem_opt_t *opt, opt0; int fd, fd2, i, c, ignore_alt = 0, no_mt_io = 0; int fixed_chunk_size = -1; @@ -372,7 +380,7 @@ int main_mem(int argc, char *argv[]) } } else update_a(opt, &opt0); bwa_fill_scmat(opt->a, opt->b, opt->mat); - + fprintf(stderr, "zzh-1\n"); aux.idx = bwa_idx_load_from_shm(argv[optind]); if (aux.idx == 0) { if ((aux.idx = bwa_idx_load(argv[optind], BWA_IDX_ALL)) == 0) return 1; // FIXME: memory leak @@ -387,6 +395,7 @@ int main_mem(int argc, char *argv[]) if (bwa_verbose >= 1) fprintf(stderr, "[E::%s] fail to open file `%s'.\n", __func__, argv[optind + 1]); return 1; } + fprintf(stderr, "zzh-2\n"); fp = gzdopen(fd, "r"); aux.ks = kseq_init(fp); if (optind + 2 < argc) { @@ -405,7 +414,6 @@ int main_mem(int argc, char *argv[]) } } - BuildBwtdFromBwt(aux.idx->bwt, &aux.idx->bwtd); bwa_print_sam_hdr(aux.idx->bns, hdr_line); aux.actual_chunk_size = fixed_chunk_size > 0? fixed_chunk_size : opt->chunk_size * opt->n_threads; @@ -423,10 +431,13 @@ int main_mem(int argc, char *argv[]) #ifdef SHOW_PERF fprintf(stderr, "\n"); fprintf(stderr, "time_bwt_smem1a: %f s\n", time_bwt_smem1a / 1000.0 / opt->n_threads); + fprintf(stderr, "time_bwt_extend: %f s\n", time_bwt_extend / 1000.0 / opt->n_threads); fprintf(stderr, "time_bwt_sa: %f s\n", time_bwt_sa / 1000.0 / opt->n_threads); fprintf(stderr, "time_ksw_extend2: %f s\n", time_ksw_extend2 / 1000.0 / opt->n_threads); fprintf(stderr, "time_bwt_sa_read: %f s\n", time_bwt_sa_read / 1000.0 / opt->n_threads); fprintf(stderr, "\n"); + + fclose(fp1); #endif return 0; diff --git a/fmt_idx.c b/fmt_idx.c index 52fe021..202caf2 100644 --- a/fmt_idx.c +++ b/fmt_idx.c @@ -9,19 +9,740 @@ Date : 2023/12/24 #include #include +#include + #include "fmt_idx.h" +#include "utils.h" +#include "bntseq.h" +#include "kvec.h" -// 创建fmt-index索引数据 -void BuildBwtdFromBwt(bwt_t *bwt, bwtd_t **bwtd_p) { - *bwtd_p = (bwtd_t *)malloc(sizeof(bwtd_t)); - bwtd_t *bwtd = *bwtd_p; +const static char BASE[4] = {'A', 'C', 'G', 'T'}; - const int baseLen = 12; - const int kmerSize = 1 << (baseLen << 1); +static bwtint_t fread_fix(FILE *fp, bwtint_t size, void *a) +{ // Mac/Darwin has a bug when reading data longer than 2GB. This function fixes this issue by reading data in small chunks + const int bufsize = 0x4000000; // 16 * 4M block + bwtint_t offset = 0; + while (size) + { + int x = bufsize < size ? bufsize : size; + if ((x = err_fread_noeof(a + offset, 1, x, fp)) == 0) + break; + size -= x; + offset += x; + } + return offset; +} - bwtd->kmer_range = (kmer_range_t *)malloc(kmerSize * sizeof(kmer_range_t)); +// 生成occ，每个字节对应一个pattern +void fmt_gen_cnt_occ(FMTIndex *fmt) +{ + // 0-8：大于a的occ，8-16：大于b的occ，16-24：b的occ + int i, a, b, ti; + uint32_t oa, ooa, ob, oob; + for (i = 0; i != 256; ++i) // 遍历单个字节的各种情况 + { + for (a = 0; a < 4; ++a) // ba格式 + { + oa = 0; + ooa = 0; + oa += ((i >> 4 & 3) == a) + ((i & 3) == a); + ooa += ((i >> 4 & 3) > a) + ((i & 3) > a); + for (b = 0; b < 4; ++b) + { + oob = ob = 0; + oob += ((i >> 6 & 3) > b && (i >> 4 & 3) == a) + ((i >> 2 & 3) > b && (i & 3) == a); + ob += ((i >> 6 & 3) == b && (i >> 4 & 3) == a) + ((i >> 2 & 3) == b && (i & 3) == a); + ti = a << 2 | b; + fmt->cnt_occ[ti][i] = ob << 24 | oob << 16 | oa << 8 | ooa; + } + } + } +} - printf("kmer size: %ld M\n", kmerSize * sizeof(kmer_range_t) / 1024 / 1024); +// fmt-index的count table，4对应着bwt碱基的累积量，0,1,2,3分别对应着bwt是A,C,G,T，pre-bwt的累积量 +void fmt_gen_cnt_table(uint32_t cnt_table[4][256]) +{ + int i, j, k; + uint32_t x = 0; + for (i = 0; i != 256; ++i) // 遍历单个字节的各种情况 + { + for (k = 0; k < 4; ++k) // bwt碱基 + { + x = 0; // for [A,C,G,T][A,C,G,T] + for (j = 0; j != 4; ++j) // pre-bwt碱基 + x |= (((i >> 6 & 3) == j && (i >> 4 & 3) == k) + ((i >> 2 & 3) == j && (i & 3) == k)) << (j << 3); + cnt_table[k][i] = x; + } + } +} - // exit(0); +// 将fmt结构数据写入到二进制文件 +void dump_fmt(const char *fn, const FMTIndex *fmt) +{ + FILE *fp; + fp = xopen(fn, "wb"); + err_fwrite(&fmt->primary, sizeof(bwtint_t), 1, fp); + err_fwrite(&fmt->sec_primary, sizeof(bwtint_t), 1, fp); + err_fwrite(&fmt->sec_bcp, sizeof(uint8_t), 1, fp); + err_fwrite(&fmt->first_base, sizeof(uint8_t), 1, fp); + err_fwrite(&fmt->last_base, sizeof(uint8_t), 1, fp); + err_fwrite(fmt->L2 + 1, sizeof(bwtint_t), 4, fp); + err_fwrite(fmt->bwt, 4, fmt->bwt_size, fp); + err_fflush(fp); + err_fclose(fp); +} + +// 从文件中读取fmt结构数据 +FMTIndex *fmt_restore_fmt(const char *fn) +{ + FMTIndex *fmt; + fmt = (FMTIndex *)calloc(1, sizeof(FMTIndex)); + FILE *fp = xopen(fn, "rb"); + + fseek(fp, 0, SEEK_END); + fmt->bwt_size = (ftell(fp) - sizeof(bwtint_t) * 6 - 3) >> 2; // 以32位word为单位计算的size + fmt->bwt = (uint32_t *)calloc(fmt->bwt_size, 4); + fseek(fp, 0, SEEK_SET); + err_fread_noeof(&fmt->primary, sizeof(bwtint_t), 1, fp); + err_fread_noeof(&fmt->sec_primary, sizeof(bwtint_t), 1, fp); + err_fread_noeof(&fmt->sec_bcp, sizeof(uint8_t), 1, fp); + err_fread_noeof(&fmt->first_base, sizeof(uint8_t), 1, fp); + err_fread_noeof(&fmt->last_base, sizeof(uint8_t), 1, fp); + err_fread_noeof(fmt->L2 + 1, sizeof(bwtint_t), 4, fp); + fread_fix(fp, fmt->bwt_size << 2, fmt->bwt); + fmt->seq_len = fmt->L2[4]; + err_fclose(fp); + fmt_gen_cnt_occ(fmt); // 字节所能表示的各种碱基组合中，各个碱基的累积数量 + return fmt; +} + +// 将kmer hash数据写入到文件 +void dump_kmer_idx(const char *fn, const KmerEntry *kmer_entry) +{ + FILE *fp; + fp = xopen(fn, "wb"); + err_fwrite(kmer_entry, 1, KMER_ARR_SIZE * sizeof(KmerEntry), fp); + err_fflush(fp); + err_fclose(fp); +} + +// 从文件中读取kmer hash信息 +KmerEntry *fmt_restore_kmer_idx(const char *fn) +{ + FILE *fp = xopen(fn, "rb"); + uint32_t kmer_bytes = KMER_ARR_SIZE * sizeof(KmerEntry); + KmerEntry *kmer_entry = (KmerEntry *)malloc(kmer_bytes); + fread_fix(fp, kmer_bytes, kmer_entry); + err_fclose(fp); + return kmer_entry; +} + +// 读取sa数据 +void fmt_restore_sa(const char *fn, FMTIndex *fmt) +{ + char skipped[256]; + FILE *fp; + bwtint_t primary; + + fp = xopen(fn, "rb"); + err_fread_noeof(&primary, sizeof(bwtint_t), 1, fp); + xassert(primary == fmt->primary, "SA-BWT inconsistency: primary is not the same."); + err_fread_noeof(skipped, sizeof(bwtint_t), 4, fp); // skip + err_fread_noeof(&fmt->sa_intv, sizeof(bwtint_t), 1, fp); + err_fread_noeof(&primary, sizeof(bwtint_t), 1, fp); + xassert(primary == fmt->seq_len, "SA-BWT inconsistency: seq_len is not the same."); + + fmt->n_sa = (fmt->seq_len + fmt->sa_intv) / fmt->sa_intv; + fmt->sa = (uint8_t *)malloc(SA_BYTES(fmt->n_sa)); + + fread_fix(fp, SA_BYTES(fmt->n_sa), fmt->sa); + err_fclose(fp); +} + +// 根据interval-bwt创建fmt-index +FMTIndex *create_fmt_from_bwt(bwt_t *bwt) +{ + // FILE *fmt_out = fopen("fmt.txt", "w"); + FMTIndex *fmt = (FMTIndex *)calloc(1, sizeof(FMTIndex)); + fmt_gen_cnt_occ(fmt); + + bwtint_t i, j, k, m, n, n_occ, cnt[4], cnt2[4]; + uint32_t c[4], c2[16]; /*c用来保存原来的bwt碱基串的累积值，c2用来保存pre-bwt和bwt碱基对的累计值，如AA..TT*/ + uint32_t *buf; /* 计算之后变成fmt结构中bwt部分 */ + +#ifdef FMT_MID_INTERVAL + // 加入中间的check point + uint32_t mc[4] = {0}; + uint32_t cnt_table[4][256]; // 4对应原来的cnt_table，0,1,2,3,分别对应该碱基的扩展 + fmt_gen_cnt_table(cnt_table); +#endif + + fmt->seq_len = bwt->seq_len; // bwt碱基序列的长度，不包含$字符，也就是该长度比bwt matrix长度少1 + for (i = 0; i < 5; ++i) + fmt->L2[i] = bwt->L2[i]; // 每个碱基的总累积值 + fmt->primary = bwt->primary; // $在末尾的行，在bwt matrix行中的排序位置 + + n_occ = (bwt->seq_len + FMT_OCC_INTERVAL - 1) / FMT_OCC_INTERVAL + 1; // check point 个数 + fmt->bwt_size = (fmt->seq_len * 2 + 15) >> 4; // 要保存最后两列碱基 + fmt->bwt_size += n_occ * 20; // A,C,G,T和AA,AC.....TG,TT共20个 + +#ifdef FMT_MID_INTERVAL + uint32_t s1; + bwtint_t mn_occ = (bwt->seq_len >> FMT_OCC_INTV_SHIFT) * (FMT_OCC_INTERVAL / FMT_MID_INTERVAL - 1); + bwtint_t last_seq_len = bwt->seq_len % FMT_OCC_INTERVAL; + mn_occ += (last_seq_len + FMT_MID_INTERVAL - 1) / FMT_MID_INTERVAL - 1; + fmt->bwt_size += mn_occ * 4; + i = 0; +#endif + + buf = (uint32_t *)calloc(fmt->bwt_size, 4); // 开辟计算fmt用到的缓存 + + c[0] = c[1] = c[2] = c[3] = 0; + // 首行的c2，应该是对应的ACGT对应的行，减去1的occ + for (i = 0; i < 4; ++i) + { + bwtint_t before_first_line = fmt->L2[i]; + bwt_occ4(bwt, before_first_line, cnt); + for (j = i * 4, k = 0; k < 4; ++j, ++k) + c2[j] = cnt[k]; + } + // k表示buf存储的偏移量 + for (i = k = 0; i < bwt->seq_len; ++i) + { + // 记录occ + if (i % FMT_OCC_INTERVAL == 0) + { + memcpy(buf + k, c, sizeof(uint32_t) * 4); // bwt str中各个碱基的occ + k += 4; + memcpy(buf + k, c2, sizeof(uint32_t) * 16); // pre-bwt:bwt碱基对的occ + k += 16; +#ifdef FMT_MID_INTERVAL + mc[0] = mc[1] = mc[2] = mc[3] = 0; +#endif + } + // 每个32位整数保存8个倒数第二列碱基（pre-bwt）和8个倒数第一列(bwt)碱基 + if (i % 16 == 0) // 每个32位整数可以包含16个碱基，每次需要处理16个碱基，也就是间隔最小可以设置为16 + { + uint32_t pre_bwt_16_seq = 0; // 16个pre-bwt碱基串 + uint32_t *bwt_addr = bwt_occ_intv(bwt, i) + 4; // 这里加4还是加8要看保存occ的是是uint32还是uint64，bwt字符串i对应的基准行，因为原始的bwt-cp（check point）包含由4个uint32_t(8个uint32_t)组成的occ信息 + int offset = (i % OCC_INTERVAL) / 16; // 每OCC_INTERVAL个碱基共享同一个基准地址，每16个碱基共用一个uint32整型，因此需要偏移量来获取当前碱基串的首地址 + uint32_t bwt_16_seq = *(bwt_addr + offset); // 待处理的当前16个碱基串的首地址 + for (j = 0; j < 16; ++j) // 对于bwt碱基串，一个一个碱基分别处理 + { + bwtint_t cur_str_line = i + j; // 当前碱基在bwt str中的行排序 + if (cur_str_line < bwt->seq_len) // 当前碱基行不应超出bwt str总碱基长度（bwt str长度比bwt matrix长度少1，因为bwt str不包含$） + { + uint8_t bwt_base = bwt_B0(bwt, cur_str_line); // 对应行的bwt的碱基 + // 先求出该碱基对应在第一列的行（对应的bwt matrix行） + bwtint_t cur_mtx_line = cur_str_line; + if (cur_str_line >= bwt->primary) // 因为bwt序列里除去了$符号，所以，超过$所在行之后，对应的seq位置应该加一，才是真正对应bwt matrix的行 + cur_mtx_line += 1; + bwt_occ4(bwt, cur_mtx_line, cnt); // 获取原来bwt-checkpoint中的occ值 + for (m = 0; m < 4; ++m) + c[m] = (uint32_t)cnt[m]; // 碱基m在cur_bwt_mtx_line(包含)之前的累积值，直接拷贝原bwt中的occ即可 + + cnt[bwt_base] -= 1; // 得到cur_bwt_mtx_line(不包含)之前的累积量，即bwt_occ4(bwt, cur_bwt_mtx_line-1, cnt) + bwtint_t bwt_base_mtx_line = bwt->L2[bwt_base] + 1 + cnt[bwt_base]; // bwt_base对应的bwt matrix行（LF变换） + + bwt_occ4(bwt, bwt_base_mtx_line, cnt2); // 计算bwt_base_mtx_line之前的occ + for (n = 0; n < 4; ++n) + { + int c2_idx = bwt_base << 2 | n; // bwt base放在前边 + c2[c2_idx] = (uint32_t)cnt2[n]; // pre-bwt:bwt 碱基对的累计值 + } + bwtint_t bwt_base_str_line = bwt_base_mtx_line; // bwt-str中对应的行排序 + if (bwt_base_str_line >= bwt->primary) // base_line表示在bwt str中的位置，所以超出$为最尾所在行之后，要减掉1 + bwt_base_str_line -= 1; // bwt碱基序列行（不包含$） + uint32_t pre_bwt_base = bwt_B0(bwt, bwt_base_str_line); // bwt列碱基对应的前一个碱基pre-bwt + // 此时，bwt_base对应的bwt matrix首行，是$排在最尾的行，说明bwt_base就是序列的第一个碱基， + // 此时计算出来的pre_bwt_base就是primary前一行的bwt base，以此来代替$字符，在后续的计算过程中需要考虑 + if (bwt_base_mtx_line == bwt->primary) + { + // 计算sec_bcp + fmt->sec_bcp = pre_bwt_base << 2 | bwt_base; // 因为把$当成A处理了 + fmt->sec_primary = cur_mtx_line; // pre-bwt base为$的行排序（bwt-matrix行） + fmt->first_base = bwt_base; // 原始序列第一个碱基 + fmt->last_base = pre_bwt_base; // 计算后替代$字符的碱基（应该是primary行上边一行对应的bwt base） + } + // 暂存 pre-bwt碱基序列 + pre_bwt_16_seq = pre_bwt_16_seq | (pre_bwt_base << (15 - j) * 2); // 序列靠前的碱基排在uint32_t数据中的高位 + } + else + break; + } + // 保存bwt和pre_bwt + uint32_t pre_and_bwt_seq = 0; + uint32_t pre_and_bwt_seq_2 = 0; + for (m = 0; m < 8; ++m) + { + int lshift_bit = 30 - 2 * m; + pre_and_bwt_seq |= (((pre_bwt_16_seq & (3 << lshift_bit)) >> (m * 2)) | ((bwt_16_seq & (3 << lshift_bit)) >> ((m * 2) + 2))); + } + buf[k++] = pre_and_bwt_seq; + + if (j > 8) + { + for (m = 8; m > 0; --m) + { + int lshift_bit = 2 * m - 2; + pre_and_bwt_seq_2 |= (((pre_bwt_16_seq & (3 << lshift_bit)) << (m * 2)) | ((bwt_16_seq & (3 << lshift_bit)) << (m * 2 - 2))); + } + +#ifdef FMT_MID_INTERVAL // 计算前边8+8个碱基的mid interval occ + s1 = pre_and_bwt_seq; + for (m = 0; m < 4; ++m) + mc[m] += cnt_table[m][s1 & 0xff] + cnt_table[m][s1 >> 8 & 0xff] + cnt_table[m][s1 >> 16 & 0xff] + cnt_table[m][s1 >> 24]; +#endif + +#if FMT_MID_INTERVAL == 8 // 如果mid interval是8的话，这里要保存一次 + for (m = 0; m < 4; ++m) + buf[k++] = mc[m]; +#endif + + buf[k++] = pre_and_bwt_seq_2; + +#ifdef FMT_MID_INTERVAL + s1 = pre_and_bwt_seq_2; + for (m = 0; m < 4; ++m) + mc[m] += cnt_table[m][s1 & 0xff] + cnt_table[m][s1 >> 8 & 0xff] + cnt_table[m][s1 >> 16 & 0xff] + cnt_table[m][s1 >> 24]; + + if ((i + 16) % FMT_OCC_INTERVAL != 0 && j == 16 && ((i + 16) & FMT_MID_INTV_MASK) == 0) + for (m = 0; m < 4; ++m) + buf[k++] = mc[m]; +#endif + } + } + } + // the last element + memcpy(buf + k, c, sizeof(uint32_t) * 4); + k += 4; + memcpy(buf + k, c2, sizeof(uint32_t) * 16); + k += 16; + xassert(k == fmt->bwt_size, "inconsistent bwt_size"); + // update fmt + fmt->bwt = buf; + return fmt; +} + +// 扩展两个个碱基，计算bwt base为b的pre-bwt str中各个碱基的occ +void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4]) +{ + uint32_t x = 0; + uint32_t *p, *q, tmp; + bwtint_t str_line = k, cp_line = k & (~FMT_OCC_INTV_MASK); + int i, ti = b1 << 2 | b2; + cnt[0] = 0; + cnt[1] = 0; + cnt[2] = 0; + if (k == (bwtint_t)(-1)) + { + p = fmt->bwt + 4 + b1 * 4; + for (i = b2 + 1; i < 4; ++i) + cnt[2] += p[i]; + cnt[3] = p[b2]; + return; + } + k -= (k >= fmt->primary); // k由bwt矩阵对应的行转换成bwt字符串对应的行（去掉了$，所以大于$的行，都减掉1） + p = fmt_occ_intv(fmt, k); + // fprintf(stderr, "k: %ld\n", k); + for (i = b1 + 1; i < 4; ++i) + cnt[0] += p[i]; // 大于b1的碱基的occ之和 + cnt[1] = p[b1]; // b1的occ + q = p + 4 + b1 * 4; + for (i = b2 + 1; i < 4; ++i) + cnt[2] += q[i]; // 大于b2的occ之和 + cnt[3] = q[b2]; // b2的occ + p += 20; + +#ifdef FMT_MID_INTERVAL // 加入了middle checkpoint + // 使用mid interval信息 + int mk = k % FMT_OCC_INTERVAL; + int n_mintv = mk >> FMT_MID_INTV_SHIFT; + if (n_mintv > 0) // 至少超过了第一个mid interval + { + p += n_mintv * (4 + (FMT_MID_INTERVAL >> 3)) - 4; // 对应的mid interval check point的首地址，即A C G T的局部累积量 + q = p + b1; + for (i = b1 + 1; i < 4; ++i) + x += p[i]; // 大于b1的碱基的occ之和 + cnt[0] += __fmt_mid_sum(x); + x = *q; + cnt[1] += __fmt_mid_sum(x); // b1的occ + for (i = 3; i > b2; --i) + cnt[2] += x >> (i << 3) & 0xff; // 大于b2的occ之和 + cnt[3] += x >> (b2 << 3) & 0xff; // b2的occ + x = 0; + p += 4; + } +#if FMT_MID_INTERVAL == 16 // middle check point interval等于16时候，只需要判断一下是不是要计算两个uint32表示的碱基序列 + if ((mk & FMT_MID_INTV_MASK) >> 3) + { + x += __fmt_occ_e2_aux2(fmt, ti, *p); + ++p; + } +#elif FMT_MID_INTERVAL > 16 // 该地址是bwt和pre_bwt字符串数据的首地址 + uint32_t *end = p + ((k >> 3) - ((k & ~FMT_MID_INTV_MASK) >> 3)); + for (; p < end; ++p) + { + x += __fmt_occ_e2_aux2(fmt, ti, *p); + } +#endif +#else // 没有加入middle check point interval +#if FMT_OCC_INTERVAL > 16 + uint32_t *end = p + ((k >> 3) - ((k & ~FMT_OCC_INTV_MASK) >> 3)); + // p = end - (end - p) / 4; + for (; p < end; ++p) + { + x += __fmt_occ_e2_aux2(fmt, ti, *p); + } +#else // FMT_OCC_INTERVAL等于16的时候，只需要判断一次就可以 + if ((k & FMT_OCC_INTV_MASK) >> 3) + { + x += __fmt_occ_e2_aux2(fmt, ti, *p); + ++p; + } +#endif +#endif + tmp = *p & ~((1U << ((~k & 7) << 2)) - 1); + x += __fmt_occ_e2_aux2(fmt, ti, tmp); + + if (b1 == 0) + { + x -= (~k & 7) << 8; + if (b2 == 0) + x -= (~k & 7) << 24; + } + // 如果跨过了second_primary,那么可能需要减掉一次累积值 + if (b1 == fmt->first_base && cp_line < fmt->sec_primary && str_line >= fmt->sec_primary) + { + if (b2 < fmt->last_base) + cnt[2] -= 1; + else if (b2 == fmt->last_base) + cnt[3] -= 1; + } + cnt[0] += x & 0xff; + cnt[1] += x >> 8 & 0xff; + cnt[2] += x >> 16 & 0xff; + cnt[3] += x >> 24 & 0xff; +} + +// 扩展两个碱基 +inline void fmt_extend2(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t *ok1, bwtintv_t *ok2, int is_back, int b1, int b2) +{ +#ifdef SHOW_PERF +#if 0 + int64_t tmp_time = realtime_msec(); +#endif +#endif + bwtint_t tk[4], tl[4]; + bwtintv_t intv; + // tk表示在k行之前所有各个碱基累积出现次数，tl表示在l行之前的累积 + fmt_e2_occ(fmt, ik->x[!is_back] - 1, b1, b2, tk); + fmt_e2_occ(fmt, ik->x[!is_back] - 1 + ik->x[2], b1, b2, tl); + // 第一次扩展 + intv.x[!is_back] = fmt->L2[b1] + 1 + tk[1]; + intv.x[is_back] = ik->x[is_back] + (ik->x[!is_back] <= fmt->primary && ik->x[!is_back] + ik->x[2] - 1 >= fmt->primary) + tl[0] - tk[0]; + intv.x[2] = tl[1] - tk[1]; + *ok1 = intv; + // 第二次扩展 + intv.x[is_back] = intv.x[is_back] + (intv.x[!is_back] <= fmt->primary && intv.x[!is_back] + intv.x[2] - 1 >= fmt->primary) + tl[2] - tk[2]; + intv.x[!is_back] = fmt->L2[b2] + 1 + tk[3]; + intv.x[2] = tl[3] - tk[3]; + *ok2 = intv; +#ifdef SHOW_PERF +#if 0 + tmp_time = realtime_msec() - tmp_time; + __sync_fetch_and_add(&time_bwt_extend, tmp_time); +#endif +#endif +} + +// 扩展一个碱基 +inline void fmt_extend1(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t *ok, int is_back, int b1) +{ +#ifdef SHOW_PERF +#if 0 + int64_t tmp_time = realtime_msec(); +#endif +#endif + bwtint_t tk[4], tl[4]; + int b2 = 3; // 如果只扩展一次，那么第二个碱基设置成T，可以减小一些计算量，如计算大于b2的累积数量 + // tk表示在k行之前所有各个碱基累积出现次数，tl表示在l行之前的累积 + fmt_e2_occ(fmt, ik->x[!is_back] - 1, b1, b2, tk); + fmt_e2_occ(fmt, ik->x[!is_back] - 1 + ik->x[2], b1, b2, tl); + // 这里是反向扩展 + ok->x[!is_back] = fmt->L2[b1] + 1 + tk[1]; + ok->x[2] = tl[1] - tk[1]; + // 第一次正向扩展 + ok->x[is_back] = ik->x[is_back] + (ik->x[!is_back] <= fmt->primary && ik->x[!is_back] + ik->x[2] - 1 >= fmt->primary) + tl[0] - tk[0]; +#ifdef SHOW_PERF +#if 0 + tmp_time = realtime_msec() - tmp_time; + __sync_fetch_and_add(&time_bwt_extend, tmp_time); +#endif +#endif +} + +// 获取kmer的fmt匹配信息 +inline void kmer_getval_at(KmerEntry *ke, bwtintv_t *ok, int pos) +{ + bwtint_t x0, x1, x2; + int byte_idx = pos * 14; + uint8_t *arr = ke->intv_arr + byte_idx; + x0 = *arr; + x0 = (x0 << 32) | *((uint32_t *)(arr + 1)); + arr += 5; + x1 = *arr; + x1 = (x1 << 32) | *((uint32_t *)(arr + 1)); + arr += 5; + x2 = *((uint32_t *)arr); + ok->x[0] = x0; + ok->x[1] = x1; + ok->x[2] = x2; +} + +// 设置kmer第pos个碱基对应的fmt匹配信息 +inline void kmer_setval_at(KmerEntry *ke, bwtintv_t ik, int pos) +{ + int byte_idx = pos * 14; + uint8_t *arr = ke->intv_arr + byte_idx; + arr[0] = (uint8_t)(ik.x[0] >> 32); + *((uint32_t *)(arr + 1)) = (uint32_t)ik.x[0]; + arr += 5; + arr[0] = (uint8_t)(ik.x[1] >> 32); + *((uint32_t *)(arr + 1)) = (uint32_t)ik.x[1]; + arr += 5; + *((uint32_t *)arr) = (uint32_t)ik.x[2]; +} + +// 查找并保存kmer中每扩展一个碱基对应的fmt位置信息 +void fmt_search_store_kmer(FMTIndex *fmt, const char *q, int qlen, KmerEntry *ke) +{ + bwtintv_t ik; + int i, c1, c2; + bwtint_t tk[4], tl[4]; + + fmt_set_intv(fmt, nst_nt4_table[(int)q[0]], ik); + kmer_setval_at(ke, ik, 0); + // 每次扩展两个碱基 + for (i = 1; i + 1 < qlen; i += 2) + { + // 默认kmer是由ACGT组成的，不含其他字符 + c1 = 3 - nst_nt4_table[(int)q[i]]; + c2 = 3 - nst_nt4_table[(int)q[i + 1]]; + fmt_e2_occ(fmt, ik.x[1] - 1, c1, c2, tk); + fmt_e2_occ(fmt, ik.x[1] - 1 + ik.x[2], c1, c2, tl); + // 第一次扩展的结果 + ik.x[0] = ik.x[0] + (ik.x[1] <= fmt->primary && ik.x[1] + ik.x[2] - 1 >= fmt->primary) + tl[0] - tk[0]; + ik.x[1] = fmt->L2[c1] + 1 + tk[1]; + ik.x[2] = tl[1] - tk[1]; + kmer_setval_at(ke, ik, i); + + // 第二次扩展的结果 + ik.x[0] = ik.x[0] + (ik.x[1] <= fmt->primary && ik.x[1] + ik.x[2] - 1 >= fmt->primary) + tl[2] - tk[2]; + ik.x[1] = fmt->L2[c2] + 1 + tk[3]; + ik.x[2] = tl[3] - tk[3]; + kmer_setval_at(ke, ik, i + 1); + } + if (i < qlen) + { // 最后一次扩展 + c1 = 3 - nst_nt4_table[(int)q[i]]; + c2 = 3; + fmt_e2_occ(fmt, ik.x[1] - 1, c1, c2, tk); + fmt_e2_occ(fmt, ik.x[1] - 1 + ik.x[2], c1, c2, tl); + // 第一次扩展的结果 + ik.x[0] = ik.x[0] + (ik.x[1] <= fmt->primary && ik.x[1] + ik.x[2] - 1 >= fmt->primary) + tl[0] - tk[0]; + ik.x[1] = fmt->L2[c1] + 1 + tk[1]; + ik.x[2] = tl[1] - tk[1]; + kmer_setval_at(ke, ik, i); + } +} + +// 生成所有KMER_LEN长度的序列，字符串表示 +void gen_all_seq(char **seq_arr, int kmer_len) +{ + uint32_t seq_up_val = (1 << (kmer_len << 1)); + for (uint32_t i = 0; i < seq_up_val; ++i) + { + seq_arr[i] = (char *)malloc(kmer_len); + for (int j = kmer_len - 1; j >= 0; --j) + { + seq_arr[i][kmer_len - 1 - j] = BASE[(i >> (j << 1)) & 3]; + } + } +} + +static void fmt_reverse_intvs(bwtintv_v *p) +{ + if (p->n > 1) + { + int j; + for (j = 0; j < p->n >> 1; ++j) + { + bwtintv_t tmp = p->a[p->n - 1 - j]; + p->a[p->n - 1 - j] = p->a[j]; + p->a[j] = tmp; + } + } +} + +// 找smem（seed） +int fmt_smem(const FMTIndex *fmt, int len, const uint8_t *q, int x, int min_intv, bwtintv_v *mem, bwtintv_v *tmpvec[2]) +{ + int i, j, ret, kmer_end = x + KMER_LEN; + bwtintv_t ik, ok1, ok2; + bwtintv_v a[2], *curr; + uint32_t qbit = 0; + //int only_forward = 0; + //if (x == 0 || q[x-1] > 3) only_forward = 1; // 只用向前扩展 + + mem->n = 0; + if (q[x] > 3) return x + 1; + if (min_intv < 1) min_intv = 1; // the interval size should be at least 1 + kv_init(a[0]); kv_init(a[1]); + curr = tmpvec && tmpvec[1] ? tmpvec[1] : &a[1]; // use the temporary vector if provided + kmer_end = kmer_end > len ? len : kmer_end; + // 计算kmer hash key + for (i = x; i < kmer_end; ++i) + { + if (q[i] > 3) // 要考虑碱基是N + break; + qbit |= q[i] << ((KMER_LEN - 1 - (i - x)) << 1); + } + KmerEntry *entry = &fmt->kmer_entry[qbit]; + kmer_getval_at(entry, &ik, 0); // 初始碱基位置 + ik.info = x + 1; + +// check change of the interval size and whether the interval size is too small to be extended further +#define CHECK_INTV_CHANGE(iv, ov) \ + if (ov.x[2] != iv.x[2]) { kv_push(bwtintv_t, *curr, iv); if (ov.x[2] < min_intv) break; } +#define PUSH_VAL_AND_SKIP(iv) \ + do { kv_push(bwtintv_t, *curr, iv); goto backward_search; } while(0) + + // 处理kmer对应的匹配信息 + for (j = 1, curr->n = 0; j < i - x; ++j) + { + kmer_getval_at(entry, &ok1, j); + CHECK_INTV_CHANGE(ik, ok1); + ik = ok1; + ik.info = x + j + 1; + } + if (i != kmer_end) // 遇到了N + PUSH_VAL_AND_SKIP(ik); + + // 扩展kmer之后的碱基 + for (; i + 1 < len; i += 2) + { // forward search + if (q[i] < 4 && q[i + 1] < 4) + { + fmt_extend2(fmt, &ik, &ok1, &ok2, 0, 3 - q[i], 3 - q[i + 1]); + CHECK_INTV_CHANGE(ik, ok1); + ok1.info = i + 1; + CHECK_INTV_CHANGE(ok1, ok2); + ik = ok2; + ik.info = i + 2; + } else if (q[i] < 4) // q[i+1] >= 4 + { + fmt_extend1(fmt, &ik, &ok1, 0, 3 - q[i]); + CHECK_INTV_CHANGE(ik, ok1); + ik = ok1; + ik.info = i + 1; + PUSH_VAL_AND_SKIP(ik); + } + else // q[i] >= 4 + { + PUSH_VAL_AND_SKIP(ik); + } + } + if (i == len - 1) // 扩展到了最后一个碱基 + { + if (q[i] < 4) { + fmt_extend1(fmt, &ik, &ok1, 0, 3 - q[i]); + if (ok1.x[2] != ik.x[2]) { + kv_push(bwtintv_t, *curr, ik); + if (ok1.x[2] < min_intv) + goto backward_search; + } + ik = ok1; + ik.info = i + 1; + } + else + PUSH_VAL_AND_SKIP(ik); + ++i; + } + if (i == len) + kv_push(bwtintv_t, *curr, ik); // push the last interval if we reach the end + + +backward_search: + fmt_reverse_intvs(curr); // s.t. smaller intervals (i.e. longer matches) visited first + ret = curr->a[0].info; // this will be the returned value，扩展到的最远的位置 + // swap = curr; + // curr = prev; + // prev = swap; + // 按照种子进行遍历，反向扩展 +#define CHECK_PUT_MEM(ok, pos, intv) \ + if (ok.x[2] < min_intv) { \ + if (mem->n == 0 || (pos) < mem->a[mem->n - 1].info >> 32) { \ + (intv).info |= (uint64_t)(pos) << 32; \ + kv_push(bwtintv_t, *mem, intv); \ + } \ + break; } + + for (j = 0; j < curr->n; ++j) + { + bwtintv_t *p = &curr->a[j]; // 前向扩展的种子 + for (i = x - 1; i > 0; i -= 2) + { + if (q[i] < 4 && q[i - 1] < 4) // 两个都可以扩展 + { + fmt_extend2(fmt, p, &ok1, &ok2, 1, q[i], q[i - 1]); + CHECK_PUT_MEM(ok1, i + 1, *p); + ok1.info = p->info; + CHECK_PUT_MEM(ok2, i, ok1); + ok2.info = p->info; + *p = ok2; + } + else if (q[i] < 4) // 只能扩展一个 + { + fmt_extend1(fmt, p, &ok1, 0, q[i]); + CHECK_PUT_MEM(ok1, i + 1, *p); + } else + { // 不能扩展 + if (mem->n == 0 || (i + 1) < mem->a[mem->n - 1].info >> 32) + { + p->info |= (uint64_t)(i + 1) << 32; + kv_push(bwtintv_t, *mem, *p); + } + goto fmt_smem_end; + } + } + if (i == 0) { // 扩展到了第一个碱基 + if (q[i] < 4) { + fmt_extend1(fmt, p, &ok1, 0, q[i]); + CHECK_PUT_MEM(ok1, i + 1, *p); + } else { + if (mem->n == 0 || (i + 1) < mem->a[mem->n - 1].info >> 32) + { + p->info |= (uint64_t)(i + 1) << 32; + kv_push(bwtintv_t, *mem, *p); + } + goto fmt_smem_end; + } + --i; + } + if (i == -1) { + if (mem->n == 0 || (i + 1) < mem->a[mem->n - 1].info >> 32) + { + p->info |= (uint64_t)(i + 1) << 32; + kv_push(bwtintv_t, *mem, *p); + } + goto fmt_smem_end; + } + } + +fmt_smem_end: + fmt_reverse_intvs(mem); // s.t. sorted by the start coordinate + if (tmpvec == 0 || tmpvec[0] == 0) + free(a[0].a); + return ret; } \ No newline at end of file diff --git a/fmt_idx.h b/fmt_idx.h index edafff7..730b355 100644 --- a/fmt_idx.h +++ b/fmt_idx.h @@ -1,44 +1,116 @@ /* -Description: 通过fmt-idx数据结构对seed过程进行加速（fm-index twice search in one time） +Description: 通过fmt-idx数据结构对seed过程进行加速 (fm-index twice extend in one search step) Copyright : All right reserved by ICT Author : Zhang Zhonghai Date : 2023/12/24 */ -#ifndef BWA_FMT_IDX_H -#define BWA_FMT_IDX_H +#ifndef FMT_INDEX_H_ +#define FMT_INDEX_H_ + #include #include #include "bwt.h" -typedef uint64_t bwtint_t; +#define FMT_OCC_INTV_SHIFT 8 +#define FMT_OCC_INTERVAL (1 << FMT_OCC_INTV_SHIFT) +#define FMT_OCC_INTV_MASK (FMT_OCC_INTERVAL - 1) -/* - ignore the first 12 bases, and give the entrence directly -*/ -typedef struct { - uint8_t range[9]; // 前36位表示起始行，后36位表示结束行（fm-index索引行） -} kmer_range_t; +#define FMT_MID_INTV_SHIFT 5 +#define FMT_MID_INTERVAL (1 << FMT_MID_INTV_SHIFT) +#define FMT_MID_INTV_MASK (FMT_MID_INTERVAL - 1) +// #undef FMT_MID_INTERVAL + +// 获取碱基c（待查找序列的首个碱基）和对应的互补碱基对应的行，以及间隔 +#define fmt_set_intv(fmt, c, ik) ((ik).x[0] = (fmt)->L2[(int)(c)] + 1, (ik).x[2] = (fmt)->L2[(int)(c) + 1] - (fmt)->L2[(int)(c)], (ik).x[1] = (fmt)->L2[3 - (c)] + 1, (ik).info = 0) +// k行（bwt str行（不包含$））对应的check point occ数据起始地址（小于k且是OCC_INTERVAL的整数倍） + +#if FMT_MID_INTERVAL == 8 +#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 144)) +#elif FMT_MID_INTERVAL == 16 +#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 80)) +#elif FMT_MID_INTERVAL == 32 +#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 48)) +#elif FMT_MID_INTERVAL == 64 +#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 32)) +#elif FMT_MID_INTERVAL == 128 +#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 24)) +#else +#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 20)) +#endif + +// 字节val中包含bwt base为b的pre-bwt中T G C A（按顺序保存在32位整数里（每个占8bit））的数量，下边应该是计算扩展的两个碱基的occ和大于碱基的occ +#define __fmt_occ_e2_aux2(fmt, b, val) \ + ((fmt)->cnt_occ[(b)][(val) & 0xff] + (fmt)->cnt_occ[b][(val) >> 8 & 0xff] + (fmt)->cnt_occ[b][(val) >> 16 & 0xff] + (fmt)->cnt_occ[b][(val) >> 24]) + +#define __fmt_mid_sum(x) \ + ((x) >> 24 & 0xff) + ((x) >> 16 & 0xff) + ((x) >> 8 & 0xff) + ((x) & 0xff) + +// sa存储的行间隔 +#define SA_INTV 2 + +#define KMER_LEN 12 +#define KMER_ARR_SIZE ((1 << (KMER_LEN << 1))) + +// 用来保存kmer对应的fmt的位置信息 +typedef struct +{ + // 40+40+32 14个byte，这样好处理 + uint8_t intv_arr[14 * KMER_LEN]; // 保存kmer中每扩展一个碱基对应的bwtintv_t数据 +} KmerEntry; + +// fm-index, extend twice in one search step (one memory access) typedef struct { - bwtint_t primary; // S^{-1}(0), or the primary index of BWT - bwtint_t L2[5]; // C(), cumulative count - bwtint_t seq_len; // sequence length - bwtint_t bwt_size; // size of bwt, about seq_len/4 - uint32_t *bwt; // BWT - // kmer entry - kmer_range_t *kmer_range; + bwtint_t primary; // S^{-1}(0), or the primary index of BWT + bwtint_t sec_primary; // second primary line + bwtint_t L2[5]; // C(), cumulative count + bwtint_t seq_len; // sequence length + bwtint_t bwt_size; // size of bwt, about seq_len/4 + uint32_t *bwt; // BWT // occurance array, separated to two parts - uint32_t cnt_table[256]; + uint32_t cnt_occ[16][256]; // 前16-24位表示b（碱基）的occ，8-16位表示大于b的occ，0-8表示大于a的occ，ba格式 + uint8_t sec_bcp; // base couple for sec primary line, AA=>0, AC=>1 ... TT=>15 + uint8_t first_base; // 序列的第一个碱基2bit的int类型，0,1,2,3 + uint8_t last_base; // dollar转换成的base + // 保存kmer对应的fmt位置信息 + KmerEntry *kmer_entry; // suffix array int sa_intv; bwtint_t n_sa; - bwtint_t *sa; -} bwtd_t; + uint8_t *sa; +} FMTIndex; -// 创建fmt-index索引数据 -void BuildBwtdFromBwt(bwt_t *bwt, bwtd_t **bwtd_p); +// 将fmt结构数据写入到二进制文件 +void dump_fmt(const char *fn, const FMTIndex *fmt); +// 从文件中读取fmt结构数据 +FMTIndex *fmt_restore_fmt(const char *fn); +// 将kmer hash数据写入到文件 +void dump_kmer_idx(const char *fn, const KmerEntry *kmer_entry); +// 从文件中读取kmer hash信息 +KmerEntry *fmt_restore_kmer_idx(const char *fn); +// 读取sa数据 +void fmt_restore_sa(const char *fn, FMTIndex *fmt); +// 根据interval-bwt创建fmt-index +FMTIndex *create_fmt_from_bwt(bwt_t *bwt); +// 扩展两个个碱基，计算bwt base为b的pre-bwt str中各个碱基的occ +void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4]); +// 扩展两个碱基 +void fmt_extend2(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t *ok1, bwtintv_t *ok2, int is_back, int b1, int b2); +// 扩展一个碱基 +void fmt_extend1(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t *ok, int is_back, int b1); +// 查找并保存kmer中每扩展一个碱基对应的fmt位置信息 +void fmt_search_store_kmer(FMTIndex *fmt, const char *q, int qlen, KmerEntry *ke); +// 生成所有KMER_LEN长度的序列，字符串表示 +void gen_all_seq(char **seq_arr, int kmer_len); +// 设置kmer第pos个碱基对应的fmt匹配信息 +void kmer_setval_at(KmerEntry *ke, bwtintv_t ik, int pos); +// 获取kmer的fmt匹配信息 +void kmer_getval_at(KmerEntry *ke, bwtintv_t *ok, int pos); + +// 找smem（seed） +int fmt_smem(const FMTIndex *fmt, int len, const uint8_t *q, int x, int min_intv, bwtintv_v *mem, bwtintv_v *tmpvec[2]); #endif \ No newline at end of file diff --git a/run.sh b/run.sh index 690b563..f61e300 100755 --- a/run.sh +++ b/run.sh @@ -1,3 +1,14 @@ +thread=1 +#n_r1=~/data/fastq/ZY2105177532213000/n_r1.fq +#n_r2=~/data/fastq/ZY2105177532213000/n_r2.fq +#n_r1=~/data/fastq/ZY2105177532213000/sn_r1.fq +#n_r2=~/data/fastq/ZY2105177532213000/sn_r2.fq +#reference=~/data/reference/human_g1k_v37_decoy.fasta +n_r1=~/fastq/sn_r1.fq +n_r2=~/fastq/sn_r2.fq +#n_r1=~/fastq/tiny_n_r1.fq +#n_r2=~/fastq/tiny_n_r2.fq +reference=~/reference/human_g1k_v37_decoy.fasta #time ./bwa mem -t 12 -M -R @RG\\tID:normal\\tSM:normal\\tPL:illumina\\tLB:normal\\tPG:bwa \ # /home/zzh/data/reference/human_g1k_v37_decoy.fasta \ # /home/zzh/data/fastq/nm1.fq \ @@ -9,10 +20,10 @@ # /mnt/d/data/fastq/ZY2105177532213000/ZY2105177532213010_L4_2.fq.gz \ # -o /dev/null -time ./bwa mem -t 12 -M -R @RG\\tID:normal\\tSM:normal\\tPL:illumina\\tLB:normal\\tPG:bwa \ - /home/zzh/data/reference/human_g1k_v37_decoy.fasta \ - /mnt/d/data/fastq/ZY2105177532213000/ZY2105177532213030_L3_1.fq.gz \ - /mnt/d/data/fastq/ZY2105177532213000/ZY2105177532213030_L3_2.fq.gz \ +time ./bwa mem -t $thread -M -R @RG\\tID:normal\\tSM:normal\\tPL:illumina\\tLB:normal\\tPG:bwa \ + $reference \ + $n_r1 \ + $n_r2 \ -o /dev/null diff --git a/utils.h b/utils.h index 1b88b4e..881f912 100644 --- a/utils.h +++ b/utils.h @@ -37,10 +37,15 @@ extern int64_t time_ksw_extend2, time_ksw_global2, time_ksw_align2, time_bwt_smem1a, + time_bwt_extend, time_bwt_occ4, time_bwt_sa, time_bwt_sa_read; +extern int64_t dn; + +extern FILE *fp1; + #endif #ifdef __GNUC__