sam_header.c

#include "sam_header.h"
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <stdlib.h>
#include <stdarg.h>
#include "khash.h"

KHASH_MAP_INIT_STR(str, const char *)

struct _HeaderList
{
	struct _HeaderList *last;   // Hack: Used and maintained only by list_append_to_end. Maintained in the root node only.
	struct _HeaderList *next;
	void *data;
};

typedef struct _HeaderList list_t;
typedef list_t HeaderDict;

typedef struct
{
	char key[2];
	char *value;
} HeaderTag;

typedef struct
{
	char type[2];
	list_t *tags;
} HeaderLine;

const char *o_hd_tags[]      = {"SO", "GO", NULL};
const char *r_hd_tags[]      = {"VN", NULL};
const char *o_sq_tags[]      = {"AS", "M5", "UR", "SP", NULL};
const char *r_sq_tags[]      = {"SN", "LN", NULL};
const char *u_sq_tags[]      = {"SN", NULL};
const char *o_rg_tags[]      = {"CN", "DS", "DT", "FO", "KS", "LB", "PG", "PI", "PL", "PU", NULL};
const char *r_rg_tags[]      = {"ID", "SM", "PO", NULL};
const char *u_rg_tags[]      = {"ID", NULL};
const char *o_pg_tags[]      = {"VN", "CL", NULL};
const char *r_pg_tags[]      = {"ID", NULL};
const char *types[]          = {"HD", "SQ", "RG", "PG", "CO", NULL};
const char **optional_tags[] = {o_hd_tags, o_sq_tags, o_rg_tags, o_pg_tags, NULL, NULL};
const char **required_tags[] = {r_hd_tags ,r_sq_tags, r_rg_tags, r_pg_tags, NULL, NULL};
const char **unique_tags[]   = {NULL, u_sq_tags, u_rg_tags, NULL, NULL, NULL};


static void debug(const char *format, ...)
{
	va_list ap;
	va_start(ap, format);
	vfprintf(stderr, format, ap);
	va_end(ap);
}

// Relies on the root->last being correct. Do not use with the other list_*
//  routines unless they are fixed to modify root->last as well.
static list_t *list_append_to_end(list_t *root, void *data)
{
	list_t *l = (list_t*)malloc(sizeof(list_t));

	l->last = l;
	l->next = NULL;
	l->data = data;

	if (!root)
		return l;

	root->last->next = l;
	root->last = l;

	return root;
}

static list_t *list_append(list_t *root, void *data)
{
	list_t *l = root;

	while (l && l->next)
		l = l->next;

	if (l)
	{
		l->next = (list_t*)malloc(sizeof(list_t));
		l = l->next;
	}
	else
	{
		l = (list_t*)malloc(sizeof(list_t));
		root = l;
	}

	l->data = data;
	l->next = NULL;

	return root;
}

static void list_free(list_t *root)
{
	list_t *l = root;

	while (root)
	{
		l = root;
		root = root->next;
		free(l);
	}
}

// Look for a tag "XY" in a predefined const char *[] array.
static int tag_exists(const char *tag, const char **tags)
{
	int itag = 0;

	if (!tags)
		return -1;

	while (tags[itag])
	{
		if ((tags[itag][0] == tag[0]) && (tags[itag][1] == tag[1]))
			return itag;

		itag++;
	}

	return -1;
}

// Mimics the behaviour of getline, except it returns pointer to the next chunk of the text
//  or NULL if everything has been read. The lineptr should be freed by the caller. The
//  newline character is stripped.
static const char *nextline(char **lineptr, size_t *n, const char *text)
{
	int len;
	const char *to = text;

	if (!*to)
		return NULL;

	while (*to && (*to != '\n') && (*to != '\r'))
		to++;

	len = to - text + 1;

	if (*to)
	{
		// Advance the pointer for the next call
		if (*to == '\n')
			to++;
		else if ((*to=='\r') && (*(to+1)=='\n'))
			to+=2;
	}

	if (!len)
		return to;

	if (!*lineptr)
	{
		*lineptr = (char*)malloc(len);
		*n = len;
	}
	else if (*n < len)
	{
		*lineptr = (char*)realloc(*lineptr, len);
		*n = len;
	}

	if (!*lineptr)
	{
		debug("[nextline] Insufficient memory!\n");
		return 0;
	}

	memcpy(*lineptr, text, len);
	(*lineptr)[len-1] = 0;

	return to;
}

// name points to "XY", value_from points to the first character of the value string and
//  value_to points to the last character of the value string.
static HeaderTag *new_tag(const char *name, const char *value_from, const char *value_to)
{
	HeaderTag *tag = (HeaderTag*)malloc(sizeof(HeaderTag));
	int len = value_to - value_from + 1;

	tag->key[0] = name[0];
	tag->key[1] = name[1];
	tag->value = (char*)malloc(len+1);
	memcpy(tag->value, value_from, len+1);
	tag->value[len] = 0;

	return tag;
}

static HeaderTag *header_line_has_tag(HeaderLine *hline, const char *key)
{
	list_t *tags = hline->tags;

	while (tags)
	{
		HeaderTag *tag = (HeaderTag*)tags->data;

		if ((tag->key[0] == key[0]) && (tag->key[1] == key[1]))
			return tag;

		tags = tags->next;
	}

	return NULL;
}

// Return codes:
//   0 .. different types or unique tags differ or conflicting tags, cannot be merged
//   1 .. all tags identical -> no need to merge, drop one
//   2 .. the unique tags match and there are some conflicting tags (same tag, different value) -> error, cannot be merged nor duplicated
//   3 .. there are some missing complementary tags and no unique conflict -> can be merged into a single line
static int sam_header_compare_lines(HeaderLine *hline1, HeaderLine *hline2)
{
	int itype;
	int itag;
	int missing;
	HeaderTag *t1, *t2;

	if ((hline1->type[0] != hline2->type[0]) || (hline1->type[1] != hline2->type[1]))
		return 0;

	itype = tag_exists(hline1->type, types);

	if (itype == -1)
	{
		debug("[sam_header_compare_lines] Unknown type [%c%c]\n", hline1->type[0], hline1->type[1]);
		return -1; // FIXME (lh3): error; I do not know how this will be handled in Petr's code
	}

	if (unique_tags[itype])
	{
		t1 = header_line_has_tag(hline1, unique_tags[itype][0]);
		t2 = header_line_has_tag(hline2, unique_tags[itype][0]);

		if (!t1 || !t2) // this should never happen, the unique tags are required
			return 2;

		if (strcmp(t1->value, t2->value))
			return 0;   // the unique tags differ, cannot be merged
	}

	if (!required_tags[itype] && !optional_tags[itype])
	{
		t1 = (HeaderTag*)hline1->tags->data;
		t2 = (HeaderTag*)hline2->tags->data;

		if (!strcmp(t1->value, t2->value))
			return 1; // identical comments

		return 0;
	}

	missing = 0;
	itag = 0;

	while (required_tags[itype] && required_tags[itype][itag])
	{
		t1 = header_line_has_tag(hline1, required_tags[itype][itag]);
		t2 = header_line_has_tag(hline2, required_tags[itype][itag]);

		if (!t1 && !t2)
			return 2;       // this should never happen
		else if (!t1 || !t2)
			missing = 1;    // there is some tag missing in one of the hlines
		else if (strcmp(t1->value, t2->value))
		{
			if (unique_tags[itype])
				return 2;   // the lines have a matching unique tag but have a conflicting tag

			return 0;    // the lines contain conflicting tags, cannot be merged
		}

		itag++;
	}

	itag = 0;

	while (optional_tags[itype] && optional_tags[itype][itag])
	{
		t1 = header_line_has_tag(hline1, optional_tags[itype][itag]);
		t2 = header_line_has_tag(hline2, optional_tags[itype][itag]);

		if (!t1 && !t2)
		{
			itag++;
			continue;
		}

		if (!t1 || !t2)
			missing = 1;    // there is some tag missing in one of the hlines
		else if (strcmp(t1->value, t2->value))
		{
			if (unique_tags[itype])
				return 2;   // the lines have a matching unique tag but have a conflicting tag

			return 0;   // the lines contain conflicting tags, cannot be merged
		}

		itag++;
	}

	if (missing)
		return 3;    // there are some missing complementary tags with no conflicts, can be merged

	return 1;
}

static HeaderLine *sam_header_line_clone(const HeaderLine *hline)
{
	list_t *tags;
	HeaderLine *out = (HeaderLine*)malloc(sizeof(HeaderLine));
	HeaderTag *old;
	HeaderTag *nw;

	out->type[0] = hline->type[0];
	out->type[1] = hline->type[1];
	out->tags = NULL;
	tags = hline->tags;

	while (tags)
	{
		old = (HeaderTag*)tags->data;
		nw = (HeaderTag*)malloc(sizeof(HeaderTag));
		nw->key[0] = old->key[0];
		nw->key[1] = old->key[1];
		nw->value  = strdup(old->value);
		out->tags = list_append(out->tags, nw);
		tags = tags->next;
		free(nw);
	}

	return out;
}

static int sam_header_line_merge_with(HeaderLine *out_hline, const HeaderLine *tmpl_hline)
{
	list_t *tmpl_tags;

	if ((out_hline->type[0] != tmpl_hline->type[0]) || (out_hline->type[1] != tmpl_hline->type[1]))
		return 0;

	tmpl_tags = tmpl_hline->tags;

	while (tmpl_tags)
	{
		HeaderTag *tmpl_tag = (HeaderTag*)tmpl_tags->data;
		HeaderTag *out_tag  = header_line_has_tag(out_hline, tmpl_tag->key);
		if (!out_tag)
		{
			HeaderTag *tag = (HeaderTag*)malloc(sizeof(HeaderTag));
			tag->key[0] = tmpl_tag->key[0];
			tag->key[1] = tmpl_tag->key[1];
			tag->value  = strdup(tmpl_tag->value);
			out_hline->tags = list_append(out_hline->tags, tag);
		}

		tmpl_tags = tmpl_tags->next;
	}

	return 1;
}

static HeaderLine *sam_header_line_parse(const char *headerLine)
{
	int itype;
	HeaderLine *hline;
	HeaderTag *tag;
	const char *from, *to;
	from = headerLine;

	if (*from != '@')
	{
		debug("[sam_header_line_parse] expected '@', got [%s]\n", headerLine);
		return 0;
	}

	to = ++from;
	while (*to && (*to != '\t'))
		to++;

	if ((to - from) != 2)
	{
		debug("[sam_header_line_parse] expected '@XY', got [%s]\nHint: The header tags must be tab-separated.\n", headerLine);
		return 0;
	}

	hline = (HeaderLine*)malloc(sizeof(HeaderLine));
	hline->type[0] = from[0];
	hline->type[1] = from[1];
	hline->tags = NULL;
	itype = tag_exists(hline->type, types);
	from = to;

	while (*to && (*to == '\t'))
		to++;

	if ((to - from) != 1)
	{
		debug("[sam_header_line_parse] multiple tabs on line [%s] (%d)\n", headerLine, (int)(to - from));
		return 0;
	}

	from = to;

	while (*from)
	{
		while (*to && (*to != '\t'))
			to++;

		if (!required_tags[itype] && !optional_tags[itype])
		{
			// CO is a special case, it can contain anything, including tabs
			if (*to)
			{
				to++;
				continue;
			}
			tag = new_tag("  ", from, to - 1);
		}
		else
			tag = new_tag(from, from + 3, to - 1);

		if (header_line_has_tag(hline, tag->key))
			debug("The tag '%c%c' present (at least) twice on line [%s]\n", tag->key[0], tag->key[1], headerLine);

		hline->tags = list_append(hline->tags, tag);
		from = to;

		while (*to && (*to == '\t'))
			to++;

		if (*to && (to-from != 1))
		{
			debug("[sam_header_line_parse] multiple tabs on line [%s] (%d)\n", headerLine,(int)(to-from));
			return 0;
		}

		from = to;
	}

	return hline;
}

// Must be of an existing type, all tags must be recognised and all required tags must be present
static int sam_header_line_validate(HeaderLine *hline)
{
	list_t *tags;
	HeaderTag *tag;
	int itype, itag;

	// Is the type correct?
	itype = tag_exists(hline->type, types);
	if (itype == -1)
	{
		debug("The type [%c%c] not recognised.\n", hline->type[0], hline->type[1]);
		return 0;
	}

	// Has all required tags?
	itag = 0;

	while (required_tags[itype] && required_tags[itype][itag])
	{
		if (!header_line_has_tag(hline,required_tags[itype][itag]))
		{
			debug("The tag [%c%c] required for [%c%c] not present.\n", required_tags[itype][itag][0], required_tags[itype][itag][1], hline->type[0], hline->type[1]);
			return 0;
		}

		itag++;
	}

	// Are all tags recognised?
	tags = hline->tags;

	while (tags)
	{
		tag = (HeaderTag*)tags->data;

		if (!tag_exists(tag->key, required_tags[itype]) && !tag_exists(tag->key, optional_tags[itype]))
		{
			// Lower case tags are user-defined values.
			if(!(islower(tag->key[0]) || islower(tag->key[1])))
			{
				// Neither is lower case, but tag was not recognized.
				debug("Unknown tag [%c%c] for [%c%c].\n", tag->key[0],tag->key[1], hline->type[0],hline->type[1]);
				// return 0; // Even unknown tags are allowed - for forward compatibility with new attributes
			}
			// else - allow user defined tag
		}

		tags = tags->next;
	}

	return 1;
}

static void print_header_line(FILE *fp, HeaderLine *hline)
{
	list_t *tags = hline->tags;
	HeaderTag *tag;

	fprintf(fp, "@%c%c", hline->type[0], hline->type[1]);

	while (tags)
	{
		tag = (HeaderTag*)tags->data;
		fprintf(fp, "\t");

		if ((tag->key[0]!=' ') || (tag->key[1]!=' '))
			fprintf(fp, "%c%c:", tag->key[0], tag->key[1]);

		fprintf(fp, "%s", tag->value);
		tags = tags->next;
	}

	fprintf(fp,"\n");
}


static void sam_header_line_free(HeaderLine *hline)
{
	list_t *tags = hline->tags;

	while (tags)
	{
		HeaderTag *tag = (HeaderTag*)tags->data;

		free(tag->value);
		free(tag);
		tags = tags->next;
	}

	list_free(hline->tags);
	free(hline);
}

void sam_header_free(void *_header)
{
	HeaderDict *header = (HeaderDict*)_header;
	list_t *hlines = header;

	while (hlines)
	{
		sam_header_line_free((HeaderLine*)hlines->data);
		hlines = hlines->next;
	}

	list_free(header);
}

HeaderDict *sam_header_clone(const HeaderDict *dict)
{
	HeaderDict *out = NULL;

	while (dict)
	{
		HeaderLine *hline = (HeaderLine*)dict->data;

		out = list_append(out, sam_header_line_clone(hline));
		dict = dict->next;
	}

	return out;
}

// Returns a newly allocated string
char *sam_header_write(const void *_header)
{
	const HeaderDict *header = (const HeaderDict*)_header;
	char *out = NULL;
	int len = 0;
	int nout = 0;
	const list_t *hlines;
	list_t *tags;
	HeaderLine *hline;
	HeaderTag *tag;

	// Calculate the length of the string to allocate
	hlines = header;

	while (hlines)
	{
		// @XY and \n
		len += 4;
		hline = (HeaderLine*)hlines->data;
		tags = hline->tags;

		while (tags)
		{
			tag = (HeaderTag*)tags->data;
			// \t
			len += strlen(tag->value) + 1;

			// XY:
			if ((tag->key[0] != ' ') || (tag->key[1] != ' '))
				len += strlen(tag->value) + 3;

			tags = tags->next;
		}

		hlines = hlines->next;
	}

	nout = 0;
	out = (char*)malloc(len + 1);
	hlines = header;

	while (hlines)
	{
		hline = (HeaderLine*)hlines->data;
		nout += sprintf(out+nout, "@%c%c", hline->type[0], hline->type[1]);
		tags = hline->tags;

		while (tags)
		{
			tag = (HeaderTag*)tags->data;
			nout += sprintf(out+nout, "\t");

			if ((tag->key[0] != ' ') || (tag->key[1] != ' '))
				nout += sprintf(out+nout, "%c%c:", tag->key[0], tag->key[1]);

			nout += sprintf(out+nout, "%s", tag->value);
			tags = tags->next;
		}

		hlines = hlines->next;
		nout += sprintf(out+nout, "\n");
	}

	out[len] = 0;

	return out;
}

void *sam_header_parse2(const char *headerText)
{
	list_t *hlines = NULL;
	HeaderLine *hline;
	const char *text;
	char *buf = NULL;
	size_t nbuf = 0;
	int tovalidate = 0;

	if (!headerText)
		return 0;

	text = headerText;

	while ((text = nextline(&buf, &nbuf, text)))
	{
		hline = sam_header_line_parse(buf);

		if (hline && (!tovalidate || sam_header_line_validate(hline)))
			// With too many (~250,000) reference sequences the header parsing was too slow with list_append.
			hlines = list_append_to_end(hlines, hline);
		else
		{
			if (hline)
				sam_header_line_free(hline);

			sam_header_free(hlines);

			if (buf)
				free(buf);

			return NULL;
		}
	}

	if (buf)
		free(buf);

	return hlines;
}

void *sam_header2tbl(const void *_dict, char type[2], char key_tag[2], char value_tag[2])
{
	const HeaderDict *dict = (const HeaderDict*)_dict;
	const list_t *l = dict;
	khash_t(str) *tbl = kh_init(str);
	khiter_t k;
	int ret;
	HeaderLine *hline;
	HeaderTag *key;
	HeaderTag *value;

	// return an empty (not null) hash table
	if (_dict == 0)
		return tbl;

	while (l)
	{
		hline = (HeaderLine*)l->data;

		if ((hline->type[0] != type[0]) || (hline->type[1]!=type[1]))
		{
			l = l->next;
			continue;
		}

		key = header_line_has_tag(hline, key_tag);
		value = header_line_has_tag(hline, value_tag);

		if (!key || !value)
		{
			l = l->next;
			continue;
		}

		k = kh_get(str, tbl, key->value);

		if (k != kh_end(tbl))
			debug("[sam_header_lookup_table] They key %s not unique.\n", key->value);

		k = kh_put(str, tbl, key->value, &ret);
		kh_value(tbl, k) = value->value;
		l = l->next;
	}

	return tbl;
}

char **sam_header2list(const void *_dict, char type[2], char key_tag[2], int *_n)
{
	const HeaderDict *dict = (const HeaderDict*)_dict;
	const list_t *l = dict;
	int max;
	int n;
	char **ret;
	HeaderLine *hline;
	HeaderTag *key;

	ret = 0;
	*_n = max = n = 0;

	while (l)
	{
		hline = (HeaderLine*)l->data;

		if ((hline->type[0] != type[0]) || (hline->type[1] != type[1]))
		{
			l = l->next;
			continue;
		}

		key   = header_line_has_tag(hline, key_tag);

		if (!key)
		{
			l = l->next;
			continue;
		}

		if (n == max)
		{
			max = max ? max << 1 : 4;
			ret = (char**)realloc(ret, (size_t)max * sizeof(char*));
		}

		ret[n++] = key->value;
		l = l->next;
	}

	*_n = n;

	return ret;
}

const char *sam_tbl_get(void *h, const char *key)
{
	khash_t(str) *tbl = (khash_t(str)*)h;
	khint_t k;

	k = kh_get(str, tbl, key);

	return k == kh_end(tbl) ? 0 : kh_val(tbl, k);
}

int sam_tbl_size(void *h)
{
	khash_t(str) *tbl = (khash_t(str)*)h;

	return h ? kh_size(tbl) : 0;
}

void sam_tbl_destroy(void *h)
{
	khash_t(str) *tbl = (khash_t(str)*)h;
	kh_destroy(str, tbl);
}

void *sam_header_merge(int n, const void **_dicts)
{
	const HeaderDict **dicts = (const HeaderDict**)_dicts;
	HeaderDict *out_dict;
	int idict;
	int status;

	if (n < 2)
		return NULL;

	out_dict = sam_header_clone(dicts[0]);

	for (idict=1; idict < n; idict++)
	{
		const list_t *tmpl_hlines = dicts[idict];

		while (tmpl_hlines)
		{
			list_t *out_hlines = out_dict;
			int inserted = 0;

			while (out_hlines)
			{
				status = sam_header_compare_lines((HeaderLine*)tmpl_hlines->data, (HeaderLine*)out_hlines->data);

				if (status == 0)
				{
					out_hlines = out_hlines->next;
					continue;
				}

				if (status == 2)
				{
					print_header_line(stderr, (HeaderLine*)tmpl_hlines->data);
					print_header_line(stderr, (HeaderLine*)out_hlines->data);
					debug("Conflicting lines, cannot merge the headers.\n");
					return 0;
				}

				if (status == 3)
					sam_header_line_merge_with((HeaderLine*)out_hlines->data, (HeaderLine*)tmpl_hlines->data);

				inserted = 1;
				break;
			}

			if (!inserted)
				out_dict = list_append(out_dict, sam_header_line_clone((HeaderLine*)tmpl_hlines->data));

			tmpl_hlines = tmpl_hlines->next;
		}
	}

	return out_dict;
}