From dc61301d9fc866e9e179b78980df656233559671 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 1 Mar 2018 21:58:12 -0500 Subject: [PATCH] sam2paf cleanup --- misc/paftools.js | 85 +++++++++++++++++++++++++++++++----------------- 1 file changed, 56 insertions(+), 29 deletions(-) diff --git a/misc/paftools.js b/misc/paftools.js index 59bfc9b..a61e47e 100755 --- a/misc/paftools.js +++ b/misc/paftools.js @@ -1031,8 +1031,12 @@ function paf_sam2paf(args) var c, pri_only = false; while ((c = getopt(args, "p")) != null) if (c == 'p') pri_only = true; + if (args.length == getopt.ind) { + print("Usage: paftools.js sam2paf [-p] "); + exit(1); + } - var file = args.length == getopt.ind || args[getopt.ind] == "-"? new File() : new File(args[getopt.ind]); + var file = args[getopt.ind] == "-"? new File() : new File(args[getopt.ind]); var buf = new Bytes(); var re = /(\d+)([MIDSHNX=])/g; @@ -1050,31 +1054,40 @@ function paf_sam2paf(args) } var t = line.split("\t"); var flag = parseInt(t[1]); - if (t[9] != '*' && t[10] != '*' && t[9].length != t[10].length) throw Error("ERROR at line " + lineno + ": inconsistent SEQ and QUAL lengths - " + t[9].length + " != " + t[10].length); - if (t[2] == '*' || (flag&4)) continue; + if (t[9] != '*' && t[10] != '*' && t[9].length != t[10].length) + throw Error("ERROR at line " + lineno + ": inconsistent SEQ and QUAL lengths - " + t[9].length + " != " + t[10].length); + if (t[2] == '*' || (flag&4) || t[5] == '*') continue; if (pri_only && (flag&0x100)) continue; var tlen = len[t[2]]; if (tlen == null) throw Error("ERROR at line " + lineno + ": can't find the length of contig " + t[2]); - var nn = (m = /\tnn:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : 0; - var NM = (m = /\tNM:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : null; - var have_NM = NM == null? false : true; - NM += nn; - var clip = [0, 0], I = [0, 0], D = [0, 0], M = 0, N = 0, ql = 0, tl = 0, mm = 0, ext_cigar = false; + // find tags + var nn = 0, NM = null, MD = null; + for (var i = 11; i < t.length; ++i) { + if (t[i].indexOf("NM:i:") == 0) NM = parseInt(t[i].substr(5)); + else if (t[i].indexOf("nn:i:") == 0) nn = parseInt(t[i].substr(5)); + else if (t[i].indexOf("MD:Z:") == 0) MD = t[i].substr(5); + } + // infer various lengths from CIGAR + var clip = [0, 0], soft_clip = 0, I = [0, 0], D = [0, 0], M = 0, N = 0, mm = 0, have_M = false, have_ext = false; while ((m = re.exec(t[5])) != null) { var l = parseInt(m[1]); - if (m[2] == 'M') M += l, ql += l, tl += l, ext_cigar = false; - else if (m[2] == 'I') ++I[0], I[1] += l, ql += l; - else if (m[2] == 'D') ++D[0], D[1] += l, tl += l; - else if (m[2] == 'N') N += l, tl += l; - else if (m[2] == 'S') clip[M == 0? 0 : 1] = l, ql += l; - else if (m[2] == 'H') clip[M == 0? 0 : 1] = l; - else if (m[2] == '=') M += l, ql += l, tl += l, ext_cigar = true; - else if (m[2] == 'X') M += l, ql += l, tl += l, mm += l, ext_cigar = true; + if (m[2] == 'M') M += l, have_M = true; + else if (m[2] == 'I') ++I[0], I[1] += l; + else if (m[2] == 'D') ++D[0], D[1] += l; + else if (m[2] == 'N') N += l; + else if (m[2] == 'S') clip[n_cigar == 0? 0 : 1] = l, soft_clip += l; + else if (m[2] == 'H') clip[n_cigar == 0? 0 : 1] = l; + else if (m[2] == '=') M += l, have_ext = true; + else if (m[2] == 'X') M += l, mm += l, have_ext = true; ++n_cigar; } + var ql = M + I[1] + soft_clip; + var tl = M + D[1] + N; + var ts = parseInt(t[3]) - 1, te = ts + tl; + // checking coordinate and length consistencies if (n_cigar > 65535) warn("WARNING at line " + lineno + ": " + n_cigar + " CIGAR operations"); - if (tl + parseInt(t[3]) - 1 > tlen) { + if (te > tlen) { warn("WARNING at line " + lineno + ": alignment end position larger than ref length; skipped"); continue; } @@ -1082,24 +1095,38 @@ function paf_sam2paf(args) warn("WARNING at line " + lineno + ": SEQ length inconsistent with CIGAR (" + t[9].length + " != " + ql + "); skipped"); continue; } - if (!have_NM || ext_cigar) NM = I[1] + D[1] + mm; - if (NM < I[1] + D[1] + mm) { - warn("WARNING at line " + lineno + ": NM is less than the total number of gaps (" + NM + " < " + (I[1]+D[1]+mm) + ")"); + // compute matching length, block length and calibrate NM + if (have_ext && !have_M) { // extended CIGAR + if (NM != null && NM != I[1] + D[1] + mm) + warn("WARNING at line " + lineno + ": NM is different from sum of gaps and mismatches"); NM = I[1] + D[1] + mm; + } else if (NM != null) { // standard CIGAR; NM present + if (NM < I[1] + D[1]) { + warn("WARNING at line " + lineno + ": NM is less than the total number of gaps (" + NM + " < " + (I[1]+D[1]) + ")"); + NM = I[1] + D[1]; + } + mm = NM - (I[1] + D[1]); + } else { // no way to compute mm + warn("WARNING at line " + lineno + ": unable to find the number of mismatches; assuming zero"); + mm = 0; } - var extra = ["mm:i:"+(NM-I[1]-D[1]), "io:i:"+I[0], "in:i:"+I[1], "do:i:"+D[0], "dn:i:"+D[1]]; - var match = M - (NM - I[1] - D[1]); + var mlen = M - mm; var blen = M + I[1] + D[1]; + // find query name, start and end var qlen = M + I[1] + clip[0] + clip[1]; - var qs, qe; - if (flag&16) qs = clip[1], qe = qlen - clip[0]; - else qs = clip[0], qe = qlen - clip[1]; - var ts = parseInt(t[3]) - 1, te = ts + M + D[1] + N; - var qname = t[0]; + var qname = t[0], qs, qe; if ((flag&1) && (flag&0x40)) qname += '/1'; if ((flag&1) && (flag&0x80)) qname += '/2'; - var a = [qname, qlen, qs, qe, flag&16? '-' : '+', t[2], tlen, ts, te, match, blen, t[4]]; - print(a.join("\t"), extra.join("\t")); + if (flag&16) qs = clip[1], qe = qlen - clip[0]; + else qs = clip[0], qe = qlen - clip[1]; + // optional tags + var type = flag&0x100? 'S' : 'P'; + var tags = ["tp:A:" + type]; + if (NM != null) tags.push("mm:i:"+mm); + tags.push("gn:i:"+(I[1]+D[1]), "go:i:"+(I[0]+D[0]), "cg:Z:" + t[5].replace(/\d+[SH]/g, '')); + // print out + var a = [qname, qlen, qs, qe, flag&16? '-' : '+', t[2], tlen, ts, te, mlen, blen, t[4]]; + print(a.join("\t"), tags.join("\t")); } buf.destroy();