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