mirror of
https://github.com/lh3/minimap2.git
synced 2026-09-30 21:18:12 +08:00
Merge branch 'master' into utec
This commit is contained in:
Heng Li
+447
-22
@@ -1,6 +1,6 @@
|
||||
#!/usr/bin/env k8
|
||||
|
||||
var paftools_version = '2.17-r982-dirty';
|
||||
var paftools_version = '2.21-r1071';
|
||||
|
||||
/*****************************
|
||||
***** Library functions *****
|
||||
@@ -962,12 +962,13 @@ function paf_asmgene(args)
|
||||
|
||||
function paf_stat(args)
|
||||
{
|
||||
var c, gap_out_len = null;
|
||||
while ((c = getopt(args, "l:")) != null)
|
||||
var c, gap_out_len = null, count_err = false;
|
||||
while ((c = getopt(args, "cl:")) != null)
|
||||
if (c == 'l') gap_out_len = parseInt(getopt.arg);
|
||||
else if (c == 'c') count_err = true;
|
||||
|
||||
if (getopt.ind == args.length) {
|
||||
print("Usage: paftools.js stat [-l gapOutLen] <in.sam>|<in.paf>");
|
||||
print("Usage: paftools.js stat [-c] [-l gapOutLen] <in.sam>|<in.paf>");
|
||||
exit(1);
|
||||
}
|
||||
|
||||
@@ -998,7 +999,7 @@ function paf_stat(args)
|
||||
if (line.charAt(0) != '@') {
|
||||
var t = line.split("\t", 12);
|
||||
var m, rs, cigar = null, is_pri = false, is_sam = false, is_rev = false, tname = null;
|
||||
var atlen = null, aqlen, qs, qe, mapq, ori_qlen;
|
||||
var atlen = null, aqlen, qs, qe, mapq, ori_qlen, NM = null;
|
||||
if (t.length < 2) continue;
|
||||
if (t[4] == '+' || t[4] == '-' || t[4] == '*') { // PAF
|
||||
if (t[4] == '*') continue; // unmapped
|
||||
@@ -1006,6 +1007,8 @@ function paf_stat(args)
|
||||
++n_2nd;
|
||||
continue;
|
||||
}
|
||||
if ((m = /\tNM:i:(\d+)/.exec(line)) != null)
|
||||
NM = parseInt(m[1]);
|
||||
if ((m = /\tcg:Z:(\S+)/.exec(line)) != null)
|
||||
cigar = m[1];
|
||||
if (cigar == null) {
|
||||
@@ -1027,6 +1030,8 @@ function paf_stat(args)
|
||||
++n_2nd;
|
||||
continue;
|
||||
}
|
||||
if ((m = /\tNM:i:(\d+)/.exec(line)) != null)
|
||||
NM = parseInt(m[1]);
|
||||
cigar = t[5];
|
||||
tname = t[2];
|
||||
rs = parseInt(t[3]) - 1;
|
||||
@@ -1045,11 +1050,13 @@ function paf_stat(args)
|
||||
++n_seq, last = t[0];
|
||||
}
|
||||
var M = 0, tl = 0, ql = 0, clip = [0, 0], n_cigar = 0, sclip = 0;
|
||||
var n_gapo = 0, n_gap_all = 0, l_match = 0;
|
||||
while ((m = re.exec(cigar)) != null) {
|
||||
var l = parseInt(m[1]);
|
||||
++n_cigar;
|
||||
if (m[2] == 'M' || m[2] == '=' || m[2] == 'X') {
|
||||
tl += l, ql += l, M += l;
|
||||
l_match += l;
|
||||
} else if (m[2] == 'I' || m[2] == 'D') {
|
||||
var type;
|
||||
if (l < 50) type = 0;
|
||||
@@ -1062,6 +1069,7 @@ function paf_stat(args)
|
||||
else tl += l, ++n_gap[1][type];
|
||||
if (gap_out_len != null && l >= gap_out_len)
|
||||
print(t[0], ql, is_rev? '-' : '+', tname, rs + tl, m[2], l);
|
||||
++n_gapo, n_gap_all += l;
|
||||
} else if (m[2] == 'N') {
|
||||
tl += l;
|
||||
} else if (m[2] == 'S') {
|
||||
@@ -1079,6 +1087,12 @@ function paf_stat(args)
|
||||
qs = clip[is_rev? 1 : 0], qe = qs + ql;
|
||||
ori_qlen = clip[0] + ql + clip[1];
|
||||
}
|
||||
if (count_err && NM != null) {
|
||||
var n_mm = NM - n_gap_all;
|
||||
if (n_mm < 0) warn("WARNING: NM is smaller than the number of gaps at line " + lineno);
|
||||
if (n_mm < 0) n_mm = 0;
|
||||
print(t[0], ori_qlen, t[11], ori_qlen - (qe - qs), NM, l_match + n_gap_all, n_mm + n_gapo, l_match + n_gapo);
|
||||
}
|
||||
regs.push([qs, qe]);
|
||||
last_qlen = ori_qlen;
|
||||
}
|
||||
@@ -1091,7 +1105,7 @@ function paf_stat(args)
|
||||
file.close();
|
||||
buf.destroy();
|
||||
|
||||
if (gap_out_len == null) {
|
||||
if (gap_out_len == null && !count_err) {
|
||||
print("Number of mapped sequences: " + n_seq);
|
||||
print("Number of primary alignments: " + n_pri);
|
||||
print("Number of secondary alignments: " + n_2nd);
|
||||
@@ -1405,7 +1419,7 @@ function paf_view(args)
|
||||
|
||||
var s_ref = new Bytes(), s_qry = new Bytes(), s_mid = new Bytes(); // these are used to show padded alignment
|
||||
var re_cs = /([:=\-\+\*])(\d+|[A-Za-z]+)/g;
|
||||
var re_cg = /(\d+)([MIDNSH])/g;
|
||||
var re_cg = /(\d+)([MIDNSHP=X])/g;
|
||||
|
||||
var buf = new Bytes();
|
||||
var file = args[getopt.ind] == "-"? new File() : new File(args[getopt.ind]);
|
||||
@@ -1885,7 +1899,7 @@ function paf_splice2bed(args)
|
||||
a.length = 0;
|
||||
}
|
||||
|
||||
var re = /(\d+)([MIDNSH])/g;
|
||||
var re = /(\d+)([MIDNSHP=X])/g;
|
||||
var c, fmt = "bed", fn_name_conv = null, keep_multi = false;
|
||||
while ((c = getopt(args, "f:n:m")) != null) {
|
||||
if (c == 'f') fmt = getopt.arg;
|
||||
@@ -2355,28 +2369,43 @@ function paf_junceval(args)
|
||||
|
||||
file = getopt.ind+1 >= args.length || args[getopt.ind+1] == '-'? new File() : new File(args[getopt.ind+1]);
|
||||
var last_qname = null;
|
||||
var re_cigar = /(\d+)([MIDNSHX=])/g;
|
||||
var re_cigar = /(\d+)([MIDNSHP=X])/g;
|
||||
while (file.readline(buf) >= 0) {
|
||||
var m, t = buf.toString().split("\t");
|
||||
var ctg_name = null, cigar = null, pos = null, qname = t[0];
|
||||
|
||||
if (t[0].charAt(0) == '@') continue;
|
||||
if (chr_only && !/^(chr)?([0-9]+|X|Y)$/.test(t[2])) continue;
|
||||
var flag = parseInt(t[1]);
|
||||
if (flag&0x100) continue;
|
||||
if (first_only && last_qname == t[0]) continue;
|
||||
if (t[2] == '*') {
|
||||
if (t[4] == '+' || t[4] == '-' || t[4] == '*') { // PAF
|
||||
ctg_name = t[5], pos = parseInt(t[7]);
|
||||
var type = 'P';
|
||||
for (i = 12; i < t.length; ++i) {
|
||||
if ((m = /^(tp:A|cg:Z):(\S+)/.exec(t[i])) != null) {
|
||||
if (m[1] == 'tp:A') type = m[2];
|
||||
else cigar = m[2];
|
||||
}
|
||||
}
|
||||
if (type == 'S') continue; // secondary
|
||||
} else { // SAM
|
||||
ctg_name = t[2], pos = parseInt(t[3]) - 1, cigar = t[5];
|
||||
var flag = parseInt(t[1]);
|
||||
if (flag&0x100) continue; // secondary
|
||||
}
|
||||
|
||||
if (chr_only && !/^(chr)?([0-9]+|X|Y)$/.test(ctg_name)) continue;
|
||||
if (first_only && last_qname == qname) continue;
|
||||
if (ctg_name == '*') { // unmapped
|
||||
++n_unmapped;
|
||||
continue;
|
||||
} else {
|
||||
++n_pri;
|
||||
if (last_qname != t[0]) {
|
||||
if (last_qname != qname) {
|
||||
++n_mapped;
|
||||
last_qname = t[0];
|
||||
last_qname = qname;
|
||||
}
|
||||
}
|
||||
|
||||
var pos = parseInt(t[3]) - 1, intron = [];
|
||||
while ((m = re_cigar.exec(t[5])) != null) {
|
||||
var intron = [];
|
||||
while ((m = re_cigar.exec(cigar)) != null) {
|
||||
var len = parseInt(m[1]), op = m[2];
|
||||
if (op == 'N') {
|
||||
intron.push([pos, pos + len]);
|
||||
@@ -2389,7 +2418,7 @@ function paf_junceval(args)
|
||||
}
|
||||
n_splice += intron.length;
|
||||
|
||||
var chr = anno[t[2]];
|
||||
var chr = anno[ctg_name];
|
||||
if (chr != null) {
|
||||
for (var i = 0; i < intron.length; ++i) {
|
||||
var o = Interval.find_ovlp(chr, intron[i][0], intron[i][1]);
|
||||
@@ -2413,12 +2442,12 @@ function paf_junceval(args)
|
||||
x += '(' + o[j][0] + "," + o[j][1] + ')';
|
||||
}
|
||||
x += ']';
|
||||
print(type, t[0], i+1, t[2], intron[i][0], intron[i][1], x);
|
||||
print(type, qname, i+1, ctg_name, intron[i][0], intron[i][1], x);
|
||||
}
|
||||
} else {
|
||||
++n_splice_novel;
|
||||
if (print_ovlp)
|
||||
print('N', t[0], i+1, t[2], intron[i][0], intron[i][1]);
|
||||
print('N', qname, i+1, ctg_name, intron[i][0], intron[i][1]);
|
||||
}
|
||||
}
|
||||
} else {
|
||||
@@ -2512,6 +2541,395 @@ function paf_ov_eval(args)
|
||||
print((100 * (1 - n_missing / n_ovlp)).toFixed(2) + "% sensitivity");
|
||||
}
|
||||
|
||||
function paf_vcfstat(args)
|
||||
{
|
||||
var c, ts = { "AG":1, "GA":1, "CT":1, "TC":1 };
|
||||
while ((c = getopt(args, "")) != null) {
|
||||
}
|
||||
var buf = new Bytes();
|
||||
var file = args.length == getopt.ind? new File() : new File(args[getopt.ind]);
|
||||
var x = { sub:0, ts:0, tv:0, ins:0, del:0, ins1:0, del1:0, ins2:0, del2:0, ins50:0, del50:0, ins1k:0, del1k:0, ins7k:0, del7k:0, insinf:0, delinf:0 };
|
||||
while (file.readline(buf) >= 0) {
|
||||
var t = buf.toString().split("\t");
|
||||
if (t[0][0] == '#') continue;
|
||||
var alt = t[4].split(",");
|
||||
var ref = t[3];
|
||||
for (var i = 0; i < alt.length; ++i) {
|
||||
var a = alt[i];
|
||||
if (a[0] == '<' || a[1] == '>') continue;
|
||||
var l = ref.length < a.length? ref.length : a.length;
|
||||
for (var j = 0; j < l; ++j) {
|
||||
if (ref[j] != a[j]) {
|
||||
++x.sub;
|
||||
if (ts[ref[j] + a[j]]) ++x.ts;
|
||||
else ++x.tv;
|
||||
}
|
||||
}
|
||||
var d = a.length - ref.length;
|
||||
if (d > 0) {
|
||||
++x.ins;
|
||||
if (d == 1) ++x.ins1;
|
||||
else if (d == 2) ++x.ins2;
|
||||
else if (d < 50) ++x.ins50;
|
||||
else if (d < 1000) ++x.ins1k;
|
||||
else if (d < 7000) ++x.ins7k;
|
||||
else ++x.insinf;
|
||||
} else if (d < 0) {
|
||||
d = -d;
|
||||
++x.del;
|
||||
if (d == 1) ++x.del1;
|
||||
else if (d == 2) ++x.del2;
|
||||
else if (d < 50) ++x.del50;
|
||||
else if (d < 1000) ++x.del1k;
|
||||
else if (d < 7000) ++x.del7k;
|
||||
else ++x.delinf;
|
||||
}
|
||||
}
|
||||
}
|
||||
file.close();
|
||||
buf.destroy();
|
||||
print("# substitutions: " + x.sub);
|
||||
print("ts/tv: " + (x.ts / x.tv).toFixed(3));
|
||||
print("# insertions: " + x.ins);
|
||||
print("# 1bp insertions: " + x.ins1);
|
||||
print("# 2bp insertions: " + x.ins2);
|
||||
print("# [3,50) insertions: " + x.ins50);
|
||||
print("# [50,1000) insertions: " + x.ins1k);
|
||||
print("# [1000,7000) insertions: " + x.ins7k);
|
||||
print("# >=7000 insertions: " + x.insinf);
|
||||
print("# deletions: " + x.del);
|
||||
print("# 1bp deletions: " + x.del1);
|
||||
print("# 2bp deletions: " + x.del2);
|
||||
print("# [3,50) deletions: " + x.del50);
|
||||
print("# [50,1000) deletions: " + x.del1k);
|
||||
print("# [1000,7000) deletions: " + x.del7k);
|
||||
print("# >=7000 deletions: " + x.delinf);
|
||||
}
|
||||
|
||||
function paf_parseNum(s) {
|
||||
var m, x = null;
|
||||
if ((m = /^(\d*\.?\d*)([mMgGkK]?)/.exec(s)) != null) {
|
||||
x = parseFloat(m[1]);
|
||||
if (m[2] == 'k' || m[2] == 'K') x *= 1000;
|
||||
else if (m[2] == 'm' || m[2] == 'M') x *= 1000000;
|
||||
else if (m[2] == 'g' || m[2] == 'G') x *= 1000000000;
|
||||
}
|
||||
return Math.floor(x + .499);
|
||||
}
|
||||
|
||||
function paf_misjoin(args)
|
||||
{
|
||||
var c, min_seg_len = 1000000, max_gap = 1000000, fn_cen = null, show_long = false, show_err = false, cen_ratio = 0.5;
|
||||
var n_diff = [0, 0], n_gap = [0, 0], n_inv = [0, 0], n_inv_end = [0, 0];
|
||||
while ((c = getopt(args, "l:g:c:per:")) != null) {
|
||||
if (c == 'l') min_seg_len = paf_parseNum(getopt.arg);
|
||||
else if (c == 'g') max_gap = paf_parseNum(getopt.arg);
|
||||
else if (c == 'c') fn_cen = getopt.arg;
|
||||
else if (c == 'r') cen_ratio = parseFloat(getopt.arg);
|
||||
else if (c == 'p') show_long = true;
|
||||
else if (c == 'e') show_err = true;
|
||||
}
|
||||
if (args.length == getopt.ind) {
|
||||
print("Usage: paftools.js misjoin [options] <in.paf>");
|
||||
print("Options:");
|
||||
print(" -c FILE BED for centromeres []");
|
||||
print(" -r FLOAT count a centromeric event if overlap ratio > FLOAT [" + cen_ratio + "]");
|
||||
print(" -l NUM min alignment block length [1m]");
|
||||
print(" -g NUM max gap size [1m]");
|
||||
print(" -e output misjoins not involving centromeres");
|
||||
print(" -p output long alignment blocks for debugging");
|
||||
return;
|
||||
}
|
||||
var cen = {};
|
||||
var file, buf = new Bytes();
|
||||
if (fn_cen != null) {
|
||||
file = new File(fn_cen);
|
||||
while (file.readline(buf) >= 0) {
|
||||
var t = buf.toString().split("\t");
|
||||
if (cen[t[0]] == null) cen[t[0]] = [];
|
||||
cen[t[0]].push([parseInt(t[1]), parseInt(t[2])]);
|
||||
}
|
||||
file.close();
|
||||
}
|
||||
|
||||
function test_cen(cen, chr, st, en) {
|
||||
var b = cen[chr], len = 0;
|
||||
if (b == null) return false;
|
||||
for (var j = 0; j < b.length; ++j)
|
||||
if (b[j][0] < en && b[j][1] > st) {
|
||||
var s = b[j][0] > st? b[j][0] : st;
|
||||
var e = b[j][1] < en? b[j][1] : en;
|
||||
len += e - s;
|
||||
}
|
||||
return len < (en - st) * cen_ratio? false : true;
|
||||
}
|
||||
|
||||
function process(a) {
|
||||
var k = 0;
|
||||
for (var i = 0; i < a.length; ++i) {
|
||||
for (var j = 1; j <= 3; ++j) a[i][j] = parseInt(a[i][j]);
|
||||
for (var j = 6; j <= 11; ++j) a[i][j] = parseInt(a[i][j]);
|
||||
if (a[i][10] >= min_seg_len) a[k++] = a[i];
|
||||
}
|
||||
a.length = k;
|
||||
if (a.length == 1) return;
|
||||
a = a.sort(function(x,y){return x[2]-y[2]});
|
||||
if (show_long) for (var i = 0; i < a.length; ++i) print(a[i].join("\t"));
|
||||
for (var i = 1; i < a.length; ++i) {
|
||||
var ov = [false, false];
|
||||
ov[0] = test_cen(cen, a[i-1][5], a[i-1][7], a[i-1][8]);
|
||||
ov[1] = test_cen(cen, a[i][5], a[i][7], a[i][8]);
|
||||
if (a[i-1][5] != a[i][5]) { // different chr
|
||||
if (ov[0] || ov[1]) ++n_diff[1];
|
||||
else if (show_err) {
|
||||
print("J", a[i-1].slice(0, 12).join("\t"));
|
||||
print("J", a[i].slice(0, 12).join("\t"));
|
||||
}
|
||||
++n_diff[0];
|
||||
} else if (a[i-1][4] == a[i][4]) { // a gap
|
||||
var dq = a[i][2] - a[i-1][3];
|
||||
var dr = a[i][4] == '+'? a[i][7] - a[i-1][8] : a[i-1][7] - a[i][8];
|
||||
var gap = dr > dq? dr - dq : dq - dr;
|
||||
if (gap > max_gap) {
|
||||
if (ov[0] || ov[1]) ++n_gap[1];
|
||||
else if (show_err) {
|
||||
print("G", a[i-1].slice(0, 12).join("\t"));
|
||||
print("G", a[i].slice(0, 12).join("\t"));
|
||||
}
|
||||
++n_gap[0];
|
||||
}
|
||||
} else if (i + 1 < a.length && a[i+1][4] == a[i-1][4]) { // bracketed inversion
|
||||
if (ov[0] || ov[1]) ++n_inv[1];
|
||||
else if (show_err) {
|
||||
print("M", a[i-1].slice(0, 12).join("\t"));
|
||||
print("M", a[i].slice(0, 12).join("\t"));
|
||||
print("M", a[i+1].slice(0, 12).join("\t"));
|
||||
}
|
||||
++n_inv[0];
|
||||
++i;
|
||||
} else { // hanging inversion
|
||||
if (ov[0] || ov[1]) ++n_inv_end[1];
|
||||
++n_inv_end[0];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
file = args[getopt.ind] == "-"? new File() : new File(args[getopt.ind]);
|
||||
var a = [];
|
||||
while (file.readline(buf) >= 0) {
|
||||
var t = buf.toString().split("\t");
|
||||
if (a.length > 0 && a[0][0] != t[0]) {
|
||||
process(a);
|
||||
a.length = 0;
|
||||
}
|
||||
a.push(t);
|
||||
}
|
||||
if (a.length > 0) process(a);
|
||||
file.close();
|
||||
buf.destroy();
|
||||
print("# inter-chromosomal misjoins: " + n_diff.join(","));
|
||||
print("# intra-chromosomal gaps: " + n_gap.join(","));
|
||||
print("# candidate inversions in the middle: " + n_inv.join(","));
|
||||
print("# candidate inversions at contig ends: " + n_inv_end.join(","));
|
||||
}
|
||||
|
||||
function _paf_get_alen(t)
|
||||
{
|
||||
var svlen = null, alen = null;
|
||||
if ((m = /(^|;)SVLEN=(-?\d+)/.exec(t[7])) != null)
|
||||
svlen = parseInt(m[2]);
|
||||
var s = t[4].split(",");
|
||||
var min_abs_diff = 1<<30, max_abs_diff = 0;
|
||||
if (svlen != null && svlen != 0)
|
||||
alen = svlen, min_abs_diff = max_abs_diff = svlen > 0? svlen : -svlen;
|
||||
var rlen = t[3].length;
|
||||
for (var i = 0; i < s.length; ++i) {
|
||||
if (/^<\S+>$/.test(s[i])) continue;
|
||||
var diff = s[i].length - rlen;
|
||||
var abs_diff = diff > 0? diff : -diff;
|
||||
min_abs_diff = min_abs_diff < abs_diff? min_abs_diff : abs_diff;
|
||||
if (max_abs_diff < abs_diff)
|
||||
max_abs_diff = abs_diff, alen = diff;
|
||||
}
|
||||
return [alen, min_abs_diff, max_abs_diff];
|
||||
}
|
||||
|
||||
function paf_sveval(args)
|
||||
{
|
||||
var c, min_flt = 30, min_size = 50, max_size = 100000, win_size = 500, print_err = false, print_match = false, bed_fn = null;
|
||||
var len_diff_ratio = 0.5;
|
||||
while ((c = getopt(args, "f:i:x:w:er:pd:")) != null) {
|
||||
if (c == 'f') min_flt = paf_parseNum(getopt.arg);
|
||||
else if (c == 'i') min_size = paf_parseNum(getopt.arg);
|
||||
else if (c == 'x') max_size = paf_parseNum(getopt.arg);
|
||||
else if (c == 'w') win_size = paf_parseNum(getopt.arg);
|
||||
else if (c == 'd') len_diff_ratio = parseFloat(getopt.arg);
|
||||
else if (c == 'r') bed_fn = getopt.arg;
|
||||
else if (c == 'e') print_err = true;
|
||||
else if (c == 'p') print_match = true;
|
||||
}
|
||||
if (args.length - getopt.ind < 2) {
|
||||
print("Usage: paftools.js sveval [options] <base.vcf> <call.vcf>");
|
||||
print("Options:");
|
||||
print(" -r FILE confident region in BED []");
|
||||
print(" -f INT min length to discard [" + min_flt + "]");
|
||||
print(" -i INT min SV length [" + min_size + "]");
|
||||
print(" -x INT max SV length [" + max_size + "]");
|
||||
print(" -w INT fuzzy windown size [" + win_size + "]");
|
||||
print(" -d FLOAT max allele diff if there is a single allele in the window [" + len_diff_ratio + "]");
|
||||
print(" -e print errors");
|
||||
return;
|
||||
}
|
||||
|
||||
function read_bed(fn) {
|
||||
var buf = new Bytes();
|
||||
var file = new File(fn);
|
||||
var bed = {};
|
||||
while (file.readline(buf) >= 0) {
|
||||
var t = buf.toString().split("\t");
|
||||
if (bed[t[0]] == null) bed[t[0]] = [];
|
||||
bed[t[0]].push([parseInt(t[1]), parseInt(t[2])]);
|
||||
}
|
||||
file.close();
|
||||
buf.destroy();
|
||||
for (var x in bed) {
|
||||
Interval.sort(bed[x]);
|
||||
Interval.merge(bed[x]);
|
||||
Interval.index_end(bed[x]);
|
||||
}
|
||||
return bed;
|
||||
}
|
||||
|
||||
var bed = bed_fn != null? read_bed(bed_fn) : null;
|
||||
|
||||
function read_vcf(fn, bed) {
|
||||
var buf = new Bytes();
|
||||
var file = new File(fn);
|
||||
var v = {};
|
||||
while (file.readline(buf) >= 0) {
|
||||
var m, t = buf.toString().split("\t");
|
||||
if (t[0][0] == '#') continue;
|
||||
if (bed != null && bed[t[0]] == null) continue;
|
||||
if (t[4] == '<INV>' || t[4] == '<INVDUP>') continue; // no inversion
|
||||
if (/[\[\]]/.test(t[4])) continue; // no break points
|
||||
var st = parseInt(t[1]) - 1, en = st + t[3].length;
|
||||
// parse svlen
|
||||
var b = _paf_get_alen(t), svlen = b[0];
|
||||
var abslen = svlen == null? 0 : svlen > 0? svlen : -svlen;
|
||||
if (abslen < min_flt || abslen > max_size) continue;
|
||||
// update end
|
||||
if ((m = /(^|;)END=(\d+)/.exec(t[7])) != null)
|
||||
en = parseInt(m[2]);
|
||||
else if (svlen != null && svlen < 0)
|
||||
en = st + (-svlen);
|
||||
if (en < st) en = st;
|
||||
if (st == en) --st, ++en;
|
||||
if (bed != null && Interval.find_ovlp(bed[t[0]], st, en).length == 0) continue;
|
||||
// insert
|
||||
if (v[t[0]] == null) v[t[0]] = [];
|
||||
v[t[0]].push([st, en, svlen, abslen]);
|
||||
}
|
||||
file.close();
|
||||
buf.destroy();
|
||||
for (var x in v) {
|
||||
Interval.sort(v[x]);
|
||||
Interval.index_end(v[x]);
|
||||
}
|
||||
return v;
|
||||
}
|
||||
|
||||
function compare_vcf(v0, v1, label) {
|
||||
var m = 0, n = 0;
|
||||
for (var x in v1) {
|
||||
var a1 = v1[x], a0 = v0[x];
|
||||
for (var i = 0; i < a1.length; ++i) {
|
||||
if (a1[i][3] < min_size) continue;
|
||||
++n;
|
||||
if (a0 == null) continue;
|
||||
var ws = win_size + (a1[i][3]>>1);
|
||||
var st = a1[i][0] > ws? a1[i][0] - ws : 0;
|
||||
b = Interval.find_ovlp(a0, st, a1[i][1] + ws);
|
||||
var n_ins = 0, n_del = 0, sv_del = null, sv_ins = null;
|
||||
for (var j = 0; j < b.length; ++j) {
|
||||
if (b[j][2] < 0) ++n_del, sv_del = -b[j][2];
|
||||
else if (b[j][2] > 0) ++n_ins, sv_ins = b[j][2];
|
||||
if (print_match)
|
||||
print("MA", x, a1[i].slice(0, 3).join("\t"), b[j].slice(0, 3).join("\t"));
|
||||
}
|
||||
var match = false;
|
||||
if (a1[i][2] > 0) { // insertion
|
||||
if (n_ins == 1) {
|
||||
var diff = sv_ins - a1[i][3];
|
||||
if (diff < 0) diff = -diff;
|
||||
if (diff < min_size || diff / a1[i][3] < len_diff_ratio)
|
||||
match = true;
|
||||
} else if (n_ins > 1) match = true; // multiple insertions; ambiguous
|
||||
} else if (a1[i][2] < 0) {
|
||||
if (n_del == 1) { // deletion
|
||||
var diff = sv_del - a1[i][3];
|
||||
if (diff < 0) diff = -diff;
|
||||
if (diff < min_size || diff / a1[i][3] < len_diff_ratio)
|
||||
match = true;
|
||||
} else if (n_del > 1) match = true; // multiple deletions; ambiguous
|
||||
}
|
||||
if (match) ++m;
|
||||
else if (print_err) {
|
||||
if ((a1[i][2] > 0 && n_ins > 0) || (a1[i][2] < 0 && n_del > 0))
|
||||
print("MM", x, a1[i].slice(0, 3).join("\t"));
|
||||
print(label, x, a1[i].slice(0, 3).join("\t"));
|
||||
}
|
||||
}
|
||||
}
|
||||
return [n, m];
|
||||
}
|
||||
|
||||
var v_base = read_vcf(args[getopt.ind+0], bed);
|
||||
var v_call = read_vcf(args[getopt.ind+1], bed);
|
||||
var fn = compare_vcf(v_call, v_base, 'FN');
|
||||
var fp = compare_vcf(v_base, v_call, 'FP');
|
||||
print('SN', fn[0], fn[1], (fn[1] / fn[0]).toFixed(6));
|
||||
print('PC', fp[0], fp[1], (fp[1] / fp[0]).toFixed(6));
|
||||
print('F1', ((fn[1] / fn[0] + fp[1] / fp[0]) / 2).toFixed(6));
|
||||
}
|
||||
|
||||
function paf_vcfsel(args)
|
||||
{
|
||||
var c, min_l = 0, max_l = 1<<30;
|
||||
while ((c = getopt(args, "l:L:")) != null) {
|
||||
if (c == 'l') min_l = parseInt(getopt.arg);
|
||||
else if (c == 'L') max_l = parseInt(getopt.arg);
|
||||
}
|
||||
|
||||
var buf = new Bytes();
|
||||
if (getopt.ind == args.length) {
|
||||
print("Usage: paftools.js vcfsel [options] <in.vcf>");
|
||||
return 1;
|
||||
}
|
||||
var file = args[getopt.ind] == "-"? new File() : new File(args[getopt.ind]);
|
||||
while (file.readline(buf) >= 0) {
|
||||
var m, line = buf.toString();
|
||||
if (line[0] == '#') {
|
||||
print(line);
|
||||
continue;
|
||||
}
|
||||
var t = line.split("\t");
|
||||
var st = parseInt(t[1]), en = st + t[3].length - 1;
|
||||
if ((m = /(^|;)END=(\d+)/.exec(t[7])) != null)
|
||||
en = parseInt(m[2]);
|
||||
if (en < st) {
|
||||
warn("END is smaller than POS: " + en + " < " + st);
|
||||
en = st;
|
||||
}
|
||||
var b = _paf_get_alen(t);
|
||||
var alen = b[0], min_abs_diff = b[1], max_abs_diff = b[2];
|
||||
if (max_abs_diff < min_l || min_abs_diff > max_l)
|
||||
continue;
|
||||
print(line);
|
||||
}
|
||||
file.close();
|
||||
buf.destroy();
|
||||
}
|
||||
|
||||
/*************************
|
||||
***** main function *****
|
||||
*************************/
|
||||
@@ -2529,10 +2947,13 @@ function main(args)
|
||||
print("");
|
||||
print(" stat collect basic mapping information in PAF/SAM");
|
||||
print(" asmstat collect basic assembly information");
|
||||
print(" asmgene evaluate gene completeness (EXPERIMENTAL)");
|
||||
print(" asmgene evaluate gene completeness");
|
||||
print(" misjoin evaluate large-scale misjoins");
|
||||
print(" liftover simplistic liftOver");
|
||||
print(" call call variants from asm-to-ref alignment with the cs tag");
|
||||
print(" bedcov compute the number of bases covered");
|
||||
print(" vcfstat VCF statistics");
|
||||
print(" sveval compare two SV callsets in VCF");
|
||||
print(" version print paftools.js version");
|
||||
print("");
|
||||
print(" mapeval evaluate mapping accuracy using mason2/PBSIM-simulated FASTQ");
|
||||
@@ -2552,6 +2973,7 @@ function main(args)
|
||||
else if (cmd == 'stat') paf_stat(args);
|
||||
else if (cmd == 'asmstat') paf_asmstat(args);
|
||||
else if (cmd == 'asmgene') paf_asmgene(args);
|
||||
else if (cmd == 'misjoin') paf_misjoin(args);
|
||||
else if (cmd == 'liftover' || cmd == 'liftOver') paf_liftover(args);
|
||||
else if (cmd == 'vcfpair') paf_vcfpair(args);
|
||||
else if (cmd == 'call') paf_call(args);
|
||||
@@ -2561,6 +2983,9 @@ function main(args)
|
||||
else if (cmd == 'pbsim2fq') paf_pbsim2fq(args);
|
||||
else if (cmd == 'junceval') paf_junceval(args);
|
||||
else if (cmd == 'ov-eval') paf_ov_eval(args);
|
||||
else if (cmd == 'vcfstat') paf_vcfstat(args);
|
||||
else if (cmd == 'sveval') paf_sveval(args);
|
||||
else if (cmd == 'vcfsel') paf_vcfsel(args);
|
||||
else if (cmd == 'version') print(paftools_version);
|
||||
else throw Error("unrecognized command: " + cmd);
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user