Skip to content

format SAM/PAF records with inline appenders #5

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
nh13 wants to merge 1 commit into
base: master
Choose a base branch
from
Draft

format SAM/PAF records with inline appenders #5

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
223 changes: 141 additions & 82 deletions format.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,43 +6,114 @@

static char mb_rg_id[256];

/******************
* string helpers *
******************/

static inline void str_enlarge(kstring_t *s, size_t l)
{
if (s->l + l + 1 > s->m) {
s->m = s->l + l + 1;
kom_roundup64(s->m);
s->s = kom_realloc(char, s->s, s->m);
}
}

static inline void str_copy(kstring_t *s, const char *st, const char *en)
{
str_enlarge(s, en - st);
memcpy(&s->s[s->l], st, en - st);
s->l += en - st;
}

// inline equivalents of the kom_sprintf_lite "%c"/"%s"/"%d"/"%u" conversions,
// without the per-field format-string parsing
static inline void str_putc(kstring_t *s, char c)
{
str_enlarge(s, 1);
s->s[s->l++] = c;
}

static inline void str_puts(kstring_t *s, const char *t)
{
str_copy(s, t, t + strlen(t));
}

static inline void str_puti(kstring_t *s, long c) // matches "%d" and "%ld"
{
char buf[24];
int i, l = 0;
unsigned long x = c >= 0? (unsigned long)c : -(unsigned long)c;
do { buf[l++] = x % 10 + '0'; x /= 10; } while (x > 0);
if (c < 0) buf[l++] = '-';
str_enlarge(s, l);
for (i = l - 1; i >= 0; --i) s->s[s->l++] = buf[i];
}

static inline void str_putu(kstring_t *s, unsigned x) // matches "%u"
{
char buf[16];
int i, l = 0;
do { buf[l++] = x % 10 + '0'; x /= 10; } while (x > 0);
str_enlarge(s, l);
for (i = l - 1; i >= 0; --i) s->s[s->l++] = buf[i];
}

/**************
* PAF output *
**************/

static inline void write_tags(kstring_t *s, const mb_hit_t *p)
{
int32_t nm = p->blen - p->mlen + p->p->n_ambi;
kom_sprintf_lite(s, "\tNM:i:%d\tAS:i:%d\tms:i:%d\tmd:i:%d", nm, p->p->dp_score, p->p->dp_max0, p->p->dp_max - p->p->dp_max2);
str_puts(s, "\tNM:i:"); str_puti(s, nm);
str_puts(s, "\tAS:i:"); str_puti(s, p->p->dp_score);
str_puts(s, "\tms:i:"); str_puti(s, p->p->dp_max0);
str_puts(s, "\tmd:i:"); str_puti(s, p->p->dp_max - p->p->dp_max2);
}

void mb_fmt_paf(kstring_t *s, const l2b_t *l2b, const mb_bseq1_t *t, const mb_hit_t *p, uint64_t opt_flag, int n_seg, int seg_idx)
{
kom_sprintf_lite(s, "%s", t->name);
if (n_seg > 1 && seg_idx >= 0)
kom_sprintf_lite(s, "/%d", seg_idx + 1);
kom_sprintf_lite(s, "\t%ld", (long)t->l_seq);
str_puts(s, t->name);
if (n_seg > 1 && seg_idx >= 0) {
str_putc(s, '/'); str_puti(s, seg_idx + 1);
}
str_putc(s, '\t'); str_puti(s, t->l_seq);
if (p == 0) { // for unmapped reads
kom_sprintf_lite(s, "\t*\t*\t*\t*\t*\t*\t*\t0\t0\t0\n");
str_puts(s, "\t*\t*\t*\t*\t*\t*\t*\t0\t0\t0\n");
s->s[s->l] = 0; // keep the always-room-for-an-extra-byte NUL invariant
return;
}
kom_sprintf_lite(s, "\t%d\t%d\t%c\t%s\t%ld\t%ld\t%ld\t%d\t%d\t%d\ttp:A:%c\ts1:i:%d\tcm:i:%d",
p->qs, p->qe, p->rev? '-' : '+', l2b->ctg[p->tid].name, (long)l2b->ctg[p->tid].len, (long)p->ts, (long)p->te,
p->mlen, p->blen, p->mapq, p->parent == p->id? 'P' : 'S', p->score, p->cnt);
if (p->parent == p->id) kom_sprintf_lite(s, "\ts2:i:%d", p->subsc >= 0? p->subsc : 0);
str_putc(s, '\t'); str_puti(s, p->qs);
str_putc(s, '\t'); str_puti(s, p->qe);
str_putc(s, '\t'); str_putc(s, p->rev? '-' : '+');
str_putc(s, '\t'); str_puts(s, l2b->ctg[p->tid].name);
str_putc(s, '\t'); str_puti(s, l2b->ctg[p->tid].len);
str_putc(s, '\t'); str_puti(s, p->ts);
str_putc(s, '\t'); str_puti(s, p->te);
str_putc(s, '\t'); str_puti(s, p->mlen);
str_putc(s, '\t'); str_puti(s, p->blen);
str_putc(s, '\t'); str_puti(s, p->mapq);
str_puts(s, "\ttp:A:"); str_putc(s, p->parent == p->id? 'P' : 'S');
str_puts(s, "\ts1:i:"); str_puti(s, p->score);
str_puts(s, "\tcm:i:"); str_puti(s, p->cnt);
if (p->parent == p->id) { str_puts(s, "\ts2:i:"); str_puti(s, p->subsc >= 0? p->subsc : 0); }
if (p->p) {
write_tags(s, p);
if (p->p->n_cigar > 0) {
int32_t i;
kom_sprintf_lite(s, "\tcg:Z:");
for (i = 0; i < p->p->n_cigar; ++i)
kom_sprintf_lite(s, "%d%c", p->p->cigar[i]>>4, MB_CIGAR_STR[p->p->cigar[i]&0xf]);
str_puts(s, "\tcg:Z:");
for (i = 0; i < p->p->n_cigar; ++i) {
str_puti(s, p->p->cigar[i]>>4); str_putc(s, MB_CIGAR_STR[p->p->cigar[i]&0xf]);
}
}
if (p->p->cs) kom_sprintf_lite(s, "\t%s", (char*)&p->p->cigar[p->p->n_cigar]);
if (p->p->cs) { str_putc(s, '\t'); str_puts(s, (char*)&p->p->cigar[p->p->n_cigar]); }
}
if ((opt_flag & MB_F_COPY_COMMENT) && t->comment) {
str_putc(s, '\t'); str_puts(s, t->comment);
}
if ((opt_flag & MB_F_COPY_COMMENT) && t->comment)
kom_sprintf_lite(s, "\t%s", t->comment);
kom_sprintf_lite(s, "\n");
str_putc(s, '\n');
s->s[s->l] = 0; // keep the always-room-for-an-extra-byte NUL invariant (matches mb_fmt_sam)
}

/**************
Expand Down Expand Up @@ -126,22 +197,6 @@ int mb_fmt_sam_hdr(kstring_t *str, const l2b_t *idx, const char *rg, const char
* SAM output *
**************/

static inline void str_enlarge(kstring_t *s, int l)
{
if (s->l + l + 1 > s->m) {
s->m = s->l + l + 1;
kom_roundup64(s->m);
s->s = kom_realloc(char, s->s, s->m);
}
}

static inline void str_copy(kstring_t *s, const char *st, const char *en)
{
str_enlarge(s, en - st);
memcpy(&s->s[s->l], st, en - st);
s->l += en - st;
}

static void sam_write_sq(kstring_t *s, char *seq, int l, int rev, int comp)
{
if (rev) {
Expand All @@ -168,27 +223,29 @@ static inline const mb_hit_t *get_sam_pri(int n_hit, const mb_hit_t *hit)
static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, const mb_hit_t *r, int64_t opt_flag)
{
if (r->p == 0) {
kom_sprintf_lite(s, "*");
str_putc(s, '*');
} else {
uint32_t k, clip_len[2];
clip_len[0] = r->rev? qlen - r->qe : r->qs;
clip_len[1] = r->rev? r->qs : qlen - r->qe;
if (in_tag) {
int clip_char = (((sam_flag&0x800) || ((sam_flag&0x100) && (opt_flag&MB_F_2ND_SEQ))) &&
!(opt_flag&MB_F_SUPP_SOFT)) ? 5 : 4;
kom_sprintf_lite(s, "\tCG:B:I");
if (clip_len[0]) kom_sprintf_lite(s, ",%u", clip_len[0]<<4|clip_char);
for (k = 0; k < r->p->n_cigar; ++k)
kom_sprintf_lite(s, ",%u", r->p->cigar[k]);
if (clip_len[1]) kom_sprintf_lite(s, ",%u", clip_len[1]<<4|clip_char);
str_puts(s, "\tCG:B:I");
if (clip_len[0]) { str_putc(s, ','); str_putu(s, clip_len[0]<<4|clip_char); }
for (k = 0; k < r->p->n_cigar; ++k) {
str_putc(s, ','); str_putu(s, r->p->cigar[k]);
}
if (clip_len[1]) { str_putc(s, ','); str_putu(s, clip_len[1]<<4|clip_char); }
} else {
int clip_char = (((sam_flag&0x800) || ((sam_flag&0x100) && (opt_flag&MB_F_2ND_SEQ))) &&
!(opt_flag&MB_F_SUPP_SOFT)) ? 'H' : 'S';
assert(clip_len[0] < qlen && clip_len[1] < qlen);
if (clip_len[0]) kom_sprintf_lite(s, "%d%c", clip_len[0], clip_char);
for (k = 0; k < r->p->n_cigar; ++k)
kom_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, MB_CIGAR_STR[r->p->cigar[k]&0xf]);
if (clip_len[1]) kom_sprintf_lite(s, "%d%c", clip_len[1], clip_char);
if (clip_len[0]) { str_puti(s, clip_len[0]); str_putc(s, clip_char); }
for (k = 0; k < r->p->n_cigar; ++k) {
str_puti(s, r->p->cigar[k]>>4); str_putc(s, MB_CIGAR_STR[r->p->cigar[k]&0xf]);
}
if (clip_len[1]) { str_puti(s, clip_len[1]); str_putc(s, clip_char); }
}
}
}
Expand All @@ -209,7 +266,7 @@ void mb_fmt_sam(void *km, kstring_t *s, const l2b_t *l2b, const mb_bseq1_t *t, i
} else r_prev = r_next = NULL;

// write QNAME and FLAG
kom_sprintf_lite(s, "%s", t->name);
str_puts(s, t->name);
flag = n_seg > 1? 0x1 : 0x0;
if (r == 0) {
flag |= 0x4;
Expand All @@ -225,17 +282,17 @@ void mb_fmt_sam(void *km, kstring_t *s, const l2b_t *l2b, const mb_bseq1_t *t, i
if (r_next == NULL) flag |= 0x8;
else if (r_next->rev) flag |= 0x20;
}
kom_sprintf_lite(s, "\t%d", flag);
str_putc(s, '\t'); str_puti(s, flag);

// write coordinate, MAPQ and CIGAR
if (r == 0) {
if (r_prev) {
this_tid = r_prev->tid, this_pos = r_prev->ts;
kom_sprintf_lite(s, "\t%s\t%d\t0\t*", l2b->ctg[this_tid].name, this_pos+1);
} else kom_sprintf_lite(s, "\t*\t0\t0\t*");
str_putc(s, '\t'); str_puts(s, l2b->ctg[this_tid].name); str_putc(s, '\t'); str_puti(s, this_pos+1); str_puts(s, "\t0\t*");
} else str_puts(s, "\t*\t0\t0\t*");
} else {
this_tid = r->tid, this_pos = r->ts;
kom_sprintf_lite(s, "\t%s\t%d\t%d\t", l2b->ctg[r->tid].name, r->ts+1, r->mapq);
str_putc(s, '\t'); str_puts(s, l2b->ctg[r->tid].name); str_putc(s, '\t'); str_puti(s, r->ts+1); str_putc(s, '\t'); str_puti(s, r->mapq); str_putc(s, '\t');
write_sam_cigar(s, flag, 0, t->l_seq, r, opt->flag);
}

Expand All @@ -249,65 +306,66 @@ void mb_fmt_sam(void *km, kstring_t *s, const l2b_t *l2b, const mb_bseq1_t *t, i
int next_pos5 = r_next->rev? r_next->te - 1 : r_next->ts;
tlen = next_pos5 - this_pos5;
}
kom_sprintf_lite(s, "\t=\t");
} else kom_sprintf_lite(s, "\t%s\t", l2b->ctg[r_next->tid].name);
kom_sprintf_lite(s, "%d\t", r_next->ts + 1);
str_puts(s, "\t=\t");
} else { str_putc(s, '\t'); str_puts(s, l2b->ctg[r_next->tid].name); str_putc(s, '\t'); }
str_puti(s, r_next->ts + 1); str_putc(s, '\t');
} else if (r_next) { // && this_tid < 0
kom_sprintf_lite(s, "\t%s\t%d\t", l2b->ctg[r_next->tid].name, r_next->ts + 1);
str_putc(s, '\t'); str_puts(s, l2b->ctg[r_next->tid].name); str_putc(s, '\t'); str_puti(s, r_next->ts + 1); str_putc(s, '\t');
} else if (this_tid >= 0) { // && r_next == NULL. In this case, mate is unmapped and placed at this segment's primary alignment
const mb_hit_t *r_pri;
r_pri = get_sam_pri(n_h, h);
if (r->tid == r_pri->tid)
kom_sprintf_lite(s, "\t=\t%d\t", r_pri->ts + 1);
else
kom_sprintf_lite(s, "\t%s\t%d\t", l2b->ctg[r_pri->tid].name, r_pri->ts + 1);
} else kom_sprintf_lite(s, "\t*\t0\t"); // neither has coordinates
if (r->tid == r_pri->tid) {
str_puts(s, "\t=\t"); str_puti(s, r_pri->ts + 1); str_putc(s, '\t');
} else {
str_putc(s, '\t'); str_puts(s, l2b->ctg[r_pri->tid].name); str_putc(s, '\t'); str_puti(s, r_pri->ts + 1); str_putc(s, '\t');
}
} else str_puts(s, "\t*\t0\t"); // neither has coordinates
if (tlen > 0) ++tlen;
else if (tlen < 0) --tlen;
kom_sprintf_lite(s, "%d\t", tlen);
} else kom_sprintf_lite(s, "\t*\t0\t0\t");
str_puti(s, tlen); str_putc(s, '\t');
} else str_puts(s, "\t*\t0\t0\t");

// write SEQ and QUAL
if (r == 0) {
sam_write_sq(s, t->seq, t->l_seq, 0, 0);
kom_sprintf_lite(s, "\t");
str_putc(s, '\t');
if (t->qual) sam_write_sq(s, t->qual, t->l_seq, 0, 0);
else kom_sprintf_lite(s, "*");
else str_putc(s, '*');
} else {
if ((flag & 0x900) == 0 || (opt->flag & MB_F_SUPP_SOFT)) {
sam_write_sq(s, t->seq, t->l_seq, r->rev, r->rev);
kom_sprintf_lite(s, "\t");
str_putc(s, '\t');
if (t->qual) sam_write_sq(s, t->qual, t->l_seq, r->rev, 0);
else kom_sprintf_lite(s, "*");
else str_putc(s, '*');
} else if ((flag & 0x100) && !(opt->flag & MB_F_2ND_SEQ)){
kom_sprintf_lite(s, "*\t*");
str_puts(s, "*\t*");
} else {
sam_write_sq(s, t->seq + r->qs, r->qe - r->qs, r->rev, r->rev);
kom_sprintf_lite(s, "\t");
str_putc(s, '\t');
if (t->qual) sam_write_sq(s, t->qual + r->qs, r->qe - r->qs, r->rev, 0);
else kom_sprintf_lite(s, "*");
else str_putc(s, '*');
}
}

// write tags
if (mb_rg_id[0]) kom_sprintf_lite(s, "\tRG:Z:%s", mb_rg_id);
if (n_seg > 2) kom_sprintf_lite(s, "\tFI:i:%d", seg_idx);
if (mb_rg_id[0]) { str_puts(s, "\tRG:Z:"); str_puts(s, mb_rg_id); }
if (n_seg > 2) { str_puts(s, "\tFI:i:"); str_puti(s, seg_idx); }
if (r) {
write_tags(s, r);
// MC:Z mate CIGAR and MQ:i mate MAPQ; r_next is the mate's primary (see above).
if (n_seg > 1 && r_next && r_next->p && r_next->p->n_cigar > 0 && mate_qlen > 0) {
kom_sprintf_lite(s, "\tMC:Z:");
str_puts(s, "\tMC:Z:");
write_sam_cigar(s, 0, 0, mate_qlen, r_next, opt->flag);
kom_sprintf_lite(s, "\tMQ:i:%d", r_next->mapq);
str_puts(s, "\tMQ:i:"); str_puti(s, r_next->mapq);
}
if (r->p->cs) kom_sprintf_lite(s, "\t%s", (char*)&r->p->cigar[r->p->n_cigar]);
if (r->p->cs) { str_putc(s, '\t'); str_puts(s, (char*)&r->p->cigar[r->p->n_cigar]); }
if (r->parent == r->id && r->p && n_h > 1 && h && r >= h && r - h < n_h) { // supplementary aln may exist
int i, n_sa = 0; // n_sa: number of SA fields
for (i = 0; i < n_h; ++i)
if (i != r - h && h[i].parent == h[i].id && h[i].p)
++n_sa;
if (n_sa > 0) {
kom_sprintf_lite(s, "\tSA:Z:");
str_puts(s, "\tSA:Z:");
for (i = 0; i < n_h; ++i) {
const mb_hit_t *q = &h[i];
int l_M, l_I = 0, l_D = 0, clip5 = 0, clip3 = 0;
Expand All @@ -316,13 +374,13 @@ void mb_fmt_sam(void *km, kstring_t *s, const l2b_t *l2b, const mb_bseq1_t *t, i
else l_M = q->te - q->ts, 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;
kom_sprintf_lite(s, "%s,%d,%c,", l2b->ctg[q->tid].name, q->ts+1, "+-"[q->rev]);
if (clip5) kom_sprintf_lite(s, "%dS", clip5);
if (l_M) kom_sprintf_lite(s, "%dM", l_M);
if (l_I) kom_sprintf_lite(s, "%dI", l_I);
if (l_D) kom_sprintf_lite(s, "%dD", l_D);
if (clip3) kom_sprintf_lite(s, "%dS", clip3);
kom_sprintf_lite(s, ",%d,%d;", q->mapq, q->blen - q->mlen + q->p->n_ambi);
str_puts(s, l2b->ctg[q->tid].name); str_putc(s, ','); str_puti(s, q->ts+1); str_putc(s, ','); str_putc(s, "+-"[q->rev]); str_putc(s, ',');
if (clip5) { str_puti(s, clip5); str_putc(s, 'S'); }
if (l_M) { str_puti(s, l_M); str_putc(s, 'M'); }
if (l_I) { str_puti(s, l_I); str_putc(s, 'I'); }
if (l_D) { str_puti(s, l_D); str_putc(s, 'D'); }
if (clip3) { str_puti(s, clip3); str_putc(s, 'S'); }
str_putc(s, ','); str_puti(s, q->mapq); str_putc(s, ','); str_puti(s, q->blen - q->mlen + q->p->n_ambi); str_putc(s, ';');
}
}
if (opt->xa_max > 0) {
Expand All @@ -346,9 +404,10 @@ void mb_fmt_sam(void *km, kstring_t *s, const l2b_t *l2b, const mb_bseq1_t *t, i
}
}

if ((opt->flag & MB_F_COPY_COMMENT) && t->comment)
kom_sprintf_lite(s, "\t%s", t->comment);
kom_sprintf_lite(s, "\n");
if ((opt->flag & MB_F_COPY_COMMENT) && t->comment) {
str_putc(s, '\t'); str_puts(s, t->comment);
}
str_putc(s, '\n');
s->s[s->l] = 0; // we always have room for an extra byte
}

Expand Down