Merge commit 'upstream/0.1.7a.dfsg'
[samtools.git] / 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