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