diff --git a/bseq.c b/bseq.c index 3179fa1..16ead16 100644 --- a/bseq.c +++ b/bseq.c @@ -48,6 +48,7 @@ bseq1_t *bseq_read(bseq_file_t *fp, int chunk_size, int *n_) s = &seqs[n]; s->name = strdup(ks->name.s); s->seq = strdup(ks->seq.s); + s->qual = ks->qual.l? strdup(ks->qual.s) : 0; s->l_seq = ks->seq.l; size += seqs[n++].l_seq; if (size >= chunk_size) break; diff --git a/bseq.h b/bseq.h index 0595d98..bb06766 100644 --- a/bseq.h +++ b/bseq.h @@ -8,7 +8,7 @@ typedef struct bseq_file_s bseq_file_t; typedef struct { int l_seq, rid; - char *name, *seq; + char *name, *seq, *qual; } bseq1_t; bseq_file_t *bseq_open(const char *fn); diff --git a/format.c b/format.c index 092aaf4..7690baf 100644 --- a/format.c +++ b/format.c @@ -110,23 +110,40 @@ void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const bseq1_t *t, const mm_r { int flag = 0; s->l = 0; - 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 (r->p) { // TODO: using hard clippings - 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); - for (k = 0; k < r->p->n_cigar; ++k) - mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MID"[r->p->cigar[k]&0xf]); - 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, "*"); - mm_sprintf_lite(s, "\t*\t0\t0\t"); - if ((flag & 0x900) == 0) sam_write_sq(s, t->seq, t->l_seq, r->rev, r->rev); - else if (flag & 0x100) mm_sprintf_lite(s, "*"); - else sam_write_sq(s, t->seq + r->qs, r->qe - r->qs, r->rev, r->rev); - mm_sprintf_lite(s, "\t*"); // quality - write_tags(s, r); + if (r == 0) { + mm_sprintf_lite(s, "%s\t0\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, "*"); + } 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 (r->p) { // TODO: using hard clippings + 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); + for (k = 0; k < r->p->n_cigar; ++k) + mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MID"[r->p->cigar[k]&0xf]); + 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, "*"); + mm_sprintf_lite(s, "\t*\t0\t0\t"); + if ((flag & 0x900) == 0) { + sam_write_sq(s, t->seq, t->l_seq, r->rev, r->rev); + mm_sprintf_lite(s, "\t"); + if (t->qual) sam_write_sq(s, t->qual, t->l_seq, r->rev, 0); + else mm_sprintf_lite(s, "*"); + } else if (flag & 0x100) { + mm_sprintf_lite(s, "*\t*"); + } else { + sam_write_sq(s, t->seq + r->qs, r->qe - r->qs, r->rev, r->rev); + mm_sprintf_lite(s, "\t"); + 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(s, r); + } } diff --git a/main.c b/main.c index 24b9e34..d24c952 100644 --- a/main.c +++ b/main.c @@ -10,7 +10,7 @@ #include "minimap.h" #include "mmpriv.h" -#define MM_VERSION "2.0-r120-pre" +#define MM_VERSION "2.0-r121-pre" void liftrlimit() { diff --git a/map.c b/map.c index c725ba8..1ec5fec 100644 --- a/map.c +++ b/map.c @@ -321,12 +321,15 @@ static void *worker_pipeline(void *shared, int step, void *in) bseq1_t *t = &s->seq[i]; for (j = 0; j < s->n_reg[i]; ++j) { mm_reg1_t *r = &s->reg[i][j]; - if (r->cnt == 0) continue; if (p->opt->flag & MM_F_OUT_SAM) mm_write_sam(&p->str, mi, t, r); 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); + puts(p->str.s); + } free(s->reg[i]); free(s->seq[i].seq); free(s->seq[i].name); }