zzh
/
fast-bwa
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

mapQ is not quite right

This commit is contained in:
Heng Li 2014-09-19 23:39:59 -04:00
parent a33ca20499
commit 5fc0f44e27
1 changed files with 17 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
bwa-postalt.js
View File

@ -351,12 +351,9 @@ function bwa_postalt(args)
for (var i = 0; i < hits.length; ++i) for (var i = 0; i < hits.length; ++i)
if (hits[i].g == reported_g) if (hits[i].g == reported_g)
++n_group0; ++n_group0;
if (n_group0 == 1) { // then keep the reported alignment and mapQ
buf2.push(t);
continue;
}
// re-estimate mapQ var mapQ, ori_mapQ = t[4];
if (n_group0 > 1) { // then re-estimate mapQ
var group_max = []; var group_max = [];
for (var i = 0; i < hits.length; ++i) { for (var i = 0; i < hits.length; ++i) {
var g = hits[i].g; var g = hits[i].g;
@ -365,13 +362,12 @@ function bwa_postalt(args)
} }
if (group_max.length > 1) if (group_max.length > 1)
group_max.sort(function(x,y) {return y[0]-x[0]}); group_max.sort(function(x,y) {return y[0]-x[0]});
var mapQ;
if (group_max[0][1] == reported_g) { // the best hit is the hit reported in SAM if (group_max[0][1] == reported_g) { // the best hit is the hit reported in SAM
mapQ = group_max.length == 1? 60 : 6 * (group_max[0][0] - group_max[1][0]); mapQ = group_max.length == 1? 60 : 6 * (group_max[0][0] - group_max[1][0]);
} else mapQ = 0; } else mapQ = 0;
mapQ = mapQ < 60? mapQ : 60; mapQ = mapQ < 60? mapQ : 60;
var ori_mapQ = t[4];
mapQ = mapQ > ori_mapQ? mapQ : ori_mapQ; mapQ = mapQ > ori_mapQ? mapQ : ori_mapQ;
} else mapQ = t[4];
// check if the reported hit overlaps a hit to the primary assembly; if so, don't reduce mapping quality // check if the reported hit overlaps a hit to the primary assembly; if so, don't reduce mapping quality
if (opt.recover_mapq && n_rpt_lifted == 1 && mapQ > 0) { if (opt.recover_mapq && n_rpt_lifted == 1 && mapQ > 0) {
@ -420,7 +416,7 @@ function bwa_postalt(args)
// print // print
t[4] = mapQ; t[4] = mapQ;
t.push("om:i:"+ori_mapQ); if (n_group0 > 1) t.push("om:i:"+ori_mapQ);
if (hits[reported_i].lifted_str) t.push("lt:Z:" + hits[reported_i].lifted_str); if (hits[reported_i].lifted_str) t.push("lt:Z:" + hits[reported_i].lifted_str);
buf2.push(t); buf2.push(t);
var cnt = 0; var cnt = 0;