diff --git a/.gitignore b/.gitignore index 9ae2235..dbfdc4d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ +*.paf +*.sam .cproject .project .*.swp diff --git a/.vscode/launch.json b/.vscode/launch.json index 22af17e..b70e95b 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -5,7 +5,7 @@ "version": "0.2.0", "configurations": [ { - "name": "Launch", + "name": "read overlap", "preLaunchTask": "Build", "type": "cppdbg", "request": "launch", @@ -14,13 +14,31 @@ "-x", "ava-ont", "-t", - "1", + "4", "/public/home/zzh/work/3gseq/TGM-2021YFF/Acinetobacter_pittii.fastq", "/public/home/zzh/work/3gseq/TGM-2021YFF/Acinetobacter_pittii.fastq", "-o", "reads.paf" ], "cwd": "${workspaceFolder}", // 当前工作路径：当前文件所在的工作空间 + }, + { + "name": "mapping", + "preLaunchTask": "Build", + "type": "cppdbg", + "request": "launch", + "program": "${workspaceRoot}/minimap2", + "args": [ + "-ax", + "map-ont", + "-t", + "1", + "/public/home/zzh/work/3gseq/TGM-2021YFF/reads.fasta", + "/public/home/zzh/work/3gseq/TGM-2021YFF/Acinetobacter_pittii.fastq", + "-o", + "aln.sam" + ], + "cwd": "${workspaceFolder}", // 当前工作路径：当前文件所在的工作空间 } ] } \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json index 8909a6f..5549e14 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -1,6 +1,7 @@ { "files.associations": { "minimap.h": "c", - "time.h": "c" + "time.h": "c", + "kalloc.h": "c" } } \ No newline at end of file diff --git a/align.c b/align.c index f13d264..d50a2eb 100644 --- a/align.c +++ b/align.c @@ -5,16 +5,21 @@ #include "minimap.h" #include "mmpriv.h" #include "ksw2.h" +#ifdef SHOW_PERF +extern int64_t get_mseconds(); +extern int64_t time_ksw_extd2_sse; +#endif static void ksw_gen_simple_mat(int m, int8_t *mat, int8_t a, int8_t b, int8_t sc_ambi) { int i, j; - a = a < 0? -a : a; - b = b > 0? -b : b; - sc_ambi = sc_ambi > 0? -sc_ambi : sc_ambi; - for (i = 0; i < m - 1; ++i) { + a = a < 0 ? -a : a; + b = b > 0 ? -b : b; + sc_ambi = sc_ambi > 0 ? -sc_ambi : sc_ambi; + for (i = 0; i < m - 1; ++i) + { for (j = 0; j < m - 1; ++j) - mat[i * m + j] = i == j? a : b; + mat[i * m + j] = i == j ? a : b; mat[i * m + m - 1] = sc_ambi; } for (j = 0; j < m; ++j) @@ -23,36 +28,40 @@ static void ksw_gen_simple_mat(int m, int8_t *mat, int8_t a, int8_t b, int8_t sc static void ksw_gen_ts_mat(int m, int8_t *mat, int8_t a, int8_t b, int8_t transition, int8_t sc_ambi) { - assert(m==5); - ksw_gen_simple_mat(m,mat,a,b,sc_ambi); - transition = transition > 0? -transition : transition; - mat[0*m+2]=transition; // A->G - mat[1*m+3]=transition; // C->T - mat[2*m+0]=transition; // G->A - mat[3*m+1]=transition; // T->C + assert(m == 5); + ksw_gen_simple_mat(m, mat, a, b, sc_ambi); + transition = transition > 0 ? -transition : transition; + mat[0 * m + 2] = transition; // A->G + mat[1 * m + 3] = transition; // C->T + mat[2 * m + 0] = transition; // G->A + mat[3 * m + 1] = transition; // T->C } static inline void mm_seq_rev(uint32_t len, uint8_t *seq) { uint32_t i; uint8_t t; - for (i = 0; i < len>>1; ++i) + for (i = 0; i < len >> 1; ++i) t = seq[i], seq[i] = seq[len - 1 - i], seq[len - 1 - i] = t; } static inline void update_max_zdrop(int32_t score, int i, int j, int32_t *max, int *max_i, int *max_j, int e, int *max_zdrop, int pos[2][2]) { - if (score < *max) { + if (score < *max) + { int li = i - *max_i; int lj = j - *max_j; - int diff = li > lj? li - lj : lj - li; + int diff = li > lj ? li - lj : lj - li; int z = *max - score - diff * e; - if (z > *max_zdrop) { + if (z > *max_zdrop) + { *max_zdrop = z; pos[0][0] = *max_i, pos[0][1] = i; pos[1][0] = *max_j, pos[1][1] = j; } - } else *max = score, *max_i = i, *max_j = j; + } + else + *max = score, *max_i = i, *max_j = j; } static int mm_test_zdrop(void *km, const mm_mapopt_t *opt, const uint8_t *qseq, const uint8_t *tseq, uint32_t n_cigar, uint32_t *cigar, const int8_t *mat) @@ -62,32 +71,41 @@ static int mm_test_zdrop(void *km, const mm_mapopt_t *opt, const uint8_t *qseq, int pos[2][2] = {{-1, -1}, {-1, -1}}, q_len, t_len; // find the score and the region where score drops most along diagonal - for (k = 0, score = 0; k < n_cigar; ++k) { - uint32_t l, op = cigar[k]&0xf, len = cigar[k]>>4; - if (op == MM_CIGAR_MATCH) { - for (l = 0; l < len; ++l) { + for (k = 0, score = 0; k < n_cigar; ++k) + { + uint32_t l, op = cigar[k] & 0xf, len = cigar[k] >> 4; + if (op == MM_CIGAR_MATCH) + { + for (l = 0; l < len; ++l) + { score += mat[tseq[i + l] * 5 + qseq[j + l]]; - update_max_zdrop(score, i+l, j+l, &max, &max_i, &max_j, opt->e, &max_zdrop, pos); + update_max_zdrop(score, i + l, j + l, &max, &max_i, &max_j, opt->e, &max_zdrop, pos); } i += len, j += len; - } else if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL || op == MM_CIGAR_N_SKIP) { + } + else if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL || op == MM_CIGAR_N_SKIP) + { score -= opt->q + opt->e * len; - if (op == MM_CIGAR_INS) j += len; - else i += len; + if (op == MM_CIGAR_INS) + j += len; + else + i += len; update_max_zdrop(score, i, j, &max, &max_i, &max_j, opt->e, &max_zdrop, pos); } } // test if there is an inversion in the most dropped region q_len = pos[1][1] - pos[1][0], t_len = pos[0][1] - pos[0][0]; - if (!(opt->flag&(MM_F_SPLICE|MM_F_SR|MM_F_FOR_ONLY|MM_F_REV_ONLY)) && max_zdrop > opt->zdrop_inv && q_len < opt->max_gap && t_len < opt->max_gap) { + if (!(opt->flag & (MM_F_SPLICE | MM_F_SR | MM_F_FOR_ONLY | MM_F_REV_ONLY)) && max_zdrop > opt->zdrop_inv && q_len < opt->max_gap && t_len < opt->max_gap) + { uint8_t *qseq2; void *qp; int q_off, t_off; - qseq2 = (uint8_t*)kmalloc(km, q_len); - for (i = 0; i < q_len; ++i) { + qseq2 = (uint8_t *)kmalloc(km, q_len); + for (i = 0; i < q_len; ++i) + { int c = qseq[pos[1][1] - i - 1]; - qseq2[i] = c >= 4? 4 : 3 - c; + qseq2[i] = c >= 4 ? 4 : 3 - c; } qp = ksw_ll_qinit(km, 2, q_len, qseq2, 5, mat); score = ksw_ll_i16(qp, t_len, tseq + pos[0][0], opt->q, opt->e, &q_off, &t_off); @@ -96,7 +114,7 @@ static int mm_test_zdrop(void *km, const mm_mapopt_t *opt, const uint8_t *qseq, if (score >= opt->min_chain_score * opt->a && score >= opt->min_dp_max) return 2; // there is a potential inversion } - return max_zdrop > opt->zdrop? 1 : 0; + return max_zdrop > opt->zdrop ? 1 : 0; } static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, int *qshift, int *tshift) @@ -105,47 +123,67 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, int32_t toff = 0, qoff = 0, to_shrink = 0; uint32_t k; *qshift = *tshift = 0; - if (p->n_cigar <= 1) return; - for (k = 0; k < p->n_cigar; ++k) { // indel left alignment - uint32_t op = p->cigar[k]&0xf, len = p->cigar[k]>>4; - if (len == 0) to_shrink = 1; - if (op == MM_CIGAR_MATCH) { + if (p->n_cigar <= 1) + return; + for (k = 0; k < p->n_cigar; ++k) + { // indel left alignment + uint32_t op = p->cigar[k] & 0xf, len = p->cigar[k] >> 4; + if (len == 0) + to_shrink = 1; + if (op == MM_CIGAR_MATCH) + { toff += len, qoff += len; - } else if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL) { - if (k > 0 && k < p->n_cigar - 1 && (p->cigar[k-1]&0xf) == 0 && (p->cigar[k+1]&0xf) == 0) { - int l, prev_len = p->cigar[k-1] >> 4; - if (op == MM_CIGAR_INS) { + } + else if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL) + { + if (k > 0 && k < p->n_cigar - 1 && (p->cigar[k - 1] & 0xf) == 0 && (p->cigar[k + 1] & 0xf) == 0) + { + int l, prev_len = p->cigar[k - 1] >> 4; + if (op == MM_CIGAR_INS) + { for (l = 0; l < prev_len; ++l) if (qseq[qoff - 1 - l] != qseq[qoff + len - 1 - l]) break; - } else { + } + else + { for (l = 0; l < prev_len; ++l) if (tseq[toff - 1 - l] != tseq[toff + len - 1 - l]) break; } if (l > 0) - p->cigar[k-1] -= l<<4, p->cigar[k+1] += l<<4, qoff -= l, toff -= l; - if (l == prev_len) to_shrink = 1; + p->cigar[k - 1] -= l << 4, p->cigar[k + 1] += l << 4, qoff -= l, toff -= l; + if (l == prev_len) + to_shrink = 1; } - if (op == MM_CIGAR_INS) qoff += len; - else toff += len; - } else if (op == MM_CIGAR_N_SKIP) { + if (op == MM_CIGAR_INS) + qoff += len; + else + toff += len; + } + else if (op == MM_CIGAR_N_SKIP) + { toff += len; } } assert(qoff == r->qe - r->qs && toff == r->re - r->rs); - for (k = 0; k < p->n_cigar - 2; ++k) { // fix CIGAR like 5I6D7I - if ((p->cigar[k]&0xf) > 0 && (p->cigar[k]&0xf) + (p->cigar[k+1]&0xf) == 3) { - uint32_t l, s[3] = {0,0,0}; - for (l = k; l < p->n_cigar; ++l) { // count number of adjacent I and D - uint32_t op = p->cigar[l]&0xf; - if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL || p->cigar[l]>>4 == 0) + for (k = 0; k < p->n_cigar - 2; ++k) + { // fix CIGAR like 5I6D7I + if ((p->cigar[k] & 0xf) > 0 && (p->cigar[k] & 0xf) + (p->cigar[k + 1] & 0xf) == 3) + { + uint32_t l, s[3] = {0, 0, 0}; + for (l = k; l < p->n_cigar; ++l) + { // count number of adjacent I and D + uint32_t op = p->cigar[l] & 0xf; + if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL || p->cigar[l] >> 4 == 0) s[op] += p->cigar[l] >> 4; - else break; + else + break; } - if (s[1] > 0 && s[2] > 0 && l - k > 2) { // turn to a single I and a single D - p->cigar[k] = s[1]<<4|MM_CIGAR_INS; - p->cigar[k+1] = s[2]<<4|MM_CIGAR_DEL; + if (s[1] > 0 && s[2] > 0 && l - k > 2) + { // turn to a single I and a single D + p->cigar[k] = s[1] << 4 | MM_CIGAR_INS; + p->cigar[k + 1] = s[2] << 4 | MM_CIGAR_DEL; for (k += 2; k < l; ++k) p->cigar[k] &= 0xf; to_shrink = 1; @@ -153,25 +191,33 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, k = l; } } - if (to_shrink) { // squeeze out zero-length operations + if (to_shrink) + { // squeeze out zero-length operations int32_t l = 0; for (k = 0; k < p->n_cigar; ++k) // squeeze out zero-length operations - if (p->cigar[k]>>4 != 0) + if (p->cigar[k] >> 4 != 0) p->cigar[l++] = p->cigar[k]; p->n_cigar = l; for (k = l = 0; k < p->n_cigar; ++k) // merge two adjacent operations if they are the same - if (k == p->n_cigar - 1 || (p->cigar[k]&0xf) != (p->cigar[k+1]&0xf)) + if (k == p->n_cigar - 1 || (p->cigar[k] & 0xf) != (p->cigar[k + 1] & 0xf)) p->cigar[l++] = p->cigar[k]; - else p->cigar[k+1] += p->cigar[k]>>4<<4; // add length to the next CIGAR operator + else + p->cigar[k + 1] += p->cigar[k] >> 4 << 4; // add length to the next CIGAR operator p->n_cigar = l; } - if ((p->cigar[0]&0xf) == MM_CIGAR_INS || (p->cigar[0]&0xf) == MM_CIGAR_DEL) { // get rid of leading I or D + if ((p->cigar[0] & 0xf) == MM_CIGAR_INS || (p->cigar[0] & 0xf) == MM_CIGAR_DEL) + { // get rid of leading I or D int32_t l = p->cigar[0] >> 4; - if ((p->cigar[0]&0xf) == MM_CIGAR_INS) { - if (r->rev) r->qe -= l; - else r->qs += l; + if ((p->cigar[0] & 0xf) == MM_CIGAR_INS) + { + if (r->rev) + r->qe -= l; + else + r->qs += l; *qshift = l; - } else r->rs += l, *tshift = l; + } + else + r->rs += l, *tshift = l; --p->n_cigar; memmove(p->cigar, p->cigar + 1, p->n_cigar * 4); } @@ -182,63 +228,107 @@ static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t uint32_t n_EQX = 0; uint32_t k, l, m, cap, toff = 0, qoff = 0, n_M = 0; mm_extra_t *p; - if (r->p == 0) return; - for (k = 0; k < r->p->n_cigar; ++k) { - uint32_t op = r->p->cigar[k]&0xf, len = r->p->cigar[k]>>4; - if (op == MM_CIGAR_MATCH) { - while (len > 0) { - for (l = 0; l < len && qseq[qoff + l] == tseq[toff + l]; ++l) {} // run of "="; TODO: N<=>N is converted to "=" - if (l > 0) { ++n_EQX; len -= l; toff += l; qoff += l; } + if (r->p == 0) + return; + for (k = 0; k < r->p->n_cigar; ++k) + { + uint32_t op = r->p->cigar[k] & 0xf, len = r->p->cigar[k] >> 4; + if (op == MM_CIGAR_MATCH) + { + while (len > 0) + { + for (l = 0; l < len && qseq[qoff + l] == tseq[toff + l]; ++l) + { + } // run of "="; TODO: N<=>N is converted to "=" + if (l > 0) + { + ++n_EQX; + len -= l; + toff += l; + qoff += l; + } - for (l = 0; l < len && qseq[qoff + l] != tseq[toff + l]; ++l) {} // run of "X" - if (l > 0) { ++n_EQX; len -= l; toff += l; qoff += l; } + for (l = 0; l < len && qseq[qoff + l] != tseq[toff + l]; ++l) + { + } // run of "X" + if (l > 0) + { + ++n_EQX; + len -= l; + toff += l; + qoff += l; + } } ++n_M; - } else if (op == MM_CIGAR_INS) { + } + else if (op == MM_CIGAR_INS) + { qoff += len; - } else if (op == MM_CIGAR_DEL) { + } + else if (op == MM_CIGAR_DEL) + { toff += len; - } else if (op == MM_CIGAR_N_SKIP) { + } + else if (op == MM_CIGAR_N_SKIP) + { toff += len; } } // update in-place if we can - if (n_EQX == n_M) { - for (k = 0; k < r->p->n_cigar; ++k) { - uint32_t op = r->p->cigar[k]&0xf, len = r->p->cigar[k]>>4; - if (op == MM_CIGAR_MATCH) r->p->cigar[k] = len << 4 | MM_CIGAR_EQ_MATCH; + if (n_EQX == n_M) + { + for (k = 0; k < r->p->n_cigar; ++k) + { + uint32_t op = r->p->cigar[k] & 0xf, len = r->p->cigar[k] >> 4; + if (op == MM_CIGAR_MATCH) + r->p->cigar[k] = len << 4 | MM_CIGAR_EQ_MATCH; } return; } // allocate new storage cap = r->p->n_cigar + (n_EQX - n_M) + sizeof(mm_extra_t); kroundup32(cap); - p = (mm_extra_t*)calloc(cap, 4); + p = (mm_extra_t *)calloc(cap, 4); memcpy(p, r->p, sizeof(mm_extra_t)); p->capacity = cap; // update cigar while copying toff = qoff = m = 0; - for (k = 0; k < r->p->n_cigar; ++k) { - uint32_t op = r->p->cigar[k]&0xf, len = r->p->cigar[k]>>4; - if (op == MM_CIGAR_MATCH) { - while (len > 0) { + for (k = 0; k < r->p->n_cigar; ++k) + { + uint32_t op = r->p->cigar[k] & 0xf, len = r->p->cigar[k] >> 4; + if (op == MM_CIGAR_MATCH) + { + while (len > 0) + { // match - for (l = 0; l < len && qseq[qoff + l] == tseq[toff + l]; ++l) {} - if (l > 0) p->cigar[m++] = l << 4 | MM_CIGAR_EQ_MATCH; + for (l = 0; l < len && qseq[qoff + l] == tseq[toff + l]; ++l) + { + } + if (l > 0) + p->cigar[m++] = l << 4 | MM_CIGAR_EQ_MATCH; len -= l; toff += l, qoff += l; // mismatch - for (l = 0; l < len && qseq[qoff + l] != tseq[toff + l]; ++l) {} - if (l > 0) p->cigar[m++] = l << 4 | MM_CIGAR_X_MISMATCH; + for (l = 0; l < len && qseq[qoff + l] != tseq[toff + l]; ++l) + { + } + if (l > 0) + p->cigar[m++] = l << 4 | MM_CIGAR_X_MISMATCH; len -= l; toff += l, qoff += l; } continue; - } else if (op == MM_CIGAR_INS) { + } + else if (op == MM_CIGAR_INS) + { qoff += len; - } else if (op == MM_CIGAR_DEL) { + } + else if (op == MM_CIGAR_DEL) + { toff += len; - } else if (op == MM_CIGAR_N_SKIP) { + } + else if (op == MM_CIGAR_N_SKIP) + { toff += len; } p->cigar[m++] = r->p->cigar[k]; @@ -254,71 +344,102 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts int32_t qshift, tshift, toff = 0, qoff = 0; double s = 0.0, max = 0.0; mm_extra_t *p = r->p; - if (p == 0) return; + if (p == 0) + return; mm_fix_cigar(r, qseq, tseq, &qshift, &tshift); qseq += qshift, tseq += tshift; // qseq and tseq may be shifted due to the removal of leading I/D r->blen = r->mlen = 0; - for (k = 0; k < p->n_cigar; ++k) { - uint32_t op = p->cigar[k]&0xf, len = p->cigar[k]>>4; - if (op == MM_CIGAR_MATCH) { + for (k = 0; k < p->n_cigar; ++k) + { + uint32_t op = p->cigar[k] & 0xf, len = p->cigar[k] >> 4; + if (op == MM_CIGAR_MATCH) + { int n_ambi = 0, n_diff = 0; - for (l = 0; l < len; ++l) { + for (l = 0; l < len; ++l) + { int cq = qseq[qoff + l], ct = tseq[toff + l]; - if (ct > 3 || cq > 3) ++n_ambi; - else if (ct != cq) ++n_diff; + if (ct > 3 || cq > 3) + ++n_ambi; + else if (ct != cq) + ++n_diff; s += mat[ct * 5 + cq]; - if (s < 0) s = 0; - else max = max > s? max : s; + if (s < 0) + s = 0; + else + max = max > s ? max : s; } r->blen += len - n_ambi, r->mlen += len - (n_ambi + n_diff), p->n_ambi += n_ambi; toff += len, qoff += len; - } else if (op == MM_CIGAR_INS) { + } + else if (op == MM_CIGAR_INS) + { int n_ambi = 0; for (l = 0; l < len; ++l) - if (qseq[qoff + l] > 3) ++n_ambi; + if (qseq[qoff + l] > 3) + ++n_ambi; r->blen += len - n_ambi, p->n_ambi += n_ambi; - if (log_gap) s -= q + (double)e * mg_log2(1.0 + len); - else s -= q + e; - if (s < 0) s = 0; + if (log_gap) + s -= q + (double)e * mg_log2(1.0 + len); + else + s -= q + e; + if (s < 0) + s = 0; qoff += len; - } else if (op == MM_CIGAR_DEL) { + } + else if (op == MM_CIGAR_DEL) + { int n_ambi = 0; for (l = 0; l < len; ++l) - if (tseq[toff + l] > 3) ++n_ambi; + if (tseq[toff + l] > 3) + ++n_ambi; r->blen += len - n_ambi, p->n_ambi += n_ambi; - if (log_gap) s -= q + (double)e * mg_log2(1.0 + len); - else s -= q + e; - if (s < 0) s = 0; + if (log_gap) + s -= q + (double)e * mg_log2(1.0 + len); + else + s -= q + e; + if (s < 0) + s = 0; toff += len; - } else if (op == MM_CIGAR_N_SKIP) { + } + else if (op == MM_CIGAR_N_SKIP) + { toff += len; } } p->dp_max = (int32_t)(max + .499); assert(qoff == r->qe - r->qs && toff == r->re - r->rs); - if (is_eqx) mm_update_cigar_eqx(r, qseq, tseq); // NB: it has to be called here as changes to qseq and tseq are not returned + if (is_eqx) + mm_update_cigar_eqx(r, qseq, tseq); // NB: it has to be called here as changes to qseq and tseq are not returned } static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) // TODO: this calls the libc realloc() { mm_extra_t *p; - if (n_cigar == 0) return; - if (r->p == 0) { - uint32_t capacity = n_cigar + sizeof(mm_extra_t)/4; + if (n_cigar == 0) + return; + if (r->p == 0) + { + uint32_t capacity = n_cigar + sizeof(mm_extra_t) / 4; kroundup32(capacity); - r->p = (mm_extra_t*)calloc(capacity, 4); + r->p = (mm_extra_t *)calloc(capacity, 4); r->p->capacity = capacity; - } else if (r->p->n_cigar + n_cigar + sizeof(mm_extra_t)/4 > r->p->capacity) { - r->p->capacity = r->p->n_cigar + n_cigar + sizeof(mm_extra_t)/4; + } + else if (r->p->n_cigar + n_cigar + sizeof(mm_extra_t) / 4 > r->p->capacity) + { + r->p->capacity = r->p->n_cigar + n_cigar + sizeof(mm_extra_t) / 4; kroundup32(r->p->capacity); - r->p = (mm_extra_t*)realloc(r->p, r->p->capacity * 4); + r->p = (mm_extra_t *)realloc(r->p, r->p->capacity * 4); } p = r->p; - if (p->n_cigar > 0 && (p->cigar[p->n_cigar-1]&0xf) == (cigar[0]&0xf)) { // same CIGAR op at the boundary - p->cigar[p->n_cigar-1] += cigar[0]>>4<<4; - if (n_cigar > 1) memcpy(p->cigar + p->n_cigar, cigar + 1, (n_cigar - 1) * 4); + if (p->n_cigar > 0 && (p->cigar[p->n_cigar - 1] & 0xf) == (cigar[0] & 0xf)) + { // same CIGAR op at the boundary + p->cigar[p->n_cigar - 1] += cigar[0] >> 4 << 4; + if (n_cigar > 1) + memcpy(p->cigar + p->n_cigar, cigar + 1, (n_cigar - 1) * 4); p->n_cigar += n_cigar - 1; - } else { + } + else + { memcpy(p->cigar + p->n_cigar, cigar, n_cigar * 4); p->n_cigar += n_cigar; } @@ -326,31 +447,53 @@ static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) // static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const uint8_t *junc, const int8_t *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez) { - if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) { +#ifdef SHOW_PERF + int64_t tmp_cur_time = 0, tmp_diff = 0; +#endif + if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) + { int i; fprintf(stderr, "===> q=(%d,%d), e=(%d,%d), bw=%d, flag=%d, zdrop=%d <===\n", opt->q, opt->q2, opt->e, opt->e2, w, flag, opt->zdrop); - for (i = 0; i < tlen; ++i) fputc("ACGTN"[tseq[i]], stderr); + for (i = 0; i < tlen; ++i) + fputc("ACGTN"[tseq[i]], stderr); fputc('\n', stderr); - for (i = 0; i < qlen; ++i) fputc("ACGTN"[qseq[i]], stderr); + for (i = 0; i < qlen; ++i) + fputc("ACGTN"[qseq[i]], stderr); fputc('\n', stderr); } - if (opt->b != opt->transition) flag |= KSW_EZ_GENERIC_SC; - if (opt->max_sw_mat > 0 && (int64_t)tlen * qlen > opt->max_sw_mat) { + if (opt->b != opt->transition) + flag |= KSW_EZ_GENERIC_SC; + if (opt->max_sw_mat > 0 && (int64_t)tlen * qlen > opt->max_sw_mat) + { ksw_reset_extz(ez); ez->zdropped = 1; - } else if (opt->flag & MM_F_SPLICE) { + } + else if (opt->flag & MM_F_SPLICE) + { int flag_tmp = flag; - if (!(opt->flag & MM_F_SPLICE_OLD)) flag_tmp |= KSW_EZ_SPLICE_CMPLX; + if (!(opt->flag & MM_F_SPLICE_OLD)) + flag_tmp |= KSW_EZ_SPLICE_CMPLX; ksw_exts2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->noncan, zdrop, opt->junc_bonus, flag_tmp, junc, ez); - } else if (opt->q == opt->q2 && opt->e == opt->e2) + } + else if (opt->q == opt->q2 && opt->e == opt->e2) ksw_extz2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, w, zdrop, end_bonus, flag, ez); else + { +#ifdef SHOW_PERF + tmp_cur_time = get_mseconds(); +#endif ksw_extd2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez); - if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) { +#ifdef SHOW_PERF + tmp_diff = get_mseconds() - tmp_cur_time; + __sync_fetch_and_add(&time_ksw_extd2_sse, tmp_diff); +#endif + } + if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) + { int i; fprintf(stderr, "score=%d, cigar=", ez->score); for (i = 0; i < ez->n_cigar; ++i) - fprintf(stderr, "%d%c", ez->cigar[i]>>4, MM_CIGAR_STR[ez->cigar[i]&0xf]); + fprintf(stderr, "%d%c", ez->cigar[i] >> 4, MM_CIGAR_STR[ez->cigar[i] & 0xf]); fprintf(stderr, "\n"); } } @@ -360,24 +503,29 @@ static inline int mm_get_hplen_back(const mm_idx_t *mi, uint32_t rid, uint32_t x int64_t i, off0 = mi->seq[rid].offset, off = off0 + x; int c = mm_seq4_get(mi->S, off); for (i = off - 1; i >= off0; --i) - if (mm_seq4_get(mi->S, i) != c) break; + if (mm_seq4_get(mi->S, i) != c) + break; return (int)(off - i); } static inline void mm_adjust_minier(const mm_idx_t *mi, uint8_t *const qseq0[2], mm128_t *a, int32_t *r, int32_t *q) { - if (mi->flag & MM_I_HPC) { - const uint8_t *qseq = qseq0[a->x>>63]; + if (mi->flag & MM_I_HPC) + { + const uint8_t *qseq = qseq0[a->x >> 63]; int i, c; *q = (int32_t)a->y; for (i = *q - 1, c = qseq[*q]; i > 0; --i) - if (qseq[i] != c) break; + if (qseq[i] != c) + break; *q = i + 1; - c = mm_get_hplen_back(mi, a->x<<1>>33, (int32_t)a->x); + c = mm_get_hplen_back(mi, a->x << 1 >> 33, (int32_t)a->x); *r = (int32_t)a->x + 1 - c; - } else { - *r = (int32_t)a->x - (mi->k>>1); - *q = (int32_t)a->y - (mi->k>>1); + } + else + { + *r = (int32_t)a->x - (mi->k >> 1); + *q = (int32_t)a->y - (mi->k >> 1); } } @@ -385,13 +533,17 @@ static int *collect_long_gaps(void *km, int as1, int cnt1, mm128_t *a, int min_g { int i, n, *K; *n_ = 0; - for (i = 1, n = 0; i < cnt1; ++i) { // count the number of gaps longer than min_gap + for (i = 1, n = 0; i < cnt1; ++i) + { // count the number of gaps longer than min_gap int gap = ((int32_t)a[as1 + i].y - a[as1 + i - 1].y) - ((int32_t)a[as1 + i].x - a[as1 + i - 1].x); - if (gap < -min_gap || gap > min_gap) ++n; + if (gap < -min_gap || gap > min_gap) + ++n; } - if (n <= 1) return 0; - K = (int*)kmalloc(km, n * sizeof(int)); - for (i = 1, n = 0; i < cnt1; ++i) { // store the positions of long gaps + if (n <= 1) + return 0; + K = (int *)kmalloc(km, n * sizeof(int)); + for (i = 1, n = 0; i < cnt1; ++i) + { // store the positions of long gaps int gap = ((int32_t)a[as1 + i].y - a[as1 + i - 1].y) - ((int32_t)a[as1 + i].x - a[as1 + i - 1].x); if (gap < -min_gap || gap > min_gap) K[n++] = i; @@ -404,29 +556,39 @@ static void mm_filter_bad_seeds(void *km, int as1, int cnt1, mm128_t *a, int min { int max_st, max_en, n, i, k, max, *K; K = collect_long_gaps(km, as1, cnt1, a, min_gap, &n); - if (K == 0) return; + if (K == 0) + return; max = 0, max_st = max_en = -1; - for (k = 0;; ++k) { // traverse long gaps + for (k = 0;; ++k) + { // traverse long gaps int gap, l, n_ins = 0, n_del = 0, qs, rs, max_diff = 0, max_diff_l = -1; - if (k == n || k >= max_en) { + if (k == n || k >= max_en) + { if (max_en > 0) for (i = K[max_st]; i < K[max_en]; ++i) a[as1 + i].y |= MM_SEED_IGNORE; max = 0, max_st = max_en = -1; - if (k == n) break; + if (k == n) + break; } i = K[k]; gap = ((int32_t)a[as1 + i].y - (int32_t)a[as1 + i - 1].y) - (int32_t)(a[as1 + i].x - a[as1 + i - 1].x); - if (gap > 0) n_ins += gap; - else n_del += -gap; + if (gap > 0) + n_ins += gap; + else + n_del += -gap; qs = (int32_t)a[as1 + i - 1].y; rs = (int32_t)a[as1 + i - 1].x; - for (l = k + 1; l < n && l <= k + max_ext_cnt; ++l) { + for (l = k + 1; l < n && l <= k + max_ext_cnt; ++l) + { int j = K[l], diff; - if ((int32_t)a[as1 + j].y - qs > max_ext_len || (int32_t)a[as1 + j].x - rs > max_ext_len) break; + if ((int32_t)a[as1 + j].y - qs > max_ext_len || (int32_t)a[as1 + j].x - rs > max_ext_len) + break; gap = ((int32_t)a[as1 + j].y - (int32_t)a[as1 + j - 1].y) - (int32_t)(a[as1 + j].x - a[as1 + j - 1].x); - if (gap > 0) n_ins += gap; - else n_del += -gap; + if (gap > 0) + n_ins += gap; + else + n_del += -gap; diff = n_ins + n_del - abs(n_ins - n_del); if (max_diff < diff) max_diff = diff, max_diff_l = l; @@ -441,28 +603,34 @@ static void mm_filter_bad_seeds_alt(void *km, int as1, int cnt1, mm128_t *a, int { int n, k, *K; K = collect_long_gaps(km, as1, cnt1, a, min_gap, &n); - if (K == 0) return; - for (k = 0; k < n;) { + if (K == 0) + return; + for (k = 0; k < n;) + { int i = K[k], l; int gap1 = ((int32_t)a[as1 + i].y - (int32_t)a[as1 + i - 1].y) - ((int32_t)a[as1 + i].x - (int32_t)a[as1 + i - 1].x); int re1 = (int32_t)a[as1 + i].x; int qe1 = (int32_t)a[as1 + i].y; - gap1 = gap1 > 0? gap1 : -gap1; - for (l = k + 1; l < n; ++l) { + gap1 = gap1 > 0 ? gap1 : -gap1; + for (l = k + 1; l < n; ++l) + { int j = K[l], gap2, q_span_pre, rs2, qs2, m; - if ((int32_t)a[as1 + j].y - qe1 > max_ext || (int32_t)a[as1 + j].x - re1 > max_ext) break; + if ((int32_t)a[as1 + j].y - qe1 > max_ext || (int32_t)a[as1 + j].x - re1 > max_ext) + break; gap2 = ((int32_t)a[as1 + j].y - (int32_t)a[as1 + j - 1].y) - (int32_t)(a[as1 + j].x - a[as1 + j - 1].x); q_span_pre = a[as1 + j - 1].y >> 32 & 0xff; rs2 = (int32_t)a[as1 + j - 1].x + q_span_pre; qs2 = (int32_t)a[as1 + j - 1].y + q_span_pre; - m = rs2 - re1 < qs2 - qe1? rs2 - re1 : qs2 - qe1; - gap2 = gap2 > 0? gap2 : -gap2; - if (m > gap1 + gap2) break; + m = rs2 - re1 < qs2 - qe1 ? rs2 - re1 : qs2 - qe1; + gap2 = gap2 > 0 ? gap2 : -gap2; + if (m > gap1 + gap2) + break; re1 = (int32_t)a[as1 + j].x; qe1 = (int32_t)a[as1 + j].y; gap1 = gap2; } - if (l > k + 1) { + if (l > k + 1) + { int j, end = K[l - 1]; for (j = K[k]; j < end; ++j) a[as1 + j].y |= MM_SEED_IGNORE; @@ -477,35 +645,44 @@ static void mm_fix_bad_ends(const mm_reg1_t *r, const mm128_t *a, int bw, int mi { int32_t i, l, m; *as = r->as, *cnt = r->cnt; - if (r->cnt < 3) return; + if (r->cnt < 3) + return; m = l = a[r->as].y >> 32 & 0xff; - for (i = r->as + 1; i < r->as + r->cnt - 1; ++i) { + for (i = r->as + 1; i < r->as + r->cnt - 1; ++i) + { int32_t lq, lr, min, max; int32_t q_span = a[i].y >> 32 & 0xff; - if (a[i].y & MM_SEED_LONG_JOIN) break; - lr = (int32_t)a[i].x - (int32_t)a[i-1].x; - lq = (int32_t)a[i].y - (int32_t)a[i-1].y; - min = lr < lq? lr : lq; - max = lr > lq? lr : lq; - if (max - min > l >> 1) *as = i; + if (a[i].y & MM_SEED_LONG_JOIN) + break; + lr = (int32_t)a[i].x - (int32_t)a[i - 1].x; + lq = (int32_t)a[i].y - (int32_t)a[i - 1].y; + min = lr < lq ? lr : lq; + max = lr > lq ? lr : lq; + if (max - min > l >> 1) + *as = i; l += min; - m += min < q_span? min : q_span; - if (l >= bw << 1 || (m >= min_match && m >= bw) || m >= r->mlen >> 1) break; + m += min < q_span ? min : q_span; + if (l >= bw << 1 || (m >= min_match && m >= bw) || m >= r->mlen >> 1) + break; } *cnt = r->as + r->cnt - *as; m = l = a[r->as + r->cnt - 1].y >> 32 & 0xff; - for (i = r->as + r->cnt - 2; i > *as; --i) { + for (i = r->as + r->cnt - 2; i > *as; --i) + { int32_t lq, lr, min, max; - int32_t q_span = a[i+1].y >> 32 & 0xff; - if (a[i+1].y & MM_SEED_LONG_JOIN) break; - lr = (int32_t)a[i+1].x - (int32_t)a[i].x; - lq = (int32_t)a[i+1].y - (int32_t)a[i].y; - min = lr < lq? lr : lq; - max = lr > lq? lr : lq; - if (max - min > l >> 1) *cnt = i + 1 - *as; + int32_t q_span = a[i + 1].y >> 32 & 0xff; + if (a[i + 1].y & MM_SEED_LONG_JOIN) + break; + lr = (int32_t)a[i + 1].x - (int32_t)a[i].x; + lq = (int32_t)a[i + 1].y - (int32_t)a[i].y; + min = lr < lq ? lr : lq; + max = lr > lq ? lr : lq; + if (max - min > l >> 1) + *cnt = i + 1 - *as; l += min; - m += min < q_span? min : q_span; - if (l >= bw << 1 || (m >= min_match && m >= bw) || m >= r->mlen >> 1) break; + m += min < q_span ? min : q_span; + if (l >= bw << 1 || (m >= min_match && m >= bw) || m >= r->mlen >> 1) + break; } } @@ -514,19 +691,24 @@ static void mm_max_stretch(const mm_reg1_t *r, const mm128_t *a, int32_t *as, in int32_t i, score, max_score, len, max_i, max_len; *as = r->as, *cnt = r->cnt; - if (r->cnt < 2) return; + if (r->cnt < 2) + return; max_score = -1, max_i = -1, max_len = 0; score = a[r->as].y >> 32 & 0xff, len = 1; - for (i = r->as + 1; i < r->as + r->cnt; ++i) { + for (i = r->as + 1; i < r->as + r->cnt; ++i) + { int32_t lq, lr, q_span; q_span = a[i].y >> 32 & 0xff; - lr = (int32_t)a[i].x - (int32_t)a[i-1].x; - lq = (int32_t)a[i].y - (int32_t)a[i-1].y; - if (lq == lr) { - score += lq < q_span? lq : q_span; + lr = (int32_t)a[i].x - (int32_t)a[i - 1].x; + lq = (int32_t)a[i].y - (int32_t)a[i - 1].y; + if (lq == lr) + { + score += lq < q_span ? lq : q_span; ++len; - } else { + } + else + { if (score > max_score) max_score = score, max_len = len, max_i = i - len; score = q_span, len = 1; @@ -540,21 +722,24 @@ static void mm_max_stretch(const mm_reg1_t *r, const mm128_t *a, int32_t *as, in static int mm_seed_ext_score(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, const int8_t mat[25], int qlen, uint8_t *qseq0[2], const mm128_t *a) { uint8_t *qseq, *tseq; - int q_span = a->y>>32&0xff, qs, qe, rs, re, rid, score, q_off, t_off, ext_len = opt->anchor_ext_len; + int q_span = a->y >> 32 & 0xff, qs, qe, rs, re, rid, score, q_off, t_off, ext_len = opt->anchor_ext_len; void *qp; - rid = a->x<<1>>33; + rid = a->x << 1 >> 33; re = (uint32_t)a->x + 1, rs = re - q_span; qe = (uint32_t)a->y + 1, qs = qe - q_span; - rs = rs - ext_len > 0? rs - ext_len : 0; - qs = qs - ext_len > 0? qs - ext_len : 0; - re = re + ext_len < (int32_t)mi->seq[rid].len? re + ext_len : mi->seq[rid].len; - qe = qe + ext_len < qlen? qe + ext_len : qlen; - tseq = (uint8_t*)kmalloc(km, re - rs); - if (opt->flag & MM_F_QSTRAND) { + rs = rs - ext_len > 0 ? rs - ext_len : 0; + qs = qs - ext_len > 0 ? qs - ext_len : 0; + re = re + ext_len < (int32_t)mi->seq[rid].len ? re + ext_len : mi->seq[rid].len; + qe = qe + ext_len < qlen ? qe + ext_len : qlen; + tseq = (uint8_t *)kmalloc(km, re - rs); + if (opt->flag & MM_F_QSTRAND) + { qseq = qseq0[0] + qs; - mm_idx_getseq2(mi, a->x>>63, rid, rs, re, tseq); - } else { - qseq = qseq0[a->x>>63] + qs; + mm_idx_getseq2(mi, a->x >> 63, rid, rs, re, tseq); + } + else + { + qseq = qseq0[a->x >> 63] + qs; mm_idx_getseq(mi, rid, rs, re, tseq); } qp = ksw_ll_qinit(km, 2, qe - qs, qseq, 5, mat); @@ -569,15 +754,18 @@ static void mm_fix_bad_ends_splice(void *km, const mm_mapopt_t *opt, const mm_id int score; double log_gap; *as1 = r->as, *cnt1 = r->cnt; - if (r->cnt < 3) return; + if (r->cnt < 3) + return; log_gap = log((int32_t)a[r->as + 1].x - (int32_t)a[r->as].x); - if ((a[r->as].y>>32&0xff) < log_gap + opt->anchor_ext_shift) { + if ((a[r->as].y >> 32 & 0xff) < log_gap + opt->anchor_ext_shift) + { score = mm_seed_ext_score(km, opt, mi, mat, qlen, qseq0, &a[r->as]); if ((double)score / mat[0] < log_gap + opt->anchor_ext_shift) // a more exact format is "score < log_4(gap) + shift" ++(*as1), --(*cnt1); } log_gap = log((int32_t)a[r->as + r->cnt - 1].x - (int32_t)a[r->as + r->cnt - 2].x); - if ((a[r->as + r->cnt - 1].y>>32&0xff) < log_gap + opt->anchor_ext_shift) { + if ((a[r->as + r->cnt - 1].y >> 32 & 0xff) < log_gap + opt->anchor_ext_shift) + { score = mm_seed_ext_score(km, opt, mi, mat, qlen, qseq0, &a[r->as + r->cnt - 1]); if ((double)score / mat[0] < log_gap + opt->anchor_ext_shift) --(*cnt1); @@ -587,46 +775,59 @@ static void mm_fix_bad_ends_splice(void *km, const mm_mapopt_t *opt, const mm_id static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], mm_reg1_t *r, mm_reg1_t *r2, int n_a, mm128_t *a, ksw_extz_t *ez, int splice_flag) { int is_sr = !!(opt->flag & MM_F_SR), is_splice = !!(opt->flag & MM_F_SPLICE); - int32_t rid = a[r->as].x<<1>>33, rev = a[r->as].x>>63, as1, cnt1; + int32_t rid = a[r->as].x << 1 >> 33, rev = a[r->as].x >> 63, as1, cnt1; uint8_t *tseq, *qseq, *junc; int32_t i, l, bw, bw_long, dropped = 0, extra_flag = 0, rs0, re0, qs0, qe0; int32_t rs, re, qs, qe; int32_t rs1, qs1, re1, qe1; int8_t mat[25]; - if (is_sr) assert(!(mi->flag & MM_I_HPC)); // HPC won't work with SR because with HPC we can't easily tell if there is a gap + if (is_sr) + assert(!(mi->flag & MM_I_HPC)); // HPC won't work with SR because with HPC we can't easily tell if there is a gap r2->cnt = 0; - if (r->cnt == 0) return; + if (r->cnt == 0) + return; ksw_gen_ts_mat(5, mat, opt->a, opt->b, opt->transition, opt->sc_ambi); bw = (int)(opt->bw * 1.5 + 1.); bw_long = (int)(opt->bw_long * 1.5 + 1.); - if (bw_long < bw) bw_long = bw; + if (bw_long < bw) + bw_long = bw; - if (is_sr && !(mi->flag & MM_I_HPC)) { + if (is_sr && !(mi->flag & MM_I_HPC)) + { mm_max_stretch(r, a, &as1, &cnt1); - rs = (int32_t)a[as1].x + 1 - (int32_t)(a[as1].y>>32&0xff); - qs = (int32_t)a[as1].y + 1 - (int32_t)(a[as1].y>>32&0xff); - re = (int32_t)a[as1+cnt1-1].x + 1; - qe = (int32_t)a[as1+cnt1-1].y + 1; - } else { - if (!(opt->flag & MM_F_NO_END_FLT)) { + rs = (int32_t)a[as1].x + 1 - (int32_t)(a[as1].y >> 32 & 0xff); + qs = (int32_t)a[as1].y + 1 - (int32_t)(a[as1].y >> 32 & 0xff); + re = (int32_t)a[as1 + cnt1 - 1].x + 1; + qe = (int32_t)a[as1 + cnt1 - 1].y + 1; + } + else + { + if (!(opt->flag & MM_F_NO_END_FLT)) + { if (is_splice) mm_fix_bad_ends_splice(km, opt, mi, r, mat, qlen, qseq0, a, &as1, &cnt1); else mm_fix_bad_ends(r, a, opt->bw, opt->min_chain_score * 2, &as1, &cnt1); - } else as1 = r->as, cnt1 = r->cnt; - mm_filter_bad_seeds(km, as1, cnt1, a, 10, 40, opt->max_gap>>1, 10); - mm_filter_bad_seeds_alt(km, as1, cnt1, a, 30, opt->max_gap>>1); + } + else + as1 = r->as, cnt1 = r->cnt; + mm_filter_bad_seeds(km, as1, cnt1, a, 10, 40, opt->max_gap >> 1, 10); + mm_filter_bad_seeds_alt(km, as1, cnt1, a, 30, opt->max_gap >> 1); mm_adjust_minier(mi, qseq0, &a[as1], &rs, &qs); mm_adjust_minier(mi, qseq0, &a[as1 + cnt1 - 1], &re, &qe); } assert(cnt1 > 0); - if (is_splice) { - if (splice_flag & MM_F_SPLICE_FOR) extra_flag |= rev? KSW_EZ_SPLICE_REV : KSW_EZ_SPLICE_FOR; - if (splice_flag & MM_F_SPLICE_REV) extra_flag |= rev? KSW_EZ_SPLICE_FOR : KSW_EZ_SPLICE_REV; - if (opt->flag & MM_F_SPLICE_FLANK) extra_flag |= KSW_EZ_SPLICE_FLANK; + if (is_splice) + { + if (splice_flag & MM_F_SPLICE_FOR) + extra_flag |= rev ? KSW_EZ_SPLICE_REV : KSW_EZ_SPLICE_FOR; + if (splice_flag & MM_F_SPLICE_REV) + extra_flag |= rev ? KSW_EZ_SPLICE_FOR : KSW_EZ_SPLICE_REV; + if (opt->flag & MM_F_SPLICE_FLANK) + extra_flag |= KSW_EZ_SPLICE_FLANK; } /* Look for the start and end of regions to perform DP. This sounds easy @@ -634,88 +835,114 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int * clippings and lose alignable sequences. Excessively large regions * occasionally lead to large overlaps between two chains and may cause * loss of alignments in corner cases. */ - if (is_sr) { + if (is_sr) + { qs0 = 0, qe0 = qlen; l = qs; - l += l * opt->a + opt->end_bonus > opt->q? (l * opt->a + opt->end_bonus - opt->q) / opt->e : 0; - rs0 = rs - l > 0? rs - l : 0; + l += l * opt->a + opt->end_bonus > opt->q ? (l * opt->a + opt->end_bonus - opt->q) / opt->e : 0; + rs0 = rs - l > 0 ? rs - l : 0; l = qlen - qe; - l += l * opt->a + opt->end_bonus > opt->q? (l * opt->a + opt->end_bonus - opt->q) / opt->e : 0; - re0 = re + l < (int32_t)mi->seq[rid].len? re + l : mi->seq[rid].len; - } else { + l += l * opt->a + opt->end_bonus > opt->q ? (l * opt->a + opt->end_bonus - opt->q) / opt->e : 0; + re0 = re + l < (int32_t)mi->seq[rid].len ? re + l : mi->seq[rid].len; + } + else + { // compute rs0 and qs0 - rs0 = (int32_t)a[r->as].x + 1 - (int32_t)(a[r->as].y>>32&0xff); - qs0 = (int32_t)a[r->as].y + 1 - (int32_t)(a[r->as].y>>32&0xff); - if (rs0 < 0) rs0 = 0; // this may happen when HPC is in use + rs0 = (int32_t)a[r->as].x + 1 - (int32_t)(a[r->as].y >> 32 & 0xff); + qs0 = (int32_t)a[r->as].y + 1 - (int32_t)(a[r->as].y >> 32 & 0xff); + if (rs0 < 0) + rs0 = 0; // this may happen when HPC is in use assert(qs0 >= 0); // this should never happen, or it is logic error rs1 = qs1 = 0; - for (i = r->as - 1, l = 0; i >= 0 && a[i].x>>32 == a[r->as].x>>32; --i) { // inspect nearby seeds - int32_t x = (int32_t)a[i].x + 1 - (int32_t)(a[i].y>>32&0xff); - int32_t y = (int32_t)a[i].y + 1 - (int32_t)(a[i].y>>32&0xff); - if (x < rs0 && y < qs0) { - if (++l > opt->min_cnt) { - l = rs0 - x > qs0 - y? rs0 - x : qs0 - y; + for (i = r->as - 1, l = 0; i >= 0 && a[i].x >> 32 == a[r->as].x >> 32; --i) + { // inspect nearby seeds + int32_t x = (int32_t)a[i].x + 1 - (int32_t)(a[i].y >> 32 & 0xff); + int32_t y = (int32_t)a[i].y + 1 - (int32_t)(a[i].y >> 32 & 0xff); + if (x < rs0 && y < qs0) + { + if (++l > opt->min_cnt) + { + l = rs0 - x > qs0 - y ? rs0 - x : qs0 - y; rs1 = rs0 - l, qs1 = qs0 - l; - if (rs1 < 0) rs1 = 0; // not strictly necessary; better have this guard for explicit + if (rs1 < 0) + rs1 = 0; // not strictly necessary; better have this guard for explicit break; } } } - if (qs > 0 && rs > 0) { - l = qs < opt->max_gap? qs : opt->max_gap; - qs1 = qs1 > qs - l? qs1 : qs - l; - qs0 = qs0 < qs1? qs0 : qs1; // at least include qs0 - l += l * opt->a > opt->q? (l * opt->a - opt->q) / opt->e : 0; - l = l < opt->max_gap? l : opt->max_gap; - l = l < rs? l : rs; - rs1 = rs1 > rs - l? rs1 : rs - l; - rs0 = rs0 < rs1? rs0 : rs1; - rs0 = rs0 < rs? rs0 : rs; - } else rs0 = rs, qs0 = qs; + if (qs > 0 && rs > 0) + { + l = qs < opt->max_gap ? qs : opt->max_gap; + qs1 = qs1 > qs - l ? qs1 : qs - l; + qs0 = qs0 < qs1 ? qs0 : qs1; // at least include qs0 + l += l * opt->a > opt->q ? (l * opt->a - opt->q) / opt->e : 0; + l = l < opt->max_gap ? l : opt->max_gap; + l = l < rs ? l : rs; + rs1 = rs1 > rs - l ? rs1 : rs - l; + rs0 = rs0 < rs1 ? rs0 : rs1; + rs0 = rs0 < rs ? rs0 : rs; + } + else + rs0 = rs, qs0 = qs; // compute re0 and qe0 re0 = (int32_t)a[r->as + r->cnt - 1].x + 1; qe0 = (int32_t)a[r->as + r->cnt - 1].y + 1; re1 = mi->seq[rid].len, qe1 = qlen; - for (i = r->as + r->cnt, l = 0; i < n_a && a[i].x>>32 == a[r->as].x>>32; ++i) { // inspect nearby seeds + for (i = r->as + r->cnt, l = 0; i < n_a && a[i].x >> 32 == a[r->as].x >> 32; ++i) + { // inspect nearby seeds int32_t x = (int32_t)a[i].x + 1; int32_t y = (int32_t)a[i].y + 1; - if (x > re0 && y > qe0) { - if (++l > opt->min_cnt) { - l = x - re0 > y - qe0? x - re0 : y - qe0; + if (x > re0 && y > qe0) + { + if (++l > opt->min_cnt) + { + l = x - re0 > y - qe0 ? x - re0 : y - qe0; re1 = re0 + l, qe1 = qe0 + l; break; } } } - if (qe < qlen && re < (int32_t)mi->seq[rid].len) { - l = qlen - qe < opt->max_gap? qlen - qe : opt->max_gap; - qe1 = qe1 < qe + l? qe1 : qe + l; - qe0 = qe0 > qe1? qe0 : qe1; // at least include qe0 - l += l * opt->a > opt->q? (l * opt->a - opt->q) / opt->e : 0; - l = l < opt->max_gap? l : opt->max_gap; - l = l < (int32_t)mi->seq[rid].len - re? l : mi->seq[rid].len - re; - re1 = re1 < re + l? re1 : re + l; - re0 = re0 > re1? re0 : re1; - } else re0 = re, qe0 = qe; + if (qe < qlen && re < (int32_t)mi->seq[rid].len) + { + l = qlen - qe < opt->max_gap ? qlen - qe : opt->max_gap; + qe1 = qe1 < qe + l ? qe1 : qe + l; + qe0 = qe0 > qe1 ? qe0 : qe1; // at least include qe0 + l += l * opt->a > opt->q ? (l * opt->a - opt->q) / opt->e : 0; + l = l < opt->max_gap ? l : opt->max_gap; + l = l < (int32_t)mi->seq[rid].len - re ? l : mi->seq[rid].len - re; + re1 = re1 < re + l ? re1 : re + l; + re0 = re0 > re1 ? re0 : re1; + } + else + re0 = re, qe0 = qe; } - if (a[r->as].y & MM_SEED_SELF) { - int max_ext = r->qs > r->rs? r->qs - r->rs : r->rs - r->qs; - if (r->rs - rs0 > max_ext) rs0 = r->rs - max_ext; - if (r->qs - qs0 > max_ext) qs0 = r->qs - max_ext; - max_ext = r->qe > r->re? r->qe - r->re : r->re - r->qe; - if (re0 - r->re > max_ext) re0 = r->re + max_ext; - if (qe0 - r->qe > max_ext) qe0 = r->qe + max_ext; + if (a[r->as].y & MM_SEED_SELF) + { + int max_ext = r->qs > r->rs ? r->qs - r->rs : r->rs - r->qs; + if (r->rs - rs0 > max_ext) + rs0 = r->rs - max_ext; + if (r->qs - qs0 > max_ext) + qs0 = r->qs - max_ext; + max_ext = r->qe > r->re ? r->qe - r->re : r->re - r->qe; + if (re0 - r->re > max_ext) + re0 = r->re + max_ext; + if (qe0 - r->qe > max_ext) + qe0 = r->qe + max_ext; } assert(re0 > rs0); - tseq = (uint8_t*)kmalloc(km, re0 - rs0); - junc = (uint8_t*)kmalloc(km, re0 - rs0); + tseq = (uint8_t *)kmalloc(km, re0 - rs0); + junc = (uint8_t *)kmalloc(km, re0 - rs0); - if (qs > 0 && rs > 0) { // left extension; probably the condition can be changed to "qs > qs0 && rs > rs0" - if (opt->flag & MM_F_QSTRAND) { + if (qs > 0 && rs > 0) + { // left extension; probably the condition can be changed to "qs > qs0 && rs > rs0" + if (opt->flag & MM_F_QSTRAND) + { qseq = &qseq0[0][qs0]; mm_idx_getseq2(mi, rev, rid, rs0, rs, tseq); - } else { + } + else + { qseq = &qseq0[rev][qs0]; mm_idx_getseq(mi, rid, rs0, rs, tseq); } @@ -723,110 +950,146 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int mm_seq_rev(qs - qs0, qseq); mm_seq_rev(rs - rs0, tseq); mm_seq_rev(rs - rs0, junc); - mm_align_pair(km, opt, qs - qs0, qseq, rs - rs0, tseq, junc, mat, bw, opt->end_bonus, r->split_inv? opt->zdrop_inv : opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY|KSW_EZ_RIGHT|KSW_EZ_REV_CIGAR, ez); - if (ez->n_cigar > 0) { + mm_align_pair(km, opt, qs - qs0, qseq, rs - rs0, tseq, junc, mat, bw, opt->end_bonus, r->split_inv ? opt->zdrop_inv : opt->zdrop, extra_flag | KSW_EZ_EXTZ_ONLY | KSW_EZ_RIGHT | KSW_EZ_REV_CIGAR, ez); + if (ez->n_cigar > 0) + { mm_append_cigar(r, ez->n_cigar, ez->cigar); r->p->dp_score += ez->max; } - rs1 = rs - (ez->reach_end? ez->mqe_t + 1 : ez->max_t + 1); - qs1 = qs - (ez->reach_end? qs - qs0 : ez->max_q + 1); + rs1 = rs - (ez->reach_end ? ez->mqe_t + 1 : ez->max_t + 1); + qs1 = qs - (ez->reach_end ? qs - qs0 : ez->max_q + 1); mm_seq_rev(qs - qs0, qseq); - } else rs1 = rs, qs1 = qs; + } + else + rs1 = rs, qs1 = qs; re1 = rs, qe1 = qs; assert(qs1 >= 0 && rs1 >= 0); - for (i = is_sr? cnt1 - 1 : 1; i < cnt1; ++i) { // gap filling - if ((a[as1+i].y & (MM_SEED_IGNORE|MM_SEED_TANDEM)) && i != cnt1 - 1) continue; - if (is_sr && !(mi->flag & MM_I_HPC)) { + for (i = is_sr ? cnt1 - 1 : 1; i < cnt1; ++i) + { // gap filling + if ((a[as1 + i].y & (MM_SEED_IGNORE | MM_SEED_TANDEM)) && i != cnt1 - 1) + continue; + if (is_sr && !(mi->flag & MM_I_HPC)) + { re = (int32_t)a[as1 + i].x + 1; qe = (int32_t)a[as1 + i].y + 1; - } else mm_adjust_minier(mi, qseq0, &a[as1 + i], &re, &qe); + } + else + mm_adjust_minier(mi, qseq0, &a[as1 + i], &re, &qe); re1 = re, qe1 = qe; - if (i == cnt1 - 1 || (a[as1+i].y&MM_SEED_LONG_JOIN) || (qe - qs >= opt->min_ksw_len && re - rs >= opt->min_ksw_len)) { + if (i == cnt1 - 1 || (a[as1 + i].y & MM_SEED_LONG_JOIN) || (qe - qs >= opt->min_ksw_len && re - rs >= opt->min_ksw_len)) + { int j, bw1 = bw_long, zdrop_code; - if (a[as1+i].y & MM_SEED_LONG_JOIN) - bw1 = qe - qs > re - rs? qe - qs : re - rs; + if (a[as1 + i].y & MM_SEED_LONG_JOIN) + bw1 = qe - qs > re - rs ? qe - qs : re - rs; // perform alignment - if (opt->flag & MM_F_QSTRAND) { + if (opt->flag & MM_F_QSTRAND) + { qseq = &qseq0[0][qs]; mm_idx_getseq2(mi, rev, rid, rs, re, tseq); - } else { + } + else + { qseq = &qseq0[rev][qs]; mm_idx_getseq(mi, rid, rs, re, tseq); } mm_idx_bed_junc(mi, rid, rs, re, junc); - if (is_sr) { // perform ungapped alignment + if (is_sr) + { // perform ungapped alignment assert(qe - qs == re - rs); ksw_reset_extz(ez); - for (j = 0, ez->score = 0; j < qe - qs; ++j) { - if (qseq[j] >= 4 || tseq[j] >= 4) ez->score += opt->e2; - else ez->score += qseq[j] == tseq[j]? opt->a : -opt->b; + for (j = 0, ez->score = 0; j < qe - qs; ++j) + { + if (qseq[j] >= 4 || tseq[j] >= 4) + ez->score += opt->e2; + else + ez->score += qseq[j] == tseq[j] ? opt->a : -opt->b; } ez->cigar = ksw_push_cigar(km, &ez->n_cigar, &ez->m_cigar, ez->cigar, MM_CIGAR_MATCH, qe - qs); - } else { // perform normal gapped alignment - mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, junc, mat, bw1, -1, opt->zdrop, extra_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop + } + else + { // perform normal gapped alignment + mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, junc, mat, bw1, -1, opt->zdrop, extra_flag | KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop } // test Z-drop and inversion Z-drop if ((zdrop_code = mm_test_zdrop(km, opt, qseq, tseq, ez->n_cigar, ez->cigar, mat)) != 0) - mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, junc, mat, bw1, -1, zdrop_code == 2? opt->zdrop_inv : opt->zdrop, extra_flag, ez); // second pass: lift approximate + mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, junc, mat, bw1, -1, zdrop_code == 2 ? opt->zdrop_inv : opt->zdrop, extra_flag, ez); // second pass: lift approximate // update CIGAR if (ez->n_cigar > 0) mm_append_cigar(r, ez->n_cigar, ez->cigar); - if (ez->zdropped) { // truncated by Z-drop; TODO: sometimes Z-drop kicks in because the next seed placement is wrong. This can be fixed in principle. - if (!r->p) { + if (ez->zdropped) + { // truncated by Z-drop; TODO: sometimes Z-drop kicks in because the next seed placement is wrong. This can be fixed in principle. + if (!r->p) + { assert(ez->n_cigar == 0); - uint32_t capacity = sizeof(mm_extra_t)/4; + uint32_t capacity = sizeof(mm_extra_t) / 4; kroundup32(capacity); - r->p = (mm_extra_t*)calloc(capacity, 4); + r->p = (mm_extra_t *)calloc(capacity, 4); r->p->capacity = capacity; } for (j = i - 1; j >= 0; --j) if ((int32_t)a[as1 + j].x <= rs + ez->max_t) break; dropped = 1; - if (j < 0) j = 0; + if (j < 0) + j = 0; r->p->dp_score += ez->max; re1 = rs + (ez->max_t + 1); qe1 = qs + (ez->max_q + 1); - if (cnt1 - (j + 1) >= opt->min_cnt) { - mm_split_reg(r, r2, as1 + j + 1 - r->as, qlen, a, !!(opt->flag&MM_F_QSTRAND)); - if (zdrop_code == 2) r2->split_inv = 1; + if (cnt1 - (j + 1) >= opt->min_cnt) + { + mm_split_reg(r, r2, as1 + j + 1 - r->as, qlen, a, !!(opt->flag & MM_F_QSTRAND)); + if (zdrop_code == 2) + r2->split_inv = 1; } break; - } else r->p->dp_score += ez->score; + } + else + r->p->dp_score += ez->score; rs = re, qs = qe; } } - if (!dropped && qe < qe0 && re < re0) { // right extension - if (opt->flag & MM_F_QSTRAND) { + if (!dropped && qe < qe0 && re < re0) + { // right extension + if (opt->flag & MM_F_QSTRAND) + { qseq = &qseq0[0][qe]; mm_idx_getseq2(mi, rev, rid, re, re0, tseq); - } else { + } + else + { qseq = &qseq0[rev][qe]; mm_idx_getseq(mi, rid, re, re0, tseq); } mm_idx_bed_junc(mi, rid, re, re0, junc); - mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, junc, mat, bw, opt->end_bonus, opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY, ez); - if (ez->n_cigar > 0) { + mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, junc, mat, bw, opt->end_bonus, opt->zdrop, extra_flag | KSW_EZ_EXTZ_ONLY, ez); + if (ez->n_cigar > 0) + { mm_append_cigar(r, ez->n_cigar, ez->cigar); r->p->dp_score += ez->max; } - re1 = re + (ez->reach_end? ez->mqe_t + 1 : ez->max_t + 1); - qe1 = qe + (ez->reach_end? qe0 - qe : ez->max_q + 1); + re1 = re + (ez->reach_end ? ez->mqe_t + 1 : ez->max_t + 1); + qe1 = qe + (ez->reach_end ? qe0 - qe : ez->max_q + 1); } assert(qe1 <= qlen); r->rs = rs1, r->re = re1; - if (!rev || (opt->flag & MM_F_QSTRAND)) r->qs = qs1, r->qe = qe1; - else r->qs = qlen - qe1, r->qe = qlen - qs1; + if (!rev || (opt->flag & MM_F_QSTRAND)) + r->qs = qs1, r->qe = qe1; + else + r->qs = qlen - qe1, r->qe = qlen - qs1; assert(re1 - rs1 <= re0 - rs0); - if (r->p) { - if (opt->flag & MM_F_QSTRAND) { + if (r->p) + { + if (opt->flag & MM_F_QSTRAND) + { mm_idx_getseq2(mi, r->rev, rid, rs1, re1, tseq); qseq = &qseq0[0][qs1]; - } else { + } + else + { mm_idx_getseq(mi, rid, rs1, re1, tseq); qseq = &qseq0[r->rev][qs1]; } @@ -847,19 +1110,25 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i void *qp; memset(r_inv, 0, sizeof(mm_reg1_t)); - if (!(r1->split&1) || !(r2->split&2)) return 0; - if (r1->id != r1->parent && r1->parent != MM_PARENT_TMP_PRI) return 0; - if (r2->id != r2->parent && r2->parent != MM_PARENT_TMP_PRI) return 0; - if (r1->rid != r2->rid || r1->rev != r2->rev) return 0; - ql = r1->rev? r1->qs - r2->qe : r2->qs - r1->qe; + if (!(r1->split & 1) || !(r2->split & 2)) + return 0; + if (r1->id != r1->parent && r1->parent != MM_PARENT_TMP_PRI) + return 0; + if (r2->id != r2->parent && r2->parent != MM_PARENT_TMP_PRI) + return 0; + if (r1->rid != r2->rid || r1->rev != r2->rev) + return 0; + ql = r1->rev ? r1->qs - r2->qe : r2->qs - r1->qe; tl = r2->rs - r1->re; - if (ql < opt->min_chain_score || ql > opt->max_gap) return 0; - if (tl < opt->min_chain_score || tl > opt->max_gap) return 0; + if (ql < opt->min_chain_score || ql > opt->max_gap) + return 0; + if (tl < opt->min_chain_score || tl > opt->max_gap) + return 0; ksw_gen_ts_mat(5, mat, opt->a, opt->b, opt->transition, opt->sc_ambi); - tseq = (uint8_t*)kmalloc(km, tl); + tseq = (uint8_t *)kmalloc(km, tl); mm_idx_getseq(mi, r1->rid, r1->re, r2->rs, tseq); - qseq = r1->rev? &qseq0[0][r2->qe] : &qseq0[1][qlen - r2->qs]; + qseq = r1->rev ? &qseq0[0][r2->qe] : &qseq0[1][qlen - r2->qs]; mm_seq_rev(ql, qseq); mm_seq_rev(tl, tseq); @@ -868,10 +1137,12 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i kfree(km, qp); mm_seq_rev(ql, qseq); mm_seq_rev(tl, tseq); - if (score < opt->min_dp_max) goto end_align1_inv; + if (score < opt->min_dp_max) + goto end_align1_inv; q_off = ql - (q_off + 1), t_off = tl - (t_off + 1); mm_align_pair(km, opt, ql - q_off, qseq + q_off, tl - t_off, tseq + t_off, 0, mat, (int)(opt->bw * 1.5), -1, opt->zdrop, KSW_EZ_EXTZ_ONLY, ez); - if (ez->n_cigar == 0) goto end_align1_inv; // should never be here + if (ez->n_cigar == 0) + goto end_align1_inv; // should never be here mm_append_cigar(r_inv, ez->n_cigar, ez->cigar); r_inv->p->dp_score = ez->max; r_inv->id = -1; @@ -880,10 +1151,13 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i r_inv->rev = !r1->rev; r_inv->rid = r1->rid; r_inv->div = -1.0f; - if (r_inv->rev == 0) { + if (r_inv->rev == 0) + { r_inv->qs = r2->qe + q_off; r_inv->qe = r_inv->qs + ez->max_q + 1; - } else { + } + else + { r_inv->qe = r2->qs - q_off; r_inv->qs = r_inv->qe - (ez->max_q + 1); } @@ -898,7 +1172,7 @@ end_align1_inv: static inline mm_reg1_t *mm_insert_reg(const mm_reg1_t *r, int i, int *n_regs, mm_reg1_t *regs) { - regs = (mm_reg1_t*)realloc(regs, (*n_regs + 1) * sizeof(mm_reg1_t)); + regs = (mm_reg1_t *)realloc(regs, (*n_regs + 1) * sizeof(mm_reg1_t)); if (i + 1 != *n_regs) memmove(®s[i + 2], ®s[i + 1], sizeof(mm_reg1_t) * (*n_regs - i - 1)); regs[i + 1] = *r; @@ -911,8 +1185,10 @@ static inline void mm_count_gaps(const mm_reg1_t *r, int32_t *n_gap_, int32_t *n uint32_t i; int32_t n_gapo = 0, n_gap = 0; *n_gap_ = *n_gapo_ = -1; - if (r->p == 0) return; - for (i = 0; i < r->p->n_cigar; ++i) { + if (r->p == 0) + return; + for (i = 0; i < r->p->n_cigar; ++i) + { int32_t op = r->p->cigar[i] & 0xf, len = r->p->cigar[i] >> 4; if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL) ++n_gapo, n_gap += len; @@ -923,7 +1199,8 @@ static inline void mm_count_gaps(const mm_reg1_t *r, int32_t *n_gap_, int32_t *n double mm_event_identity(const mm_reg1_t *r) { int32_t n_gap, n_gapo; - if (r->p == 0) return -1.0f; + if (r->p == 0) + return -1.0f; mm_count_gaps(r, &n_gap, &n_gapo); return (double)r->mlen / (r->blen + r->p->n_ambi - n_gap + n_gapo); } @@ -933,10 +1210,13 @@ static int32_t mm_recal_max_dp(const mm_reg1_t *r, double b2, int32_t match_sc) uint32_t i; int32_t n_gap = 0, n_gapo = 0, n_mis; double gap_cost = 0.0; - if (r->p == 0) return -1; - for (i = 0; i < r->p->n_cigar; ++i) { + if (r->p == 0) + return -1; + for (i = 0; i < r->p->n_cigar; ++i) + { int32_t op = r->p->cigar[i] & 0xf, len = r->p->cigar[i] >> 4; - if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL) { + if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL) + { gap_cost += b2 + (double)mg_log2(1.0 + len); ++n_gapo, n_gap += len; } @@ -949,25 +1229,38 @@ void mm_update_dp_max(int qlen, int n_regs, mm_reg1_t *regs, float frac, int a, { int32_t max = -1, max2 = -1, i, max_i = -1; double div, b2; - if (n_regs < 2) return; - for (i = 0; i < n_regs; ++i) { + if (n_regs < 2) + return; + for (i = 0; i < n_regs; ++i) + { mm_reg1_t *r = ®s[i]; - if (r->p == 0) continue; - if (r->p->dp_max > max) max2 = max, max = r->p->dp_max, max_i = i; - else if (r->p->dp_max > max2) max2 = r->p->dp_max; + if (r->p == 0) + continue; + if (r->p->dp_max > max) + max2 = max, max = r->p->dp_max, max_i = i; + else if (r->p->dp_max > max2) + max2 = r->p->dp_max; } - if (max_i < 0 || max < 0 || max2 < 0) return; - if (regs[max_i].qe - regs[max_i].qs < (double)qlen * frac) return; - if (max2 < (double)max * frac) return; + if (max_i < 0 || max < 0 || max2 < 0) + return; + if (regs[max_i].qe - regs[max_i].qs < (double)qlen * frac) + return; + if (max2 < (double)max * frac) + return; div = 1. - mm_event_identity(®s[max_i]); - if (div < 0.02) div = 0.02; + if (div < 0.02) + div = 0.02; b2 = 0.5 / div; // max value: 25 - if (b2 * a < b) b2 = (double)a / b; - for (i = 0; i < n_regs; ++i) { + if (b2 * a < b) + b2 = (double)a / b; + for (i = 0; i < n_regs; ++i) + { mm_reg1_t *r = ®s[i]; - if (r->p == 0) continue; + if (r->p == 0) + continue; r->p->dp_max = mm_recal_max_dp(r, b2, a); - if (r->p->dp_max < 0) r->p->dp_max = 0; + if (r->p->dp_max < 0) + r->p->dp_max = 0; } } @@ -979,43 +1272,57 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m ksw_extz_t ez; // encode the query sequence - qseq0[0] = (uint8_t*)kmalloc(km, qlen * 2); + qseq0[0] = (uint8_t *)kmalloc(km, qlen * 2); qseq0[1] = qseq0[0] + qlen; - for (i = 0; i < qlen; ++i) { + for (i = 0; i < qlen; ++i) + { qseq0[0][i] = seq_nt4_table[(uint8_t)qstr[i]]; - qseq0[1][qlen - 1 - i] = qseq0[0][i] < 4? 3 - qseq0[0][i] : 4; + qseq0[1][qlen - 1 - i] = qseq0[0][i] < 4 ? 3 - qseq0[0][i] : 4; } // align through seed hits n_a = mm_squeeze_a(km, n_regs, regs, a); memset(&ez, 0, sizeof(ksw_extz_t)); - for (i = 0; i < n_regs; ++i) { + for (i = 0; i < n_regs; ++i) + { mm_reg1_t r2; - if ((opt->flag&MM_F_SPLICE) && (opt->flag&MM_F_SPLICE_FOR) && (opt->flag&MM_F_SPLICE_REV)) { // then do two rounds of alignments for both strands + if ((opt->flag & MM_F_SPLICE) && (opt->flag & MM_F_SPLICE_FOR) && (opt->flag & MM_F_SPLICE_REV)) + { // then do two rounds of alignments for both strands mm_reg1_t s[2], s2[2]; int which, trans_strand; s[0] = s[1] = regs[i]; mm_align1(km, opt, mi, qlen, qseq0, &s[0], &s2[0], n_a, a, &ez, MM_F_SPLICE_FOR); mm_align1(km, opt, mi, qlen, qseq0, &s[1], &s2[1], n_a, a, &ez, MM_F_SPLICE_REV); - if (s[0].p->dp_score > s[1].p->dp_score) which = 0, trans_strand = 1; - else if (s[0].p->dp_score < s[1].p->dp_score) which = 1, trans_strand = 2; - else trans_strand = 3, which = (qlen + s[0].p->dp_score) & 1; // randomly choose a strand, effectively - if (which == 0) { + if (s[0].p->dp_score > s[1].p->dp_score) + which = 0, trans_strand = 1; + else if (s[0].p->dp_score < s[1].p->dp_score) + which = 1, trans_strand = 2; + else + trans_strand = 3, which = (qlen + s[0].p->dp_score) & 1; // randomly choose a strand, effectively + if (which == 0) + { regs[i] = s[0], r2 = s2[0]; free(s[1].p); - } else { + } + else + { regs[i] = s[1], r2 = s2[1]; free(s[0].p); } regs[i].p->trans_strand = trans_strand; - } else { // one round of alignment - mm_align1(km, opt, mi, qlen, qseq0, ®s[i], &r2, n_a, a, &ez, opt->flag); - if (opt->flag&MM_F_SPLICE) - regs[i].p->trans_strand = opt->flag&MM_F_SPLICE_FOR? 1 : 2; } - if (r2.cnt > 0) regs = mm_insert_reg(&r2, i, &n_regs, regs); - if (i > 0 && regs[i].split_inv && !(opt->flag & MM_F_NO_INV)) { - if (mm_align1_inv(km, opt, mi, qlen, qseq0, ®s[i-1], ®s[i], &r2, &ez)) { + else + { // one round of alignment + mm_align1(km, opt, mi, qlen, qseq0, ®s[i], &r2, n_a, a, &ez, opt->flag); + if (opt->flag & MM_F_SPLICE) + regs[i].p->trans_strand = opt->flag & MM_F_SPLICE_FOR ? 1 : 2; + } + if (r2.cnt > 0) + regs = mm_insert_reg(&r2, i, &n_regs, regs); + if (i > 0 && regs[i].split_inv && !(opt->flag & MM_F_NO_INV)) + { + if (mm_align1_inv(km, opt, mi, qlen, qseq0, ®s[i - 1], ®s[i], &r2, &ez)) + { regs = mm_insert_reg(&r2, i, &n_regs, regs); ++i; // skip the inserted INV alignment } @@ -1025,7 +1332,8 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m kfree(km, qseq0[0]); kfree(km, ez.cigar); mm_filter_regs(opt, qlen, n_regs_, regs); - if (!(opt->flag&MM_F_SR) && !opt->split_prefix && qlen >= opt->rank_min_len) { + if (!(opt->flag & MM_F_SR) && !opt->split_prefix && qlen >= opt->rank_min_len) + { mm_update_dp_max(qlen, *n_regs_, regs, opt->rank_frac, opt->a, opt->b); mm_filter_regs(opt, qlen, n_regs_, regs); } diff --git a/index.c b/index.c index 8a3e309..37e7673 100644 --- a/index.c +++ b/index.c @@ -15,33 +15,36 @@ #include "kvec.h" #include "khash.h" -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF extern int64_t get_mseconds(); extern int64_t time_mm_idx_reader_read; #endif -#define idx_hash(a) ((a)>>1) -#define idx_eq(a, b) ((a)>>1 == (b)>>1) +#define idx_hash(a) ((a) >> 1) +#define idx_eq(a, b) ((a) >> 1 == (b) >> 1) KHASH_INIT(idx, uint64_t, uint64_t, 1, idx_hash, idx_eq) typedef khash_t(idx) idxhash_t; KHASH_MAP_INIT_STR(str, uint32_t) -#define kroundup64(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, (x)|=(x)>>32, ++(x)) +#define kroundup64(x) (--(x), (x) |= (x) >> 1, (x) |= (x) >> 2, (x) |= (x) >> 4, (x) |= (x) >> 8, (x) |= (x) >> 16, (x) |= (x) >> 32, ++(x)) -typedef struct mm_idx_bucket_s { - mm128_v a; // (minimizer, position) array - int32_t n; // size of the _p_ array +typedef struct mm_idx_bucket_s +{ + mm128_v a; // (minimizer, position) array + int32_t n; // size of the _p_ array uint64_t *p; // position array for minimizers appearing >1 times - void *h; // hash table indexing _p_ and minimizers appearing once + void *h; // hash table indexing _p_ and minimizers appearing once } mm_idx_bucket_t; -typedef struct { +typedef struct +{ int32_t st, en, max; // max is not used for now - int32_t score:30, strand:2; + int32_t score : 30, strand : 2; } mm_idx_intv1_t; -typedef struct mm_idx_intv_s { +typedef struct mm_idx_intv_s +{ int32_t n, m; mm_idx_intv1_t *a; } mm_idx_intv_t; @@ -49,56 +52,74 @@ typedef struct mm_idx_intv_s { mm_idx_t *mm_idx_init(int w, int k, int b, int flag) { mm_idx_t *mi; - if (k*2 < b) b = k * 2; - if (w < 1) w = 1; - mi = (mm_idx_t*)calloc(1, sizeof(mm_idx_t)); + if (k * 2 < b) + b = k * 2; + if (w < 1) + w = 1; + mi = (mm_idx_t *)calloc(1, sizeof(mm_idx_t)); mi->w = w, mi->k = k, mi->b = b, mi->flag = flag; - mi->B = (mm_idx_bucket_t*)calloc(1<km = km_init(); + mi->B = (mm_idx_bucket_t *)calloc(1 << b, sizeof(mm_idx_bucket_t)); + if (!(mm_dbg_flag & 1)) + mi->km = km_init(); return mi; } void mm_idx_destroy(mm_idx_t *mi) { uint32_t i; - if (mi == 0) return; - if (mi->h) kh_destroy(str, (khash_t(str)*)mi->h); - if (mi->B) { - for (i = 0; i < 1U<b; ++i) { + if (mi == 0) + return; + if (mi->h) + kh_destroy(str, (khash_t(str) *)mi->h); + if (mi->B) + { + for (i = 0; i < 1U << mi->b; ++i) + { free(mi->B[i].p); free(mi->B[i].a.a); - kh_destroy(idx, (idxhash_t*)mi->B[i].h); + kh_destroy(idx, (idxhash_t *)mi->B[i].h); } } - if (mi->I) { + if (mi->I) + { for (i = 0; i < mi->n_seq; ++i) free(mi->I[i].a); free(mi->I); } - if (!mi->km) { + if (!mi->km) + { for (i = 0; i < mi->n_seq; ++i) free(mi->seq[i].name); free(mi->seq); - } else km_destroy(mi->km); - free(mi->B); free(mi->S); free(mi); + } + else + km_destroy(mi->km); + free(mi->B); + free(mi->S); + free(mi); } const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n) { - int mask = (1<b) - 1; + int mask = (1 << mi->b) - 1; khint_t k; - mm_idx_bucket_t *b = &mi->B[minier&mask]; - idxhash_t *h = (idxhash_t*)b->h; + mm_idx_bucket_t *b = &mi->B[minier & mask]; + idxhash_t *h = (idxhash_t *)b->h; *n = 0; - if (h == 0) return 0; - k = kh_get(idx, h, minier>>mi->b<<1); - if (k == kh_end(h)) return 0; - if (kh_key(h, k)&1) { // special casing when there is only one k-mer + if (h == 0) + return 0; + k = kh_get(idx, h, minier >> mi->b << 1); + if (k == kh_end(h)) + return 0; + if (kh_key(h, k) & 1) + { // special casing when there is only one k-mer *n = 1; return &kh_val(h, k); - } else { + } + else + { *n = (uint32_t)kh_val(h, k); - return &b->p[kh_val(h, k)>>32]; + return &b->p[kh_val(h, k) >> 32]; } } @@ -107,37 +128,46 @@ void mm_idx_stat(const mm_idx_t *mi) int n = 0, n1 = 0; uint32_t i; uint64_t sum = 0, len = 0; - fprintf(stderr, "[M::%s] kmer size: %d; skip: %d; is_hpc: %d; #seq: %d\n", __func__, mi->k, mi->w, mi->flag&MM_I_HPC, mi->n_seq); + fprintf(stderr, "[M::%s] kmer size: %d; skip: %d; is_hpc: %d; #seq: %d\n", __func__, mi->k, mi->w, mi->flag & MM_I_HPC, mi->n_seq); for (i = 0; i < mi->n_seq; ++i) len += mi->seq[i].len; - for (i = 0; i < 1U<b; ++i) - if (mi->B[i].h) n += kh_size((idxhash_t*)mi->B[i].h); - for (i = 0; i < 1U<b; ++i) { - idxhash_t *h = (idxhash_t*)mi->B[i].h; + for (i = 0; i < 1U << mi->b; ++i) + if (mi->B[i].h) + n += kh_size((idxhash_t *)mi->B[i].h); + for (i = 0; i < 1U << mi->b; ++i) + { + idxhash_t *h = (idxhash_t *)mi->B[i].h; khint_t k; - if (h == 0) continue; + if (h == 0) + continue; for (k = 0; k < kh_end(h); ++k) - if (kh_exist(h, k)) { - sum += kh_key(h, k)&1? 1 : (uint32_t)kh_val(h, k); - if (kh_key(h, k)&1) ++n1; + if (kh_exist(h, k)) + { + sum += kh_key(h, k) & 1 ? 1 : (uint32_t)kh_val(h, k); + if (kh_key(h, k) & 1) + ++n1; } } fprintf(stderr, "[M::%s::%.3f*%.2f] distinct minimizers: %d (%.2f%% are singletons); average occurrences: %.3lf; average spacing: %.3lf; total length: %ld\n", - __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), n, 100.0*n1/n, (double)sum / n, (double)len / sum, (long)len); + __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), n, 100.0 * n1 / n, (double)sum / n, (double)len / sum, (long)len); } int mm_idx_index_name(mm_idx_t *mi) { - khash_t(str) *h; + khash_t(str) * h; uint32_t i; int has_dup = 0, absent; - if (mi->h) return 0; + if (mi->h) + return 0; h = kh_init(str); - for (i = 0; i < mi->n_seq; ++i) { + for (i = 0; i < mi->n_seq; ++i) + { khint_t k; k = kh_put(str, h, mi->seq[i].name, &absent); - if (absent) kh_val(h, k) = i; - else has_dup = 1; + if (absent) + kh_val(h, k) = i; + else + has_dup = 1; } mi->h = h; if (has_dup && mm_verbose >= 2) @@ -147,18 +177,21 @@ int mm_idx_index_name(mm_idx_t *mi) int mm_idx_name2id(const mm_idx_t *mi, const char *name) { - khash_t(str) *h = (khash_t(str)*)mi->h; + khash_t(str) *h = (khash_t(str) *)mi->h; khint_t k; - if (h == 0) return -2; + if (h == 0) + return -2; k = kh_get(str, h, name); - return k == kh_end(h)? -1 : kh_val(h, k); + return k == kh_end(h) ? -1 : kh_val(h, k); } int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq) { uint64_t i, st1, en1; - if (rid >= mi->n_seq || st >= mi->seq[rid].len) return -1; - if (en > mi->seq[rid].len) en = mi->seq[rid].len; + if (rid >= mi->n_seq || st >= mi->seq[rid].len) + return -1; + if (en > mi->seq[rid].len) + en = mi->seq[rid].len; st1 = mi->seq[rid].offset + st; en1 = mi->seq[rid].offset + en; for (i = st1; i < en1; ++i) @@ -170,22 +203,27 @@ int mm_idx_getseq_rev(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en { uint64_t i, st1, en1; const mm_idx_seq_t *s; - if (rid >= mi->n_seq || st >= mi->seq[rid].len) return -1; + if (rid >= mi->n_seq || st >= mi->seq[rid].len) + return -1; s = &mi->seq[rid]; - if (en > s->len) en = s->len; + if (en > s->len) + en = s->len; st1 = s->offset + (s->len - en); en1 = s->offset + (s->len - st); - for (i = st1; i < en1; ++i) { + for (i = st1; i < en1; ++i) + { uint8_t c = mm_seq4_get(mi->S, i); - seq[en1 - i - 1] = c < 4? 3 - c : c; + seq[en1 - i - 1] = c < 4 ? 3 - c : c; } return en - st; } int mm_idx_getseq2(const mm_idx_t *mi, int is_rev, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq) { - if (is_rev) return mm_idx_getseq_rev(mi, rid, st, en, seq); - else return mm_idx_getseq(mi, rid, st, en, seq); + if (is_rev) + return mm_idx_getseq_rev(mi, rid, st, en, seq); + else + return mm_idx_getseq(mi, rid, st, en, seq); } int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f) @@ -194,16 +232,22 @@ int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f) size_t n = 0; uint32_t thres; khint_t *a, k; - if (f <= 0.) return INT32_MAX; - for (i = 0; i < 1<b; ++i) - if (mi->B[i].h) n += kh_size((idxhash_t*)mi->B[i].h); - a = (uint32_t*)malloc(n * 4); - for (i = n = 0; i < 1<b; ++i) { - idxhash_t *h = (idxhash_t*)mi->B[i].h; - if (h == 0) continue; - for (k = 0; k < kh_end(h); ++k) { - if (!kh_exist(h, k)) continue; - a[n++] = kh_key(h, k)&1? 1 : (uint32_t)kh_val(h, k); + if (f <= 0.) + return INT32_MAX; + for (i = 0; i < 1 << mi->b; ++i) + if (mi->B[i].h) + n += kh_size((idxhash_t *)mi->B[i].h); + a = (uint32_t *)malloc(n * 4); + for (i = n = 0; i < 1 << mi->b; ++i) + { + idxhash_t *h = (idxhash_t *)mi->B[i].h; + if (h == 0) + continue; + for (k = 0; k < kh_end(h); ++k) + { + if (!kh_exist(h, k)) + continue; + a[n++] = kh_key(h, k) & 1 ? 1 : (uint32_t)kh_val(h, k); } } thres = ks_ksmall_uint32_t(n, a, (uint32_t)((1. - f) * n)) + 1; @@ -220,46 +264,59 @@ static void worker_post(void *g, long i, int tid) int n, n_keys; size_t j, start_a, start_p; idxhash_t *h; - mm_idx_t *mi = (mm_idx_t*)g; + mm_idx_t *mi = (mm_idx_t *)g; mm_idx_bucket_t *b = &mi->B[i]; - if (b->a.n == 0) return; + if (b->a.n == 0) + return; // sort by minimizer radix_sort_128x(b->a.a, b->a.a + b->a.n); // count and preallocate - for (j = 1, n = 1, n_keys = 0, b->n = 0; j <= b->a.n; ++j) { - if (j == b->a.n || b->a.a[j].x>>8 != b->a.a[j-1].x>>8) { + for (j = 1, n = 1, n_keys = 0, b->n = 0; j <= b->a.n; ++j) + { + if (j == b->a.n || b->a.a[j].x >> 8 != b->a.a[j - 1].x >> 8) + { ++n_keys; - if (n > 1) b->n += n; + if (n > 1) + b->n += n; n = 1; - } else ++n; + } + else + ++n; } h = kh_init(idx); kh_resize(idx, h, n_keys); - b->p = (uint64_t*)calloc(b->n, 8); + b->p = (uint64_t *)calloc(b->n, 8); // create the hash table - for (j = 1, n = 1, start_a = start_p = 0; j <= b->a.n; ++j) { - if (j == b->a.n || b->a.a[j].x>>8 != b->a.a[j-1].x>>8) { + for (j = 1, n = 1, start_a = start_p = 0; j <= b->a.n; ++j) + { + if (j == b->a.n || b->a.a[j].x >> 8 != b->a.a[j - 1].x >> 8) + { khint_t itr; int absent; - mm128_t *p = &b->a.a[j-1]; - itr = kh_put(idx, h, p->x>>8>>mi->b<<1, &absent); + mm128_t *p = &b->a.a[j - 1]; + itr = kh_put(idx, h, p->x >> 8 >> mi->b << 1, &absent); assert(absent && j == start_a + n); - if (n == 1) { + if (n == 1) + { kh_key(h, itr) |= 1; kh_val(h, itr) = p->y; - } else { + } + else + { int k; for (k = 0; k < n; ++k) b->p[start_p + k] = b->a.a[start_a + k].y; radix_sort_64(&b->p[start_p], &b->p[start_p + n]); // sort by position; needed as in-place radix_sort_128x() is not stable - kh_val(h, itr) = (uint64_t)start_p<<32 | n; + kh_val(h, itr) = (uint64_t)start_p << 32 | n; start_p += n; } start_a = j, n = 1; - } else ++n; + } + else + ++n; } b->h = h; assert(b->n == (int32_t)start_p); @@ -268,10 +325,10 @@ static void worker_post(void *g, long i, int tid) kfree(0, b->a.a); b->a.n = b->a.m = 0, b->a.a = 0; } - + static void mm_idx_post(mm_idx_t *mi, int n_threads) { - kt_for(n_threads, worker_post, mi, 1<b); + kt_for(n_threads, worker_post, mi, 1 << mi->b); } /****************** @@ -282,24 +339,27 @@ static void mm_idx_post(mm_idx_t *mi, int n_threads) #include #include "bseq.h" -typedef struct { +typedef struct +{ int mini_batch_size; uint64_t batch_size, sum_len; mm_bseq_file_t *fp; mm_idx_t *mi; } pipeline_t; -typedef struct { - int n_seq; +typedef struct +{ + int n_seq; mm_bseq1_t *seq; mm128_v a; } step_t; static void mm_idx_add(mm_idx_t *mi, int n, const mm128_t *a) { - int i, mask = (1<b) - 1; - for (i = 0; i < n; ++i) { - mm128_v *p = &mi->B[a[i].x>>8&mask].a; + int i, mask = (1 << mi->b) - 1; + for (i = 0; i < n; ++i) + { + mm128_v *p = &mi->B[a[i].x >> 8 & mask].a; kv_push(mm128_t, 0, *p, a[i]); } } @@ -307,46 +367,60 @@ static void mm_idx_add(mm_idx_t *mi, int n, const mm128_t *a) static void *worker_pipeline(void *shared, int step, void *in) { int i; - pipeline_t *p = (pipeline_t*)shared; - if (step == 0) { // step 0: read sequences - step_t *s; - if (p->sum_len > p->batch_size) return 0; - s = (step_t*)calloc(1, sizeof(step_t)); + pipeline_t *p = (pipeline_t *)shared; + if (step == 0) + { // step 0: read sequences + step_t *s; + if (p->sum_len > p->batch_size) + return 0; + s = (step_t *)calloc(1, sizeof(step_t)); s->seq = mm_bseq_read(p->fp, p->mini_batch_size, 0, &s->n_seq); // read a mini-batch - if (s->seq) { + if (s->seq) + { uint32_t old_m, m; assert((uint64_t)p->mi->n_seq + s->n_seq <= UINT32_MAX); // to prevent integer overflow // make room for p->mi->seq old_m = p->mi->n_seq, m = p->mi->n_seq + s->n_seq; - kroundup32(m); kroundup32(old_m); + kroundup32(m); + kroundup32(old_m); if (old_m != m) - p->mi->seq = (mm_idx_seq_t*)krealloc(p->mi->km, p->mi->seq, m * sizeof(mm_idx_seq_t)); + p->mi->seq = (mm_idx_seq_t *)krealloc(p->mi->km, p->mi->seq, m * sizeof(mm_idx_seq_t)); // make room for p->mi->S - if (!(p->mi->flag & MM_I_NO_SEQ)) { + if (!(p->mi->flag & MM_I_NO_SEQ)) + { uint64_t sum_len, old_max_len, max_len; - for (i = 0, sum_len = 0; i < s->n_seq; ++i) sum_len += s->seq[i].l_seq; + for (i = 0, sum_len = 0; i < s->n_seq; ++i) + sum_len += s->seq[i].l_seq; old_max_len = (p->sum_len + 7) / 8; max_len = (p->sum_len + sum_len + 7) / 8; - kroundup64(old_max_len); kroundup64(max_len); - if (old_max_len != max_len) { - p->mi->S = (uint32_t*)realloc(p->mi->S, max_len * 4); + kroundup64(old_max_len); + kroundup64(max_len); + if (old_max_len != max_len) + { + p->mi->S = (uint32_t *)realloc(p->mi->S, max_len * 4); memset(&p->mi->S[old_max_len], 0, 4 * (max_len - old_max_len)); } } // populate p->mi->seq - for (i = 0; i < s->n_seq; ++i) { + for (i = 0; i < s->n_seq; ++i) + { mm_idx_seq_t *seq = &p->mi->seq[p->mi->n_seq]; uint32_t j; - if (!(p->mi->flag & MM_I_NO_NAME)) { - seq->name = (char*)kmalloc(p->mi->km, strlen(s->seq[i].name) + 1); + if (!(p->mi->flag & MM_I_NO_NAME)) + { + seq->name = (char *)kmalloc(p->mi->km, strlen(s->seq[i].name) + 1); strcpy(seq->name, s->seq[i].name); - } else seq->name = 0; + } + else + seq->name = 0; seq->len = s->seq[i].l_seq; seq->offset = p->sum_len; seq->is_alt = 0; // copy the sequence - if (!(p->mi->flag & MM_I_NO_SEQ)) { - for (j = 0; j < seq->len; ++j) { // TODO: this is not the fastest way, but let's first see if speed matters here + if (!(p->mi->flag & MM_I_NO_SEQ)) + { + for (j = 0; j < seq->len; ++j) + { // TODO: this is not the fastest way, but let's first see if speed matters here uint64_t o = p->sum_len + j; int c = seq_nt4_table[(uint8_t)s->seq[i].seq[j]]; mm_seq4_set(p->mi->S, o, c); @@ -357,38 +431,50 @@ static void *worker_pipeline(void *shared, int step, void *in) s->seq[i].rid = p->mi->n_seq++; } return s; - } else free(s); - } else if (step == 1) { // step 1: compute sketch - step_t *s = (step_t*)in; - for (i = 0; i < s->n_seq; ++i) { + } + else + free(s); + } + else if (step == 1) + { // step 1: compute sketch + step_t *s = (step_t *)in; + for (i = 0; i < s->n_seq; ++i) + { mm_bseq1_t *t = &s->seq[i]; if (t->l_seq > 0) - mm_sketch(0, t->seq, t->l_seq, p->mi->w, p->mi->k, t->rid, p->mi->flag&MM_I_HPC, &s->a); + mm_sketch(0, t->seq, t->l_seq, p->mi->w, p->mi->k, t->rid, p->mi->flag & MM_I_HPC, &s->a); else if (mm_verbose >= 2) fprintf(stderr, "[WARNING] the length database sequence '%s' is 0\n", t->name); - free(t->seq); free(t->name); + free(t->seq); + free(t->name); } - free(s->seq); s->seq = 0; + free(s->seq); + s->seq = 0; return s; - } else if (step == 2) { // dispatch sketch to buckets - step_t *s = (step_t*)in; - mm_idx_add(p->mi, s->a.n, s->a.a); - kfree(0, s->a.a); free(s); } - return 0; + else if (step == 2) + { // dispatch sketch to buckets + step_t *s = (step_t *)in; + mm_idx_add(p->mi, s->a.n, s->a.a); + kfree(0, s->a.a); + free(s); + } + return 0; } +// 生成index mm_idx_t *mm_idx_gen(mm_bseq_file_t *fp, int w, int k, int b, int flag, int mini_batch_size, int n_threads, uint64_t batch_size) { pipeline_t pl; - if (fp == 0 || mm_bseq_eof(fp)) return 0; + if (fp == 0 || mm_bseq_eof(fp)) + return 0; memset(&pl, 0, sizeof(pipeline_t)); - pl.mini_batch_size = (uint64_t)mini_batch_size < batch_size? mini_batch_size : batch_size; + pl.mini_batch_size = (uint64_t)mini_batch_size < batch_size ? mini_batch_size : batch_size; pl.batch_size = batch_size; pl.fp = fp; pl.mi = mm_idx_init(w, k, b, flag); - kt_pipeline(n_threads < 3? n_threads : 3, worker_pipeline, &pl, 3); + kt_pipeline(n_threads < 3 ? n_threads : 3, worker_pipeline, &pl, 3); if (mm_verbose >= 3) fprintf(stderr, "[M::%s::%.3f*%.2f] collected minimizers\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0)); @@ -404,8 +490,9 @@ mm_idx_t *mm_idx_build(const char *fn, int w, int k, int flag, int n_threads) // mm_bseq_file_t *fp; mm_idx_t *mi; fp = mm_bseq_open(fn); - if (fp == 0) return 0; - mi = mm_idx_gen(fp, w, k, 14, flag, 1<<18, n_threads, UINT64_MAX); + if (fp == 0) + return 0; + mi = mm_idx_gen(fp, w, k, 14, flag, 1 << 18, n_threads, UINT64_MAX); mm_bseq_close(fp); return mi; } @@ -413,29 +500,35 @@ mm_idx_t *mm_idx_build(const char *fn, int w, int k, int flag, int n_threads) // mm_idx_t *mm_idx_str(int w, int k, int is_hpc, int bucket_bits, int n, const char **seq, const char **name) { uint64_t sum_len = 0; - mm128_v a = {0,0,0}; + mm128_v a = {0, 0, 0}; mm_idx_t *mi; - khash_t(str) *h; + khash_t(str) * h; int i, flag = 0; - if (n <= 0) return 0; + if (n <= 0) + return 0; for (i = 0; i < n; ++i) // get the total length sum_len += strlen(seq[i]); - if (is_hpc) flag |= MM_I_HPC; - if (name == 0) flag |= MM_I_NO_NAME; - if (bucket_bits < 0) bucket_bits = 14; + if (is_hpc) + flag |= MM_I_HPC; + if (name == 0) + flag |= MM_I_NO_NAME; + if (bucket_bits < 0) + bucket_bits = 14; mi = mm_idx_init(w, k, bucket_bits, flag); mi->n_seq = n; - mi->seq = (mm_idx_seq_t*)kcalloc(mi->km, n, sizeof(mm_idx_seq_t)); // ->seq is allocated from km - mi->S = (uint32_t*)calloc((sum_len + 7) / 8, 4); + mi->seq = (mm_idx_seq_t *)kcalloc(mi->km, n, sizeof(mm_idx_seq_t)); // ->seq is allocated from km + mi->S = (uint32_t *)calloc((sum_len + 7) / 8, 4); mi->h = h = kh_init(str); - for (i = 0, sum_len = 0; i < n; ++i) { + for (i = 0, sum_len = 0; i < n; ++i) + { const char *s = seq[i]; mm_idx_seq_t *p = &mi->seq[i]; uint32_t j; - if (name && name[i]) { + if (name && name[i]) + { int absent; - p->name = (char*)kmalloc(mi->km, strlen(name[i]) + 1); + p->name = (char *)kmalloc(mi->km, strlen(name[i]) + 1); strcpy(p->name, name[i]); kh_put(str, h, p->name, &absent); assert(absent); @@ -443,13 +536,15 @@ mm_idx_t *mm_idx_str(int w, int k, int is_hpc, int bucket_bits, int n, const cha p->offset = sum_len; p->len = strlen(s); p->is_alt = 0; - for (j = 0; j < p->len; ++j) { + for (j = 0; j < p->len; ++j) + { int c = seq_nt4_table[(uint8_t)s[j]]; uint64_t o = sum_len + j; mm_seq4_set(mi->S, o, c); } sum_len += p->len; - if (p->len > 0) { + if (p->len > 0) + { a.n = 0; mm_sketch(0, s, p->len, w, k, i, is_hpc, &a); mm_idx_add(mi, a.n, a.a); @@ -472,30 +567,38 @@ void mm_idx_dump(FILE *fp, const mm_idx_t *mi) x[0] = mi->w, x[1] = mi->k, x[2] = mi->b, x[3] = mi->n_seq, x[4] = mi->flag; fwrite(MM_IDX_MAGIC, 1, 4, fp); fwrite(x, 4, 5, fp); - for (i = 0; i < mi->n_seq; ++i) { - if (mi->seq[i].name) { + for (i = 0; i < mi->n_seq; ++i) + { + if (mi->seq[i].name) + { uint8_t l = strlen(mi->seq[i].name); fwrite(&l, 1, 1, fp); fwrite(mi->seq[i].name, 1, l, fp); - } else { + } + else + { uint8_t l = 0; fwrite(&l, 1, 1, fp); } fwrite(&mi->seq[i].len, 4, 1, fp); sum_len += mi->seq[i].len; } - for (i = 0; i < 1<b; ++i) { + for (i = 0; i < 1 << mi->b; ++i) + { mm_idx_bucket_t *b = &mi->B[i]; khint_t k; - idxhash_t *h = (idxhash_t*)b->h; - uint32_t size = h? h->size : 0; + idxhash_t *h = (idxhash_t *)b->h; + uint32_t size = h ? h->size : 0; fwrite(&b->n, 4, 1, fp); fwrite(b->p, 8, b->n, fp); fwrite(&size, 4, 1, fp); - if (size == 0) continue; - for (k = 0; k < kh_end(h); ++k) { + if (size == 0) + continue; + for (k = 0; k < kh_end(h); ++k) + { uint64_t x[2]; - if (!kh_exist(h, k)) continue; + if (!kh_exist(h, k)) + continue; x[0] = kh_key(h, k), x[1] = kh_val(h, k); fwrite(x, 8, 2, fp); } @@ -512,18 +615,23 @@ mm_idx_t *mm_idx_load(FILE *fp) uint64_t sum_len = 0; mm_idx_t *mi; - if (fread(magic, 1, 4, fp) != 4) return 0; - if (strncmp(magic, MM_IDX_MAGIC, 4) != 0) return 0; - if (fread(x, 4, 5, fp) != 5) return 0; + if (fread(magic, 1, 4, fp) != 4) + return 0; + if (strncmp(magic, MM_IDX_MAGIC, 4) != 0) + return 0; + if (fread(x, 4, 5, fp) != 5) + return 0; mi = mm_idx_init(x[0], x[1], x[2], x[4]); mi->n_seq = x[3]; - mi->seq = (mm_idx_seq_t*)kcalloc(mi->km, mi->n_seq, sizeof(mm_idx_seq_t)); - for (i = 0; i < mi->n_seq; ++i) { + mi->seq = (mm_idx_seq_t *)kcalloc(mi->km, mi->n_seq, sizeof(mm_idx_seq_t)); + for (i = 0; i < mi->n_seq; ++i) + { uint8_t l; mm_idx_seq_t *s = &mi->seq[i]; fread(&l, 1, 1, fp); - if (l) { - s->name = (char*)kmalloc(mi->km, l + 1); + if (l) + { + s->name = (char *)kmalloc(mi->km, l + 1); fread(s->name, 1, l, fp); s->name[l] = 0; } @@ -532,19 +640,22 @@ mm_idx_t *mm_idx_load(FILE *fp) s->is_alt = 0; sum_len += s->len; } - for (i = 0; i < 1<b; ++i) { + for (i = 0; i < 1 << mi->b; ++i) + { mm_idx_bucket_t *b = &mi->B[i]; uint32_t j, size; khint_t k; idxhash_t *h; fread(&b->n, 4, 1, fp); - b->p = (uint64_t*)malloc(b->n * 8); + b->p = (uint64_t *)malloc(b->n * 8); fread(b->p, 8, b->n, fp); fread(&size, 4, 1, fp); - if (size == 0) continue; + if (size == 0) + continue; b->h = h = kh_init(idx); kh_resize(idx, h, size); - for (j = 0; j < size; ++j) { + for (j = 0; j < size; ++j) + { uint64_t x[2]; int absent; fread(x, 8, 2, fp); @@ -553,8 +664,9 @@ mm_idx_t *mm_idx_load(FILE *fp) kh_val(h, k) = x[1]; } } - if (!(mi->flag & MM_I_NO_SEQ)) { - mi->S = (uint32_t*)malloc((sum_len + 7) / 8 * 4); + if (!(mi->flag & MM_I_NO_SEQ)) + { + mi->S = (uint32_t *)malloc((sum_len + 7) / 8 * 4); fread(mi->S, 4, (sum_len + 7) / 8, fp); } return mi; @@ -566,14 +678,18 @@ int64_t mm_idx_is_idx(const char *fn) int64_t ret, off_end; char magic[4]; - if (strcmp(fn, "-") == 0) return 0; // read from pipe; not an index + if (strcmp(fn, "-") == 0) + return 0; // read from pipe; not an index fd = open(fn, O_RDONLY); - if (fd < 0) return -1; // error + if (fd < 0) + return -1; // error #ifdef WIN32 - if ((off_end = _lseeki64(fd, 0, SEEK_END)) >= 4) { + if ((off_end = _lseeki64(fd, 0, SEEK_END)) >= 4) + { _lseeki64(fd, 0, SEEK_SET); #else - if ((off_end = lseek(fd, 0, SEEK_END)) >= 4) { + if ((off_end = lseek(fd, 0, SEEK_END)) >= 4) + { lseek(fd, 0, SEEK_SET); #endif // WIN32 ret = read(fd, magic, 4); @@ -581,7 +697,7 @@ int64_t mm_idx_is_idx(const char *fn) is_idx = 1; } close(fd); - return is_idx? off_end : 0; + return is_idx ? off_end : 0; } mm_idx_reader_t *mm_idx_reader_open(const char *fn, const mm_idxopt_t *opt, const char *fn_out) @@ -589,44 +705,58 @@ mm_idx_reader_t *mm_idx_reader_open(const char *fn, const mm_idxopt_t *opt, cons int64_t is_idx; mm_idx_reader_t *r; is_idx = mm_idx_is_idx(fn); - if (is_idx < 0) return 0; // failed to open the index - r = (mm_idx_reader_t*)calloc(1, sizeof(mm_idx_reader_t)); + if (is_idx < 0) + return 0; // failed to open the index + r = (mm_idx_reader_t *)calloc(1, sizeof(mm_idx_reader_t)); r->is_idx = is_idx; - if (opt) r->opt = *opt; - else mm_idxopt_init(&r->opt); - if (r->is_idx) { + if (opt) + r->opt = *opt; + else + mm_idxopt_init(&r->opt); + if (r->is_idx) + { r->fp.idx = fopen(fn, "rb"); r->idx_size = is_idx; - } else r->fp.seq = mm_bseq_open(fn); - if (fn_out) r->fp_out = fopen(fn_out, "wb"); + } + else + r->fp.seq = mm_bseq_open(fn); + if (fn_out) + r->fp_out = fopen(fn_out, "wb"); return r; } void mm_idx_reader_close(mm_idx_reader_t *r) { - if (r->is_idx) fclose(r->fp.idx); - else mm_bseq_close(r->fp.seq); - if (r->fp_out) fclose(r->fp_out); + if (r->is_idx) + fclose(r->fp.idx); + else + mm_bseq_close(r->fp.seq); + if (r->fp_out) + fclose(r->fp_out); free(r); } mm_idx_t *mm_idx_reader_read(mm_idx_reader_t *r, int n_threads) { mm_idx_t *mi; -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF int64_t tmp_cur_time = get_mseconds(); #endif - if (r->is_idx) { + if (r->is_idx) + { mi = mm_idx_load(r->fp.idx); - if (mi && mm_verbose >= 2 && (mi->k != r->opt.k || mi->w != r->opt.w || (mi->flag&MM_I_HPC) != (r->opt.flag&MM_I_HPC))) + if (mi && mm_verbose >= 2 && (mi->k != r->opt.k || mi->w != r->opt.w || (mi->flag & MM_I_HPC) != (r->opt.flag & MM_I_HPC))) fprintf(stderr, "[WARNING]\033[1;31m Indexing parameters (-k, -w or -H) overridden by parameters used in the prebuilt index.\033[0m\n"); - } else + } + else mi = mm_idx_gen(r->fp.seq, r->opt.w, r->opt.k, r->opt.bucket_bits, r->opt.flag, r->opt.mini_batch_size, n_threads, r->opt.batch_size); - if (mi) { - if (r->fp_out) mm_idx_dump(r->fp_out, mi); + if (mi) + { + if (r->fp_out) + mm_idx_dump(r->fp_out, mi); mi->index = r->n_parts++; } -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF time_mm_idx_reader_read += get_mseconds() - tmp_cur_time; #endif return mi; @@ -634,7 +764,7 @@ mm_idx_t *mm_idx_reader_read(mm_idx_reader_t *r, int n_threads) int mm_idx_reader_eof(const mm_idx_reader_t *r) // TODO: in extremely rare cases, mm_bseq_eof() might not work { - return r->is_idx? (feof(r->fp.idx) || ftell(r->fp.idx) == r->idx_size) : mm_bseq_eof(r->fp.seq); + return r->is_idx ? (feof(r->fp.idx) || ftell(r->fp.idx) == r->idx_size) : mm_bseq_eof(r->fp.seq); } #include @@ -648,18 +778,24 @@ int mm_idx_alt_read(mm_idx_t *mi, const char *fn) int n_alt = 0; gzFile fp; kstream_t *ks; - kstring_t str = {0,0,0}; - fp = fn && strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r"); - if (fp == 0) return -1; + kstring_t str = {0, 0, 0}; + fp = fn && strcmp(fn, "-") ? gzopen(fn, "r") : gzdopen(fileno(stdin), "r"); + if (fp == 0) + return -1; ks = ks_init(fp); - if (mi->h == 0) mm_idx_index_name(mi); - while (ks_getuntil(ks, KS_SEP_LINE, &str, 0) >= 0) { + if (mi->h == 0) + mm_idx_index_name(mi); + while (ks_getuntil(ks, KS_SEP_LINE, &str, 0) >= 0) + { char *p; int id; - for (p = str.s; *p && !isspace(*p); ++p) { } + for (p = str.s; *p && !isspace(*p); ++p) + { + } *p = 0; id = mm_idx_name2id(mi, str.s); - if (id >= 0) mi->seq[id].is_alt = 1, ++n_alt; + if (id >= 0) + mi->seq[id].is_alt = 1, ++n_alt; } mi->n_alt = n_alt; if (mm_verbose >= 3) @@ -674,71 +810,108 @@ mm_idx_intv_t *mm_idx_read_bed(const mm_idx_t *mi, const char *fn, int read_junc { gzFile fp; kstream_t *ks; - kstring_t str = {0,0,0}; + kstring_t str = {0, 0, 0}; mm_idx_intv_t *I; - fp = fn && strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r"); - if (fp == 0) return 0; - I = (mm_idx_intv_t*)calloc(mi->n_seq, sizeof(*I)); + fp = fn && strcmp(fn, "-") ? gzopen(fn, "r") : gzdopen(fileno(stdin), "r"); + if (fp == 0) + return 0; + I = (mm_idx_intv_t *)calloc(mi->n_seq, sizeof(*I)); ks = ks_init(fp); - while (ks_getuntil(ks, KS_SEP_LINE, &str, 0) >= 0) { + while (ks_getuntil(ks, KS_SEP_LINE, &str, 0) >= 0) + { mm_idx_intv_t *r; - mm_idx_intv1_t t = {-1,-1,-1,-1,0}; + mm_idx_intv1_t t = {-1, -1, -1, -1, 0}; char *p, *q, *bl, *bs; int32_t i, id = -1, n_blk = 0; - for (p = q = str.s, i = 0;; ++p) { - if (*p == 0 || *p == '\t') { + for (p = q = str.s, i = 0;; ++p) + { + if (*p == 0 || *p == '\t') + { int32_t c = *p; *p = 0; - if (i == 0) { // chr + if (i == 0) + { // chr id = mm_idx_name2id(mi, q); - if (id < 0) break; // unknown name; TODO: throw a warning - } else if (i == 1) { // start + if (id < 0) + break; // unknown name; TODO: throw a warning + } + else if (i == 1) + { // start t.st = atol(q); // TODO: watch out integer overflow! - if (t.st < 0) break; - } else if (i == 2) { // end + if (t.st < 0) + break; + } + else if (i == 2) + { // end t.en = atol(q); - if (t.en < 0) break; - } else if (i == 4) { // BED score + if (t.en < 0) + break; + } + else if (i == 4) + { // BED score t.score = atol(q); - } else if (i == 5) { // strand - t.strand = *q == '+'? 1 : *q == '-'? -1 : 0; - } else if (i == 9) { - if (!isdigit(*q)) break; + } + else if (i == 5) + { // strand + t.strand = *q == '+' ? 1 : *q == '-' ? -1 + : 0; + } + else if (i == 9) + { + if (!isdigit(*q)) + break; n_blk = atol(q); - } else if (i == 10) { + } + else if (i == 10) + { bl = q; - } else if (i == 11) { + } + else if (i == 11) + { bs = q; break; } - if (c == 0) break; + if (c == 0) + break; ++i, q = p + 1; } } - if (id < 0 || t.st < 0 || t.st >= t.en) continue; + if (id < 0 || t.st < 0 || t.st >= t.en) + continue; r = &I[id]; - if (i >= 11 && read_junc) { // BED12 + if (i >= 11 && read_junc) + { // BED12 int32_t st, sz, en; - st = strtol(bs, &bs, 10); ++bs; - sz = strtol(bl, &bl, 10); ++bl; + st = strtol(bs, &bs, 10); + ++bs; + sz = strtol(bl, &bl, 10); + ++bl; en = t.st + st + sz; - for (i = 1; i < n_blk; ++i) { + for (i = 1; i < n_blk; ++i) + { mm_idx_intv1_t s = t; - if (r->n == r->m) { - r->m = r->m? r->m + (r->m>>1) : 16; - r->a = (mm_idx_intv1_t*)realloc(r->a, sizeof(*r->a) * r->m); + if (r->n == r->m) + { + r->m = r->m ? r->m + (r->m >> 1) : 16; + r->a = (mm_idx_intv1_t *)realloc(r->a, sizeof(*r->a) * r->m); } - st = strtol(bs, &bs, 10); ++bs; - sz = strtol(bl, &bl, 10); ++bl; + st = strtol(bs, &bs, 10); + ++bs; + sz = strtol(bl, &bl, 10); + ++bl; s.st = en, s.en = t.st + st; en = t.st + st + sz; - if (s.en > s.st) r->a[r->n++] = s; + if (s.en > s.st) + r->a[r->n++] = s; } - } else { - if (r->n == r->m) { - r->m = r->m? r->m + (r->m>>1) : 16; - r->a = (mm_idx_intv1_t*)realloc(r->a, sizeof(*r->a) * r->m); + } + else + { + if (r->n == r->m) + { + r->m = r->m ? r->m + (r->m >> 1) : 16; + r->a = (mm_idx_intv1_t *)realloc(r->a, sizeof(*r->a) * r->m); } r->a[r->n++] = t; } @@ -752,9 +925,11 @@ mm_idx_intv_t *mm_idx_read_bed(const mm_idx_t *mi, const char *fn, int read_junc int mm_idx_bed_read(mm_idx_t *mi, const char *fn, int read_junc) { int32_t i; - if (mi->h == 0) mm_idx_index_name(mi); + if (mi->h == 0) + mm_idx_index_name(mi); mi->I = mm_idx_read_bed(mi, fn, read_junc); - if (mi->I == 0) return -1; + if (mi->I == 0) + return -1; for (i = 0; i < mi->n_seq; ++i) // TODO: eliminate redundant intervals radix_sort_bed(mi->I[i].a, mi->I[i].a + mi->I[i].n); return 0; @@ -765,19 +940,28 @@ int mm_idx_bed_junc(const mm_idx_t *mi, int32_t ctg, int32_t st, int32_t en, uin int32_t i, left, right; mm_idx_intv_t *r; memset(s, 0, en - st); - if (mi->I == 0 || ctg < 0 || ctg >= mi->n_seq) return -1; + if (mi->I == 0 || ctg < 0 || ctg >= mi->n_seq) + return -1; r = &mi->I[ctg]; left = 0, right = r->n; - while (right > left) { + while (right > left) + { int32_t mid = left + ((right - left) >> 1); - if (r->a[mid].st >= st) right = mid; - else left = mid + 1; + if (r->a[mid].st >= st) + right = mid; + else + left = mid + 1; } - for (i = left; i < r->n; ++i) { - if (st <= r->a[i].st && en >= r->a[i].en && r->a[i].strand != 0) { - if (r->a[i].strand > 0) { + for (i = left; i < r->n; ++i) + { + if (st <= r->a[i].st && en >= r->a[i].en && r->a[i].strand != 0) + { + if (r->a[i].strand > 0) + { s[r->a[i].st - st] |= 1, s[r->a[i].en - 1 - st] |= 2; - } else { + } + else + { s[r->a[i].st - st] |= 8, s[r->a[i].en - 1 - st] |= 4; } } diff --git a/lchain.c b/lchain.c index ee72d68..5498649 100644 --- a/lchain.c +++ b/lchain.c @@ -5,7 +5,7 @@ #include "mmpriv.h" #include "kalloc.h" #include "krmq.h" -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF extern int64_t get_mseconds(); extern int64_t time_mg_lchain_dp, time_mg_chain_backtrack; @@ -155,7 +155,7 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int int32_t *f, *t, *v, n_u, n_v, mmax_f = 0, max_drop = bw; int64_t *p, i, j, max_ii, st = 0, n_iter = 0; uint64_t *u; -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF int64_t tmp_cur_time = get_mseconds(), tmp_diff = 0; #endif if (_u) *_u = 0, *n_u_ = 0; @@ -211,17 +211,17 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_ii = i; if (mmax_f < max_f) mmax_f = max_f; } -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF int64_t tmp_inner_time = get_mseconds(); #endif u = mg_chain_backtrack(km, n, f, p, v, t, min_cnt, min_sc, max_drop, &n_u, &n_v); -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF tmp_diff = get_mseconds() - tmp_inner_time; __sync_fetch_and_add(&time_mg_chain_backtrack, tmp_diff); #endif *n_u_ = n_u, *_u = u; // NB: note that u[] may not be sorted by score here kfree(km, p); kfree(km, f); kfree(km, t); -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF tmp_diff = get_mseconds() - tmp_cur_time; __sync_fetch_and_add(&time_mg_lchain_dp, tmp_diff); #endif diff --git a/main.c b/main.c index 6f899b6..76c75d2 100644 --- a/main.c +++ b/main.c @@ -18,7 +18,7 @@ int64_t get_mseconds() } // 记录运行时间的变量 -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF int64_t time_mm_idx_reader_read, time_mm_map_file_frag, @@ -34,7 +34,8 @@ int64_t time_mm_idx_reader_read, time_mg_lchain_dp = 0, time_collect_seed_hits_heap = 0, time_collect_seed_hits = 0, - time_mg_chain_backtrack = 0; + time_mg_chain_backtrack = 0, + time_ksw_extd2_sse = 0; #endif ////////////////////////////////// @@ -181,7 +182,7 @@ int main(int argc, char *argv[]) mm_realtime0 = realtime(); mm_set_opt(0, &ipt, &opt); -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF time_mm_idx_reader_read = 0; time_mm_map_file_frag = 0; @@ -699,7 +700,7 @@ int main(int argc, char *argv[]) fprintf(stderr, " %s", argv[i]); fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB\n", __func__, realtime() - mm_realtime0, cputime(), peakrss() / 1024.0 / 1024.0 / 1024.0); } -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF fprintf(stderr, "\n"); fprintf(stderr, "time_mm_idx_reader_read: %f s\n", time_mm_idx_reader_read / 1000.0); @@ -717,7 +718,8 @@ int main(int argc, char *argv[]) fprintf(stderr, "time_collect_seed_hits: %f s\n", time_collect_seed_hits / 1000.0 / n_threads); fprintf(stderr, "time_mg_lchain_dp: %f s\n", time_mg_lchain_dp / 1000.0 / n_threads); fprintf(stderr, "time_mg_chain_backtrack: %f s\n", time_mg_chain_backtrack / 1000.0 / n_threads); - + fprintf(stderr, "time_ksw_extd2_sse: %f s\n", time_ksw_extd2_sse / 1000.0 / n_threads); + fprintf(stderr, "\n"); #endif return 0; diff --git a/map.c b/map.c index e68567b..5809a5c 100644 --- a/map.c +++ b/map.c @@ -10,7 +10,7 @@ #include "bseq.h" #include "khash.h" -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF extern int64_t get_mseconds(); extern int64_t time_mm_map_file_frag, time_map_work_for_block_1, @@ -145,7 +145,7 @@ static mm128_t *collect_seed_hits_heap(void *km, const mm_mapopt_t *opt, int max int64_t j, n_for = 0, n_rev = 0; mm_seed_t *m; mm128_t *a, *heap; -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF int64_t tmp_cur_time = get_mseconds(), tmp_diff = 0; #endif m = mm_collect_matches(km, &n_m, qlen, max_occ, opt->max_max_occ, opt->occ_dist, mi, mv, n_a, rep_len, n_mini_pos, mini_pos); @@ -217,7 +217,7 @@ static mm128_t *collect_seed_hits_heap(void *km, const mm_mapopt_t *opt, int max memmove(a + n_for, a + (*n_a) - n_rev, n_rev * sizeof(mm128_t)); *n_a = n_for + n_rev; } -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF tmp_diff = get_mseconds() - tmp_cur_time; __sync_fetch_and_add(&time_collect_seed_hits_heap, tmp_diff); #endif @@ -230,14 +230,10 @@ static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ, int i, n_m; mm_seed_t *m; mm128_t *a; -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF int64_t tmp_cur_time = get_mseconds(), tmp_diff = 0; #endif m = mm_collect_matches(km, &n_m, qlen, max_occ, opt->max_max_occ, opt->occ_dist, mi, mv, n_a, rep_len, n_mini_pos, mini_pos); -#ifdef ANALYSIS_PERF - tmp_diff = get_mseconds() - tmp_cur_time; - __sync_fetch_and_add(&time_collect_seed_hits, tmp_diff); -#endif a = (mm128_t *)kmalloc(km, *n_a * sizeof(mm128_t)); for (i = 0, *n_a = 0; i < n_m; ++i) { @@ -276,6 +272,10 @@ static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ, } kfree(km, m); radix_sort_128x(a, a + (*n_a)); +#ifdef SHOW_PERF + tmp_diff = get_mseconds() - tmp_cur_time; + __sync_fetch_and_add(&time_collect_seed_hits, tmp_diff); +#endif return a; } @@ -317,7 +317,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** mm_reg1_t *regs0; km_stat_t kmst; float chn_pen_gap, chn_pen_skip; -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF int64_t tmp_cur_time = get_mseconds(), tmp_diff = 0; #endif for (i = 0, qlen_sum = 0; i < n_segs; ++i) @@ -333,13 +333,13 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** hash = __ac_Wang_hash(hash); collect_minimizers(b->km, opt, mi, n_segs, qlens, seqs, &mv); -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF tmp_diff = get_mseconds() - tmp_cur_time; __sync_fetch_and_add(&time_mm_map_frag_b1, tmp_diff); tmp_cur_time = get_mseconds(); #endif if (opt->q_occ_frac > 0.0f) - mm_seed_mz_flt(b->km, &mv, opt->mid_occ, opt->q_occ_frac); + mm_seed_mz_flt(b->km, &mv, opt->mid_occ, opt->q_occ_frac); // 过滤掉出现次数太多的minimizer if (opt->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos); else @@ -352,7 +352,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** fprintf(stderr, "SD\t%s\t%d\t%c\t%d\t%d\t%d\n", mi->seq[a[i].x << 1 >> 33].name, (int32_t)a[i].x, "+-"[a[i].x >> 63], (int32_t)a[i].y, (int32_t)(a[i].y >> 32 & 0xff), i == 0 ? 0 : ((int32_t)a[i].y - (int32_t)a[i - 1].y) - ((int32_t)a[i].x - (int32_t)a[i - 1].x)); } -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF tmp_diff = get_mseconds() - tmp_cur_time; __sync_fetch_and_add(&time_mm_map_frag_b2, tmp_diff); tmp_cur_time = get_mseconds(); @@ -387,7 +387,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** a = mg_lchain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score, chn_pen_gap, chn_pen_skip, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km); } -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF tmp_diff = get_mseconds() - tmp_cur_time; __sync_fetch_and_add(&time_mm_map_frag_b3, tmp_diff); tmp_cur_time = get_mseconds(); @@ -441,7 +441,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** } b->frag_gap = max_chain_gap_ref; b->rep_len = rep_len; -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF tmp_diff = get_mseconds() - tmp_cur_time; __sync_fetch_and_add(&time_mm_map_frag_b4, tmp_diff); tmp_cur_time = get_mseconds(); @@ -465,7 +465,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** mm_est_err(mi, qlen_sum, n_regs0, regs0, a, n_mini_pos, mini_pos); n_regs0 = mm_filter_strand_retained(n_regs0, regs0); } -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF tmp_diff = get_mseconds() - tmp_cur_time; __sync_fetch_and_add(&time_mm_map_frag_b5, tmp_diff); tmp_cur_time = get_mseconds(); @@ -492,7 +492,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** if (n_segs == 2 && opt->pe_ori >= 0 && (opt->flag & MM_F_CIGAR)) mm_pair(b->km, max_chain_gap_ref, opt->pe_bonus, opt->a * 2 + opt->b, opt->a, qlens, n_regs, regs); // pairing } -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF tmp_diff = get_mseconds() - tmp_cur_time; __sync_fetch_and_add(&time_mm_map_frag_b6, tmp_diff); tmp_cur_time = get_mseconds(); @@ -566,7 +566,7 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback fprintf(stderr, "QR\t%s\t%d\t%d\n", s->seq[off].name, tid, s->seq[off].l_seq); t = realtime(); } -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF int64_t tmp_cur_time = get_mseconds(), tmp_diff = 0; #endif for (j = 0; j < s->n_seg[i]; ++j) @@ -576,7 +576,7 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback qlens[j] = s->seq[off + j].l_seq; qseqs[j] = s->seq[off + j].seq; } -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF tmp_diff = get_mseconds() - tmp_cur_time; __sync_fetch_and_add(&time_map_work_for_block_1, tmp_diff); tmp_cur_time = get_mseconds(); @@ -599,7 +599,7 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback s->frag_gap[off + j] = b->frag_gap; } } -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF tmp_diff = get_mseconds() - tmp_cur_time; __sync_fetch_and_add(&time_map_work_for_block_2, tmp_diff); tmp_cur_time = get_mseconds(); @@ -618,7 +618,7 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback r->rev = !r->rev; } } -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF tmp_diff = get_mseconds() - tmp_cur_time; __sync_fetch_and_add(&time_map_work_for_block_3, tmp_diff); #endif @@ -851,7 +851,7 @@ static mm_bseq_file_t **open_bseqs(int n, const char **fn) int mm_map_file_frag(const mm_idx_t *idx, int n_segs, const char **fn, const mm_mapopt_t *opt, int n_threads) { -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF int64_t tmp_cur_time = get_mseconds(); #endif int i, pl_threads; @@ -878,7 +878,7 @@ int mm_map_file_frag(const mm_idx_t *idx, int n_segs, const char **fn, const mm_ for (i = 0; i < pl.n_fp; ++i) mm_bseq_close(pl.fp[i]); free(pl.fp); -#ifdef ANALYSIS_PERF +#ifdef SHOW_PERF time_mm_map_file_frag += get_mseconds() - tmp_cur_time; #endif return 0; diff --git a/minimap.h b/minimap.h index 18f5b8f..2d3fa98 100644 --- a/minimap.h +++ b/minimap.h @@ -8,413 +8,430 @@ #define MM_VERSION "2.26-r1175" // 用来开关调试性能分析，运行时间等信息 -#define ANALYSIS_PERF 1 +#define SHOW_PERF 1 -#define MM_F_NO_DIAG (0x001LL) // no exact diagonal hit -#define MM_F_NO_DUAL (0x002LL) // skip pairs where query name is lexicographically larger than target name -#define MM_F_CIGAR (0x004LL) -#define MM_F_OUT_SAM (0x008LL) -#define MM_F_NO_QUAL (0x010LL) -#define MM_F_OUT_CG (0x020LL) -#define MM_F_OUT_CS (0x040LL) -#define MM_F_SPLICE (0x080LL) // splice mode -#define MM_F_SPLICE_FOR (0x100LL) // match GT-AG -#define MM_F_SPLICE_REV (0x200LL) // match CT-AC, the reverse complement of GT-AG -#define MM_F_NO_LJOIN (0x400LL) -#define MM_F_OUT_CS_LONG (0x800LL) -#define MM_F_SR (0x1000LL) -#define MM_F_FRAG_MODE (0x2000LL) -#define MM_F_NO_PRINT_2ND (0x4000LL) -#define MM_F_2_IO_THREADS (0x8000LL) -#define MM_F_LONG_CIGAR (0x10000LL) -#define MM_F_INDEPEND_SEG (0x20000LL) -#define MM_F_SPLICE_FLANK (0x40000LL) -#define MM_F_SOFTCLIP (0x80000LL) -#define MM_F_FOR_ONLY (0x100000LL) -#define MM_F_REV_ONLY (0x200000LL) -#define MM_F_HEAP_SORT (0x400000LL) -#define MM_F_ALL_CHAINS (0x800000LL) -#define MM_F_OUT_MD (0x1000000LL) -#define MM_F_COPY_COMMENT (0x2000000LL) -#define MM_F_EQX (0x4000000LL) // use =/X instead of M -#define MM_F_PAF_NO_HIT (0x8000000LL) // output unmapped reads to PAF -#define MM_F_NO_END_FLT (0x10000000LL) -#define MM_F_HARD_MLEVEL (0x20000000LL) -#define MM_F_SAM_HIT_ONLY (0x40000000LL) -#define MM_F_RMQ (0x80000000LL) -#define MM_F_QSTRAND (0x100000000LL) -#define MM_F_NO_INV (0x200000000LL) -#define MM_F_NO_HASH_NAME (0x400000000LL) -#define MM_F_SPLICE_OLD (0x800000000LL) -#define MM_F_SECONDARY_SEQ (0x1000000000LL) //output SEQ field for seqondary alignments using hard clipping +#define MM_F_NO_DIAG (0x001LL) // no exact diagonal hit +#define MM_F_NO_DUAL (0x002LL) // skip pairs where query name is lexicographically larger than target name +#define MM_F_CIGAR (0x004LL) +#define MM_F_OUT_SAM (0x008LL) +#define MM_F_NO_QUAL (0x010LL) +#define MM_F_OUT_CG (0x020LL) +#define MM_F_OUT_CS (0x040LL) +#define MM_F_SPLICE (0x080LL) // splice mode +#define MM_F_SPLICE_FOR (0x100LL) // match GT-AG +#define MM_F_SPLICE_REV (0x200LL) // match CT-AC, the reverse complement of GT-AG +#define MM_F_NO_LJOIN (0x400LL) +#define MM_F_OUT_CS_LONG (0x800LL) +#define MM_F_SR (0x1000LL) +#define MM_F_FRAG_MODE (0x2000LL) +#define MM_F_NO_PRINT_2ND (0x4000LL) +#define MM_F_2_IO_THREADS (0x8000LL) +#define MM_F_LONG_CIGAR (0x10000LL) +#define MM_F_INDEPEND_SEG (0x20000LL) +#define MM_F_SPLICE_FLANK (0x40000LL) +#define MM_F_SOFTCLIP (0x80000LL) +#define MM_F_FOR_ONLY (0x100000LL) +#define MM_F_REV_ONLY (0x200000LL) +#define MM_F_HEAP_SORT (0x400000LL) +#define MM_F_ALL_CHAINS (0x800000LL) +#define MM_F_OUT_MD (0x1000000LL) +#define MM_F_COPY_COMMENT (0x2000000LL) +#define MM_F_EQX (0x4000000LL) // use =/X instead of M +#define MM_F_PAF_NO_HIT (0x8000000LL) // output unmapped reads to PAF +#define MM_F_NO_END_FLT (0x10000000LL) +#define MM_F_HARD_MLEVEL (0x20000000LL) +#define MM_F_SAM_HIT_ONLY (0x40000000LL) +#define MM_F_RMQ (0x80000000LL) +#define MM_F_QSTRAND (0x100000000LL) +#define MM_F_NO_INV (0x200000000LL) +#define MM_F_NO_HASH_NAME (0x400000000LL) +#define MM_F_SPLICE_OLD (0x800000000LL) +#define MM_F_SECONDARY_SEQ (0x1000000000LL) // output SEQ field for seqondary alignments using hard clipping -#define MM_I_HPC 0x1 -#define MM_I_NO_SEQ 0x2 -#define MM_I_NO_NAME 0x4 +#define MM_I_HPC 0x1 +#define MM_I_NO_SEQ 0x2 +#define MM_I_NO_NAME 0x4 -#define MM_IDX_MAGIC "MMI\2" +#define MM_IDX_MAGIC "MMI\2" -#define MM_MAX_SEG 255 +#define MM_MAX_SEG 255 -#define MM_CIGAR_MATCH 0 -#define MM_CIGAR_INS 1 -#define MM_CIGAR_DEL 2 -#define MM_CIGAR_N_SKIP 3 -#define MM_CIGAR_SOFTCLIP 4 -#define MM_CIGAR_HARDCLIP 5 -#define MM_CIGAR_PADDING 6 -#define MM_CIGAR_EQ_MATCH 7 +#define MM_CIGAR_MATCH 0 +#define MM_CIGAR_INS 1 +#define MM_CIGAR_DEL 2 +#define MM_CIGAR_N_SKIP 3 +#define MM_CIGAR_SOFTCLIP 4 +#define MM_CIGAR_HARDCLIP 5 +#define MM_CIGAR_PADDING 6 +#define MM_CIGAR_EQ_MATCH 7 #define MM_CIGAR_X_MISMATCH 8 -#define MM_CIGAR_STR "MIDNSHP=XB" +#define MM_CIGAR_STR "MIDNSHP=XB" #ifdef __cplusplus -extern "C" { +extern "C" +{ #endif -// emulate 128-bit integers and arrays -typedef struct { uint64_t x, y; } mm128_t; -typedef struct { size_t n, m; mm128_t *a; } mm128_v; + // emulate 128-bit integers and arrays + typedef struct + { + uint64_t x, y; + } mm128_t; + typedef struct + { + size_t n, m; + mm128_t *a; + } mm128_v; -// minimap2 index -typedef struct { - char *name; // name of the db sequence - uint64_t offset; // offset in mm_idx_t::S - uint32_t len; // length - uint32_t is_alt; -} mm_idx_seq_t; + // minimap2 index + typedef struct + { + char *name; // name of the db sequence + uint64_t offset; // offset in mm_idx_t::S + uint32_t len; // length + uint32_t is_alt; + } mm_idx_seq_t; -typedef struct { - int32_t b, w, k, flag; - uint32_t n_seq; // number of reference sequences - int32_t index; - int32_t n_alt; - mm_idx_seq_t *seq; // sequence name, length and offset - uint32_t *S; // 4-bit packed sequence - struct mm_idx_bucket_s *B; // index (hidden) - struct mm_idx_intv_s *I; // intervals (hidden) - void *km, *h; -} mm_idx_t; + typedef struct + { + int32_t b, w, k, flag; + uint32_t n_seq; // number of reference sequences + int32_t index; + int32_t n_alt; + mm_idx_seq_t *seq; // sequence name, length and offset + uint32_t *S; // 4-bit packed sequence + struct mm_idx_bucket_s *B; // index (hidden) + struct mm_idx_intv_s *I; // intervals (hidden) + void *km, *h; + } mm_idx_t; -// minimap2 alignment -typedef struct { - uint32_t capacity; // the capacity of cigar[] - int32_t dp_score, dp_max, dp_max2; // DP score; score of the max-scoring segment; score of the best alternate mappings - uint32_t n_ambi:30, trans_strand:2; // number of ambiguous bases; transcript strand: 0 for unknown, 1 for +, 2 for - - uint32_t n_cigar; // number of cigar operations in cigar[] - uint32_t cigar[]; -} mm_extra_t; + // minimap2 alignment + typedef struct + { + uint32_t capacity; // the capacity of cigar[] + int32_t dp_score, dp_max, dp_max2; // DP score; score of the max-scoring segment; score of the best alternate mappings + uint32_t n_ambi : 30, trans_strand : 2; // number of ambiguous bases; transcript strand: 0 for unknown, 1 for +, 2 for - + uint32_t n_cigar; // number of cigar operations in cigar[] + uint32_t cigar[]; + } mm_extra_t; -typedef struct { - int32_t id; // ID for internal uses (see also parent below) - int32_t cnt; // number of minimizers; if on the reverse strand - int32_t rid; // reference index; if this is an alignment from inversion rescue - int32_t score; // DP alignment score - int32_t qs, qe, rs, re; // query start and end; reference start and end - int32_t parent, subsc; // parent==id if primary; best alternate mapping score - int32_t as; // offset in the a[] array (for internal uses only) - int32_t mlen, blen; // seeded exact match length; seeded alignment block length - int32_t n_sub; // number of suboptimal mappings - int32_t score0; // initial chaining score (before chain merging/spliting) - uint32_t mapq:8, split:2, rev:1, inv:1, sam_pri:1, proper_frag:1, pe_thru:1, seg_split:1, seg_id:8, split_inv:1, is_alt:1, strand_retained:1, dummy:5; - uint32_t hash; - float div; - mm_extra_t *p; -} mm_reg1_t; + typedef struct + { + int32_t id; // ID for internal uses (see also parent below) + int32_t cnt; // number of minimizers; if on the reverse strand + int32_t rid; // reference index; if this is an alignment from inversion rescue + int32_t score; // DP alignment score + int32_t qs, qe, rs, re; // query start and end; reference start and end + int32_t parent, subsc; // parent==id if primary; best alternate mapping score + int32_t as; // offset in the a[] array (for internal uses only) + int32_t mlen, blen; // seeded exact match length; seeded alignment block length + int32_t n_sub; // number of suboptimal mappings + int32_t score0; // initial chaining score (before chain merging/spliting) + uint32_t mapq : 8, split : 2, rev : 1, inv : 1, sam_pri : 1, proper_frag : 1, pe_thru : 1, seg_split : 1, seg_id : 8, split_inv : 1, is_alt : 1, strand_retained : 1, dummy : 5; + uint32_t hash; + float div; + mm_extra_t *p; + } mm_reg1_t; -// indexing and mapping options -typedef struct { - short k, w, flag, bucket_bits; - int64_t mini_batch_size; - uint64_t batch_size; -} mm_idxopt_t; + // indexing and mapping options + typedef struct + { + short k, w, flag, bucket_bits; + int64_t mini_batch_size; + uint64_t batch_size; + } mm_idxopt_t; -typedef struct { - int64_t flag; // see MM_F_* macros - int seed; - int sdust_thres; // score threshold for SDUST; 0 to disable + typedef struct + { + int64_t flag; // see MM_F_* macros + int seed; + int sdust_thres; // score threshold for SDUST; 0 to disable - int max_qlen; // max query length + int max_qlen; // max query length - int bw, bw_long; // bandwidth - int max_gap, max_gap_ref; // break a chain if there are no minimizers in a max_gap window - int max_frag_len; - int max_chain_skip, max_chain_iter; - int min_cnt; // min number of minimizers on each chain - int min_chain_score; // min chaining score - float chain_gap_scale; - float chain_skip_scale; - int rmq_size_cap, rmq_inner_dist; - int rmq_rescue_size; - float rmq_rescue_ratio; + int bw, bw_long; // bandwidth + int max_gap, max_gap_ref; // break a chain if there are no minimizers in a max_gap window + int max_frag_len; + int max_chain_skip, max_chain_iter; + int min_cnt; // min number of minimizers on each chain + int min_chain_score; // min chaining score + float chain_gap_scale; + float chain_skip_scale; + int rmq_size_cap, rmq_inner_dist; + int rmq_rescue_size; + float rmq_rescue_ratio; - float mask_level; - int mask_len; - float pri_ratio; - int best_n; // top best_n chains are subjected to DP alignment + float mask_level; + int mask_len; + float pri_ratio; + int best_n; // top best_n chains are subjected to DP alignment - float alt_drop; + float alt_drop; - int a, b, q, e, q2, e2; // matching score, mismatch, gap-open and gap-ext penalties - int transition; // transition mismatch score (A:G, C:T) - int sc_ambi; // score when one or both bases are "N" - int noncan; // cost of non-canonical splicing sites - int junc_bonus; - int zdrop, zdrop_inv; // break alignment if alignment score drops too fast along the diagonal - int end_bonus; - int min_dp_max; // drop an alignment if the score of the max scoring segment is below this threshold - int min_ksw_len; - int anchor_ext_len, anchor_ext_shift; - float max_clip_ratio; // drop an alignment if BOTH ends are clipped above this ratio + int a, b, q, e, q2, e2; // matching score, mismatch, gap-open and gap-ext penalties + int transition; // transition mismatch score (A:G, C:T) + int sc_ambi; // score when one or both bases are "N" + int noncan; // cost of non-canonical splicing sites + int junc_bonus; + int zdrop, zdrop_inv; // break alignment if alignment score drops too fast along the diagonal + int end_bonus; + int min_dp_max; // drop an alignment if the score of the max scoring segment is below this threshold + int min_ksw_len; + int anchor_ext_len, anchor_ext_shift; + float max_clip_ratio; // drop an alignment if BOTH ends are clipped above this ratio - int rank_min_len; - float rank_frac; + int rank_min_len; + float rank_frac; - int pe_ori, pe_bonus; + int pe_ori, pe_bonus; - float mid_occ_frac; // only used by mm_mapopt_update(); see below - float q_occ_frac; - int32_t min_mid_occ, max_mid_occ; - int32_t mid_occ; // ignore seeds with occurrences above this threshold - int32_t max_occ, max_max_occ, occ_dist; - int64_t mini_batch_size; // size of a batch of query bases to process in parallel - int64_t max_sw_mat; - int64_t cap_kalloc; + float mid_occ_frac; // only used by mm_mapopt_update(); see below + float q_occ_frac; + int32_t min_mid_occ, max_mid_occ; + int32_t mid_occ; // ignore seeds with occurrences above this threshold + int32_t max_occ, max_max_occ, occ_dist; + int64_t mini_batch_size; // size of a batch of query bases to process in parallel + int64_t max_sw_mat; + int64_t cap_kalloc; - const char *split_prefix; -} mm_mapopt_t; + const char *split_prefix; + } mm_mapopt_t; -// index reader -typedef struct { - int is_idx, n_parts; - int64_t idx_size; - mm_idxopt_t opt; - FILE *fp_out; - union { - struct mm_bseq_file_s *seq; - FILE *idx; - } fp; -} mm_idx_reader_t; + // index reader + typedef struct + { + int is_idx, n_parts; + int64_t idx_size; + mm_idxopt_t opt; + FILE *fp_out; + union + { + struct mm_bseq_file_s *seq; + FILE *idx; + } fp; + } mm_idx_reader_t; -// memory buffer for thread-local storage during mapping -struct mm_tbuf_s { - void *km; - int rep_len, frag_gap; -}; + // memory buffer for thread-local storage during mapping + struct mm_tbuf_s + { + void *km; + int rep_len, frag_gap; + }; -typedef struct mm_tbuf_s mm_tbuf_t; + typedef struct mm_tbuf_s mm_tbuf_t; -// global variables -extern int mm_verbose, mm_dbg_flag; // verbose level: 0 for no info, 1 for error, 2 for warning, 3 for message (default); debugging flag -extern double mm_realtime0; // wall-clock timer + // global variables + extern int mm_verbose, mm_dbg_flag; // verbose level: 0 for no info, 1 for error, 2 for warning, 3 for message (default); debugging flag + extern double mm_realtime0; // wall-clock timer -/** - * Set default or preset parameters - * - * @param preset NULL to set all parameters as default; otherwise apply preset to affected parameters - * @param io pointer to indexing parameters - * @param mo pointer to mapping parameters - * - * @return 0 if success; -1 if _present_ unknown - */ -int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo); -int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo); + /** + * Set default or preset parameters + * + * @param preset NULL to set all parameters as default; otherwise apply preset to affected parameters + * @param io pointer to indexing parameters + * @param mo pointer to mapping parameters + * + * @return 0 if success; -1 if _present_ unknown + */ + int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo); + int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo); -/** - * Update mm_mapopt_t::mid_occ via mm_mapopt_t::mid_occ_frac - * - * If mm_mapopt_t::mid_occ is 0, this function sets it to a number such that no - * more than mm_mapopt_t::mid_occ_frac of minimizers in the index have a higher - * occurrence. - * - * @param opt mapping parameters - * @param mi minimap2 index - */ -void mm_mapopt_update(mm_mapopt_t *opt, const mm_idx_t *mi); + /** + * Update mm_mapopt_t::mid_occ via mm_mapopt_t::mid_occ_frac + * + * If mm_mapopt_t::mid_occ is 0, this function sets it to a number such that no + * more than mm_mapopt_t::mid_occ_frac of minimizers in the index have a higher + * occurrence. + * + * @param opt mapping parameters + * @param mi minimap2 index + */ + void mm_mapopt_update(mm_mapopt_t *opt, const mm_idx_t *mi); -void mm_mapopt_max_intron_len(mm_mapopt_t *opt, int max_intron_len); + void mm_mapopt_max_intron_len(mm_mapopt_t *opt, int max_intron_len); -/** - * Initialize an index reader - * - * @param fn index or fasta/fastq file name (this function tests the file type) - * @param opt indexing parameters - * @param fn_out if not NULL, write built index to this file - * - * @return an index reader on success; NULL if fail to open _fn_ - */ -mm_idx_reader_t *mm_idx_reader_open(const char *fn, const mm_idxopt_t *opt, const char *fn_out); + /** + * Initialize an index reader + * + * @param fn index or fasta/fastq file name (this function tests the file type) + * @param opt indexing parameters + * @param fn_out if not NULL, write built index to this file + * + * @return an index reader on success; NULL if fail to open _fn_ + */ + mm_idx_reader_t *mm_idx_reader_open(const char *fn, const mm_idxopt_t *opt, const char *fn_out); -/** - * Read/build an index - * - * If the input file is an index file, this function reads one part of the - * index and returns. If the input file is a sequence file (fasta or fastq), - * this function constructs the index for about mm_idxopt_t::batch_size bases. - * Importantly, for a huge collection of sequences, this function may only - * return an index for part of sequences. It needs to be repeatedly called - * to traverse the entire index/sequence file. - * - * @param r index reader - * @param n_threads number of threads for constructing index - * - * @return an index on success; NULL if reaching the end of the input file - */ -mm_idx_t *mm_idx_reader_read(mm_idx_reader_t *r, int n_threads); + /** + * Read/build an index + * + * If the input file is an index file, this function reads one part of the + * index and returns. If the input file is a sequence file (fasta or fastq), + * this function constructs the index for about mm_idxopt_t::batch_size bases. + * Importantly, for a huge collection of sequences, this function may only + * return an index for part of sequences. It needs to be repeatedly called + * to traverse the entire index/sequence file. + * + * @param r index reader + * @param n_threads number of threads for constructing index + * + * @return an index on success; NULL if reaching the end of the input file + */ + mm_idx_t *mm_idx_reader_read(mm_idx_reader_t *r, int n_threads); -/** - * Destroy/deallocate an index reader - * - * @param r index reader - */ -void mm_idx_reader_close(mm_idx_reader_t *r); + /** + * Destroy/deallocate an index reader + * + * @param r index reader + */ + void mm_idx_reader_close(mm_idx_reader_t *r); -int mm_idx_reader_eof(const mm_idx_reader_t *r); + int mm_idx_reader_eof(const mm_idx_reader_t *r); -/** - * Check whether the file contains a minimap2 index - * - * @param fn file name - * - * @return the file size if fn is an index file; 0 if fn is not. - */ -int64_t mm_idx_is_idx(const char *fn); + /** + * Check whether the file contains a minimap2 index + * + * @param fn file name + * + * @return the file size if fn is an index file; 0 if fn is not. + */ + int64_t mm_idx_is_idx(const char *fn); -/** - * Load a part of an index - * - * Given a uni-part index, this function loads the entire index into memory. - * Given a multi-part index, it loads one part only and places the file pointer - * at the end of that part. - * - * @param fp pointer to FILE object - * - * @return minimap2 index read from fp - */ -mm_idx_t *mm_idx_load(FILE *fp); + /** + * Load a part of an index + * + * Given a uni-part index, this function loads the entire index into memory. + * Given a multi-part index, it loads one part only and places the file pointer + * at the end of that part. + * + * @param fp pointer to FILE object + * + * @return minimap2 index read from fp + */ + mm_idx_t *mm_idx_load(FILE *fp); -/** - * Append an index (or one part of a full index) to file - * - * @param fp pointer to FILE object - * @param mi minimap2 index - */ -void mm_idx_dump(FILE *fp, const mm_idx_t *mi); + /** + * Append an index (or one part of a full index) to file + * + * @param fp pointer to FILE object + * @param mi minimap2 index + */ + void mm_idx_dump(FILE *fp, const mm_idx_t *mi); -/** - * Create an index from strings in memory - * - * @param w minimizer window size - * @param k minimizer k-mer size - * @param is_hpc use HPC k-mer if true - * @param bucket_bits number of bits for the first level of the hash table - * @param n number of sequences - * @param seq sequences in A/C/G/T - * @param name sequence names; could be NULL - * - * @return minimap2 index - */ -mm_idx_t *mm_idx_str(int w, int k, int is_hpc, int bucket_bits, int n, const char **seq, const char **name); + /** + * Create an index from strings in memory + * + * @param w minimizer window size + * @param k minimizer k-mer size + * @param is_hpc use HPC k-mer if true + * @param bucket_bits number of bits for the first level of the hash table + * @param n number of sequences + * @param seq sequences in A/C/G/T + * @param name sequence names; could be NULL + * + * @return minimap2 index + */ + mm_idx_t *mm_idx_str(int w, int k, int is_hpc, int bucket_bits, int n, const char **seq, const char **name); -/** - * Print index statistics to stderr - * - * @param mi minimap2 index - */ -void mm_idx_stat(const mm_idx_t *idx); + /** + * Print index statistics to stderr + * + * @param mi minimap2 index + */ + void mm_idx_stat(const mm_idx_t *idx); -/** - * Destroy/deallocate an index - * - * @param r minimap2 index - */ -void mm_idx_destroy(mm_idx_t *mi); + /** + * Destroy/deallocate an index + * + * @param r minimap2 index + */ + void mm_idx_destroy(mm_idx_t *mi); -/** - * Initialize a thread-local buffer for mapping - * - * Each mapping thread requires a buffer specific to the thread (see mm_map() - * below). The primary purpose of this buffer is to reduce frequent heap - * allocations across threads. A buffer shall not be used by two or more - * threads. - * - * @return pointer to a thread-local buffer - */ -mm_tbuf_t *mm_tbuf_init(void); + /** + * Initialize a thread-local buffer for mapping + * + * Each mapping thread requires a buffer specific to the thread (see mm_map() + * below). The primary purpose of this buffer is to reduce frequent heap + * allocations across threads. A buffer shall not be used by two or more + * threads. + * + * @return pointer to a thread-local buffer + */ + mm_tbuf_t *mm_tbuf_init(void); -/** - * Destroy/deallocate a thread-local buffer for mapping - * - * @param b the buffer - */ -void mm_tbuf_destroy(mm_tbuf_t *b); + /** + * Destroy/deallocate a thread-local buffer for mapping + * + * @param b the buffer + */ + void mm_tbuf_destroy(mm_tbuf_t *b); -void *mm_tbuf_get_km(mm_tbuf_t *b); + void *mm_tbuf_get_km(mm_tbuf_t *b); -/** - * Align a query sequence against an index - * - * This function possibly finds multiple alignments of the query sequence. - * The returned array and the mm_reg1_t::p field of each element are allocated - * with malloc(). - * - * @param mi minimap2 index - * @param l_seq length of the query sequence - * @param seq the query sequence - * @param n_regs number of hits (out) - * @param b thread-local buffer; two mm_map() calls shall not use one buffer at the same time! - * @param opt mapping parameters - * @param name query name, used for all-vs-all overlapping and debugging - * - * @return an array of hits which need to be deallocated with free() together - * with mm_reg1_t::p of each element. The size is written to _n_regs_. - */ -mm_reg1_t *mm_map(const mm_idx_t *mi, int l_seq, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *name); + /** + * Align a query sequence against an index + * + * This function possibly finds multiple alignments of the query sequence. + * The returned array and the mm_reg1_t::p field of each element are allocated + * with malloc(). + * + * @param mi minimap2 index + * @param l_seq length of the query sequence + * @param seq the query sequence + * @param n_regs number of hits (out) + * @param b thread-local buffer; two mm_map() calls shall not use one buffer at the same time! + * @param opt mapping parameters + * @param name query name, used for all-vs-all overlapping and debugging + * + * @return an array of hits which need to be deallocated with free() together + * with mm_reg1_t::p of each element. The size is written to _n_regs_. + */ + mm_reg1_t *mm_map(const mm_idx_t *mi, int l_seq, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *name); -void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname); + void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname); -/** - * Align a fasta/fastq file and print alignments to stdout - * - * @param idx minimap2 index - * @param fn fasta/fastq file name - * @param opt mapping parameters - * @param n_threads number of threads - * - * @return 0 on success; -1 if _fn_ can't be read - */ -int mm_map_file(const mm_idx_t *idx, const char *fn, const mm_mapopt_t *opt, int n_threads); + /** + * Align a fasta/fastq file and print alignments to stdout + * + * @param idx minimap2 index + * @param fn fasta/fastq file name + * @param opt mapping parameters + * @param n_threads number of threads + * + * @return 0 on success; -1 if _fn_ can't be read + */ + int mm_map_file(const mm_idx_t *idx, const char *fn, const mm_mapopt_t *opt, int n_threads); -int mm_map_file_frag(const mm_idx_t *idx, int n_segs, const char **fn, const mm_mapopt_t *opt, int n_threads); + int mm_map_file_frag(const mm_idx_t *idx, int n_segs, const char **fn, const mm_mapopt_t *opt, int n_threads); -/** - * Generate the cs tag (new in 2.12) - * - * @param km memory blocks; set to NULL if unsure - * @param buf buffer to write the cs/MD tag; typicall NULL on the first call - * @param max_len max length of the buffer; typically set to 0 on the first call - * @param mi index - * @param r alignment - * @param seq query sequence - * @param no_iden true to use : instead of = - * - * @return the length of cs - */ -int mm_gen_cs(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int no_iden); -int mm_gen_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq); + /** + * Generate the cs tag (new in 2.12) + * + * @param km memory blocks; set to NULL if unsure + * @param buf buffer to write the cs/MD tag; typicall NULL on the first call + * @param max_len max length of the buffer; typically set to 0 on the first call + * @param mi index + * @param r alignment + * @param seq query sequence + * @param no_iden true to use : instead of = + * + * @return the length of cs + */ + int mm_gen_cs(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int no_iden); + int mm_gen_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq); -// query sequence name and sequence in the minimap2 index -int mm_idx_index_name(mm_idx_t *mi); -int mm_idx_name2id(const mm_idx_t *mi, const char *name); -int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq); + // query sequence name and sequence in the minimap2 index + int mm_idx_index_name(mm_idx_t *mi); + int mm_idx_name2id(const mm_idx_t *mi, const char *name); + int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq); -int mm_idx_alt_read(mm_idx_t *mi, const char *fn); -int mm_idx_bed_read(mm_idx_t *mi, const char *fn, int read_junc); -int mm_idx_bed_junc(const mm_idx_t *mi, int32_t ctg, int32_t st, int32_t en, uint8_t *s); + int mm_idx_alt_read(mm_idx_t *mi, const char *fn); + int mm_idx_bed_read(mm_idx_t *mi, const char *fn, int read_junc); + int mm_idx_bed_junc(const mm_idx_t *mi, int32_t ctg, int32_t st, int32_t en, uint8_t *s); -// deprecated APIs for backward compatibility -void mm_mapopt_init(mm_mapopt_t *opt); -mm_idx_t *mm_idx_build(const char *fn, int w, int k, int flag, int n_threads); + // deprecated APIs for backward compatibility + void mm_mapopt_init(mm_mapopt_t *opt); + mm_idx_t *mm_idx_build(const char *fn, int w, int k, int flag, int n_threads); #ifdef __cplusplus } diff --git a/run.sh b/run.sh new file mode 100644 index 0000000..e69de29 diff --git a/seed.c b/seed.c index 76a67ae..4b64c47 100644 --- a/seed.c +++ b/seed.c @@ -6,13 +6,16 @@ void mm_seed_mz_flt(void *km, mm128_v *mv, int32_t q_occ_max, float q_occ_frac) { mm128_t *a; size_t i, j, st; - if (mv->n <= q_occ_max || q_occ_frac <= 0.0f || q_occ_max <= 0) return; + if (mv->n <= q_occ_max || q_occ_frac <= 0.0f || q_occ_max <= 0) + return; a = Kmalloc(km, mm128_t, mv->n); for (i = 0; i < mv->n; ++i) a[i].x = mv->a[i].x, a[i].y = i; radix_sort_128x(a, a + mv->n); - for (st = 0, i = 1; i <= mv->n; ++i) { - if (i == mv->n || a[i].x != a[st].x) { + for (st = 0, i = 1; i <= mv->n; ++i) + { + if (i == mv->n || a[i].x != a[st].x) + { int32_t cnt = i - st; if (cnt > q_occ_max && cnt > mv->n * q_occ_frac) for (j = st; j < i; ++j) @@ -32,20 +35,24 @@ mm_seed_t *mm_seed_collect_all(void *km, const mm_idx_t *mi, const mm128_v *mv, mm_seed_t *m; size_t i; int32_t k; - m = (mm_seed_t*)kmalloc(km, mv->n * sizeof(mm_seed_t)); - for (i = k = 0; i < mv->n; ++i) { + m = (mm_seed_t *)kmalloc(km, mv->n * sizeof(mm_seed_t)); // 为每一个minimizer开辟一个mm_seed_t + for (i = k = 0; i < mv->n; ++i) + { const uint64_t *cr; mm_seed_t *q; mm128_t *p = &mv->a[i]; uint32_t q_pos = (uint32_t)p->y, q_span = p->x & 0xff; - int t; - cr = mm_idx_get(mi, p->x>>8, &t); - if (t == 0) continue; + int t; // t表示hash值的低32位，表示啥？ + cr = mm_idx_get(mi, p->x >> 8, &t); // cr是hash值的高32位，代表位置 + if (t == 0) + continue; q = &m[k++]; q->q_pos = q_pos, q->q_span = q_span, q->cr = cr, q->n = t, q->seg_id = p->y >> 32; q->is_tandem = q->flt = 0; - if (i > 0 && p->x>>8 == mv->a[i - 1].x>>8) q->is_tandem = 1; - if (i < mv->n - 1 && p->x>>8 == mv->a[i + 1].x>>8) q->is_tandem = 1; + if (i > 0 && p->x >> 8 == mv->a[i - 1].x >> 8) + q->is_tandem = 1; + if (i < mv->n - 1 && p->x >> 8 == mv->a[i + 1].x >> 8) + q->is_tandem = 1; } *n_m_ = k; return m; @@ -55,37 +62,48 @@ mm_seed_t *mm_seed_collect_all(void *km, const mm_idx_t *mi, const mm128_v *mv, void mm_seed_select(int32_t n, mm_seed_t *a, int len, int max_occ, int max_max_occ, int dist) { // for high-occ minimizers, choose up to max_high_occ in each high-occ streak - extern void ks_heapdown_uint64_t(size_t i, size_t n, uint64_t*); - extern void ks_heapmake_uint64_t(size_t n, uint64_t*); + extern void ks_heapdown_uint64_t(size_t i, size_t n, uint64_t *); + extern void ks_heapmake_uint64_t(size_t n, uint64_t *); int32_t i, last0, m; uint64_t b[MAX_MAX_HIGH_OCC]; // this is to avoid a heap allocation - if (n == 0 || n == 1) return; + if (n == 0 || n == 1) + return; for (i = m = 0; i < n; ++i) - if (a[i].n > max_occ) ++m; - if (m == 0) return; // no high-frequency k-mers; do nothing - for (i = 0, last0 = -1; i <= n; ++i) { - if (i == n || a[i].n <= max_occ) { - if (i - last0 > 1) { - int32_t ps = last0 < 0? 0 : (uint32_t)a[last0].q_pos>>1; - int32_t pe = i == n? len : (uint32_t)a[i].q_pos>>1; + if (a[i].n > max_occ) + ++m; + if (m == 0) + return; // no high-frequency k-mers; do nothing + for (i = 0, last0 = -1; i <= n; ++i) + { + if (i == n || a[i].n <= max_occ) + { + if (i - last0 > 1) + { + int32_t ps = last0 < 0 ? 0 : (uint32_t)a[last0].q_pos >> 1; + int32_t pe = i == n ? len : (uint32_t)a[i].q_pos >> 1; int32_t j, k, st = last0 + 1, en = i; int32_t max_high_occ = (int32_t)((double)(pe - ps) / dist + .499); - if (max_high_occ > 0) { + if (max_high_occ > 0) + { if (max_high_occ > MAX_MAX_HIGH_OCC) max_high_occ = MAX_MAX_HIGH_OCC; for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k) - b[k] = (uint64_t)a[j].n<<32 | j; + b[k] = (uint64_t)a[j].n << 32 | j; ks_heapmake_uint64_t(k, b); // initialize the binomial heap - for (; j < en; ++j) { // if there are more, choose top max_high_occ - if (a[j].n < (int32_t)(b[0]>>32)) { // then update the heap - b[0] = (uint64_t)a[j].n<<32 | j; + for (; j < en; ++j) + { // if there are more, choose top max_high_occ + if (a[j].n < (int32_t)(b[0] >> 32)) + { // then update the heap + b[0] = (uint64_t)a[j].n << 32 | j; ks_heapdown_uint64_t(0, k, b); } } - for (j = 0; j < k; ++j) a[(uint32_t)b[j]].flt = 1; + for (j = 0; j < k; ++j) + a[(uint32_t)b[j]].flt = 1; } - for (j = st; j < en; ++j) a[j].flt ^= 1; + for (j = st; j < en; ++j) + a[j].flt ^= 1; for (j = st; j < en; ++j) if (a[j].n > max_max_occ) a[j].flt = 1; @@ -101,27 +119,37 @@ mm_seed_t *mm_collect_matches(void *km, int *_n_m, int qlen, int max_occ, int ma size_t i; mm_seed_t *m; *n_mini_pos = 0; - *mini_pos = (uint64_t*)kmalloc(km, mv->n * sizeof(uint64_t)); + *mini_pos = (uint64_t *)kmalloc(km, mv->n * sizeof(uint64_t)); m = mm_seed_collect_all(km, mi, mv, &n_m0); - if (dist > 0 && max_max_occ > max_occ) { + if (dist > 0 && max_max_occ > max_occ) + { mm_seed_select(n_m0, m, qlen, max_occ, max_max_occ, dist); - } else { + } + else + { for (i = 0; i < n_m0; ++i) if (m[i].n > max_occ) m[i].flt = 1; } - for (i = 0, n_m = 0, *rep_len = 0, *n_a = 0; i < n_m0; ++i) { + for (i = 0, n_m = 0, *rep_len = 0, *n_a = 0; i < n_m0; ++i) + { mm_seed_t *q = &m[i]; - //fprintf(stderr, "X\t%d\t%d\t%d\n", q->q_pos>>1, q->n, q->flt); - if (q->flt) { + // fprintf(stderr, "X\t%d\t%d\t%d\n", q->q_pos>>1, q->n, q->flt); + if (q->flt) + { int en = (q->q_pos >> 1) + 1, st = en - q->q_span; - if (st > rep_en) { + if (st > rep_en) + { *rep_len += rep_en - rep_st; rep_st = st, rep_en = en; - } else rep_en = en; - } else { + } + else + rep_en = en; + } + else + { *n_a += q->n; - (*mini_pos)[(*n_mini_pos)++] = (uint64_t)q->q_span<<32 | q->q_pos>>1; + (*mini_pos)[(*n_mini_pos)++] = (uint64_t)q->q_span << 32 | q->q_pos >> 1; m[n_m++] = *q; } } diff --git a/sketch.c b/sketch.c index 30230a2..dae2951 100644 --- a/sketch.c +++ b/sketch.c @@ -7,23 +7,22 @@ #include "mmpriv.h" unsigned char seq_nt4_table[256] = { - 0, 1, 2, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 -}; + 0, 1, 2, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4}; static inline uint64_t hash64(uint64_t key, uint64_t mask) { @@ -37,7 +36,8 @@ static inline uint64_t hash64(uint64_t key, uint64_t mask) return key; } -typedef struct { // a simplified version of kdq +typedef struct +{ // a simplified version of kdq int front, count; int a[32]; } tiny_queue_t; @@ -50,7 +50,8 @@ static inline void tq_push(tiny_queue_t *q, int x) static inline int tq_shift(tiny_queue_t *q) { int x; - if (q->count == 0) return -1; + if (q->count == 0) + return -1; x = q->a[q->front++]; q->front &= 0x1f; --q->count; @@ -76,24 +77,28 @@ static inline int tq_shift(tiny_queue_t *q) */ void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, int is_hpc, mm128_v *p) { - uint64_t shift1 = 2 * (k - 1), mask = (1ULL<<2*k) - 1, kmer[2] = {0,0}; + uint64_t shift1 = 2 * (k - 1), mask = (1ULL << 2 * k) - 1, kmer[2] = {0, 0}; int i, j, l, buf_pos, min_pos, kmer_span = 0; - mm128_t buf[256], min = { UINT64_MAX, UINT64_MAX }; + mm128_t buf[256], min = {UINT64_MAX, UINT64_MAX}; tiny_queue_t tq; assert(len > 0 && (w > 0 && w < 256) && (k > 0 && k <= 28)); // 56 bits for k-mer; could use long k-mers, but 28 enough in practice memset(buf, 0xff, w * 16); memset(&tq, 0, sizeof(tiny_queue_t)); - kv_resize(mm128_t, km, *p, p->n + len/w); + kv_resize(mm128_t, km, *p, p->n + len / w); // 扩充p，将新生成len/w个minimizer - for (i = l = buf_pos = min_pos = 0; i < len; ++i) { + for (i = l = buf_pos = min_pos = 0; i < len; ++i) + { int c = seq_nt4_table[(uint8_t)str[i]]; - mm128_t info = { UINT64_MAX, UINT64_MAX }; - if (c < 4) { // not an ambiguous base + mm128_t info = {UINT64_MAX, UINT64_MAX}; + if (c < 4) + { // not an ambiguous base int z; - if (is_hpc) { + if (is_hpc) + { int skip_len = 1; - if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) { + if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) + { for (skip_len = 2; i + skip_len < len; ++skip_len) if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c) break; @@ -101,42 +106,63 @@ void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, i } tq_push(&tq, skip_len); kmer_span += skip_len; - if (tq.count > k) kmer_span -= tq_shift(&tq); - } else kmer_span = l + 1 < k? l + 1 : k; - kmer[0] = (kmer[0] << 2 | c) & mask; // forward k-mer - kmer[1] = (kmer[1] >> 2) | (3ULL^c) << shift1; // reverse k-mer - if (kmer[0] == kmer[1]) continue; // skip "symmetric k-mers" as we don't know it strand - z = kmer[0] < kmer[1]? 0 : 1; // strand // kmer的strand到底是什么意思？为什么通过比较就能确定正反？ + if (tq.count > k) + kmer_span -= tq_shift(&tq); + } + else + kmer_span = l + 1 < k ? l + 1 : k; + kmer[0] = (kmer[0] << 2 | c) & mask; // forward k-mer + kmer[1] = (kmer[1] >> 2) | (3ULL ^ c) << shift1; // reverse k-mer + if (kmer[0] == kmer[1]) + continue; // skip "symmetric k-mers" as we don't know it strand + z = kmer[0] < kmer[1] ? 0 : 1; // strand // 选取小的那个kmer，kmer的strand到底是什么意思？为什么通过比较就能确定正反？ ++l; - if (l >= k && kmer_span < 256) { + if (l >= k && kmer_span < 256) + { info.x = hash64(kmer[z], mask) << 8 | kmer_span; - info.y = (uint64_t)rid<<32 | (uint32_t)i<<1 | z; + info.y = (uint64_t)rid << 32 | (uint32_t)i << 1 | z; } - } else l = 0, tq.count = tq.front = 0, kmer_span = 0; + } + else + l = 0, tq.count = tq.front = 0, kmer_span = 0; buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below - if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet + if (l == w + k - 1 && min.x != UINT64_MAX) + { // special case for the first window - because identical k-mers are not stored yet for (j = buf_pos + 1; j < w; ++j) - if (min.x == buf[j].x && buf[j].y != min.y) kv_push(mm128_t, km, *p, buf[j]); + if (min.x == buf[j].x && buf[j].y != min.y) + kv_push(mm128_t, km, *p, buf[j]); for (j = 0; j < buf_pos; ++j) - if (min.x == buf[j].x && buf[j].y != min.y) kv_push(mm128_t, km, *p, buf[j]); + if (min.x == buf[j].x && buf[j].y != min.y) + kv_push(mm128_t, km, *p, buf[j]); } - if (info.x <= min.x) { // a new minimum; then write the old min - if (l >= w + k && min.x != UINT64_MAX) kv_push(mm128_t, km, *p, min); + if (info.x <= min.x) + { // a new minimum; then write the old min + if (l >= w + k && min.x != UINT64_MAX) + kv_push(mm128_t, km, *p, min); min = info, min_pos = buf_pos; - } else if (buf_pos == min_pos) { // old min has moved outside the window - if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(mm128_t, km, *p, min); + } + else if (buf_pos == min_pos) + { // old min has moved outside the window + if (l >= w + k - 1 && min.x != UINT64_MAX) + kv_push(mm128_t, km, *p, min); for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers - if (min.x >= buf[j].x) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer + if (min.x >= buf[j].x) + min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer for (j = 0; j <= buf_pos; ++j) - if (min.x >= buf[j].x) min = buf[j], min_pos = j; - if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers - for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted - if (min.x == buf[j].x && min.y != buf[j].y) kv_push(mm128_t, km, *p, buf[j]); + if (min.x >= buf[j].x) + min = buf[j], min_pos = j; // 如果有多个min相同，取离当前位置最近的 + if (l >= w + k - 1 && min.x != UINT64_MAX) // 往回找相同值的kmer，放进p里 + { // write identical k-mers + for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted + if (min.x == buf[j].x && min.y != buf[j].y) + kv_push(mm128_t, km, *p, buf[j]); for (j = 0; j <= buf_pos; ++j) - if (min.x == buf[j].x && min.y != buf[j].y) kv_push(mm128_t, km, *p, buf[j]); + if (min.x == buf[j].x && min.y != buf[j].y) + kv_push(mm128_t, km, *p, buf[j]); } } - if (++buf_pos == w) buf_pos = 0; + if (++buf_pos == w) + buf_pos = 0; } if (min.x != UINT64_MAX) kv_push(mm128_t, km, *p, min);