sam_header.c

   1 #include "sam_header.h"
   2 #include <stdio.h>
   3 #include <string.h>
   4 #include <ctype.h>
   5 #include <stdlib.h>
   6 #include <stdarg.h>
   7
   8 #include "khash.h"
   9 KHASH_MAP_INIT_STR(str, const char *)
  10
  11 struct _HeaderList
  12 {
  13     struct _HeaderList *next;
  14     void *data;
  15 };
  16 typedef struct _HeaderList list_t;
  17 typedef list_t HeaderDict;
  18
  19 typedef struct
  20 {
  21     char key[2];
  22     char *value;
  23 }
  24 HeaderTag;
  25
  26 typedef struct
  27 {
  28     char type[2];
  29     list_t *tags;
  30 }
  31 HeaderLine;
  32
  33 const char *o_hd_tags[] = {"SO","GO",NULL};
  34 const char *r_hd_tags[] = {"VN",NULL};
  35
  36 const char *o_sq_tags[] = {"AS","M5","UR","SP",NULL};
  37 const char *r_sq_tags[] = {"SN","LN",NULL};
  38 const char *u_sq_tags[] = {"SN",NULL};
  39
  40 const char *o_rg_tags[] = {"LB","DS","PU","PI","CN","DT","PL",NULL};
  41 const char *r_rg_tags[] = {"ID",NULL};
  42 const char *u_rg_tags[] = {"ID",NULL};
  43
  44 const char *o_pg_tags[] = {"VN","CL",NULL};
  45 const char *r_pg_tags[] = {"ID",NULL};
  46
  47 const char *types[]          = {"HD","SQ","RG","PG","CO",NULL};
  48 const char **optional_tags[] = {o_hd_tags,o_sq_tags,o_rg_tags,o_pg_tags,NULL,NULL};
  49 const char **required_tags[] = {r_hd_tags,r_sq_tags,r_rg_tags,r_pg_tags,NULL,NULL};
  50 const char **unique_tags[]   = {NULL,     u_sq_tags,u_rg_tags,NULL,NULL,NULL};
  51
  52
  53 static void debug(const char *format, ...)
  54 {
  55     va_list ap;
  56     va_start(ap, format);
  57     vfprintf(stderr, format, ap);
  58     va_end(ap);
  59 }
  60
  61 static list_t *list_append(list_t *root, void *data)
  62 {
  63     list_t *l = root;
  64     while (l && l->next)
  65         l = l->next;
  66     if ( l )
  67     {
  68         l->next = malloc(sizeof(list_t));
  69         l = l->next;
  70     }
  71     else
  72     {
  73         l = malloc(sizeof(list_t));
  74         root = l;
  75     }
  76     l->data = data;
  77     l->next = NULL;
  78     return root;
  79 }
  80
  81 static void list_free(list_t *root)
  82 {
  83     list_t *l = root;
  84     while (root)
  85     {
  86         l = root;
  87         root = root->next;
  88         free(l);
  89     }
  90 }
  91
  92
  93
  94 // Look for a tag "XY" in a predefined const char *[] array.
  95 static int tag_exists(const char *tag, const char **tags)
  96 {
  97     int itag=0;
  98     if ( !tags ) return -1;
  99     while ( tags[itag] )
 100     {
 101         if ( tags[itag][0]==tag[0] && tags[itag][1]==tag[1] ) return itag;
 102         itag++;
 103     }
 104     return -1;
 105 }
 106
 107
 108
 109 // Mimics the behaviour of getline, except it returns pointer to the next chunk of the text
 110 //  or NULL if everything has been read. The lineptr should be freed by the caller. The
 111 //  newline character is stripped.
 112 static const char *nextline(char **lineptr, size_t *n, const char *text)
 113 {
 114     int len;
 115     const char *to = text;
 116
 117     if ( !*to ) return NULL;
 118
 119     while ( *to && *to!='\n' && *to!='\r' ) to++;
 120     len = to - text + 1;
 121
 122     if ( *to )
 123     {
 124         // Advance the pointer for the next call
 125         if ( *to=='\n' ) to++;
 126         else if ( *to=='\r' && *(to+1)=='\n' ) to+=2;
 127     }
 128     if ( !len )
 129         return to;
 130
 131     if ( !*lineptr )
 132     {
 133         *lineptr = malloc(len);
 134         *n = len;
 135     }
 136     else if ( *n<len )
 137     {
 138         *lineptr = realloc(*lineptr, len);
 139         *n = len;
 140     }
 141     if ( !*lineptr ) {
 142                 debug("[nextline] Insufficient memory!\n");
 143                 return 0;
 144         }
 145
 146     memcpy(*lineptr,text,len);
 147     (*lineptr)[len-1] = 0;
 148
 149     return to;
 150 }
 151
 152 // name points to "XY", value_from points to the first character of the value string and
 153 //  value_to points to the last character of the value string.
 154 static HeaderTag *new_tag(const char *name, const char *value_from, const char *value_to)
 155 {
 156     HeaderTag *tag = malloc(sizeof(HeaderTag));
 157     int len = value_to-value_from+1;
 158
 159     tag->key[0] = name[0];
 160     tag->key[1] = name[1];
 161     tag->value = malloc(len+1);
 162     memcpy(tag->value,value_from,len+1);
 163     tag->value[len] = 0;
 164     return tag;
 165 }
 166
 167 static HeaderTag *header_line_has_tag(HeaderLine *hline, const char *key)
 168 {
 169     list_t *tags = hline->tags;
 170     while (tags)
 171     {
 172         HeaderTag *tag = tags->data;
 173         if ( tag->key[0]==key[0] && tag->key[1]==key[1] ) return tag;
 174         tags = tags->next;
 175     }
 176     return NULL;
 177 }
 178
 179
 180 // Return codes:
 181 //   0 .. different types or unique tags differ or conflicting tags, cannot be merged
 182 //   1 .. all tags identical -> no need to merge, drop one
 183 //   2 .. the unique tags match and there are some conflicting tags (same tag, different value) -> error, cannot be merged nor duplicated
 184 //   3 .. there are some missing complementary tags and no unique conflict -> can be merged into a single line
 185 static int sam_header_compare_lines(HeaderLine *hline1, HeaderLine *hline2)
 186 {
 187     HeaderTag *t1, *t2;
 188
 189     if ( hline1->type[0]!=hline2->type[0] || hline1->type[1]!=hline2->type[1] )
 190         return 0;
 191
 192     int itype = tag_exists(hline1->type,types);
 193     if ( itype==-1 ) {
 194                 debug("[sam_header_compare_lines] Unknown type [%c%c]\n", hline1->type[0],hline1->type[1]);
 195                 return -1; // FIXME (lh3): error; I do not know how this will be handled in Petr's code
 196         }
 197
 198     if ( unique_tags[itype] )
 199     {
 200         t1 = header_line_has_tag(hline1,unique_tags[itype][0]);
 201         t2 = header_line_has_tag(hline2,unique_tags[itype][0]);
 202         if ( !t1 || !t2 ) // this should never happen, the unique tags are required
 203             return 2;
 204
 205         if ( strcmp(t1->value,t2->value) )
 206             return 0;   // the unique tags differ, cannot be merged
 207     }
 208     if ( !required_tags[itype] && !optional_tags[itype] )
 209     {
 210         t1 = hline1->tags->data;
 211         t2 = hline2->tags->data;
 212         if ( !strcmp(t1->value,t2->value) ) return 1; // identical comments
 213         return 0;
 214     }
 215
 216     int missing=0, itag=0;
 217     while ( required_tags[itype] && required_tags[itype][itag] )
 218     {
 219         t1 = header_line_has_tag(hline1,required_tags[itype][itag]);
 220         t2 = header_line_has_tag(hline2,required_tags[itype][itag]);
 221         if ( !t1 && !t2 )
 222             return 2;       // this should never happen
 223         else if ( !t1 || !t2 )
 224             missing = 1;    // there is some tag missing in one of the hlines
 225         else if ( strcmp(t1->value,t2->value) )
 226         {
 227             if ( unique_tags[itype] )
 228                 return 2;   // the lines have a matching unique tag but have a conflicting tag
 229
 230             return 0;    // the lines contain conflicting tags, cannot be merged
 231         }
 232         itag++;
 233     }
 234     itag = 0;
 235     while ( optional_tags[itype] && optional_tags[itype][itag] )
 236     {
 237         t1 = header_line_has_tag(hline1,optional_tags[itype][itag]);
 238         t2 = header_line_has_tag(hline2,optional_tags[itype][itag]);
 239         if ( !t1 && !t2 )
 240         {
 241             itag++;
 242             continue;
 243         }
 244         if ( !t1 || !t2 )
 245             missing = 1;    // there is some tag missing in one of the hlines
 246         else if ( strcmp(t1->value,t2->value) )
 247         {
 248             if ( unique_tags[itype] )
 249                 return 2;   // the lines have a matching unique tag but have a conflicting tag
 250
 251             return 0;   // the lines contain conflicting tags, cannot be merged
 252         }
 253         itag++;
 254     }
 255     if ( missing ) return 3;    // there are some missing complementary tags with no conflicts, can be merged
 256     return 1;
 257 }
 258
 259
 260 static HeaderLine *sam_header_line_clone(const HeaderLine *hline)
 261 {
 262     list_t *tags;
 263     HeaderLine *out = malloc(sizeof(HeaderLine));
 264     out->type[0] = hline->type[0];
 265     out->type[1] = hline->type[1];
 266     out->tags = NULL;
 267
 268     tags = hline->tags;
 269     while (tags)
 270     {
 271         HeaderTag *old = tags->data;
 272
 273         HeaderTag *new = malloc(sizeof(HeaderTag));
 274         new->key[0] = old->key[0];
 275         new->key[1] = old->key[1];
 276         new->value  = strdup(old->value);
 277         out->tags = list_append(out->tags, new);
 278
 279         tags = tags->next;
 280     }
 281     return out;
 282 }
 283
 284 static int sam_header_line_merge_with(HeaderLine *out_hline, const HeaderLine *tmpl_hline)
 285 {
 286     list_t *tmpl_tags;
 287
 288     if ( out_hline->type[0]!=tmpl_hline->type[0] || out_hline->type[1]!=tmpl_hline->type[1] )
 289         return 0;
 290
 291     tmpl_tags = tmpl_hline->tags;
 292     while (tmpl_tags)
 293     {
 294         HeaderTag *tmpl_tag = tmpl_tags->data;
 295         HeaderTag *out_tag  = header_line_has_tag(out_hline, tmpl_tag->key);
 296         if ( !out_tag )
 297         {
 298             HeaderTag *tag = malloc(sizeof(HeaderTag));
 299             tag->key[0] = tmpl_tag->key[0];
 300             tag->key[1] = tmpl_tag->key[1];
 301             tag->value  = strdup(tmpl_tag->value);
 302             out_hline->tags = list_append(out_hline->tags,tag);
 303         }
 304         tmpl_tags = tmpl_tags->next;
 305     }
 306     return 1;
 307 }
 308
 309
 310 static HeaderLine *sam_header_line_parse(const char *headerLine)
 311 {
 312     HeaderLine *hline;
 313     HeaderTag *tag;
 314     const char *from, *to;
 315     from = headerLine;
 316
 317     if ( *from != '@' ) {
 318                 debug("[sam_header_line_parse] expected '@', got [%s]\n", headerLine);
 319                 return 0;
 320         }
 321     to = ++from;
 322
 323     while (*to && *to!='\t') to++;
 324     if ( to-from != 2 ) {
 325                 debug("[sam_header_line_parse] expected '@XY', got [%s]\n", headerLine);
 326                 return 0;
 327         }
 328
 329     hline = malloc(sizeof(HeaderLine));
 330     hline->type[0] = from[0];
 331     hline->type[1] = from[1];
 332     hline->tags = NULL;
 333
 334     int itype = tag_exists(hline->type, types);
 335
 336     from = to;
 337     while (*to && *to=='\t') to++;
 338     if ( to-from != 1 ) {
 339         debug("[sam_header_line_parse] multiple tabs on line [%s] (%d)\n", headerLine,(int)(to-from));
 340                 return 0;
 341         }
 342     from = to;
 343     while (*from)
 344     {
 345         while (*to && *to!='\t') to++;
 346
 347         if ( !required_tags[itype] && !optional_tags[itype] )
 348         {
 349             // CO is a special case, it can contain anything, including tabs
 350             if ( *to ) { to++; continue; }
 351             tag = new_tag("  ",from,to-1);
 352         }
 353         else
 354             tag = new_tag(from,from+3,to-1);
 355
 356         if ( header_line_has_tag(hline,tag->key) )
 357                 debug("The tag '%c%c' present (at least) twice on line [%s]\n", tag->key[0],tag->key[1], headerLine);
 358         hline->tags = list_append(hline->tags, tag);
 359
 360         from = to;
 361         while (*to && *to=='\t') to++;
 362         if ( *to && to-from != 1 ) {
 363                         debug("[sam_header_line_parse] multiple tabs on line [%s] (%d)\n", headerLine,(int)(to-from));
 364                         return 0;
 365                 }
 366
 367         from = to;
 368     }
 369     return hline;
 370 }
 371
 372
 373 // Must be of an existing type, all tags must be recognised and all required tags must be present
 374 static int sam_header_line_validate(HeaderLine *hline)
 375 {
 376     list_t *tags;
 377     HeaderTag *tag;
 378     int itype, itag;
 379
 380     // Is the type correct?
 381     itype = tag_exists(hline->type, types);
 382     if ( itype==-1 )
 383     {
 384         debug("The type [%c%c] not recognised.\n", hline->type[0],hline->type[1]);
 385         return 0;
 386     }
 387
 388     // Has all required tags?
 389     itag = 0;
 390     while ( required_tags[itype] && required_tags[itype][itag] )
 391     {
 392         if ( !header_line_has_tag(hline,required_tags[itype][itag]) )
 393         {
 394             debug("The tag [%c%c] required for [%c%c] not present.\n", required_tags[itype][itag][0],required_tags[itype][itag][1],
 395                 hline->type[0],hline->type[1]);
 396             return 0;
 397         }
 398         itag++;
 399     }
 400
 401     // Are all tags recognised?
 402     tags = hline->tags;
 403     while ( tags )
 404     {
 405         tag = tags->data;
 406         if ( !tag_exists(tag->key,required_tags[itype]) && !tag_exists(tag->key,optional_tags[itype]) )
 407         {
 408             debug("Unknown tag [%c%c] for [%c%c].\n", tag->key[0],tag->key[1], hline->type[0],hline->type[1]);
 409             return 0;
 410         }
 411         tags = tags->next;
 412     }
 413
 414     return 1;
 415 }
 416
 417
 418 static void print_header_line(FILE *fp, HeaderLine *hline)
 419 {
 420     list_t *tags = hline->tags;
 421     HeaderTag *tag;
 422
 423     fprintf(fp, "@%c%c", hline->type[0],hline->type[1]);
 424     while (tags)
 425     {
 426         tag = tags->data;
 427
 428         fprintf(fp, "\t");
 429         if ( tag->key[0]!=' ' || tag->key[1]!=' ' )
 430             fprintf(fp, "%c%c:", tag->key[0],tag->key[1]);
 431         fprintf(fp, "%s", tag->value);
 432
 433         tags = tags->next;
 434     }
 435     fprintf(fp,"\n");
 436 }
 437
 438
 439 static void sam_header_line_free(HeaderLine *hline)
 440 {
 441     list_t *tags = hline->tags;
 442     while (tags)
 443     {
 444         HeaderTag *tag = tags->data;
 445         free(tag->value);
 446         free(tag);
 447         tags = tags->next;
 448     }
 449     list_free(hline->tags);
 450     free(hline);
 451 }
 452
 453 void sam_header_free(void *_header)
 454 {
 455         HeaderDict *header = (HeaderDict*)_header;
 456     list_t *hlines = header;
 457     while (hlines)
 458     {
 459         sam_header_line_free(hlines->data);
 460         hlines = hlines->next;
 461     }
 462     list_free(header);
 463 }
 464
 465 HeaderDict *sam_header_clone(const HeaderDict *dict)
 466 {
 467     HeaderDict *out = NULL;
 468     while (dict)
 469     {
 470         HeaderLine *hline = dict->data;
 471         out = list_append(out, sam_header_line_clone(hline));
 472         dict = dict->next;
 473     }
 474     return out;
 475 }
 476
 477 // Returns a newly allocated string
 478 char *sam_header_write(const void *_header)
 479 {
 480         const HeaderDict *header = (const HeaderDict*)_header;
 481     char *out = NULL;
 482     int len=0, nout=0;
 483     const list_t *hlines;
 484
 485     // Calculate the length of the string to allocate
 486     hlines = header;
 487     while (hlines)
 488     {
 489         len += 4;   // @XY and \n
 490
 491         HeaderLine *hline = hlines->data;
 492         list_t *tags = hline->tags;
 493         while (tags)
 494         {
 495             HeaderTag *tag = tags->data;
 496             len += strlen(tag->value) + 1;                  // \t
 497             if ( tag->key[0]!=' ' || tag->key[1]!=' ' )
 498                 len += strlen(tag->value) + 3;              // XY:
 499             tags = tags->next;
 500         }
 501         hlines = hlines->next;
 502     }
 503
 504     nout = 0;
 505     out  = malloc(len+1);
 506     hlines = header;
 507     while (hlines)
 508     {
 509         HeaderLine *hline = hlines->data;
 510
 511         nout += sprintf(out+nout,"@%c%c",hline->type[0],hline->type[1]);
 512
 513         list_t *tags = hline->tags;
 514         while (tags)
 515         {
 516             HeaderTag *tag = tags->data;
 517             nout += sprintf(out+nout,"\t");
 518             if ( tag->key[0]!=' ' || tag->key[1]!=' ' )
 519                 nout += sprintf(out+nout,"%c%c:", tag->key[0],tag->key[1]);
 520             nout += sprintf(out+nout,"%s", tag->value);
 521             tags = tags->next;
 522         }
 523         hlines = hlines->next;
 524         nout += sprintf(out+nout,"\n");
 525     }
 526     out[len] = 0;
 527     return out;
 528 }
 529
 530 void *sam_header_parse2(const char *headerText)
 531 {
 532     list_t *hlines = NULL;
 533     HeaderLine *hline;
 534     const char *text;
 535     char *buf=NULL;
 536     size_t nbuf = 0;
 537
 538     if ( !headerText )
 539                 return 0;
 540
 541     text = headerText;
 542     while ( (text=nextline(&buf, &nbuf, text)) )
 543     {
 544         hline = sam_header_line_parse(buf);
 545         if ( hline && sam_header_line_validate(hline) )
 546             hlines = list_append(hlines, hline);
 547         else
 548         {
 549                         if (hline) sam_header_line_free(hline);
 550                         sam_header_free(hlines);
 551             if ( buf ) free(buf);
 552             return NULL;
 553         }
 554     }
 555     if ( buf ) free(buf);
 556
 557     return hlines;
 558 }
 559
 560 void *sam_header2tbl(const void *_dict, char type[2], char key_tag[2], char value_tag[2])
 561 {
 562         const HeaderDict *dict = (const HeaderDict*)_dict;
 563     const list_t *l   = dict;
 564     khash_t(str) *tbl = kh_init(str);
 565     khiter_t k;
 566     int ret;
 567
 568         if (_dict == 0) return tbl; // return an empty (not null) hash table
 569     while (l)
 570     {
 571         HeaderLine *hline = l->data;
 572         if ( hline->type[0]!=type[0] || hline->type[1]!=type[1] )
 573         {
 574             l = l->next;
 575             continue;
 576         }
 577
 578         HeaderTag *key, *value;
 579         key   = header_line_has_tag(hline,key_tag);
 580         value = header_line_has_tag(hline,value_tag);
 581         if ( !key || !value )
 582         {
 583             l = l->next;
 584             continue;
 585         }
 586
 587         k = kh_get(str, tbl, key->value);
 588         if ( k != kh_end(tbl) )
 589             debug("[sam_header_lookup_table] They key %s not unique.\n", key->value);
 590         k = kh_put(str, tbl, key->value, &ret);
 591         kh_value(tbl, k) = value->value;
 592
 593         l = l->next;
 594     }
 595     return tbl;
 596 }
 597
 598 char **sam_header2list(const void *_dict, char type[2], char key_tag[2], int *_n)
 599 {
 600         const HeaderDict *dict = (const HeaderDict*)_dict;
 601     const list_t *l   = dict;
 602     int max, n;
 603         char **ret;
 604
 605         ret = 0; *_n = max = n = 0;
 606     while (l)
 607     {
 608         HeaderLine *hline = l->data;
 609         if ( hline->type[0]!=type[0] || hline->type[1]!=type[1] )
 610         {
 611             l = l->next;
 612             continue;
 613         }
 614
 615         HeaderTag *key;
 616         key   = header_line_has_tag(hline,key_tag);
 617         if ( !key )
 618         {
 619             l = l->next;
 620             continue;
 621         }
 622
 623                 if (n == max) {
 624                         max = max? max<<1 : 4;
 625                         ret = realloc(ret, max * sizeof(void*));
 626                 }
 627                 ret[n++] = key->value;
 628
 629         l = l->next;
 630     }
 631         *_n = n;
 632     return ret;
 633 }
 634
 635 const char *sam_tbl_get(void *h, const char *key)
 636 {
 637         khash_t(str) *tbl = (khash_t(str)*)h;
 638         khint_t k;
 639         k = kh_get(str, tbl, key);
 640         return k == kh_end(tbl)? 0 : kh_val(tbl, k);
 641 }
 642
 643 int sam_tbl_size(void *h)
 644 {
 645         khash_t(str) *tbl = (khash_t(str)*)h;
 646         return h? kh_size(tbl) : 0;
 647 }
 648
 649 void sam_tbl_destroy(void *h)
 650 {
 651         khash_t(str) *tbl = (khash_t(str)*)h;
 652         kh_destroy(str, tbl);
 653 }
 654
 655 void *sam_header_merge(int n, const void **_dicts)
 656 {
 657         const HeaderDict **dicts = (const HeaderDict**)_dicts;
 658     HeaderDict *out_dict;
 659     int idict, status;
 660
 661     if ( n<2 ) return NULL;
 662
 663     out_dict = sam_header_clone(dicts[0]);
 664
 665     for (idict=1; idict<n; idict++)
 666     {
 667         const list_t *tmpl_hlines = dicts[idict];
 668
 669         while ( tmpl_hlines )
 670         {
 671             list_t *out_hlines = out_dict;
 672             int inserted = 0;
 673             while ( out_hlines )
 674             {
 675                 status = sam_header_compare_lines(tmpl_hlines->data, out_hlines->data);
 676                 if ( status==0 )
 677                 {
 678                     out_hlines = out_hlines->next;
 679                     continue;
 680                 }
 681
 682                 if ( status==2 )
 683                 {
 684                     print_header_line(stderr,tmpl_hlines->data);
 685                     print_header_line(stderr,out_hlines->data);
 686                     debug("Conflicting lines, cannot merge the headers.\n");
 687                                         return 0;
 688                 }
 689                 if ( status==3 )
 690                     sam_header_line_merge_with(out_hlines->data, tmpl_hlines->data);
 691
 692                 inserted = 1;
 693                 break;
 694             }
 695             if ( !inserted )
 696                 out_dict = list_append(out_dict, sam_header_line_clone(tmpl_hlines->data));
 697
 698             tmpl_hlines = tmpl_hlines->next;
 699         }
 700     }
 701
 702     return out_dict;
 703 }
 704
 705