samtools/sam_header.c.pysam.c

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