diff --git a/bwa-helper.js b/bwa-helper.js index daa1dd8..ea9c334 100644 --- a/bwa-helper.js +++ b/bwa-helper.js @@ -520,7 +520,7 @@ function read_ALT_sam(fn) return idx; } -function bwa_genalt(args) +function bwa_altgen(args) { var c, opt = { a:1, b:4, o:6, e:1, verbose:3 }; @@ -529,7 +529,7 @@ function bwa_genalt(args) } if (args.length == getopt.ind) { - print("Usage: k8 bwa-helper.js genalt [aln.sam]"); + print("Usage: k8 bwa-helper.js altgen [aln.sam]"); exit(1); } @@ -690,22 +690,16 @@ function bwa_genalt(args) } // This is in effect a simplified version of bwa_genalt(). -function bwa_postalt(args) +function bwa_altlift(args) { - var c, idx = {}, opt = { verbose:3, min_pa:0.8 }; - while ((c = getopt(args, 'v:r:a:')) != null) { - if (c == 'v') opt.verbose = parseInt(getopt.arg); - else if (c == 'r') opt.min_pa = parseFloat(getopt.arg); - else if (c == 'a') idx = read_ALT_sam(getopt.arg); - } - - if (args.length == getopt.ind) { - print("Usage: k8 bwa-helper.js postalt [-r "+opt.min_pa+"] [-a alt.sam] [aln.sam]"); + if (args.length == 0) { + print("Usage: k8 bwa-helper.js altlift [aln.sam]"); exit(1); } + var idx = read_ALT_sam(args[0]); // process SAM - var file = args.length - getopt.ind >= 1? new File(args[getopt.ind]) : new File(); + var file = args.length >= 2? new File(args[1]) : new File(); var buf = new Bytes(); while (file.readline(buf) >= 0) { var m, line = buf.toString(); @@ -720,36 +714,27 @@ function bwa_postalt(args) var h = parse_hit([t[2], ((flag&16)?'-':'+') + t[3], t[5], NM], opt); // lift mapping positions to coordinates on the primary assembly - { - var a = null; - if (idx[h.ctg] != null) - a = idx[h.ctg](h.start, h.end); - if (a == null) a = []; + var a = null; + if (idx[h.ctg] != null) + a = idx[h.ctg](h.start, h.end); + if (a == null) a = []; - // find the approximate position on the primary assembly - var lifted = []; - for (var j = 0; j < a.length; ++j) { - var s, e; - if (!a[j][4]) { // ALT is mapped to the forward strand of the primary assembly - s = cigar2pos(a[j][6], h.start); - e = cigar2pos(a[j][6], h.end - 1) + 1; - } else { - s = cigar2pos(a[j][6], a[j][2] - h.end); - e = cigar2pos(a[j][6], a[j][2] - h.start - 1) + 1; - } - if (s < 0 || e < 0) continue; // read is mapped to clippings in the ALT-to-chr alignment - s += a[j][5]; e += a[j][5]; - lifted.push([a[j][3], (h.rev!=a[j][4]), s, e]); + // find the approximate position on the primary assembly + var lifted = []; + for (var j = 0; j < a.length; ++j) { + var s, e; + if (!a[j][4]) { // ALT is mapped to the forward strand of the primary assembly + s = cigar2pos(a[j][6], h.start); + e = cigar2pos(a[j][6], h.end - 1) + 1; + } else { + s = cigar2pos(a[j][6], a[j][2] - h.end); + e = cigar2pos(a[j][6], a[j][2] - h.start - 1) + 1; } - h.lifted = lifted; - } - - // update mapQ if necessary - var ori_mapQ = null; - if ((m = /\tpa:f:(\d+\.\d+)/.exec(line)) != null) { - if (parseFloat(m[1]) < opt.min_pa) - ori_mapQ = t[4], t[4] = 0; + if (s < 0 || e < 0) continue; // read is mapped to clippings in the ALT-to-chr alignment + s += a[j][5]; e += a[j][5]; + lifted.push([a[j][3], (h.rev!=a[j][4]), s, e]); } + h.lifted = lifted; // generate lifted_str if (h.lifted && h.lifted.length) { @@ -761,7 +746,6 @@ function bwa_postalt(args) } else h.lifted_str = null; // print - if (ori_mapQ != null) t.push("om:i:"+ori_mapQ); if (h.lifted_str) t.push("lt:Z:" + h.lifted_str); print(t.join("\t")); } @@ -777,8 +761,8 @@ function main(args) { if (args.length == 0) { print("\nUsage: k8 bwa-helper.js [arguments]\n"); - print("Commands: postalt post process ALT-aware alignments"); - print(" genalt generate ALT alignments for ALT-unaware alignments"); + print("Commands: altlift add lt tag to show lifted position on the primary assembly"); + print(" altgen generate ALT alignments for ALT-unaware alignments\n"); print(" sam2pas convert SAM to pairwise alignment summary format (PAS)"); print(" pas2reg extract covered regions"); print(" reg2cut regions to extract for the 2nd round bwa-mem"); @@ -795,8 +779,8 @@ function main(args) else if (cmd == 'markovlp') bwa_markOvlp(args); else if (cmd == 'pas2reg') bwa_pas2reg(args); else if (cmd == 'reg2cut') bwa_reg2cut(args); - else if (cmd == 'genalt') bwa_genalt(args); - else if (cmd == 'postalt') bwa_postalt(args); + else if (cmd == 'altgen' || cmd == 'genalt') bwa_alt(args); + else if (cmd == 'altlift') bwa_altlift(args); else if (cmd == 'shortname') bwa_shortname(args); else warn("Unrecognized command"); } diff --git a/bwamem.c b/bwamem.c index 7b6b1f4..7ba6c57 100644 --- a/bwamem.c +++ b/bwamem.c @@ -527,7 +527,7 @@ int mem_mark_primary_se(const mem_opt_t *opt, int n, mem_alnreg_t *a, int64_t id if (!p->is_alt && p->secondary >= 0 && a[p->secondary].is_alt) p->alt_sc = a[p->secondary].score; } - if (n_pri > 0 || n_pri != n) { + if (n_pri > 0 && n_pri < n) { ks_introsort(mem_ars_hash2, n, a); for (i = 0; i < n_pri; ++i) a[i].sub = 0, a[i].secondary = -1; mem_mark_primary_se_core(opt, n_pri, a, &z); diff --git a/bwamem_pair.c b/bwamem_pair.c index f754743..3882784 100644 --- a/bwamem_pair.c +++ b/bwamem_pair.c @@ -151,6 +151,7 @@ int mem_matesw(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, co memset(&b, 0, sizeof(mem_alnreg_t)); if (aln.score >= opt->min_seed_len && aln.qb >= 0) { // something goes wrong if aln.qb < 0 b.rid = a->rid; + b.is_alt = a->is_alt; b.qb = is_rev? l_ms - (aln.qe + 1) : aln.qb; b.qe = is_rev? l_ms - aln.qb : aln.qe + 1; b.rb = is_rev? (l_pac<<1) - (rb + aln.te + 1) : rb + aln.tb; diff --git a/main.c b/main.c index bedbe43..063ddda 100644 --- a/main.c +++ b/main.c @@ -4,7 +4,7 @@ #include "utils.h" #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.7.10-r855-dirty" +#define PACKAGE_VERSION "0.7.10-r858-dirty" #endif int bwa_fa2pac(int argc, char *argv[]);