From a8c861864486f4c30dc664e3378ffd6e021b6364 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Tue, 11 Aug 2009 08:29:11 +0000 Subject: [PATCH] * samtools-0.1.5-27 (r443) * SEQ and QUAL can be "*" * parse CIGAR "=" and "X" as "M" --- bam.c | 10 ++++++---- bam_import.c | 25 +++++++++++++++---------- bamtk.c | 2 +- 3 files changed, 22 insertions(+), 15 deletions(-) diff --git a/bam.c b/bam.c index 52d97a6..619b46a 100644 --- a/bam.c +++ b/bam.c @@ -271,10 +271,12 @@ char *bam_format1_core(const bam_header_t *header, const bam1_t *b, int of) else if (c->mtid == c->tid) kputs("=\t", &str); else ksprintf(&str, "%s\t", header->target_name[c->mtid]); ksprintf(&str, "%d\t%d\t", c->mpos + 1, c->isize); - for (i = 0; i < c->l_qseq; ++i) kputc(bam_nt16_rev_table[bam1_seqi(s, i)], &str); - kputc('\t', &str); - if (t[0] == 0xff) kputc('*', &str); - else for (i = 0; i < c->l_qseq; ++i) kputc(t[i] + 33, &str); + if (c->l_qseq) { + for (i = 0; i < c->l_qseq; ++i) kputc(bam_nt16_rev_table[bam1_seqi(s, i)], &str); + kputc('\t', &str); + if (t[0] == 0xff) kputc('*', &str); + else for (i = 0; i < c->l_qseq; ++i) kputc(t[i] + 33, &str); + } else ksprintf(&str, "*\t*"); s = bam1_aux(b); while (s < b->data + b->data_len) { uint8_t type, key[2]; diff --git a/bam_import.c b/bam_import.c index 4c9395b..19cd75f 100644 --- a/bam_import.c +++ b/bam_import.c @@ -359,7 +359,7 @@ int sam_read1(tamFile fp, bam_header_t *header, bam1_t *b) for (i = 0, s = str->s; i != c->n_cigar; ++i) { x = strtol(s, &t, 10); op = toupper(*t); - if (op == 'M') op = BAM_CMATCH; + if (op == 'M' || op == '=' || op == 'X') op = BAM_CMATCH; else if (op == 'I') op = BAM_CINS; else if (op == 'D') op = BAM_CDEL; else if (op == 'N') op = BAM_CREF_SKIP; @@ -373,7 +373,10 @@ int sam_read1(tamFile fp, bam_header_t *header, bam1_t *b) if (*s) parse_error(fp->n_lines, "unmatched CIGAR operation"); c->bin = bam_reg2bin(c->pos, bam_calend(c, bam1_cigar(b))); doff += c->n_cigar * 4; - } else c->bin = bam_reg2bin(c->pos, c->pos + 1); + } else { + if (!(c->flag&BAM_FUNMAP)) parse_error(fp->n_lines, "mapped sequence without CIGAR"); + c->bin = bam_reg2bin(c->pos, c->pos + 1); + } } { // mtid, mpos, isize ret = ks_getuntil(ks, KS_SEP_TAB, str, &dret); z += str->l + 1; @@ -386,16 +389,18 @@ int sam_read1(tamFile fp, bam_header_t *header, bam1_t *b) } { // seq and qual int i; - uint8_t *p; + uint8_t *p = 0; if (ks_getuntil(ks, KS_SEP_TAB, str, &dret) < 0) return -5; // seq z += str->l + 1; - c->l_qseq = strlen(str->s); - if (c->n_cigar && c->l_qseq != (int32_t)bam_cigar2qlen(c, bam1_cigar(b))) - parse_error(fp->n_lines, "CIGAR and sequence length are inconsistent"); - p = (uint8_t*)alloc_data(b, doff + c->l_qseq + (c->l_qseq+1)/2) + doff; - memset(p, 0, (c->l_qseq+1)/2); - for (i = 0; i < c->l_qseq; ++i) - p[i/2] |= bam_nt16_table[(int)str->s[i]] << 4*(1-i%2); + if (strcmp(str->s, "*")) { + c->l_qseq = strlen(str->s); + if (c->n_cigar && c->l_qseq != (int32_t)bam_cigar2qlen(c, bam1_cigar(b))) + parse_error(fp->n_lines, "CIGAR and sequence length are inconsistent"); + p = (uint8_t*)alloc_data(b, doff + c->l_qseq + (c->l_qseq+1)/2) + doff; + memset(p, 0, (c->l_qseq+1)/2); + for (i = 0; i < c->l_qseq; ++i) + p[i/2] |= bam_nt16_table[(int)str->s[i]] << 4*(1-i%2); + } else c->l_qseq = 0; if (ks_getuntil(ks, KS_SEP_TAB, str, &dret) < 0) return -6; // qual z += str->l + 1; if (strcmp(str->s, "*") && c->l_qseq != strlen(str->s)) diff --git a/bamtk.c b/bamtk.c index 56fa73a..2df9032 100644 --- a/bamtk.c +++ b/bamtk.c @@ -9,7 +9,7 @@ #endif #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.1.5-26 (r442)" +#define PACKAGE_VERSION "0.1.5-27 (r443)" #endif int bam_taf2baf(int argc, char *argv[]); -- 2.39.2