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             tag = new_tag("  ",from,to-1);
 349         else
 350             tag = new_tag(from,from+3,to-1);
 351
 352         if ( header_line_has_tag(hline,tag->key) )
 353                 debug("The tag '%c%c' present (at least) twice on line [%s]\n", tag->key[0],tag->key[1], headerLine);
 354         hline->tags = list_append(hline->tags, tag);
 355
 356         from = to;
 357         while (*to && *to=='\t') to++;
 358         if ( *to && to-from != 1 ) {
 359                         debug("[sam_header_line_parse] multiple tabs on line [%s] (%d)\n", headerLine,(int)(to-from));
 360                         return 0;
 361                 }
 362
 363         from = to;
 364     }
 365     return hline;
 366 }
 367
 368
 369 // Must be of an existing type, all tags must be recognised and all required tags must be present
 370 static int sam_header_line_validate(HeaderLine *hline)
 371 {
 372     list_t *tags;
 373     HeaderTag *tag;
 374     int itype, itag;
 375
 376     // Is the type correct?
 377     itype = tag_exists(hline->type, types);
 378     if ( itype==-1 )
 379     {
 380         debug("The type [%c%c] not recognised.\n", hline->type[0],hline->type[1]);
 381         return 0;
 382     }
 383
 384     // Has all required tags?
 385     itag = 0;
 386     while ( required_tags[itype] && required_tags[itype][itag] )
 387     {
 388         if ( !header_line_has_tag(hline,required_tags[itype][itag]) )
 389         {
 390             debug("The tag [%c%c] required for [%c%c] not present.\n", required_tags[itype][itag][0],required_tags[itype][itag][1],
 391                 hline->type[0],hline->type[1]);
 392             return 0;
 393         }
 394         itag++;
 395     }
 396
 397     // Are all tags recognised?
 398     tags = hline->tags;
 399     while ( tags )
 400     {
 401         tag = tags->data;
 402         if ( !tag_exists(tag->key,required_tags[itype]) && !tag_exists(tag->key,optional_tags[itype]) )
 403         {
 404             debug("Unknown tag [%c%c] for [%c%c].\n", tag->key[0],tag->key[1], hline->type[0],hline->type[1]);
 405             return 0;
 406         }
 407         tags = tags->next;
 408     }
 409
 410     return 1;
 411 }
 412
 413
 414 static void print_header_line(FILE *fp, HeaderLine *hline)
 415 {
 416     list_t *tags = hline->tags;
 417     HeaderTag *tag;
 418
 419     fprintf(fp, "@%c%c", hline->type[0],hline->type[1]);
 420     while (tags)
 421     {
 422         tag = tags->data;
 423
 424         fprintf(fp, "\t");
 425         if ( tag->key[0]!=' ' || tag->key[1]!=' ' )
 426             fprintf(fp, "%c%c:", tag->key[0],tag->key[1]);
 427         fprintf(fp, "%s", tag->value);
 428
 429         tags = tags->next;
 430     }
 431     fprintf(fp,"\n");
 432 }
 433
 434
 435 static void sam_header_line_free(HeaderLine *hline)
 436 {
 437     list_t *tags = hline->tags;
 438     while (tags)
 439     {
 440         HeaderTag *tag = tags->data;
 441         free(tag->value);
 442         free(tag);
 443         tags = tags->next;
 444     }
 445     list_free(hline->tags);
 446     free(hline);
 447 }
 448
 449 void sam_header_free(void *_header)
 450 {
 451         HeaderDict *header = (HeaderDict*)_header;
 452     list_t *hlines = header;
 453     while (hlines)
 454     {
 455         sam_header_line_free(hlines->data);
 456         hlines = hlines->next;
 457     }
 458     list_free(header);
 459 }
 460
 461 HeaderDict *sam_header_clone(const HeaderDict *dict)
 462 {
 463     HeaderDict *out = NULL;
 464     while (dict)
 465     {
 466         HeaderLine *hline = dict->data;
 467         out = list_append(out, sam_header_line_clone(hline));
 468         dict = dict->next;
 469     }
 470     return out;
 471 }
 472
 473 // Returns a newly allocated string
 474 char *sam_header_write(const void *_header)
 475 {
 476         const HeaderDict *header = (const HeaderDict*)_header;
 477     char *out = NULL;
 478     int len=0, nout=0;
 479     const list_t *hlines;
 480
 481     // Calculate the length of the string to allocate
 482     hlines = header;
 483     while (hlines)
 484     {
 485         len += 4;   // @XY and \n
 486
 487         HeaderLine *hline = hlines->data;
 488         list_t *tags = hline->tags;
 489         while (tags)
 490         {
 491             HeaderTag *tag = tags->data;
 492             len += strlen(tag->value) + 1;                  // \t
 493             if ( tag->key[0]!=' ' || tag->key[1]!=' ' )
 494                 len += strlen(tag->value) + 3;              // XY:
 495             tags = tags->next;
 496         }
 497         hlines = hlines->next;
 498     }
 499
 500     nout = 0;
 501     out  = malloc(len+1);
 502     hlines = header;
 503     while (hlines)
 504     {
 505         HeaderLine *hline = hlines->data;
 506
 507         nout += sprintf(out+nout,"@%c%c",hline->type[0],hline->type[1]);
 508
 509         list_t *tags = hline->tags;
 510         while (tags)
 511         {
 512             HeaderTag *tag = tags->data;
 513             nout += sprintf(out+nout,"\t");
 514             if ( tag->key[0]!=' ' || tag->key[1]!=' ' )
 515                 nout += sprintf(out+nout,"%c%c:", tag->key[0],tag->key[1]);
 516             nout += sprintf(out+nout,"%s", tag->value);
 517             tags = tags->next;
 518         }
 519         hlines = hlines->next;
 520         nout += sprintf(out+nout,"\n");
 521     }
 522     out[len] = 0;
 523     return out;
 524 }
 525
 526 void *sam_header_parse2(const char *headerText)
 527 {
 528     list_t *hlines = NULL;
 529     HeaderLine *hline;
 530     const char *text;
 531     char *buf=NULL;
 532     size_t nbuf = 0;
 533
 534     if ( !headerText )
 535                 return 0;
 536
 537     text = headerText;
 538     while ( (text=nextline(&buf, &nbuf, text)) )
 539     {
 540         hline = sam_header_line_parse(buf);
 541         if ( hline && sam_header_line_validate(hline) )
 542             hlines = list_append(hlines, hline);
 543         else
 544         {
 545                         if (hline) sam_header_line_free(hline);
 546                         sam_header_free(hlines);
 547             if ( buf ) free(buf);
 548             return NULL;
 549         }
 550     }
 551     if ( buf ) free(buf);
 552
 553     return hlines;
 554 }
 555
 556 void *sam_header2tbl(const void *_dict, char type[2], char key_tag[2], char value_tag[2])
 557 {
 558         const HeaderDict *dict = (const HeaderDict*)_dict;
 559     const list_t *l   = dict;
 560     khash_t(str) *tbl = kh_init(str);
 561     khiter_t k;
 562     int ret;
 563
 564         if (_dict == 0) return tbl; // return an empty (not null) hash table
 565     while (l)
 566     {
 567         HeaderLine *hline = l->data;
 568         if ( hline->type[0]!=type[0] || hline->type[1]!=type[1] )
 569         {
 570             l = l->next;
 571             continue;
 572         }
 573
 574         HeaderTag *key, *value;
 575         key   = header_line_has_tag(hline,key_tag);
 576         value = header_line_has_tag(hline,value_tag);
 577         if ( !key || !value )
 578         {
 579             l = l->next;
 580             continue;
 581         }
 582
 583         k = kh_get(str, tbl, key->value);
 584         if ( k != kh_end(tbl) )
 585             debug("[sam_header_lookup_table] They key %s not unique.\n", key->value);
 586         k = kh_put(str, tbl, key->value, &ret);
 587         kh_value(tbl, k) = value->value;
 588
 589         l = l->next;
 590     }
 591     return tbl;
 592 }
 593
 594 char **sam_header2list(const void *_dict, char type[2], char key_tag[2], int *_n)
 595 {
 596         const HeaderDict *dict = (const HeaderDict*)_dict;
 597     const list_t *l   = dict;
 598     int max, n;
 599         char **ret;
 600
 601         ret = 0; *_n = max = n = 0;
 602     while (l)
 603     {
 604         HeaderLine *hline = l->data;
 605         if ( hline->type[0]!=type[0] || hline->type[1]!=type[1] )
 606         {
 607             l = l->next;
 608             continue;
 609         }
 610
 611         HeaderTag *key;
 612         key   = header_line_has_tag(hline,key_tag);
 613         if ( !key )
 614         {
 615             l = l->next;
 616             continue;
 617         }
 618
 619                 if (n == max) {
 620                         max = max? max<<1 : 4;
 621                         ret = realloc(ret, max * sizeof(void*));
 622                 }
 623                 ret[n++] = key->value;
 624
 625         l = l->next;
 626     }
 627         *_n = n;
 628     return ret;
 629 }
 630
 631 const char *sam_tbl_get(void *h, const char *key)
 632 {
 633         khash_t(str) *tbl = (khash_t(str)*)h;
 634         khint_t k;
 635         k = kh_get(str, tbl, key);
 636         return k == kh_end(tbl)? 0 : kh_val(tbl, k);
 637 }
 638
 639 int sam_tbl_size(void *h)
 640 {
 641         khash_t(str) *tbl = (khash_t(str)*)h;
 642         return h? kh_size(tbl) : 0;
 643 }
 644
 645 void sam_tbl_destroy(void *h)
 646 {
 647         khash_t(str) *tbl = (khash_t(str)*)h;
 648         kh_destroy(str, tbl);
 649 }
 650
 651 void *sam_header_merge(int n, const void **_dicts)
 652 {
 653         const HeaderDict **dicts = (const HeaderDict**)_dicts;
 654     HeaderDict *out_dict;
 655     int idict, status;
 656
 657     if ( n<2 ) return NULL;
 658
 659     out_dict = sam_header_clone(dicts[0]);
 660
 661     for (idict=1; idict<n; idict++)
 662     {
 663         const list_t *tmpl_hlines = dicts[idict];
 664
 665         while ( tmpl_hlines )
 666         {
 667             list_t *out_hlines = out_dict;
 668             int inserted = 0;
 669             while ( out_hlines )
 670             {
 671                 status = sam_header_compare_lines(tmpl_hlines->data, out_hlines->data);
 672                 if ( status==0 )
 673                 {
 674                     out_hlines = out_hlines->next;
 675                     continue;
 676                 }
 677
 678                 if ( status==2 )
 679                 {
 680                     print_header_line(stderr,tmpl_hlines->data);
 681                     print_header_line(stderr,out_hlines->data);
 682                     debug("Conflicting lines, cannot merge the headers.\n");
 683                                         return 0;
 684                 }
 685                 if ( status==3 )
 686                     sam_header_line_merge_with(out_hlines->data, tmpl_hlines->data);
 687
 688                 inserted = 1;
 689                 break;
 690             }
 691             if ( !inserted )
 692                 out_dict = list_append(out_dict, sam_header_line_clone(tmpl_hlines->data));
 693
 694             tmpl_hlines = tmpl_hlines->next;
 695         }
 696     }
 697
 698     return out_dict;
 699 }
 700
 701