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