X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=samtools.git;a=blobdiff_plain;f=bam_import.c;h=9d84328bdbdafb5f4e76a19e7edf176121aca5b4;hp=1dc906eb49f36428460be547b1582ab52432f651;hb=049705238f6c2d4cdd79022e7bcf2cc7d83ba9c7;hpb=d363084f0412f3bcdeb0304aeb0974c9a10c7649 diff --git a/bam_import.c b/bam_import.c index 1dc906e..9d84328 100644 --- a/bam_import.c +++ b/bam_import.c @@ -10,6 +10,7 @@ #endif #include "kstring.h" #include "bam.h" +#include "sam_header.h" #include "kseq.h" #include "khash.h" @@ -115,7 +116,7 @@ static bam_header_t *hash2header(const kh_ref_t *hash) bam_header_t *sam_header_read2(const char *fn) { bam_header_t *header; - int c, dret, ret; + int c, dret, ret, error = 0; gzFile fp; kstream_t *ks; kstring_t *str; @@ -134,6 +135,10 @@ bam_header_t *sam_header_read2(const char *fn) ks_getuntil(ks, 0, str, &dret); len = atoi(str->s); k = kh_put(ref, hash, s, &ret); + if (ret == 0) { + fprintf(stderr, "[sam_header_read2] duplicated sequence name: %s\n", s); + error = 1; + } kh_value(hash, k) = (uint64_t)len<<32 | i; if (dret != '\n') while ((c = ks_getc(ks)) != '\n' && c != -1); @@ -142,6 +147,7 @@ bam_header_t *sam_header_read2(const char *fn) gzclose(fp); free(str->s); free(str); fprintf(stderr, "[sam_header_read2] %d sequences loaded.\n", kh_size(hash)); + if (error) return 0; header = hash2header(hash); kh_destroy(ref, hash); return header; @@ -162,115 +168,49 @@ static inline void parse_error(int64_t n_lines, const char * __restrict msg) } static inline void append_text(bam_header_t *header, kstring_t *str) { - int x = header->l_text, y = header->l_text + str->l + 2; // 2 = 1 byte dret + 1 byte null + size_t x = header->l_text, y = header->l_text + str->l + 2; // 2 = 1 byte dret + 1 byte null kroundup32(x); kroundup32(y); - if (x < y) header->text = (char*)realloc(header->text, y); + if (x < y) + { + header->n_text = y; + header->text = (char*)realloc(header->text, y); + if ( !header->text ) + { + fprintf(stderr,"realloc failed to alloc %ld bytes\n", y); + abort(); + } + } + // Sanity check + if ( header->l_text+str->l+1 >= header->n_text ) + { + fprintf(stderr,"append_text FIXME: %ld>=%ld, x=%ld,y=%ld\n", header->l_text+str->l+1,header->n_text,x,y); + abort(); + } strncpy(header->text + header->l_text, str->s, str->l+1); // we cannot use strcpy() here. header->l_text += str->l + 1; header->text[header->l_text] = 0; } -int sam_header_parse_rg(bam_header_t *h) -{ - kstring_t *rgid, *rglib; - char *p, *q, *s, *r; - int n = 0; - - // free - if (h == 0) return 0; - bam_strmap_destroy(h->rg2lib); h->rg2lib = 0; - if (h->l_text < 3) return 0; - // parse @RG lines - h->rg2lib = bam_strmap_init(); - rgid = calloc(1, sizeof(kstring_t)); - rglib = calloc(1, sizeof(kstring_t)); - s = h->text; - while ((s = strstr(s, "@RG")) != 0) { - if (rgid->l && rglib->l) { - bam_strmap_put(h->rg2lib, rgid->s, rglib->s); - ++n; - } - rgid->l = rglib->l = 0; - s += 3; - r = s; - if ((p = strstr(s, "ID:")) != 0) { - q = p + 3; - for (p = q; *p && *p != '\t' && *p != '\r' && *p != '\n'; ++p); - kputsn(q, p - q, rgid); - } else { - fprintf(stderr, "[bam_header_parse] missing ID tag in @RG lines.\n"); - break; - } - if (r < p) r = p; - if ((p = strstr(s, "LB:")) != 0) { - q = p + 3; - for (p = q; *p && *p != '\t' && *p != '\r' && *p != '\n'; ++p); - kputsn(q, p - q, rglib); - } else { - fprintf(stderr, "[bam_header_parse] missing LB tag in @RG lines.\n"); - break; - } - if (r < p) r = p; - s = r + 3; - } - if (rgid->l && rglib->l) { - bam_strmap_put(h->rg2lib, rgid->s, rglib->s); - ++n; - } - free(rgid->s); free(rgid); - free(rglib->s); free(rglib); - if (n == 0) { - bam_strmap_destroy(h->rg2lib); - h->rg2lib = 0; - } - return n; -} - int sam_header_parse(bam_header_t *h) { + char **tmp; int i; - char *s, *p, *q, *r; - - // free free(h->target_len); free(h->target_name); h->n_targets = 0; h->target_len = 0; h->target_name = 0; if (h->l_text < 3) return 0; - // count number of @SQ - s = h->text; - while ((s = strstr(s, "@SQ")) != 0) { - ++h->n_targets; - s += 3; - } + if (h->dict == 0) h->dict = sam_header_parse2(h->text); + tmp = sam_header2list(h->dict, "SQ", "SN", &h->n_targets); if (h->n_targets == 0) return 0; - h->target_len = (uint32_t*)calloc(h->n_targets, 4); - h->target_name = (char**)calloc(h->n_targets, sizeof(void*)); - // parse @SQ lines - i = 0; - s = h->text; - while ((s = strstr(s, "@SQ")) != 0) { - s += 3; - r = s; - if ((p = strstr(s, "SN:")) != 0) { - q = p + 3; - for (p = q; *p && *p != '\t' && *p != '\r' && *p != '\n'; ++p); - h->target_name[i] = (char*)calloc(p - q + 1, 1); - strncpy(h->target_name[i], q, p - q); - } else goto header_err_ret; - if (r < p) r = p; - if ((p = strstr(s, "LN:")) != 0) h->target_len[i] = strtol(p + 3, 0, 10); - else goto header_err_ret; - if (r < p) r = p; - s = r + 3; - ++i; - } - sam_header_parse_rg(h); + h->target_name = calloc(h->n_targets, sizeof(void*)); + for (i = 0; i < h->n_targets; ++i) + h->target_name[i] = strdup(tmp[i]); + free(tmp); + tmp = sam_header2list(h->dict, "SQ", "LN", &h->n_targets); + h->target_len = calloc(h->n_targets, 4); + for (i = 0; i < h->n_targets; ++i) + h->target_len[i] = atoi(tmp[i]); + free(tmp); return h->n_targets; - -header_err_ret: - fprintf(stderr, "[bam_header_parse] missing SN or LN tag in @SQ lines.\n"); - free(h->target_len); free(h->target_name); - h->n_targets = 0; h->target_len = 0; h->target_name = 0; - return 0; } bam_header_t *sam_header_read(tamFile fp)