diff --git a/bseq.h b/bseq.h index 5dfe095..015dca3 100644 --- a/bseq.h +++ b/bseq.h @@ -26,15 +26,19 @@ int mm_bseq_eof(mm_bseq_file_t *fp); extern unsigned char seq_nt4_table[256]; extern unsigned char seq_comp_table[256]; +static inline int mm_qname_len(const char *s) +{ + int l; + l = strlen(s); + return l >= 3 && s[l-1] >= '0' && s[l-1] <= '9' && s[l-2] == '/'? l - 2 : l; +} + static inline int mm_qname_same(const char *s1, const char *s2) { int l1, l2; - l1 = strlen(s1); - l2 = strlen(s2); - if (l1 != l2 || l1 < 3) return 0; - if (!(s1[l1-1] >= '0' && s1[l1-1] <= '9' && s1[l1-2] == '/')) return 0; - if (!(s2[l2-1] >= '0' && s2[l2-1] <= '9' && s2[l2-2] == '/')) return 0; - return (strncmp(s1, s2, l1 - 2) == 0); + l1 = mm_qname_len(s1); + l2 = mm_qname_len(s2); + return (l1 == l2 && strncmp(s1, s2, l1) == 0); } static inline void mm_revcomp_bseq(mm_bseq1_t *s) diff --git a/format.c b/format.c index 989e4ff..201b1c4 100644 --- a/format.c +++ b/format.c @@ -240,21 +240,30 @@ 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, 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_regss, const mm_reg1_t *const* regss, void *km) { - int flag = 0; + int flag = n_seg > 1? 0x1 : 0x0, n_regs = n_regss[seg_idx]; + 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; + + // write QNAME s->l = 0; + 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 if (r == 0) { - mm_sprintf_lite(s, "%s\t4\t*\t0\t0\t*\t*\t0\t0\t", t->name); - sam_write_sq(s, t->seq, t->l_seq, 0, 0); - mm_sprintf_lite(s, "\t"); - if (t->qual) sam_write_sq(s, t->qual, t->l_seq, 0, 0); - else mm_sprintf_lite(s, "*"); + mm_sprintf_lite(s, "\t4\t*\t0\t0\t*\t"); } else { if (r->rev) flag |= 0x10; 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 (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); 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'; @@ -264,7 +273,20 @@ void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const m clip_len = r->rev? r->qs : t->l_seq - r->qe; if (clip_len) mm_sprintf_lite(s, "%d%c", clip_len, clip_char); } else mm_sprintf_lite(s, "*"); + } + + // write mate positions + if (n_seg > 1) { mm_sprintf_lite(s, "\t*\t0\t0\t"); + } else mm_sprintf_lite(s, "\t*\t0\t0\t"); + + // write SEQ and QUAL + if (r == 0) { + sam_write_sq(s, t->seq, t->l_seq, 0, 0); + mm_sprintf_lite(s, "\t"); + if (t->qual) sam_write_sq(s, t->qual, t->l_seq, 0, 0); + else mm_sprintf_lite(s, "*"); + } else { if ((flag & 0x900) == 0) { sam_write_sq(s, t->seq, t->l_seq, r->rev, r->rev); mm_sprintf_lite(s, "\t"); @@ -278,8 +300,13 @@ void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const m if (t->qual) sam_write_sq(s, t->qual + r->qs, r->qe - r->qs, r->rev, 0); else mm_sprintf_lite(s, "*"); } + } + + // write tags + if (mm_rg_id[0]) mm_sprintf_lite(s, "\tRG:Z:%s", mm_rg_id); + if (n_seg > 2) mm_sprintf_lite(s, "\tFI:i:%d", seg_idx); + if (r) { write_tags(s, r); - if (mm_rg_id[0]) mm_sprintf_lite(s, "\tRG:Z:%s", mm_rg_id); 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) @@ -306,5 +333,14 @@ void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const m } } } + s->s[s->l] = 0; // we always have room for an extra byte (see str_enlarge) } + +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 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); +} diff --git a/main.c b/main.c index a1949b3..aefa7f4 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "getopt.h" -#define MM_VERSION "2.2-r468-dirty" +#define MM_VERSION "2.2-r469-dirty" #ifdef __linux__ #include diff --git a/map.c b/map.c index 8b5a9aa..0b85473 100644 --- a/map.c +++ b/map.c @@ -410,7 +410,7 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback static void *worker_pipeline(void *shared, int step, void *in) { - int i, j; + int i, j, k; pipeline_t *p = (pipeline_t*)shared; if (step == 0) { // step 0: read sequences int with_qual = (!!(p->opt->flag & MM_F_OUT_SAM) && !(p->opt->flag & MM_F_NO_QUAL)); @@ -448,29 +448,34 @@ static void *worker_pipeline(void *shared, int step, void *in) for (i = 0; i < p->n_threads; ++i) mm_tbuf_destroy(s->buf[i]); free(s->buf); if ((p->opt->flag & MM_F_OUT_CS) && !(mm_dbg_flag & MM_DBG_NO_KALLOC)) km = km_init(); - for (i = 0; i < s->n_seq; ++i) { - mm_bseq1_t *t = &s->seq[i]; - for (j = 0; j < s->n_reg[i]; ++j) { - mm_reg1_t *r = &s->reg[i][j]; - assert(!r->sam_pri || r->id == r->parent); - if ((p->opt->flag & MM_F_NO_PRINT_2ND) && r->id != r->parent) - continue; - 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, km, p->opt->flag); - puts(p->str.s); + for (k = 0; k < s->n_frag; ++k) { + int seg_st = s->seg_off[k], seg_en = s->seg_off[k] + s->n_seg[k]; + for (i = seg_st; i < seg_en; ++i) { + mm_bseq1_t *t = &s->seq[i]; + for (j = 0; j < s->n_reg[i]; ++j) { + mm_reg1_t *r = &s->reg[i][j]; + assert(!r->sam_pri || r->id == r->parent); + 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); + 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); + puts(p->str.s); + } } - if (s->n_reg[i] == 0 && (p->opt->flag & MM_F_OUT_SAM)) { - mm_write_sam(&p->str, 0, t, 0, 0, 0); - puts(p->str.s); + for (i = seg_st; i < seg_en; ++i) { + 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); } - 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); } - free(s->reg); free(s->n_reg); free(s->seq); + free(s->reg); free(s->n_reg); free(s->seq); // seg_off and n_seg were allocated with reg; no memory leak here km_destroy(km); if (mm_verbose >= 3) fprintf(stderr, "[M::%s::%.3f*%.2f] mapped %d sequences\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), s->n_seq); diff --git a/mmpriv.h b/mmpriv.h index a48a480..cbaef08 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -57,6 +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_idxopt_init(mm_idxopt_t *opt); const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n);