From 5ab99eb26e667ca71b6ea86188823bf17a3d4460 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 5 Oct 2017 10:59:38 -0400 Subject: [PATCH] more accurate SAM flag --- format.c | 34 +++++++++++++++++++++++----------- map.c | 4 ++-- minimap.h | 5 +++-- mmpriv.h | 2 +- pe.c | 1 + 5 files changed, 30 insertions(+), 16 deletions(-) diff --git a/format.c b/format.c index 201b1c4..5b293d8 100644 --- a/format.c +++ b/format.c @@ -240,9 +240,10 @@ static void sam_write_sq(kstring_t *s, char *seq, int l, int rev, int comp) } else str_copy(s, seq, seq + l); } -void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km) +void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int opt_flag) { - int flag = n_seg > 1? 0x1 : 0x0, n_regs = n_regss[seg_idx]; + int flag, n_regs = n_regss[seg_idx]; + int next_sid = n_seg > 1? (seg_idx + 1) % n_seg : -1; const mm_reg1_t *regs = regss[seg_idx]; const mm_reg1_t *r = n_regs > 0 && reg_idx < n_regs && reg_idx >= 0? ®s[reg_idx] : NULL; @@ -251,19 +252,28 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se mm_sprintf_lite(s, "%s", t->name); if (n_seg > 1) s->l = mm_qname_len(t->name); // trim the suffix like /1 or /2 - // write up to CIGAR + // write flag + flag = n_seg > 1? 0x1 : 0x0; + if (n_seg > 1) { + if (r && r->proper_frag) flag |= 0x2; // FIXME: this doesn't work when there are more than 2 segments + if (seg_idx == 0) flag |= 0x40; + else if (seg_idx == n_seg - 1) flag |= 0x80; + if (n_regss[next_sid] == 0) flag |= 0x8; + } if (r == 0) { - mm_sprintf_lite(s, "\t4\t*\t0\t0\t*\t"); + flag |= 0x4; } else { if (r->rev) flag |= 0x10; if (r->parent != r->id) flag |= 0x100; else if (!r->sam_pri) flag |= 0x800; - if (n_seg > 1) { - if (seg_idx == 0) flag |= 0x40; - else if (seg_idx == n_seg - 1) flag |= 0x80; - // TODO: set 0x2!!! - } - mm_sprintf_lite(s, "\t%d\t%s\t%d\t%d\t", flag, mi->seq[r->rid].name, r->rs+1, r->mapq); + } + mm_sprintf_lite(s, "\t%d", flag); + + // write up to CIGAR + if (r == 0) { + mm_sprintf_lite(s, "\t*\t0\t0\t*\t"); + } else { + mm_sprintf_lite(s, "\t%s\t%d\t%d\t", mi->seq[r->rid].name, r->rs+1, r->mapq); if (r->p) { // actually this should always be true for SAM output uint32_t k, clip_len = r->rev? t->l_seq - r->qe : r->qs; int clip_char = (flag&0x800)? 'H' : 'S'; @@ -332,6 +342,8 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se } } } + if (r->p && (opt_flag & MM_F_OUT_CS)) + write_cs(km, s, mi, t, r, MM_F_CS_NO_EQUAL); } s->s[s->l] = 0; // we always have room for an extra byte (see str_enlarge) @@ -342,5 +354,5 @@ void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const m int i; for (i = 0; i < n_regs; ++i) if (r == ®s[i]) break; - mm_write_sam2(s, mi, t, 0, i, 1, &n_regs, ®s, NULL); + mm_write_sam2(s, mi, t, 0, i, 1, &n_regs, ®s, NULL, 0); } diff --git a/map.c b/map.c index 0b85473..a8a1ed0 100644 --- a/map.c +++ b/map.c @@ -458,13 +458,13 @@ static void *worker_pipeline(void *shared, int step, void *in) if ((p->opt->flag & MM_F_NO_PRINT_2ND) && r->id != r->parent) continue; if (p->opt->flag & MM_F_OUT_SAM) - mm_write_sam2(&p->str, mi, t, i - seg_st, j, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km); + mm_write_sam2(&p->str, mi, t, i - seg_st, j, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag); else mm_write_paf(&p->str, mi, t, r, km, p->opt->flag); puts(p->str.s); } if (s->n_reg[i] == 0 && (p->opt->flag & MM_F_OUT_SAM)) { - mm_write_sam2(&p->str, mi, t, i - seg_st, -1, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km); + mm_write_sam2(&p->str, mi, t, i - seg_st, -1, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag); puts(p->str.s); } } diff --git a/minimap.h b/minimap.h index c0a645b..2d4b281 100644 --- a/minimap.h +++ b/minimap.h @@ -69,8 +69,9 @@ typedef struct { int32_t parent, subsc; // parent==id if primary; best alternate mapping score int32_t as; // offset in the a[] array (for internal uses only) int32_t fuzzy_mlen, fuzzy_blen; // seeded exact match length; seeded alignment block length (approximate) - uint32_t mapq:8, split:2, sam_pri:1, n_sub:21; // mapQ; split pattern; if SAM primary; number of suboptimal mappings - uint32_t hash, dummy; + uint32_t mapq:8, split:2, n_sub:22; // mapQ; split pattern; number of suboptimal mappings + uint32_t sam_pri:1, proper_frag:1, dummy:30; + uint32_t hash; mm_extra_t *p; } mm_reg1_t; diff --git a/mmpriv.h b/mmpriv.h index cbaef08..57e2648 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -57,7 +57,7 @@ void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, i void mm_write_sam_hdr(const mm_idx_t *mi, const char *rg, const char *ver, int argc, char *argv[]); void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag); void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int n_regs, const mm_reg1_t *regs); -void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regs, const mm_reg1_t *const* regs, void *km); +void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regs, const mm_reg1_t *const* regs, void *km, int opt_flag); void mm_idxopt_init(mm_idxopt_t *opt); const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n); diff --git a/pe.c b/pe.c index ffef9ec..70affcc 100644 --- a/pe.c +++ b/pe.c @@ -115,6 +115,7 @@ void mm_pair(void *km, int max_gap_ref, int pe_bonus, int sub_diff, int match_sc int n_sub = 0, mapq_pe; mm_reg1_t *r[2]; r[0] = a[max_idx[0]].r, r[1] = a[max_idx[1]].r; + r[0]->proper_frag = r[1]->proper_frag = 1; for (s = 0; s < 2; ++s) { if (r[s]->id != r[s]->parent) { mm_reg1_t *p = ®s[s][r[s]->parent];