From 2179e9e24bd6d05038da0e50690de549744a5c17 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Fri, 28 Jul 2017 23:08:39 -0400 Subject: [PATCH] r221: output SA in the SAM output --- format.c | 29 +++++++++++++++++++++++++++-- main.c | 2 +- map.c | 6 +++--- mmpriv.h | 2 +- 4 files changed, 32 insertions(+), 7 deletions(-) diff --git a/format.c b/format.c index c927adb..5098aac 100644 --- a/format.c +++ b/format.c @@ -105,7 +105,7 @@ 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_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r) +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) { int flag = 0; s->l = 0; @@ -120,7 +120,7 @@ void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const m if (r->parent != r->id) flag |= 0x100; else if (!r->sam_pri) flag |= 0x800; mm_sprintf_lite(s, "%s\t%d\t%s\t%d\t%d\t", t->name, flag, mi->seq[r->rid].name, r->rs+1, r->mapq); - if (r->p) { // TODO: using hard clippings + 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'; if (clip_len) mm_sprintf_lite(s, "%d%c", clip_len, clip_char); @@ -144,6 +144,31 @@ void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const m else mm_sprintf_lite(s, "*"); } write_tags(s, r); + if (r->parent == r->id && r->p && n_regs > 1 && regs && r >= regs && r - regs < n_regs) { // supplementary aln may exist + int i, n_sa = 0; // n_sa: number of SA fields + for (i = 0; i < n_regs; ++i) + if (i != r - regs && regs[i].parent == regs[i].id && regs[i].p) + ++n_sa; + if (n_sa > 0) { + mm_sprintf_lite(s, "\tSA:Z:"); + for (i = 0; i < n_regs; ++i) { + const mm_reg1_t *q = ®s[i]; + int l_M, l_I = 0, l_D = 0, clip5 = 0, clip3 = 0; + if (r == q || q->parent != q->id || q->p == 0) continue; + if (q->qe - q->qs < q->re - q->rs) l_M = q->qe - q->qs, l_D = (q->re - q->rs) - l_M; + else l_M = q->re - q->rs, l_I = (q->qe - q->qs) - l_M; + clip5 = q->rev? t->l_seq - q->qe : q->qs; + clip3 = q->rev? q->qs : t->l_seq - q->qe; + mm_sprintf_lite(s, "%s,%d,%c,", mi->seq[q->rid].name, q->rs+1, "+-"[q->rev]); + if (clip5) mm_sprintf_lite(s, "%dS", clip5); + if (l_M) mm_sprintf_lite(s, "%dM", l_M); + if (l_I) mm_sprintf_lite(s, "%dI", l_I); + if (l_D) mm_sprintf_lite(s, "%dD", l_D); + if (clip3) mm_sprintf_lite(s, "%dS", clip3); + mm_sprintf_lite(s, ",%d,%d;", q->mapq, q->p->n_diff); + } + } + } } s->s[s->l] = 0; // we always have room for an extra byte (see str_enlarge) } diff --git a/main.c b/main.c index 9f69faf..f7397ae 100644 --- a/main.c +++ b/main.c @@ -8,7 +8,7 @@ #include "minimap.h" #include "mmpriv.h" -#define MM_VERSION "2.0-r219-dirty" +#define MM_VERSION "2.0-r221-dirty" void liftrlimit() { diff --git a/map.c b/map.c index b3f3924..e9d894e 100644 --- a/map.c +++ b/map.c @@ -346,15 +346,15 @@ static void *worker_pipeline(void *shared, int step, void *in) mm_bseq1_t *t = &s->seq[i]; for (j = 0; j < s->n_reg[i]; ++j) { mm_reg1_t *r = &s->reg[i][j]; - if (p->opt->flag & MM_F_OUT_SAM) mm_write_sam(&p->str, mi, t, r); + if (p->opt->flag & MM_F_OUT_SAM) mm_write_sam(&p->str, mi, t, r, s->n_reg[i], s->reg[i]); else mm_write_paf(&p->str, mi, t, r); puts(p->str.s); - free(r->p); } if (s->n_reg[i] == 0 && (p->opt->flag & MM_F_OUT_SAM)) { - mm_write_sam(&p->str, 0, t, 0); + mm_write_sam(&p->str, 0, t, 0, 0, 0); puts(p->str.s); } + for (j = 0; j < s->n_reg[i]; ++j) free(s->reg[i][j].p); free(s->reg[i]); free(s->seq[i].seq); free(s->seq[i].name); if (s->seq[i].qual) free(s->seq[i].qual); diff --git a/mmpriv.h b/mmpriv.h index dcf49ea..e8a4156 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -40,7 +40,7 @@ void radix_sort_64(uint64_t *beg, uint64_t *end); uint32_t ks_ksmall_uint32_t(size_t n, uint32_t arr[], size_t kk); void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r); -void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r); +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); int mm_chain_dp(int max_dist, int bw, int max_skip, int min_cnt, int min_sc, int64_t n, mm128_t *a, uint64_t **_u, void *km); mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a);