diff --git a/bwa-postalt.js b/bwa-postalt.js index 2132a66..1580182 100644 --- a/bwa-postalt.js +++ b/bwa-postalt.js @@ -207,13 +207,13 @@ function bwa_postalt(args) { var c, opt = { a:1, b:4, o:6, e:1, verbose:3, show_pri:false, recover_mapq:true, min_mapq:10, min_sc:90, max_nm_sc:10, show_ev:false, min_pa_ratio:0.8 }; - while ((c = getopt(args, 'Pqev:p:r')) != null) { + while ((c = getopt(args, 'Pqev:p:r:')) != null) { if (c == 'v') opt.verbose = parseInt(getopt.arg); else if (c == 'p') opt.pre = getopt.arg; else if (c == 'P') opt.show_pri = true; else if (c == 'q') opt.recover_maq = false; else if (c == 'e') opt.show_ev = true; - else if (c == 'r') opt.min_pa_ratio = parseFloat(optarg.arg); + else if (c == 'r') opt.min_pa_ratio = parseFloat(getopt.arg); } if (opt.min_pa_ratio > 1.) opt.min_pa_ratio = 1.; @@ -485,20 +485,21 @@ function bwa_postalt(args) if (opt.recover_mapq && n_rpt_lifted == 1 && mapQ > 0) { var l = rpt_lifted; for (var i = 0; i < buf2.length; ++i) { - var s = buf2[i]; - if (l[0] != s[2]) continue; // different chr - if (((s[1]&16) != 0) != l[1]) continue; // different strand + var s = buf2[i], is_ovlp = true; + if (l[0] != s[2]) is_ovlp = false; // different chr + if (((s[1]&16) != 0) != l[1]) is_ovlp = false; // different strand var start = s[3] - 1, end = start; while ((m = re_cigar.exec(t[5])) != null) if (m[2] == 'M' || m[2] == 'D' || m[2] == 'N') end += parseInt(m[1]); + if (!(start < l[3] && l[2] < end)) is_ovlp = false; // no overlap var om = -1, pa = 10.; for (var j = 11; j < s.length; ++j) if ((m = /^om:i:(\d+)/.exec(s[j])) != null) om = parseInt(m[1]); else if ((m = /^pa:f:(\S+)/.exec(s[j])) != null) pa = parseFloat(m[1]); - if (start < l[3] && l[2] < end) { // overlapping the lifted hit + if (is_ovlp) { // overlapping the lifted hit if (om > 0) s[4] = om; s[4] = s[4] < mapQ? s[4] : mapQ; } else if (pa < opt.min_pa_ratio) { // not overlapping; has a small pa diff --git a/bwt.c b/bwt.c index a719f53..833e37d 100644 --- a/bwt.c +++ b/bwt.c @@ -305,10 +305,10 @@ int bwt_smem1a(const bwt_t *bwt, int len, const uint8_t *q, int x, int min_intv, if (q[i] < 4) { // an A/C/G/T base c = 3 - q[i]; // complement of q[i] bwt_extend(bwt, &ik, ok, 0); + if (ik.x[2] < max_intv && (i - x > max_len || ik.x[2] == 1)) break; if (ok[c].x[2] != ik.x[2]) { // change of the interval size kv_push(bwtintv_t, *curr, ik); if (ok[c].x[2] < min_intv) break; // the interval size is too small to be extended further - if (i - x > max_len && ik.x[2] < max_intv) break; } ik = ok[c]; ik.info = i + 1; } else { // an ambiguous base