zzh
/
minimap2
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

r121: output QUAL and unmapped to SAM

This commit is contained in:
Heng Li 2017-06-30 14:40:54 -04:00
parent d11049eb32
commit fce87ce7bd
5 changed files with 43 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
bseq.c
View File

@ -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;

2
bseq.h
View File

@ -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);

55
format.c
View File

@ -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);
}
}

2
main.c
View File

@ -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()
{

5
map.c
View File

@ -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);
}