header.c

/*
Copyright (c) 2018-2020, 2023, 2025-2026 Genome Research Ltd.
Authors: James Bonfield <jkb@sanger.ac.uk>, Valeriu Ohan <vo2@sanger.ac.uk>

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

   1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.

   2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.

   3. Neither the names Genome Research Ltd and Wellcome Trust Sanger
Institute nor the names of its contributors may be used to endorse or promote
products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY GENOME RESEARCH LTD AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL GENOME RESEARCH LTD OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#define HTS_BUILDING_LIBRARY // Enables HTSLIB_EXPORT, see htslib/hts_defs.h
#include <config.h>

#include <string.h>
#include <assert.h>
#include <errno.h>
#include "textutils_internal.h"
#include "header.h"
#include "htslib/bgzf.h"
#include "htslib/hfile.h"
#include "htslib/hts_alloc.h"
#include "htslib/kseq.h"

// Hash table for removing multiple lines from the header
KHASH_SET_INIT_STR(rm)
// Used for long refs in SAM files
KHASH_DECLARE(s2i, kh_cstr_t, int64_t)

typedef khash_t(rm) rmhash_t;

static int sam_hdr_link_pg(sam_hdr_t *bh);

static int sam_hrecs_vupdate(sam_hrecs_t *hrecs, sam_hrec_type_t *type, va_list ap);
static int sam_hrecs_update(sam_hrecs_t *hrecs, sam_hrec_type_t *type, ...);


#define MAX_ERROR_QUOTE 320 // Prevent over-long error messages
static void sam_hrecs_error(const char *msg, const char *line, size_t len, size_t lno) {
    int j;

    if (len > MAX_ERROR_QUOTE)
        len = MAX_ERROR_QUOTE;
    for (j = 0; j < len && line[j] != '\n'; j++)
        ;
    hts_log_error("%s at line %zd: \"%.*s\"", msg, lno, j, line);
}

/* ==== Static methods ==== */

static int sam_hrecs_init_type_order(sam_hrecs_t *hrecs, char *type_list) {
    if (!hrecs)
        return -1;

    if (!type_list) {
        hrecs->type_count = 5;
        hrecs->type_order = calloc(hrecs->type_count, 3);
        if (!hrecs->type_order)
            return -1;
        memcpy(hrecs->type_order[0], "HD", 2);
        memcpy(hrecs->type_order[1], "SQ", 2);
        memcpy(hrecs->type_order[2], "RG", 2);
        memcpy(hrecs->type_order[3], "PG", 2);
        memcpy(hrecs->type_order[4], "CO", 2);
    }

    return 0;
}

static int sam_hrecs_add_ref_altnames(sam_hrecs_t *hrecs, int nref, const char *list) {
    const char *token;
    ks_tokaux_t aux;

    if (!list)
        return 0;

    for (token = kstrtok(list, ",", &aux); token; token = kstrtok(NULL, NULL, &aux)) {
        if (aux.p == token)
            continue;

        char *name = string_ndup(hrecs->str_pool, token, aux.p - token);
        if (!name)
            return -1;
        int r;
        khint_t k = kh_put(m_s2i, hrecs->ref_hash, name, &r);
        if (r < 0) return -1;

        if (r > 0)
            kh_val(hrecs->ref_hash, k) = nref;
        else if (kh_val(hrecs->ref_hash, k) != nref)
            hts_log_warning("Duplicate entry AN:\"%s\" in sam header", name);
    }

    return 0;
}

static void sam_hrecs_remove_ref_altnames(sam_hrecs_t *hrecs, int expected, const char *list) {
    const char *token, *sn;
    ks_tokaux_t aux;
    kstring_t str = KS_INITIALIZE;

    if (expected < 0 || expected >= hrecs->nref)
        return;
    sn = hrecs->ref[expected].name;

    for (token = kstrtok(list, ",", &aux); token; token = kstrtok(NULL, NULL, &aux)) {
        kputsn(token, aux.p - token, ks_clear(&str));
        khint_t k = kh_get(m_s2i, hrecs->ref_hash, str.s);
        if (k != kh_end(hrecs->ref_hash)
            && kh_val(hrecs->ref_hash, k) == expected
            && strcmp(sn, str.s) != 0)
            kh_del(m_s2i, hrecs->ref_hash, k);
    }

    free(str.s);
}

/* Updates the hash tables in the sam_hrecs_t structure.
 *
 * Returns 0 on success;
 *        -1 on failure
 */
static int sam_hrecs_update_hashes(sam_hrecs_t *hrecs,
                                   khint32_t type,
                                   sam_hrec_type_t *h_type) {
    /* Add to reference hash? */
    if (type == TYPEKEY("SQ")) {
        sam_hrec_tag_t *tag = h_type->tag;
        int nref = hrecs->nref;
        const char *name = NULL;
        const char *altnames = NULL;
        hts_pos_t len = -1;
        int r, invLN = 0;
        khint_t k;

        while (tag) {
            if (tag->str[0] == 'S' && tag->str[1] == 'N') {
                assert(tag->len >= 3);
                name = tag->str+3;
            } else if (tag->str[0] == 'L' && tag->str[1] == 'N') {
                assert(tag->len >= 3);
                hts_pos_t tmp = len;
                len = strtoll(tag->str+3, NULL, 10);
                if (tmp != -1 && tmp != len) {  //duplicate and different LN
                    invLN = 1;
                }
            } else if (tag->str[0] == 'A' && tag->str[1] == 'N') {
                assert(tag->len >= 3);
                altnames = tag->str+3;
            }
            tag = tag->next;
        }

        if (!name) {
            hts_log_error("Header includes @SQ line with no SN: tag");
            return -1; // SN should be present, according to spec.
        }

        if (len <= 0) {
            hts_log_error("Header includes @SQ line \"%s\" with invalid or no LN: tag",
                          name);
            return -1; // LN should be present, according to spec.
        }

        if (invLN) {
            hts_log_error("Header includes @SQ line \"%s\" with multiple LN:"
            " tag with different values.", name);
            return -1; // LN should not be duplicated or be same
        }

        // Seen already?
        k = kh_get(m_s2i, hrecs->ref_hash, name);
        if (k < kh_end(hrecs->ref_hash)) {
            nref = kh_val(hrecs->ref_hash, k);
            int ref_changed_flag = 0;

            // Check for hash entry added by sam_hrecs_refs_from_targets_array()
            if (hrecs->ref[nref].ty == NULL) {
                // Attach header line to existing stub entry.
                hrecs->ref[nref].ty = h_type;
                // Check lengths match; correct if not.
                if (len != hrecs->ref[nref].len) {
                    char tmp[32];
                    snprintf(tmp, sizeof(tmp), "%" PRIhts_pos,
                             hrecs->ref[nref].len);
                    if (sam_hrecs_update(hrecs, h_type, "LN", tmp, NULL) < 0)
                        return -1;
                    ref_changed_flag = 1;
                }
                if (sam_hrecs_add_ref_altnames(hrecs, nref, altnames) < 0)
                    return -1;

                if (ref_changed_flag && (hrecs->refs_changed < 0 || hrecs->refs_changed > nref))
                    hrecs->refs_changed = nref;
                return 0;
            }

            // Check to see if an existing entry is being updated
            if (hrecs->ref[nref].ty == h_type) {
                if (hrecs->ref[nref].len != len) {
                    hrecs->ref[nref].len = len;
                    ref_changed_flag = 1;
                }
                if (!hrecs->ref[nref].name || strcmp(hrecs->ref[nref].name, name)) {
                    hrecs->ref[nref].name = name;
                    ref_changed_flag = 1;
                }
                if (sam_hrecs_add_ref_altnames(hrecs, nref, altnames) < 0)
                    return -1;

                if (ref_changed_flag && (hrecs->refs_changed < 0 || hrecs->refs_changed > nref))
                    hrecs->refs_changed = nref;
                return 0;
            }

            // If here, the name is a duplicate.
            // Check to see if it matches the SN: tag from the earlier record.
            if (strcmp(hrecs->ref[nref].name, name) == 0) {
                hts_log_error("Duplicate entry \"%s\" in sam header",
                                name);
                return -1;
            }

            // Clash with an already-seen altname
            // As SN: should be preferred to AN: add this as a new
            // record and update the hash entry to point to it.
            hts_log_warning("Ref name SN:\"%s\" is a duplicate of an existing AN key", name);
            nref = hrecs->nref;
        }

        if (nref == hrecs->ref_sz) {
            size_t new_sz = hrecs->ref_sz >= 4 ? hrecs->ref_sz + (hrecs->ref_sz / 4) : 32;
            sam_hrec_sq_t *new_ref = hts_realloc_p(hrecs->ref, sizeof(*hrecs->ref), new_sz);
            if (!new_ref)
                return -1;
            hrecs->ref = new_ref;
            hrecs->ref_sz = new_sz;
        }

        hrecs->ref[nref].name = name;
        hrecs->ref[nref].len  = len;
        hrecs->ref[nref].ty = h_type;

        k = kh_put(m_s2i, hrecs->ref_hash, hrecs->ref[nref].name, &r);
        if (-1 == r) return -1;
        kh_val(hrecs->ref_hash, k) = nref;

        if (sam_hrecs_add_ref_altnames(hrecs, nref, altnames) < 0)
            return -1;

        if (hrecs->refs_changed < 0 || hrecs->refs_changed > hrecs->nref)
            hrecs->refs_changed = hrecs->nref;
        hrecs->nref++;
    }

    /* Add to read-group hash? */
    if (type == TYPEKEY("RG")) {
        sam_hrec_tag_t *tag = sam_hrecs_find_key(h_type, "ID", NULL);
        int nrg = hrecs->nrg, r;
        khint_t k;

        if (!tag) {
            hts_log_error("Header includes @RG line with no ID: tag");
            return -1;  // ID should be present, according to spec.
        }
        assert(tag->str && tag->len >= 3);

        // Seen already?
        k = kh_get(m_s2i, hrecs->rg_hash, tag->str + 3);
        if (k < kh_end(hrecs->rg_hash)) {
            nrg = kh_val(hrecs->rg_hash, k);
            assert(hrecs->rg[nrg].ty != NULL);
            if (hrecs->rg[nrg].ty != h_type) {
                hts_log_warning("Duplicate entry \"%s\" in sam header",
                                tag->str + 3);
            } else {
                hrecs->rg[nrg].name = tag->str + 3;
                hrecs->rg[nrg].name_len = tag->len - 3;
            }
            return 0;
        }

        if (nrg == hrecs->rg_sz) {
            size_t new_sz = hrecs->rg_sz >= 4 ? hrecs->rg_sz + hrecs->rg_sz / 4 : 4;
            sam_hrec_rg_t *new_rg = hts_realloc_p(hrecs->rg, sizeof(*hrecs->rg), new_sz);
            if (!new_rg)
                return -1;
            hrecs->rg = new_rg;
            hrecs->rg_sz = new_sz;
        }

        hrecs->rg[nrg].name = tag->str + 3;
        hrecs->rg[nrg].name_len = tag->len - 3;
        hrecs->rg[nrg].ty   = h_type;
        hrecs->rg[nrg].id   = nrg;

        k = kh_put(m_s2i, hrecs->rg_hash, hrecs->rg[nrg].name, &r);
        if (-1 == r) return -1;
        kh_val(hrecs->rg_hash, k) = nrg;

        hrecs->nrg++;
    }

    /* Add to program hash? */
    if (type == TYPEKEY("PG")) {
        sam_hrec_tag_t *tag;
        sam_hrec_pg_t *new_pg;
        int npg = hrecs->npg;

        if (npg == hrecs->pg_sz) {
            size_t new_sz = hrecs->pg_sz >= 4 ? hrecs->pg_sz + hrecs->pg_sz / 4 : 4;
            new_pg = hts_realloc_p(hrecs->pg, sizeof(*hrecs->pg), new_sz);
            if (!new_pg)
                return -1;
            hrecs->pg = new_pg;
            hrecs->pg_sz = new_sz;
        }

        tag = h_type->tag;
        hrecs->pg[npg].name = NULL;
        hrecs->pg[npg].name_len = 0;
        hrecs->pg[npg].ty  = h_type;
        hrecs->pg[npg].id   = npg;
        hrecs->pg[npg].prev_id = -1;

        while (tag) {
            if (tag->str[0] == 'I' && tag->str[1] == 'D') {
                /* Avoid duplicate ID tags coming from other applications */
                if (!hrecs->pg[npg].name) {
                    assert(tag->len >= 3);
                    hrecs->pg[npg].name = tag->str + 3;
                    hrecs->pg[npg].name_len = tag->len - 3;
                } else {
                    hts_log_warning("PG line with multiple ID tags. The first encountered was preferred - ID:%s", hrecs->pg[npg].name);
                }
            } else if (tag->str[0] == 'P' && tag->str[1] == 'P') {
                // Resolve later if needed
                khint_t k;
                k = kh_get(m_s2i, hrecs->pg_hash, tag->str+3);

                if (k != kh_end(hrecs->pg_hash)) {
                    int p_id = kh_val(hrecs->pg_hash, k);
                    hrecs->pg[npg].prev_id = hrecs->pg[p_id].id;

                    /* Unmark previous entry as a PG termination */
                    if (hrecs->npg_end > 0 &&
                        hrecs->pg_end[hrecs->npg_end-1] == p_id) {
                        hrecs->npg_end--;
                    } else {
                        int i;
                        for (i = 0; i < hrecs->npg_end; i++) {
                            if (hrecs->pg_end[i] == p_id) {
                                memmove(&hrecs->pg_end[i], &hrecs->pg_end[i+1],
                                        (hrecs->npg_end-i-1)*sizeof(*hrecs->pg_end));
                                hrecs->npg_end--;
                            }
                        }
                    }
                } else {
                    hrecs->pg[npg].prev_id = -1;
                }
            }
            tag = tag->next;
        }

        if (hrecs->pg[npg].name) {
            khint_t k;
            int r;
            k = kh_put(m_s2i, hrecs->pg_hash, hrecs->pg[npg].name, &r);
            if (-1 == r) return -1;
            kh_val(hrecs->pg_hash, k) = npg;
        } else {
            return -1; // ID should be present, according to spec.
        }

        /* Add to npg_end[] array. Remove later if we find a PP line */
        if (hrecs->npg_end >= hrecs->npg_end_alloc) {
            int *new_pg_end;
            int  new_alloc = hrecs->npg_end_alloc ? hrecs->npg_end_alloc*2 : 4;

            new_pg_end = hts_realloc_p(hrecs->pg_end, sizeof(int), new_alloc);
            if (!new_pg_end)
                return -1;
            hrecs->npg_end_alloc = new_alloc;
            hrecs->pg_end = new_pg_end;
        }
        hrecs->pg_end[hrecs->npg_end++] = npg;

        hrecs->npg++;
    }

    return 0;
}

static int sam_hrecs_remove_hash_entry(sam_hrecs_t *hrecs, khint32_t type, sam_hrec_type_t *h_type) {
    if (!hrecs || !h_type)
        return -1;

    sam_hrec_tag_t *tag;
    const char *key = NULL;
    khint_t k;

    /* Remove name and any alternative names from reference hash */
    if (type == TYPEKEY("SQ")) {
        const char *altnames = NULL;

        tag = h_type->tag;

        while (tag) {
            if (tag->str[0] == 'S' && tag->str[1] == 'N') {
                assert(tag->len >= 3);
                key = tag->str + 3;
            } else if (tag->str[0] == 'A' && tag->str[1] == 'N') {
                assert(tag->len >= 3);
                altnames = tag->str + 3;
            }
            tag = tag->next;
        }

        if (key) {
            k = kh_get(m_s2i, hrecs->ref_hash, key);
            if (k != kh_end(hrecs->ref_hash)) {
                int idx = kh_val(hrecs->ref_hash, k);
                if (idx + 1 < hrecs->nref)
                    memmove(&hrecs->ref[idx], &hrecs->ref[idx+1],
                            sizeof(sam_hrec_sq_t)*(hrecs->nref - idx - 1));
                if (altnames)
                    sam_hrecs_remove_ref_altnames(hrecs, idx, altnames);
                kh_del(m_s2i, hrecs->ref_hash, k);
                hrecs->nref--;
                if (hrecs->refs_changed < 0 || hrecs->refs_changed > idx)
                    hrecs->refs_changed = idx;
                for (k = 0; k < kh_end(hrecs->ref_hash); k++) {
                    if (kh_exist(hrecs->ref_hash, k)
                        && kh_value(hrecs->ref_hash, k) > idx) {
                        kh_value(hrecs->ref_hash, k)--;
                    }
                }
            }
        }
    }

    /* Remove from read-group hash */
    if (type == TYPEKEY("RG")) {
        tag = h_type->tag;

        while (tag) {
            if (tag->str[0] == 'I' && tag->str[1] == 'D') {
                assert(tag->len >= 3);
                key = tag->str + 3;
                k = kh_get(m_s2i, hrecs->rg_hash, key);
                if (k != kh_end(hrecs->rg_hash)) {
                    int idx = kh_val(hrecs->rg_hash, k);
                    if (idx + 1 < hrecs->nrg)
                        memmove(&hrecs->rg[idx], &hrecs->rg[idx+1], sizeof(sam_hrec_rg_t)*(hrecs->nrg - idx - 1));
                    kh_del(m_s2i, hrecs->rg_hash, k);
                    hrecs->nrg--;
                    for (k = 0; k < kh_end(hrecs->rg_hash); k++) {
                        if (kh_exist(hrecs->rg_hash, k)
                            && kh_value(hrecs->rg_hash, k) > idx) {
                            kh_value(hrecs->rg_hash, k)--;
                        }
                    }
                }
                break;
            }
            tag = tag->next;
        }
    }

    return 0;
}

/** Add a header record to the global line ordering
 *
 * If @p after is not NULL, the new record will be inserted after this one,
 * otherwise it will go at the end.
 *
 * An exception is an HD record, which will always be put first unless
 * one is already present.
 */
static void sam_hrecs_global_list_add(sam_hrecs_t *hrecs,
                                      sam_hrec_type_t *h_type,
                                      sam_hrec_type_t *after) {
    const khint32_t hd_type = TYPEKEY("HD");
    int update_first_line = 0;

    // First line seen
    if (!hrecs->first_line) {
        hrecs->first_line = h_type->global_next = h_type->global_prev = h_type;
        return;
    }

    // @HD goes at the top (unless there's one already)
    if (h_type->type == hd_type && hrecs->first_line->type != hd_type) {
        after = hrecs->first_line->global_prev;
        update_first_line = 1;
    }

    // If no instructions given, put it at the end
    if (!after)
        after = hrecs->first_line->global_prev;

    h_type->global_prev = after;
    h_type->global_next = after->global_next;
    h_type->global_prev->global_next = h_type;
    h_type->global_next->global_prev = h_type;

    if (update_first_line)
        hrecs->first_line = h_type;
}

/*! Add header record with a va_list interface.
 *
 * Adds a single record to a SAM header.
 *
 * This takes a header record type, a va_list argument and one or more
 * key,value pairs, ending with the NULL key.
 *
 * Eg. sam_hrecs_vadd(h, "SQ", args, "ID", "foo", "LN", "100", NULL).
 *
 * The purpose of the additional va_list parameter is to permit other
 * varargs functions to call this while including their own additional
 * parameters; an example is in sam_hdr_add_pg().
 *
 * Note: this function invokes va_arg at least once, making the value
 * of ap indeterminate after the return. The caller should call
 * va_start/va_end before/after calling this function or use va_copy.
 *
 * @return
 * Returns >= 0 on success;
 *        -1 on failure
 */
static int sam_hrecs_vadd(sam_hrecs_t *hrecs, const char *type, va_list ap, ...) {
    va_list args;
    sam_hrec_type_t *h_type;
    sam_hrec_tag_t *h_tag, *last=NULL;
    int new;
    khint32_t type_i = TYPEKEY(type), k;

    if (!strncmp(type, "HD", 2) && (h_type = sam_hrecs_find_type_id(hrecs, "HD", NULL, NULL)))
        return sam_hrecs_vupdate(hrecs, h_type, ap);

    if (!(h_type = pool_alloc(hrecs->type_pool)))
        return -1;
    k = kh_put(sam_hrecs_t, hrecs->h, type_i, &new);
    if (new < 0)
        return -1;

    h_type->type = type_i;

    // Form the ring, either with self or other lines of this type
    if (!new) {
        sam_hrec_type_t *t = kh_val(hrecs->h, k), *p;
        p = t->prev;

        assert(p->next == t);
        p->next = h_type;
        h_type->prev = p;

        t->prev = h_type;
        h_type->next = t;
    } else {
        kh_val(hrecs->h, k) = h_type;
        h_type->prev = h_type->next = h_type;
    }
    h_type->tag = NULL;

    // Add to global line ordering after any existing line of the same type,
    // or at the end if no line of this type exists yet.
    sam_hrecs_global_list_add(hrecs, h_type, !new ? h_type->prev : NULL);

    // Check linked-list invariants
    assert(h_type->prev->next == h_type);
    assert(h_type->next->prev == h_type);
    assert(h_type->global_prev->global_next == h_type);
    assert(h_type->global_next->global_prev == h_type);

    // Any ... varargs
    va_start(args, ap);
    for (;;) {
        char *key, *val = NULL, *str;

        if (!(key = (char *)va_arg(args, char *)))
            break;
        if (strncmp(type, "CO", 2) && !(val = (char *)va_arg(args, char *)))
            break;
        if (*val == '\0')
            continue;

        if (!(h_tag = pool_alloc(hrecs->tag_pool)))
            return -1;

        if (strncmp(type, "CO", 2)) {
            h_tag->len = 3 + strlen(val);
            str = string_alloc(hrecs->str_pool, hts_add_sat2(h_tag->len, 1));
            if (!str || snprintf(str, h_tag->len+1, "%2.2s:%s", key, val) < 0)
                return -1;
            h_tag->str = str;
        } else {
            h_tag->len = strlen(key);
            h_tag->str = string_ndup(hrecs->str_pool, key, h_tag->len);
            if (!h_tag->str)
                return -1;
        }

        h_tag->next = NULL;
        if (last)
            last->next = h_tag;
        else
            h_type->tag = h_tag;

        last = h_tag;
    }
    va_end(args);

    // Plus the specified va_list params
    for (;;) {
        char *key, *val = NULL, *str;

        if (!(key = (char *)va_arg(ap, char *)))
            break;
        if (strncmp(type, "CO", 2) && !(val = (char *)va_arg(ap, char *)))
            break;

        if (!(h_tag = pool_alloc(hrecs->tag_pool)))
            return -1;

        if (strncmp(type, "CO", 2)) {
            h_tag->len = 3 + strlen(val);
            str = string_alloc(hrecs->str_pool, hts_add_sat2(h_tag->len, 1));
            if (!str || snprintf(str, h_tag->len+1, "%2.2s:%s", key, val) < 0)
                return -1;
            h_tag->str = str;
        } else {
            h_tag->len = strlen(key);
            h_tag->str = string_ndup(hrecs->str_pool, key, h_tag->len);
            if (!h_tag->str)
                return -1;
        }

        h_tag->next = NULL;
        if (last)
            last->next = h_tag;
        else
            h_type->tag = h_tag;

        last = h_tag;
    }

    if (-1 == sam_hrecs_update_hashes(hrecs, TYPEKEY(type), h_type))
        return -1;

    if (!strncmp(type, "PG", 2))
        hrecs->pgs_changed = 1;

    hrecs->dirty = 1;

    return 0;
}

// As sam_hrecs_vadd(), but without the extra va_list parameter
static int sam_hrecs_add(sam_hrecs_t *hrecs, const char *type, ...) {
    va_list args;
    int res;
    va_start(args, type);
    res = sam_hrecs_vadd(hrecs, type, args, NULL);
    va_end(args);
    return res;
}

/*
 * Function for deallocating a list of tags
 */

static void sam_hrecs_free_tags(sam_hrecs_t *hrecs, sam_hrec_tag_t *tag) {
    if (!hrecs || !tag)
        return;
    if (tag->next)
        sam_hrecs_free_tags(hrecs, tag->next);

    pool_free(hrecs->tag_pool, tag);
}

static int sam_hrecs_remove_line(sam_hrecs_t *hrecs, const char *type_name, sam_hrec_type_t *type_found, int remove_hash) {
    if (!hrecs || !type_name || !type_found)
        return -1;

    khint32_t itype = TYPEKEY(type_name);
    khint_t k = kh_get(sam_hrecs_t, hrecs->h, itype);
    if (k == kh_end(hrecs->h))
        return -1;

    // Remove from global list (remembering it could be the only line)
    if (hrecs->first_line == type_found) {
        hrecs->first_line = (type_found->global_next != type_found
                             ? type_found->global_next : NULL);
    }
    type_found->global_next->global_prev = type_found->global_prev;
    type_found->global_prev->global_next = type_found->global_next;

    /* single element in the list */
    if (type_found->prev == type_found || type_found->next == type_found) {
        kh_del(sam_hrecs_t, hrecs->h, k);
    } else {
        type_found->prev->next = type_found->next;
        type_found->next->prev = type_found->prev;
        if (kh_val(hrecs->h, k) == type_found) { //first element
            kh_val(hrecs->h, k) = type_found->next;
        }
    }

    if (remove_hash && (!strncmp(type_name, "SQ", 2) || !strncmp(type_name, "RG", 2)))
        sam_hrecs_remove_hash_entry(hrecs, itype, type_found);

    sam_hrecs_free_tags(hrecs, type_found->tag);
    pool_free(hrecs->type_pool, type_found);

    hrecs->dirty = 1;

    return 0;
}

// Paste together a line from the parsed data structures
static int build_header_line(const sam_hrec_type_t *ty, kstring_t *ks) {
    sam_hrec_tag_t *tag;
    int r = 0;
    char c[2]= { ty->type >> 8, ty->type & 0xff };

    r |= (kputc_('@', ks) == EOF);
    r |= (kputsn(c, 2, ks) == EOF);
    for (tag = ty->tag; tag; tag = tag->next) {
        r |= (kputc_('\t', ks) == EOF);
        r |= (kputsn(tag->str, tag->len, ks) == EOF);
    }

    return r;
}

static int sam_hrecs_rebuild_lines(const sam_hrecs_t *hrecs, kstring_t *ks) {
    const sam_hrec_type_t *t1, *t2;

    if (!hrecs->first_line)
        return kputsn("", 0, ks) >= 0 ? 0 : -1;

    t1 = t2 = hrecs->first_line;
    do {
        if (build_header_line(t1, ks) != 0)
            return -1;
        if (kputc('\n', ks) < 0)
            return -1;

        t1 = t1->global_next;
    } while (t1 != t2);

    return 0;
}

// Print warnings about known log messages that may have been inadvertently
// merged into SAM files.

static void known_stderr(const char *tool, const char *advice) {
    hts_log_warning("SAM file corrupted by embedded %s error/log message", tool);
    hts_log_warning("%s", advice);
}

// Look for known log messages that may have been inadvertently merged into SAM
// files, and print advice on how to avoid the problem if one is spotted.

static void warn_if_known_stderr(const char *line, size_t len) {
    int ilen = len < INT_MAX ? len : INT_MAX;
    if (kmemmem(line, ilen, "M::bwa_idx_load_from_disk", 25, NULL) != NULL)
        known_stderr("bwa", "Use `bwa mem -o file.sam ...` or `bwa sampe -f file.sam ...` instead of `bwa ... > file.sam`");
    else if (kmemmem(line, ilen, "M::mem_pestat", 13, NULL) != NULL)
        known_stderr("bwa", "Use `bwa mem -o file.sam ...` instead of `bwa mem ... > file.sam`");
    else if (kmemmem(line, ilen, "loaded/built the index", 22, NULL) != NULL)
        known_stderr("minimap2", "Use `minimap2 -o file.sam ...` instead of `minimap2 ... > file.sam`");
}

// Parse a CO line
// Returns 0 on failure; line length (>= 1) on success
static size_t parse_comment_line(const char *hdr, size_t len, size_t lno,
                                 sam_hrecs_t *hrecs, sam_hrec_type_t *h_type) {
    size_t j;
    sam_hrec_tag_t *h_tag;

    if (len == 0 || hdr[0] != '\t') {
        sam_hrecs_error("Missing tab", hdr, len, lno);
        warn_if_known_stderr(hdr, len);
        return 0;
    }

    for (j = 1; j < len && hdr[j] != '\0' && hdr[j] != '\n'; j++)
        ;

    if (!(h_type->tag = h_tag = pool_alloc(hrecs->tag_pool)))
        return -1;
    h_tag->str = string_ndup(hrecs->str_pool, hdr + 1, j - 1);
    h_tag->len = j - 1;
    h_tag->next = NULL;
    if (!h_tag->str)
        return 0;
    return j;
}

// Parse header lines other than CO
// Returns 0 on failure; line length (>= 1) on success
static size_t parse_noncomment_line(const char *hdr, size_t len, size_t lno,
                                    sam_hrecs_t *hrecs, sam_hrec_type_t *h_type,
                                    khint_t *sq_count, khint_t *an_count) {
    size_t i = 0;
    sam_hrec_tag_t *h_tag, *last = NULL;
    int is_sq = h_type->type == TYPEKEY("SQ");

    // Count reference name
    if (is_sq)
        (*sq_count)++;

    do {
        size_t j;

        if (i == len || hdr[i] != '\t') {
            sam_hrecs_error("Missing tab", hdr, len, lno);
            warn_if_known_stderr(hdr, len);
            return 0;
        }

        for (j = ++i;
             j < len && hdr[j] != '\0' && hdr[j] != '\n' && hdr[j] != '\t';
             j++)
            ;

        if (j - i < 3 || hdr[i + 2] != ':') {
            sam_hrecs_error("Malformed key:value pair", hdr, len, lno);
            warn_if_known_stderr(hdr, len);
            return 0;
        }

        // Count alternate names
        if (is_sq && hdr[i] == 'A' && hdr[i + 1] == 'N') {
            size_t k;
            (*an_count)++;
            for (k = i + 2; k < j; k++)
                (*an_count) += hdr[k] == ',';
        }

        if (!(h_tag = pool_alloc(hrecs->tag_pool)))
            return 0;
        h_tag->str = string_ndup(hrecs->str_pool, &hdr[i], j-i);
        h_tag->len = j-i;
        h_tag->next = NULL;
        if (!h_tag->str)
            return 0;

        if (last)
            last->next = h_tag;
        else
            h_type->tag = h_tag;

        last = h_tag;
        i = j;
    } while (i < len && hdr[i] != '\0' && hdr[i] != '\n');

    return i;
}

/** Parse a single header line.
 * @param         hrecs        Destination sam_hrecs_t struct
 * @param         hdr          Data to parse
 * @param         len          Size of @p hdr
 * @param         lno          Line number, for error / warning messages
 * @param[out]    line_len_out Returned length of parsed line
 * @param[in,out] sq_count     Count of SQ lines seen so far
 * @param[in,out] an_count     Count of alternate names seen so far
 *
 * Used by sam_hrecs_parse_lines() and sam_hdr_build_from_sam_file(),
 * so has to be able to deal with lines that are not NUL terminated.
 * A single line pointed to by @ph hdr will be parsed into @p hrecs.
 * @p len gives the size of @p hdr, which may include more than one line
 * if called from sam_hrecs_parse_lines(). The length of the parsed line
 * is returned in @p line_len_out.
 * @p sq_count and @p an_count are pointers to the number of SQ lines
 * (i.e. references), and the number of alternate names.  These are used
 * by the caller to set the sizes of the sam_hrecs_t::ref array and
 * sam_hrecs_t::ref_hash table once all the lines have been parsed.
 */

static sam_hrec_type_t * sam_hrecs_parse_single_line(sam_hrecs_t *hrecs,
                                                     const char *hdr,
                                                     size_t len, size_t lno,
                                                     size_t *line_len_out,
                                                     khint_t *sq_count,
                                                     khint_t *an_count) {
    int new;
    khint32_t type;
    khint_t k;
    sam_hrec_type_t *h_type;
    size_t line_len = 0;

    if (hdr[0] != '@') {
        sam_hrecs_error("Header line does not start with '@'",
                        hdr, len, lno);
        return NULL;
    }

    if (!isalpha_c(hdr[1]) || !isalpha_c(hdr[2])) {
        sam_hrecs_error("Header line does not have a two character key",
                        hdr, len, lno);
        return NULL;
    }
    type = TYPEKEY(&hdr[1]);

    if (len < 3 || hdr[3] == '\n')
        return NULL;

    // Add the header line type
    if (!(h_type = pool_alloc(hrecs->type_pool)))
        return NULL;

    h_type->tag = NULL;

    k = kh_put(sam_hrecs_t, hrecs->h, type, &new);
    if (new < 0)
        goto error;

    h_type->type = type;

    // Parse the tags on this line
    if (type == TYPEKEY("CO")) {

        line_len = parse_comment_line(hdr + 3, len - 3, lno,
                                      hrecs, h_type);
        if (line_len == 0) // Failed
            goto error;

    } else {

        line_len = parse_noncomment_line(hdr + 3, len - 3, lno,
                                         hrecs, h_type,
                                         sq_count, an_count);
        if (line_len == 0) // Failed
            goto error;
    }

    // Add to end of global list
    sam_hrecs_global_list_add(hrecs, h_type, NULL);

    // Form the ring, either with self or other lines of this type
    if (!new) {
        sam_hrec_type_t *t = kh_val(hrecs->h, k), *p;
        p = t->prev;

        assert(p->next == t);
        p->next = h_type;
        h_type->prev = p;

        t->prev = h_type;
        h_type->next = t;
    } else {
        kh_val(hrecs->h, k) = h_type;
        h_type->prev = h_type->next = h_type;
    }

    if (line_len_out)
        *line_len_out = line_len + 3;
    return h_type;

 error:
    // Tidy up as best we can
    sam_hrecs_free_tags(hrecs, h_type->tag);
    pool_free(hrecs->type_pool, h_type);
    return NULL;
}

// Parse complete header into a sam_hrecs_t structure

static int sam_hrecs_parse_lines(sam_hrecs_t *hrecs, const char *hdr, size_t len) {
    size_t i, lno;

    if (!hrecs || len > SSIZE_MAX)
        return -1;

    if (!len)
        len = strlen(hdr);

    if (len < 3) {
        if (len == 0 || *hdr == '\0') return 0;
        sam_hrecs_error("Header line too short", hdr, len, 1);
        return -1;
    }

    sam_hrec_type_t *first_h_type = NULL;
    sam_hrec_type_t *last_h_type = NULL;
    khint_t sq_count = 0, an_count = 0;
    size_t line_len = 0;

    for (i = 0, lno = 1; i < len - 3 && hdr[i] != '\0'; lno++) {
        sam_hrec_type_t *h_type;
        h_type = sam_hrecs_parse_single_line(hrecs, &hdr[i], len - i,
                                             lno, &line_len,
                                             &sq_count, &an_count);
        if (!h_type)
            return -1;

        i += line_len + 1;

        // Keep track of first / last new record so they can be passed to
        // sam_hrecs_update_hashes() below.  Ignore HD as it may be stored
        // out of order (it always goes at the start) and
        // sam_hrecs_update_hashes() does nothing with HD.
        if (h_type->type != TYPEKEY("HD")) {
            if (!first_h_type)
                first_h_type = h_type;
            last_h_type = h_type;
        }
    }

    // Reserve space in hash tables and refs array
    if (sq_count) {
        // It's possible that some entries are already in the hash table,
        // e.g. via sam_hrecs_refs_from_targets_array(), so only reserve
        // space beyond what's already there.
        khint_t ref_name_count = sq_count + an_count;
        khint_t sz = kh_size(hrecs->ref_hash);
        if (ref_name_count > sz
            && kh_grow_to_fit(m_s2i, hrecs->ref_hash, ref_name_count - sz) < 0)
            return -1;

        if (hrecs->ref_sz < sq_count) {
            sam_hrec_sq_t *new_ref = hts_realloc_p(hrecs->ref,
                                                   sizeof(*hrecs->ref), sq_count);
            if (!new_ref)
                return -1;
            hrecs->ref = new_ref;
            hrecs->ref_sz = sq_count;
        }
    }

    // Update RG/SQ/PG hashes
    if (first_h_type) {
        sam_hrec_type_t *h_type = first_h_type->global_prev;
        do {
            h_type = h_type->global_next;
            if (-1 == sam_hrecs_update_hashes(hrecs, h_type->type, h_type))
                return -1;
        } while (h_type != last_h_type);
    }
    return 0;
}

/*! Update sam_hdr_t target_name and target_len arrays
 *
 *  @return 0 on success; -1 on failure
 */
int sam_hdr_update_target_arrays(sam_hdr_t *bh, const sam_hrecs_t *hrecs,
                                 int refs_changed) {
    if (!bh || !hrecs)
        return -1;

    if (refs_changed < 0)
        return 0;

    // Grow arrays if necessary
    if (bh->n_targets < hrecs->nref) {
        char **new_names = hts_realloc_p(bh->target_name, sizeof(*new_names),
                                         hrecs->nref);
        if (!new_names)
            return -1;
        bh->target_name = new_names;
        uint32_t *new_lens = hts_realloc_p(bh->target_len, sizeof(*new_lens),
                                           hrecs->nref);
        if (!new_lens)
            return -1;
        bh->target_len = new_lens;
    }

    // Update names and lengths where changed
    // hrecs->refs_changed is the first ref that has been updated, so ones
    // before that can be skipped.
    int i;
    khint_t k;
    khash_t(s2i) *long_refs = (khash_t(s2i) *) bh->sdict;
    for (i = refs_changed; i < hrecs->nref; i++) {
        if (i >= bh->n_targets
            || strcmp(bh->target_name[i], hrecs->ref[i].name) != 0) {
            if (i < bh->n_targets)
                free(bh->target_name[i]);
            bh->target_name[i] = strdup(hrecs->ref[i].name);
            if (!bh->target_name[i])
                return -1;
        }
        if (hrecs->ref[i].len < UINT32_MAX) {
            bh->target_len[i] = hrecs->ref[i].len;

            if (!long_refs)
                continue;

            // Check if we have an old length, if so remove it.
            k = kh_get(s2i, long_refs, bh->target_name[i]);
            if (k < kh_end(long_refs))
                kh_del(s2i, long_refs, k);
        } else {
            bh->target_len[i] = UINT32_MAX;
            if (bh->hrecs != hrecs) {
                // Called from sam_hdr_dup; need to add sdict entries
                if (!long_refs) {
                    if (!(bh->sdict = long_refs = kh_init(s2i)))
                        return -1;
                }

                // Add / update length
                int absent;
                k = kh_put(s2i, long_refs, bh->target_name[i], &absent);
                if (absent < 0)
                    return -1;
                kh_val(long_refs, k) = hrecs->ref[i].len;
            }
        }
    }

    // Free up any names that have been removed
    for (; i < bh->n_targets; i++) {
        if (long_refs) {
            k = kh_get(s2i, long_refs, bh->target_name[i]);
            if (k < kh_end(long_refs))
                kh_del(s2i, long_refs, k);
        }
        free(bh->target_name[i]);
    }

    bh->n_targets = hrecs->nref;
    return 0;
}

static int rebuild_target_arrays(sam_hdr_t *bh) {
    if (!bh || !bh->hrecs)
        return -1;

    sam_hrecs_t *hrecs = bh->hrecs;
    if (hrecs->refs_changed < 0)
        return 0;

    if (sam_hdr_update_target_arrays(bh, hrecs, hrecs->refs_changed) != 0)
        return -1;

    hrecs->refs_changed = -1;
    return 0;
}

/// Populate hrecs refs array from header target_name, target_len arrays
/**
 * @return 0 on success; -1 on failure
 *
 * Pre-fills the refs hash from the target arrays.  For BAM files this
 * will ensure that they are in the correct order as the target arrays
 * are the canonical source for converting target ids to names and lengths.
 *
 * The added entries do not link to a header line. sam_hrecs_update_hashes()
 * will add the links later for lines found in the text header.
 *
 * This should be called before the text header is parsed.
 */
static int sam_hrecs_refs_from_targets_array(sam_hrecs_t *hrecs,
                                             const sam_hdr_t *bh) {
    int32_t tid = 0;

    if (!hrecs || !bh)
        return -1;

    // This should always be called before parsing the text header
    // so the ref array should start off empty, and we don't have to try
    // to reconcile any existing data.
    if (hrecs->nref > 0) {
        hts_log_error("Called with non-empty ref array");
        return -1;
    }

    if (hrecs->ref_sz < bh->n_targets) {
        sam_hrec_sq_t *new_ref = hts_realloc_p(hrecs->ref, sizeof(*new_ref),
                                               bh->n_targets);
        if (!new_ref)
            return -1;

        hrecs->ref = new_ref;
        hrecs->ref_sz = bh->n_targets;
    }

    if (kh_grow_to_fit(m_s2i, hrecs->ref_hash, bh->n_targets) != 0)
        goto fail;

    for (tid = 0; tid < bh->n_targets; tid++) {
        khint_t k;
        int r;
        hrecs->ref[tid].name = string_dup(hrecs->str_pool, bh->target_name[tid]);
        if (!hrecs->ref[tid].name) goto fail;
        if (bh->target_len[tid] < UINT32_MAX || !bh->sdict) {
            hrecs->ref[tid].len  = bh->target_len[tid];
        } else {
            khash_t(s2i) *long_refs = (khash_t(s2i) *) bh->sdict;
            k = kh_get(s2i, long_refs, hrecs->ref[tid].name);
            if (k < kh_end(long_refs)) {
                hrecs->ref[tid].len = kh_val(long_refs, k);
            } else {
                hrecs->ref[tid].len = UINT32_MAX;
            }
        }
        hrecs->ref[tid].ty   = NULL;
        k = kh_put(m_s2i, hrecs->ref_hash, hrecs->ref[tid].name, &r);
        if (r < 0) goto fail;
        if (r == 0) {
            hts_log_error("Duplicate entry \"%s\" in target list",
                            hrecs->ref[tid].name);
            return -1;
        } else {
            kh_val(hrecs->ref_hash, k) = tid;
        }
    }
    hrecs->nref = bh->n_targets;
    return 0;

 fail: {
        int32_t i;
        hts_log_error("%s", strerror(errno));
        for (i = 0; i < tid; i++) {
            khint_t k;
            if (!hrecs->ref[i].name) continue;
            k = kh_get(m_s2i, hrecs->ref_hash, hrecs->ref[tid].name);
            if (k < kh_end(hrecs->ref_hash)) kh_del(m_s2i, hrecs->ref_hash, k);
        }
        hrecs->nref = 0;
        return -1;
    }
}

/*
 * Add SQ header records for any references in the hrecs->ref array that
 * were added by sam_hrecs_refs_from_targets_array() but have not
 * been linked to an @SQ line by sam_hrecs_update_hashes() yet.
 *
 * This may be needed either because:
 *
 *   - A bam file was read that had entries in its refs list with no
 *     corresponding @SQ line.
 *
 *   - A program constructed a sam_hdr_t which has target_name and target_len
 *     array entries with no corresponding @SQ line in text.
 */
static int add_stub_ref_sq_lines(sam_hrecs_t *hrecs) {
    int tid;
    char len[32];

    for (tid = 0; tid < hrecs->nref; tid++) {
        if (hrecs->ref[tid].ty == NULL) {
            snprintf(len, sizeof(len), "%"PRIhts_pos, hrecs->ref[tid].len);
            if (sam_hrecs_add(hrecs, "SQ",
                              "SN", hrecs->ref[tid].name,
                              "LN", len, NULL) != 0)
                return -1;

            // Check that the stub has actually been filled
            if(hrecs->ref[tid].ty == NULL) {
                hts_log_error("Reference stub with tid=%d, name=\"%s\", len=%"PRIhts_pos" could not be filled",
                        tid, hrecs->ref[tid].name, hrecs->ref[tid].len);
                return -1;
            }
        }
    }
    return 0;
}

int sam_hdr_fill_hrecs(sam_hdr_t *bh) {
    sam_hrecs_t *hrecs = sam_hrecs_new();

    if (!hrecs)
        return -1;

    if (bh->target_name && bh->target_len && bh->n_targets > 0) {
        if (sam_hrecs_refs_from_targets_array(hrecs, bh) != 0) {
            sam_hrecs_free(hrecs);
            return -1;
        }
    }

    // Parse existing header text
    if (bh->text && bh->l_text > 0) {
        if (sam_hrecs_parse_lines(hrecs, bh->text, bh->l_text) != 0) {
            sam_hrecs_free(hrecs);
            return -1;
        }
    }

    if (add_stub_ref_sq_lines(hrecs) < 0) {
        sam_hrecs_free(hrecs);
        return -1;
    }

    bh->hrecs = hrecs;

    if (hrecs->refs_changed >= 0 && rebuild_target_arrays(bh) != 0)
        return -1;

    return 0;
}

// Check header type is valid according to the SAM specification

static int valid_sam_header_type(const char *s) {
    if (s[0] != '@') return 0;
    switch (s[1]) {
    case 'H':
        return s[2] == 'D' && s[3] == '\t';
    case 'S':
        return s[2] == 'Q' && s[3] == '\t';
    case 'R':
    case 'P':
        return s[2] == 'G' && s[3] == '\t';
    case 'C':
        return s[2] == 'O';
    }
    return 0;
}

/** Populate sam_hdr_t from SAM file content
 *
 * @param hdr empty header struct to be filled
 * @param fp  File to read from
 * @return 0 on success
 *        -1 on failure
 *
 * This function is used to build the header structure when reading SAM files.
 * The lines read from the file are used to create sam_hrecs_t structures.
 * The function also populates sam_hdr_t::text, sam_hdr_t::l_text,
 * sam_hdr_t::target_name and sam_hdr_t::target_len.
 */
int sam_hdr_build_from_sam_file(sam_hdr_t *hdr, htsFile* fp) {
    kstring_t str = { 0, 0, NULL };
    int ret;
    int next_c = '@';
    sam_hrec_type_t *first_h_type = NULL;
    sam_hrec_type_t *last_h_type = NULL;
    sam_hrecs_t *hrecs = sam_hrecs_new();
    size_t lno = 0;
    khint_t sq_count = 0, an_count = 0;

    if (!hrecs)
        return -1;

    hdr->hrecs = hrecs;

    while (next_c == '@' && (ret = hts_getline(fp, KS_SEP_LINE, &fp->line)) >= 0) {
        sam_hrec_type_t *h_type;

        lno++;
        if (!fp->line.s || fp->line.s[0] != '@')
            break;

        if (fp->line.l < 3
            || !isalpha_c(fp->line.s[1]) || !isalpha_c(fp->line.s[2])) {
            sam_hrecs_error("Header line does not have a two character key",
                            fp->line.s, fp->line.l, lno);
            goto error;
        }

        if (fp->line.l == 3 && strncmp(fp->line.s, "@CO", 3) == 0) {
            // HTSlib's original SAM reader let you get away with an untabbed
            // @CO line
            kputc('\t', &fp->line);
        }

        if (!valid_sam_header_type(fp->line.s)) {
            hts_log_error("Invalid header line: must start with @HD/@SQ/@RG/@PG/@CO");
            warn_if_known_stderr(fp->line.s, fp->line.l);
            goto error;
        }

        h_type = sam_hrecs_parse_single_line(hrecs, fp->line.s, fp->line.l,
                                             lno, NULL, &sq_count, &an_count);
        if (!h_type)
            goto error;

        if (h_type->type != TYPEKEY("HD")) {
            if (!first_h_type)
                first_h_type = h_type;
            last_h_type = h_type;
        }

        if (kputsn(fp->line.s, fp->line.l, &str) < 0)
            goto error;

        if (kputc('\n', &str) < 0)
            goto error;

        if (fp->is_bgzf) {
            next_c = bgzf_peek(fp->fp.bgzf);
        } else {
            unsigned char nc;
            ssize_t pret = hpeek(fp->fp.hfile, &nc, 1);
            next_c = pret > 0 ? nc : pret - 1;
        }
        if (next_c < -1)
            goto error;
    }
    if (next_c != '@')
        fp->line.l = 0;

    if (ret < -1)
        goto error;

    if (sq_count == 0 && fp->fn_aux) {
        // No @SQ lines found.  Try to get them from the fasta index.
        kstring_t line = { 0, 0, NULL };

        /* The reference index (.fai) is actually needed here */
        char *fai_fn = fp->fn_aux;
        char *fn_delim = strstr(fp->fn_aux, HTS_IDX_DELIM);
        if (fn_delim)
            fai_fn = fn_delim + strlen(HTS_IDX_DELIM);

        hFILE* f = hopen(fai_fn, "r");
        int e = 0;
        if (f == NULL)
            goto error;

        while (line.l = 0, khgetline(&line, f) >= 0) {
            char* tab = strchr(line.s, '\t');
            size_t l_start = str.l;
            sam_hrec_type_t *h_type;
            hts_pos_t ln;

            if (tab == NULL)
                continue;

            ln = strtoll(tab, NULL, 10);

            e |= kputs("@SQ\tSN:", &str) < 0;
            e |= kputsn(line.s, tab - line.s, &str) < 0;
            e |= kputs("\tLN:", &str) < 0;
            e |= kputll(ln, &str) < 0;
            e |= kputc('\n', &str) < 0;
            if (e)
                break;

            h_type = sam_hrecs_parse_single_line(hrecs, &str.s[l_start],
                                                 str.l - l_start,
                                                 lno, NULL,
                                                 &sq_count, &an_count);
            if (!h_type)
                goto error;

            if (!first_h_type)
                first_h_type = h_type;
            last_h_type = h_type;
        }

        ks_free(&line);
        if (hclose(f) != 0) {
            hts_log_error("Error on closing %s", fai_fn);
            e = 1;
        }
        if (e)
            goto error;
    }

    // Reserve space in hash tables
    if (sq_count) {
        // Hash table should be empty at the moment
        if (kh_grow_to_fit(m_s2i, hrecs->ref_hash, sq_count + an_count) < 0)
            goto error;

        if (hrecs->ref_sz < sq_count) {
            sam_hrec_sq_t *new_ref = hts_realloc_p(hrecs->ref,
                                                   sizeof(*hrecs->ref), sq_count);
            if (!new_ref)
                goto error;
            hrecs->ref = new_ref;
            hrecs->ref_sz = sq_count;
        }
    }

    // Update RG/SQ/PG hashes
    if (first_h_type) {
        sam_hrec_type_t *h_type = first_h_type->global_prev;
        do {
            h_type = h_type->global_next;
            if (-1 == sam_hrecs_update_hashes(hrecs, h_type->type, h_type))
                goto error;
        } while (h_type != last_h_type);
    }

    if (str.l == 0)
        kputsn("", 0, &str);
    hdr->l_text = str.l;
    hdr->text = ks_release(&str);

    if (sam_hdr_update_target_arrays(hdr, hrecs, 0) < 0)
        goto error;

    return 0;

 error:
    ks_free(&str);
    return -1;
}

/** Remove outdated header text

    @param bh     BAM header

    This is called when API functions have changed the header so that the
    text version is no longer valid.
 */
static void redact_header_text(sam_hdr_t *bh) {
    assert(bh->hrecs && bh->hrecs->dirty);
    bh->l_text = 0;
    free(bh->text);
    bh->text = NULL;
}

/** Find nth header record of a given type

    @param type   Header type (SQ, RG etc.)
    @param idx    0-based index

    @return sam_hrec_type_t pointer to the record on success
            NULL if no record exists with the given type and index
 */

static sam_hrec_type_t *sam_hrecs_find_type_pos(sam_hrecs_t *hrecs,
                                                const char *type, int idx) {
    sam_hrec_type_t *first, *itr;

    if (idx < 0)
        return NULL;

    if (type[0] == 'S' && type[1] == 'Q')
        return idx < hrecs->nref ? hrecs->ref[idx].ty : NULL;

    if (type[0] == 'R' && type[1] == 'G')
        return idx < hrecs->nrg ? hrecs->rg[idx].ty : NULL;

    if (type[0] == 'P' && type[1] == 'G')
        return idx < hrecs->npg ? hrecs->pg[idx].ty : NULL;

    first = itr = sam_hrecs_find_type_id(hrecs, type, NULL, NULL);
    if (!first)
        return NULL;

    while (idx > 0) {
        itr = itr->next;
        if (itr == first)
            break;
        --idx;
    }

    return idx == 0 ? itr : NULL;
}

/* ==== Public methods ==== */

size_t sam_hdr_length(sam_hdr_t *bh) {
    if (!bh || -1 == sam_hdr_rebuild(bh))
        return SIZE_MAX;

    return bh->l_text;
}

const char *sam_hdr_str(sam_hdr_t *bh) {
    if (!bh || -1 == sam_hdr_rebuild(bh))
        return NULL;

    return bh->text;
}

int sam_hdr_nref(const sam_hdr_t *bh) {
    if (!bh)
        return -1;

    return bh->hrecs ? bh->hrecs->nref : bh->n_targets;
}

/*
 * Reconstructs the text representation from the header hash table.
 * Returns 0 on success
 *        -1 on failure
 */
int sam_hdr_rebuild(sam_hdr_t *bh) {
    sam_hrecs_t *hrecs;
    if (!bh)
        return -1;

    if (!(hrecs = bh->hrecs))
        return bh->text ? 0 : -1;

    if (hrecs->refs_changed >= 0) {
        if (rebuild_target_arrays(bh) < 0) {
            hts_log_error("Header target array rebuild has failed");
            return -1;
        }
    }

    /* If header text wasn't changed or header is empty, don't rebuild it. */
    if (!hrecs->dirty)
        return 0;

    if (hrecs->pgs_changed && sam_hdr_link_pg(bh) < 0) {
        hts_log_error("Linking @PG lines has failed");
        return -1;
    }

    kstring_t ks = KS_INITIALIZE;
    if (sam_hrecs_rebuild_text(hrecs, &ks) != 0) {
        ks_free(&ks);
        hts_log_error("Header text rebuild has failed");
        return -1;
    }

    hrecs->dirty = 0;

    /* Sync */
    free(bh->text);
    bh->l_text = ks_len(&ks);
    bh->text = ks_release(&ks);

    return 0;
}

/*
 * Appends a formatted line to an existing SAM header.
 * Line is a full SAM header record, eg "@SQ\tSN:foo\tLN:100", with
 * optional new-line. If it contains more than 1 line then multiple lines
 * will be added in order.
 *
 * Input text is of maximum length len or as terminated earlier by a NUL.
 * len may be 0 if unknown, in which case lines must be NUL-terminated.
 *
 * Returns 0 on success
 *        -1 on failure
 */
int sam_hdr_add_lines(sam_hdr_t *bh, const char *lines, size_t len) {
    sam_hrecs_t *hrecs;

    if (!bh || !lines)
        return -1;

    if (len == 0 && *lines == '\0')
        return 0;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -1;
        hrecs = bh->hrecs;
    }

    if (sam_hrecs_parse_lines(hrecs, lines, len) != 0)
        return -1;

    if (hrecs->refs_changed >= 0 && rebuild_target_arrays(bh) != 0)
        return -1;

    hrecs->dirty = 1;
    redact_header_text(bh);

    return 0;
}

/*
 * Adds a single line to a SAM header.
 * Specify type and one or more key,value pairs, ending with the NULL key.
 * Eg. sam_hdr_add_line(h, "SQ", "ID", "foo", "LN", "100", NULL).
 *
 * Returns 0 on success
 *        -1 on failure
 */
int sam_hdr_add_line(sam_hdr_t *bh, const char *type, ...) {
    va_list args;
    sam_hrecs_t *hrecs;

    if (!bh || !type)
        return -1;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -1;
        hrecs = bh->hrecs;
    }

    va_start(args, type);
    int ret = sam_hrecs_vadd(hrecs, type, args, NULL);
    va_end(args);

    if (ret == 0) {
        if (hrecs->refs_changed >= 0 && rebuild_target_arrays(bh) != 0)
            return -1;

        if (hrecs->dirty)
            redact_header_text(bh);
    }

    return ret;
}

/*
 * Returns a complete line of formatted text for a specific head type/ID
 * combination. If ID_key is NULL then it returns the first line of the specified
 * type.
 */
int sam_hdr_find_line_id(sam_hdr_t *bh, const char *type,
                      const char *ID_key, const char *ID_val, kstring_t *ks) {
    sam_hrecs_t *hrecs;
    if (!bh || !type)
        return -2;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -2;
        hrecs = bh->hrecs;
    }

    sam_hrec_type_t *ty = sam_hrecs_find_type_id(hrecs, type, ID_key, ID_val);
    if (!ty)
        return -1;

    ks->l = 0;
    if (build_header_line(ty, ks) < 0) {
        return -2;
    }

    return 0;
}

int sam_hdr_find_line_pos(sam_hdr_t *bh, const char *type,
                          int pos, kstring_t *ks) {
    sam_hrecs_t *hrecs;
    if (!bh || !type)
        return -2;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -2;
        hrecs = bh->hrecs;
    }

    sam_hrec_type_t *ty = sam_hrecs_find_type_pos(hrecs, type, pos);
    if (!ty)
        return -1;

    ks->l = 0;
    if (build_header_line(ty, ks) < 0) {
        return -2;
    }

    return 0;
}

/*
 * Remove a line from the header by specifying a tag:value that uniquely
 * identifies a line, i.e. the @SQ line containing "SN:ref1".
 * @SQ line is uniquely identified by SN tag.
 * @RG line is uniquely identified by ID tag.
 * @PG line is uniquely identified by ID tag.
 *
 * Returns 0 on success and -1 on error
 */

int sam_hdr_remove_line_id(sam_hdr_t *bh, const char *type, const char *ID_key, const char *ID_value) {
    sam_hrecs_t *hrecs;
    if (!bh || !type)
        return -1;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -1;
        hrecs = bh->hrecs;
    }

    if (!strncmp(type, "PG", 2)) {
        hts_log_warning("Removing PG lines is not supported!");
        return -1;
    }

    sam_hrec_type_t *type_found = sam_hrecs_find_type_id(hrecs, type, ID_key, ID_value);
    if (!type_found)
        return 0;

    int ret = sam_hrecs_remove_line(hrecs, type, type_found, 1);
    if (ret == 0) {
        if (hrecs->refs_changed >= 0 && rebuild_target_arrays(bh) != 0)
            return -1;

        if (hrecs->dirty)
            redact_header_text(bh);
    }

    return ret;
}

/*
 * Remove a line from the header by specifying the position in the type
 * group, i.e. 3rd @SQ line.
 *
 * Returns 0 on success and -1 on error
 */

int sam_hdr_remove_line_pos(sam_hdr_t *bh, const char *type, int position) {
    sam_hrecs_t *hrecs;
    if (!bh || !type || position < 0)
        return -1;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -1;
        hrecs = bh->hrecs;
    }

    if (!strncmp(type, "PG", 2)) {
        hts_log_warning("Removing PG lines is not supported!");
        return -1;
    }

    sam_hrec_type_t *type_found = sam_hrecs_find_type_pos(hrecs, type,
                                                          position);
    if (!type_found)
        return -1;

    int ret = sam_hrecs_remove_line(hrecs, type, type_found, 1);
    if (ret == 0) {
        if (hrecs->refs_changed >= 0 && rebuild_target_arrays(bh) != 0)
            return -1;

        if (hrecs->dirty)
            redact_header_text(bh);
    }

    return ret;
}

/*
 * Check if sam_hdr_update_line() is being used to change the name of
 * a record, and if the new name is going to clash with an existing one.
 *
 * If ap includes repeated keys, we go with the last one as sam_hrecs_vupdate()
 * will go through them all and leave the final one in place.
 *
 * Returns 0 if the name does not change
 *         1 if the name changes but does not clash
 *        -1 if the name changes and the new one is already in use
 */
static int check_for_name_update(sam_hrecs_t *hrecs, sam_hrec_type_t *rec,
                                 va_list ap, const char **old_name,
                                 const char **new_name,
                                 char id_tag_out[3],
                                 khash_t(m_s2i) **hash_out) {
    char *key, *val;
    const char *id_tag;
    sam_hrec_tag_t *tag, *prev;
    khash_t(m_s2i) *hash;
    khint_t k;
    int ret = 0;

    if        (rec->type == TYPEKEY("SQ")) {
        id_tag = "SN"; hash = hrecs->ref_hash;
    } else if (rec->type == TYPEKEY("RG")) {
        id_tag = "ID"; hash = hrecs->rg_hash;
    } else if (rec->type == TYPEKEY("PG")) {
        id_tag = "ID"; hash = hrecs->pg_hash;
    } else {
        return 0;
    }

    memcpy(id_tag_out, id_tag, 3);
    *hash_out = hash;

    tag = sam_hrecs_find_key(rec, id_tag, &prev);
    if (!tag)
        return 0;
    assert(tag->len >= 3);
    *old_name = tag->str + 3;

    while ((key = va_arg(ap, char *)) != NULL) {
        val = va_arg(ap, char *);
        if (!val) val = "";
        if (strcmp(key, id_tag) != 0) continue;
        if (strcmp(val, tag->str + 3) == 0) { ret = 0; continue; }
        k = kh_get(m_s2i, hash, val);
        ret = k < kh_end(hash) ? -1 : 1;
        *new_name = val;
    }
    return ret;
}

int sam_hdr_update_line(sam_hdr_t *bh, const char *type,
        const char *ID_key, const char *ID_value, ...) {
    sam_hrecs_t *hrecs;
    if (!bh)
        return -1;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -1;
        hrecs = bh->hrecs;
    }

    int ret, rename;
    sam_hrec_type_t *ty = sam_hrecs_find_type_id(hrecs, type, ID_key, ID_value);
    if (!ty)
        return -1;

    va_list args;
    const char *old_name = "?", *new_name = "?";
    char id_tag[3];
    khash_t(m_s2i) *hash = NULL;
    va_start(args, ID_value);
    rename = check_for_name_update(hrecs, ty, args,
                                   &old_name, &new_name, id_tag, &hash);
    va_end(args);
    if (rename < 0) {
        hts_log_error("Cannot rename @%s \"%s\" to \"%s\" : already exists",
                      type, old_name, new_name);
        return -1;
    }
    if (rename > 0 && TYPEKEY(type) == TYPEKEY("PG")) {
        // This is just too complicated
        hts_log_error("Renaming @PG records is not supported");
        return -1;
    }
    va_start(args, ID_value);
    ret = sam_hrecs_vupdate(hrecs, ty, args);
    va_end(args);

    if (ret)
        return ret;

    // TODO Account for @SQ-AN altnames

    if (rename) {
        // Adjust the hash table to point to the new name
        // sam_hrecs_update_hashes() should sort out everything else
        khint_t k = kh_get(m_s2i, hash, old_name);
        sam_hrec_tag_t *new_tag = sam_hrecs_find_key(ty, id_tag, NULL);
        int r, pos;
        assert(k < kh_end(hash));        // Or we wouldn't have found it earlier
        assert(new_tag && new_tag->str); // id_tag should exist
        assert(new_tag->len > 3);
        pos = kh_val(hash, k);
        kh_del(m_s2i, hash, k);
        k = kh_put(m_s2i, hash, new_tag->str + 3, &r);
        if (r < 1) {
            hts_log_error("Failed to rename item in hash table");
            return -1;
        }
        kh_val(hash, k) = pos;
    }

    ret = sam_hrecs_update_hashes(hrecs, TYPEKEY(type), ty);

    if (!ret && hrecs->refs_changed >= 0)
        ret = rebuild_target_arrays(bh);

    if (!ret && hrecs->dirty)
        redact_header_text(bh);

    return ret;
}

static int rebuild_hash(sam_hrecs_t *hrecs, const char *type) {
    sam_hrec_type_t *head, *step;
    khiter_t k;

    if (strncmp(type, "SQ", 2) == 0) {
        hrecs->nref = 0;
        kh_clear(m_s2i, hrecs->ref_hash);
    } else if (strncmp(type, "RG", 2) == 0) {
        hrecs->nrg = 0;
        kh_clear(m_s2i, hrecs->rg_hash);
    }

    k = kh_get(sam_hrecs_t, hrecs->h, TYPEKEY(type));

    if (k != kh_end(hrecs->h)) { // something to rebuild
        head = kh_val(hrecs->h, k);
        step = head;

        do {
            if (sam_hrecs_update_hashes(hrecs, TYPEKEY(type), step) == -1) {
                hts_log_error("Unable to rebuild hashes");
                return -1;
            }

            step = step->next;
        } while (step != head);
    }

    return 0;
}

int sam_hdr_remove_except(sam_hdr_t *bh, const char *type, const char *ID_key, const char *ID_value) {
    sam_hrecs_t *hrecs;
    if (!bh || !type)
        return -1;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -1;
        hrecs = bh->hrecs;
    }

    sam_hrec_type_t *step;
    int ret = 1, remove_all = (ID_key == NULL);

    if (!strncmp(type, "PG", 2) || !strncmp(type, "CO", 2)) {
        hts_log_warning("Removing PG or CO lines is not supported!");
        return -1;
    }

    sam_hrec_type_t *type_found = sam_hrecs_find_type_id(hrecs, type, ID_key, ID_value);
    if (!type_found) { // remove all line of this type
        khint_t k = kh_get(sam_hrecs_t, hrecs->h, TYPEKEY(type));
        if (k == kh_end(hrecs->h))
            return 0;
        type_found =  kh_val(hrecs->h, k);
        if (!type_found)
            return 0;
        remove_all = 1;
    }

    step = type_found->next;
    while (step != type_found) {
        sam_hrec_type_t *to_remove = step;
        step = step->next;
        ret &= sam_hrecs_remove_line(hrecs, type, to_remove, 0);
    }

    if (remove_all)
        ret &= sam_hrecs_remove_line(hrecs, type, type_found, 0);

    /* if RG or SQ, delete then rebuild the hashes (as it is faster
       to rebuild than delete one by one).
    */

    if ((strncmp(type, "SQ", 2) == 0) || (strncmp(type, "RG", 2) == 0)) {
        if (rebuild_hash(hrecs, type)) {
            return -1;
        }
    }

    if (!ret && hrecs->dirty)
        redact_header_text(bh);

    return 0;
}

int sam_hdr_remove_lines(sam_hdr_t *bh, const char *type, const char *id, void *vrh) {
    sam_hrecs_t *hrecs;
    rmhash_t *rh = (rmhash_t *)vrh;

    if (!bh || !type)
        return -1;
    if (!rh) // remove all lines
        return sam_hdr_remove_except(bh, type, NULL, NULL);
    if (!id)
        return -1;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -1;
        hrecs = bh->hrecs;
    }

    khint_t k = kh_get(sam_hrecs_t, hrecs->h, TYPEKEY(type));
    if (k == kh_end(hrecs->h)) // nothing to remove from
        return 0;

    sam_hrec_type_t *head = kh_val(hrecs->h, k);
    if (!head) {
        hts_log_error("Header inconsistency");
        return -1;
    }

    int ret = 0;
    sam_hrec_type_t *step = head->next;
    while (step != head) {
        sam_hrec_tag_t *tag = sam_hrecs_find_key(step, id, NULL);
        if (tag && tag->str && tag->len >= 3) {
           k = kh_get(rm, rh, tag->str+3);
           if (k == kh_end(rh)) { // value is not in the hash table, so remove
               sam_hrec_type_t *to_remove = step;
               step = step->next;
               ret |= sam_hrecs_remove_line(hrecs, type, to_remove, 0);
           } else {
               step = step->next;
           }
        } else { // tag is not on the line, so skip to next line
            step = step->next;
        }
    }

    // process the first line
    sam_hrec_tag_t * tag = sam_hrecs_find_key(head, id, NULL);
    if (tag && tag->str && tag->len >= 3) {
       k = kh_get(rm, rh, tag->str+3);
       if (k == kh_end(rh)) { // value is not in the hash table, so remove
           sam_hrec_type_t *to_remove = head;
           head = head->next;
           ret |= sam_hrecs_remove_line(hrecs, type, to_remove, 0);
       }
    }

    /* if RG or SQ, delete then rebuild the hashes (as it is faster
       to rebuild than delete one by one).
    */

    if ((strncmp(type, "SQ", 2) == 0) || (strncmp(type, "RG", 2) == 0)) {
        if (rebuild_hash(hrecs, type)) {
            return -1;
        }
    }

    if (!ret && hrecs->dirty)
        redact_header_text(bh);

    return ret;
}

int sam_hdr_count_lines(sam_hdr_t *bh, const char *type) {
    int count;
    sam_hrec_type_t *first_ty, *itr_ty;

    if (!bh || !type)
        return -1;

    if (!bh->hrecs) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -1;
    }

    // Deal with types that have counts
    switch (type[0]) {
    case 'S':
        if (type[1] == 'Q')
            return bh->hrecs->nref;
        break;
    case 'R':
        if (type[1] == 'G')
            return bh->hrecs->nrg;
        break;
    case 'P':
        if (type[1] == 'G')
            return bh->hrecs->npg;
        break;
    default:
        break;
    }

    first_ty = sam_hrecs_find_type_id(bh->hrecs, type, NULL, NULL);
    if (!first_ty)
        return 0;

    count = 1;
    for (itr_ty = first_ty->next;
         itr_ty && itr_ty != first_ty; itr_ty = itr_ty->next) {
        count++;
    }

    return count;
}

int sam_hdr_line_index(sam_hdr_t *bh,
                       const char *type,
                       const char *key) {
    sam_hrecs_t *hrecs;
    if (!bh || !type || !key)
        return -2;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -2;
        hrecs = bh->hrecs;
    }

    khint_t k;
    int idx = -1;
    switch (type[0]) {
    case 'S':
        if (type[1] == 'Q') {
            k = kh_get(m_s2i, hrecs->ref_hash, key);
            if (k != kh_end(hrecs->ref_hash))
                idx = kh_val(hrecs->ref_hash, k);
        } else {
            hts_log_warning("Type '%s' not supported. Only @SQ, @RG and @PG lines are indexed", type);
        }
        break;
    case 'R':
        if (type[1] == 'G') {
            k = kh_get(m_s2i, hrecs->rg_hash, key);
            if (k != kh_end(hrecs->rg_hash))
                idx = kh_val(hrecs->rg_hash, k);
        } else {
            hts_log_warning("Type '%s' not supported. Only @SQ, @RG and @PG lines are indexed", type);
        }
        break;
    case 'P':
        if (type[1] == 'G') {
            k = kh_get(m_s2i, hrecs->pg_hash, key);
            if (k != kh_end(hrecs->pg_hash))
                idx = kh_val(hrecs->pg_hash, k);
        } else {
            hts_log_warning("Type '%s' not supported. Only @SQ, @RG and @PG lines are indexed", type);
        }
        break;
    default:
        hts_log_warning("Type '%s' not supported. Only @SQ, @RG and @PG lines are indexed", type);
    }

    return idx;
}

const char *sam_hdr_line_name(sam_hdr_t *bh,
                              const char *type,
                              int pos) {
    sam_hrecs_t *hrecs;
    if (!bh || !type || pos < 0)
        return NULL;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return NULL;
        hrecs = bh->hrecs;
    }

    switch (type[0]) {
    case 'S':
        if (type[1] == 'Q') {
            if (pos < hrecs->nref)
                return hrecs->ref[pos].name;
        } else {
            hts_log_warning("Type '%s' not supported. Only @SQ, @RG and @PG lines are indexed", type);
        }
        break;
    case 'R':
        if (type[1] == 'G') {
            if (pos < hrecs->nrg)
                return hrecs->rg[pos].name;
        } else {
            hts_log_warning("Type '%s' not supported. Only @SQ, @RG and @PG lines are indexed", type);
        }
        break;
    case 'P':
        if (type[1] == 'G') {
            if (pos < hrecs->npg)
                return hrecs->pg[pos].name;
        } else {
            hts_log_warning("Type '%s' not supported. Only @SQ, @RG and @PG lines are indexed", type);
        }
        break;
    default:
        hts_log_warning("Type '%s' not supported. Only @SQ, @RG and @PG lines are indexed", type);
    }

    return NULL;
}

/* ==== Key:val level methods ==== */

int sam_hdr_find_tag_id(sam_hdr_t *bh,
                     const char *type,
                     const char *ID_key,
                     const char *ID_value,
                     const char *key,
                     kstring_t *ks) {
    sam_hrecs_t *hrecs;
    if (!bh || !type || !key)
        return -2;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -2;
        hrecs = bh->hrecs;
    }

    sam_hrec_type_t *ty = sam_hrecs_find_type_id(hrecs, type, ID_key, ID_value);
    if (!ty)
        return -1;

    sam_hrec_tag_t *tag = sam_hrecs_find_key(ty, key, NULL);
    if (!tag || !tag->str || tag->len < 4)
        return -1;

    ks->l = 0;
    if (kputsn(tag->str+3, tag->len-3, ks) == EOF) {
        return -2;
    }

    return 0;
}

int sam_hdr_find_tag_pos(sam_hdr_t *bh,
                     const char *type,
                     int pos,
                     const char *key,
                     kstring_t *ks) {
    sam_hrecs_t *hrecs;
    if (!bh || !type || !key)
        return -2;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -2;
        hrecs = bh->hrecs;
    }

    sam_hrec_type_t *ty = sam_hrecs_find_type_pos(hrecs, type, pos);
    if (!ty)
        return -1;

    sam_hrec_tag_t *tag = sam_hrecs_find_key(ty, key, NULL);
    if (!tag || !tag->str || tag->len < 4)
        return -1;

    ks->l = 0;
    if (kputsn(tag->str+3, tag->len-3, ks) == EOF) {
        return -2;
    }

    return 0;
}

int sam_hdr_remove_tag_id(sam_hdr_t *bh,
        const char *type,
        const char *ID_key,
        const char *ID_value,
        const char *key) {
    sam_hrecs_t *hrecs;
    if (!bh || !type || !key)
        return -1;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -1;
        hrecs = bh->hrecs;
    }

    sam_hrec_type_t *ty = sam_hrecs_find_type_id(hrecs, type, ID_key, ID_value);
    if (!ty)
        return -1;

    int ret = sam_hrecs_remove_key(hrecs, ty, key);
    if (!ret && hrecs->dirty)
        redact_header_text(bh);

    return ret;
}

/*
 * Reconstructs a kstring from the header hash table.
 * Returns 0 on success
 *        -1 on failure
 */
int sam_hrecs_rebuild_text(const sam_hrecs_t *hrecs, kstring_t *ks) {
    ks->l = 0;

    if (!hrecs->h || !hrecs->h->size) {
        return kputsn("", 0, ks) >= 0 ? 0 : -1;
    }
    if (sam_hrecs_rebuild_lines(hrecs, ks) != 0)
        return -1;

    return 0;
}

/*
 * Looks up a reference sequence by name and returns the numerical ID.
 * Returns -1 if unknown reference; -2 if header could not be parsed.
 */
int sam_hdr_name2tid(sam_hdr_t *bh, const char *ref) {
    sam_hrecs_t *hrecs;
    khint_t k;

    if (!bh)
        return -1;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -2;
        hrecs = bh->hrecs;
    }

    if (!hrecs->ref_hash)
        return -1;

    k = kh_get(m_s2i, hrecs->ref_hash, ref);
    return k == kh_end(hrecs->ref_hash) ? -1 : kh_val(hrecs->ref_hash, k);
}

const char *sam_hdr_tid2name(const sam_hdr_t *h, int tid) {
    sam_hrecs_t *hrecs;

    if (!h || tid < 0)
        return NULL;

    if ((hrecs = h->hrecs) != NULL && tid < hrecs->nref) {
        return hrecs->ref[tid].name;
    } else {
        if (tid < h->n_targets)
            return h->target_name[tid];
    }

    return NULL;
}

hts_pos_t sam_hdr_tid2len(const sam_hdr_t *h, int tid) {
    sam_hrecs_t *hrecs;

    if (!h || tid < 0)
        return 0;

    if ((hrecs = h->hrecs) != NULL && tid < hrecs->nref) {
        return hrecs->ref[tid].len;
    } else {
        if (tid < h->n_targets) {
            if (h->target_len[tid] < UINT32_MAX || !h->sdict) {
                return h->target_len[tid];
            } else {
                khash_t(s2i) *long_refs = (khash_t(s2i) *) h->sdict;
                khint_t k = kh_get(s2i, long_refs, h->target_name[tid]);
                if (k < kh_end(long_refs)) {
                    return kh_val(long_refs, k);
                } else {
                    return UINT32_MAX;
                }
            }
        }
    }

    return 0;
}

/*
 * Fixes any PP links in @PG headers.
 * If the entries are in order then this doesn't need doing, but in case
 * our header is out of order this goes through the hrecs->pg[] array
 * setting the prev_id field.
 *
 * Note we can have multiple complete chains. This code should identify the
 * tails of these chains as these are the entries we have to link to in
 * subsequent PP records.
 *
 * Returns 0 on success
 *        -1 on failure (indicating broken PG/PP records)
 */
static int sam_hdr_link_pg(sam_hdr_t *bh) {
    sam_hrecs_t *hrecs;
    int i, j, ret = 0, *new_pg_end;

    if (!bh)
        return -1;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -1;
        hrecs = bh->hrecs;
    }

    if (!hrecs->pgs_changed || !hrecs->npg)
        return 0;

    hrecs->npg_end_alloc = hrecs->npg;
    new_pg_end = hts_realloc_p(hrecs->pg_end, sizeof(*new_pg_end), hrecs->npg);
    if (!new_pg_end)
        return -1;
    hrecs->pg_end = new_pg_end;
    int *chain_size = calloc(hrecs->npg, sizeof(int));
    if (!chain_size)
        return -1;

    for (i = 0; i < hrecs->npg; i++)
        hrecs->pg_end[i] = i;

    for (i = 0; i < hrecs->npg; i++) {
        khint_t k;
        sam_hrec_tag_t *tag;

        assert(hrecs->pg[i].ty != NULL);
        for (tag = hrecs->pg[i].ty->tag; tag; tag = tag->next) {
            if (tag->str[0] == 'P' && tag->str[1] == 'P')
                break;
        }
        if (!tag) {
            /* Chain start points */
            continue;
        }

        k = kh_get(m_s2i, hrecs->pg_hash, tag->str+3);

        if (k == kh_end(hrecs->pg_hash)) {
            hts_log_warning("PG line with ID:%s has a PP link to missing program '%s'",
                    hrecs->pg[i].name, tag->str+3);
            continue;
        }

        int pp_idx = kh_val(hrecs->pg_hash, k);
        if (pp_idx == i) {
            hts_log_warning("PG line with ID:%s has a PP link to itself",
                            hrecs->pg[i].name);
            continue;
        }

        hrecs->pg[i].prev_id = hrecs->pg[pp_idx].id;
        hrecs->pg_end[pp_idx] = -1;
        chain_size[i] = chain_size[pp_idx]+1;
    }

    int last_end = -1;
    for (i = j = 0; i < hrecs->npg; i++) {
        if (hrecs->pg_end[i] != -1) {
            last_end = hrecs->pg_end[i];
            if (chain_size[i] > 0)
                hrecs->pg_end[j++] = hrecs->pg_end[i];
        }
    }
    /* Only leafs? Choose the last one! */
    if (!j && hrecs->npg_end > 0 && last_end >= 0) {
        hrecs->pg_end[0] = last_end;
        j = 1;
    }

    hrecs->npg_end = j;
    hrecs->pgs_changed = 0;

    /* mark as dirty or empty for rebuild */
    hrecs->dirty = 1;
    redact_header_text(bh);
    free(chain_size);

    return ret;
}

/*
 * Returns a unique ID from a base name.
 *
 * The value returned is valid until the next call to
 * this function.
 */
const char *sam_hdr_pg_id(sam_hdr_t *bh, const char *name) {
    sam_hrecs_t *hrecs;
    size_t name_len;
    const size_t name_extra = 17;
    if (!bh || !name)
        return NULL;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return NULL;
        hrecs = bh->hrecs;
    }

    khint_t k = kh_get(m_s2i, hrecs->pg_hash, name);
    if (k == kh_end(hrecs->pg_hash))
        return name;

    name_len = strlen(name);
    if (name_len > 1000) name_len = 1000;
    if (hrecs->ID_buf_sz < name_len + name_extra) {
        char *new_ID_buf = realloc(hrecs->ID_buf,
                                   hts_add_sat2(name_len, name_extra));
        if (new_ID_buf == NULL)
            return NULL;
        hrecs->ID_buf = new_ID_buf;
        hrecs->ID_buf_sz = name_len + name_extra;
    }

    do {
        snprintf(hrecs->ID_buf, hrecs->ID_buf_sz, "%.1000s.%d", name, hrecs->ID_cnt++);
        k = kh_get(m_s2i, hrecs->pg_hash, hrecs->ID_buf);
    } while (k != kh_end(hrecs->pg_hash));

    return hrecs->ID_buf;
}

/*
 * Add an @PG line.
 *
 * If we wish complete control over this use sam_hdr_add_line() directly. This
 * function uses that, but attempts to do a lot of tedious house work for
 * you too.
 *
 * - It will generate a suitable ID if the supplied one clashes.
 * - It will generate multiple @PG records if we have multiple PG chains.
 *
 * Call it as per sam_hdr_add_line() with a series of key,value pairs ending
 * in NULL.
 *
 * Returns 0 on success
 *        -1 on failure
 */
int sam_hdr_add_pg(sam_hdr_t *bh, const char *name, ...) {
    sam_hrecs_t *hrecs;
    const char *specified_id = NULL, *specified_pn = NULL, *specified_pp = NULL;
    const char *key, *val;
    if (!bh)
        return -1;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -1;
        hrecs = bh->hrecs;
    }

    bh->hrecs->pgs_changed = 1;
    if (sam_hdr_link_pg(bh) < 0) {
        hts_log_error("Error linking @PG lines");
        return -1;
    }

    va_list args;
    // Check for ID / PN / PP tags in varargs list
    va_start(args, name);
    while ((key = va_arg(args, const char *)) != NULL) {
        val = va_arg(args, const char *);
        if (!val) break;
        if (strcmp(key, "PN") == 0 && *val != '\0')
            specified_pn = val;
        else if (strcmp(key, "PP") == 0 && *val != '\0')
            specified_pp = val;
        else if (strcmp(key, "ID") == 0 && *val != '\0')
            specified_id = val;
    }
    va_end(args);

    if (specified_id && hrecs->pg_hash) {
        khint_t k = kh_get(m_s2i, hrecs->pg_hash, specified_id);
        if (k != kh_end(hrecs->pg_hash)) {
            hts_log_error("Header @PG ID:%s already present", specified_id);
            return -1;
        }
    }

    if (specified_pp && hrecs->pg_hash) {
        khint_t k = kh_get(m_s2i, hrecs->pg_hash, specified_pp);
        if (k == kh_end(hrecs->pg_hash)) {
            hts_log_error("Header @PG ID:%s referred to by PP tag not present",
                          specified_pp);
            return -1;
        }
    }

    if (!specified_pp && hrecs->npg_end) {
        /* Copy ends array to avoid us looping while modifying it */
        int *end = hts_malloc_p(sizeof(int), hrecs->npg_end);
        int i, nends = hrecs->npg_end;

        if (!end)
            return -1;

        memcpy(end, hrecs->pg_end, nends * sizeof(*end));

        for (i = 0; i < nends; i++) {
            const char *id = !specified_id ? sam_hdr_pg_id(bh, name) : "";
            if (!id) {
                free(end);
                return -1;
            }
            assert(end[i] >= 0 && end[i] < hrecs->npg);
            va_start(args, name);
            if (-1 == sam_hrecs_vadd(hrecs, "PG", args,
                                     "ID", id,
                                     "PN", !specified_pn ? name : "",
                                     "PP", hrecs->pg[end[i]].name,
                                     NULL)) {
                free(end);
                return  -1;
            }
            va_end(args);
        }

        free(end);
    } else {
        const char *id = !specified_id ? sam_hdr_pg_id(bh, name) : "";
        if (!id)
            return -1;
        va_start(args, name);
        if (-1 == sam_hrecs_vadd(hrecs, "PG", args,
                                 "ID", id,
                                 "PN", !specified_pn ? name : "",
                                 NULL))
            return -1;
        va_end(args);
    }

    hrecs->dirty = 1;
    redact_header_text(bh);

    return 0;
}

/*! Increments a reference count on bh.
 *
 * This permits multiple files to share the same header, all calling
 * sam_hdr_destroy when done, without causing errors for other open files.
 */
void sam_hdr_incr_ref(sam_hdr_t *bh) {
    if (!bh)
        return;
    bh->ref_count++;
}

/* ==== Internal methods ==== */

/*
 * Creates an empty SAM header.  Allocates space for the SAM header
 * structures (hash tables) ready to be populated.
 *
 * Returns a sam_hrecs_t struct on success (free with sam_hrecs_free())
 *         NULL on failure
 */
sam_hrecs_t *sam_hrecs_new(void) {
    sam_hrecs_t *hrecs = calloc(1, sizeof(*hrecs));

    if (!hrecs)
        return NULL;

    hrecs->h = kh_init(sam_hrecs_t);
    if (!hrecs->h)
        goto err;

    hrecs->ID_cnt = 1;

    hrecs->nref = 0;
    hrecs->ref_sz = 0;
    hrecs->ref  = NULL;
    if (!(hrecs->ref_hash = kh_init(m_s2i)))
        goto err;
    hrecs->refs_changed = -1;

    hrecs->nrg = 0;
    hrecs->rg_sz = 0;
    hrecs->rg  = NULL;
    if (!(hrecs->rg_hash = kh_init(m_s2i)))
        goto err;

    hrecs->npg = 0;
    hrecs->pg_sz = 0;
    hrecs->pg  = NULL;
    hrecs->npg_end = hrecs->npg_end_alloc = 0;
    hrecs->pg_end = NULL;
    if (!(hrecs->pg_hash = kh_init(m_s2i)))
        goto err;

    if (!(hrecs->tag_pool = pool_create(sizeof(sam_hrec_tag_t))))
        goto err;

    if (!(hrecs->type_pool = pool_create(sizeof(sam_hrec_type_t))))
        goto err;

    if (!(hrecs->str_pool = string_pool_create(65536)))
        goto err;

    if (sam_hrecs_init_type_order(hrecs, NULL))
        goto err;

    return hrecs;

err:
    if (hrecs->h)
        kh_destroy(sam_hrecs_t, hrecs->h);

    if (hrecs->tag_pool)
        pool_destroy(hrecs->tag_pool);

    if (hrecs->type_pool)
        pool_destroy(hrecs->type_pool);

    if (hrecs->str_pool)
        string_pool_destroy(hrecs->str_pool);

    free(hrecs);

    return NULL;
}
#if 0
/*
 * Produces a duplicate copy of source and returns it.
 * Returns NULL on failure
 */
sam_hrecs_t *sam_hrecs_dup(sam_hrecs_t *source) {
        return NULL;
}
#endif
/*! Deallocates all storage used by a sam_hrecs_t struct.
 *
 * This also decrements the header reference count. If after decrementing
 * it is still non-zero then the header is assumed to be in use by another
 * caller and the free is not done.
 *
 */
void sam_hrecs_free(sam_hrecs_t *hrecs) {
    if (!hrecs)
        return;

    if (hrecs->h)
        kh_destroy(sam_hrecs_t, hrecs->h);

    if (hrecs->ref_hash)
        kh_destroy(m_s2i, hrecs->ref_hash);

    if (hrecs->ref)
        free(hrecs->ref);

    if (hrecs->rg_hash)
        kh_destroy(m_s2i, hrecs->rg_hash);

    if (hrecs->rg)
        free(hrecs->rg);

    if (hrecs->pg_hash)
        kh_destroy(m_s2i, hrecs->pg_hash);

    if (hrecs->pg)
        free(hrecs->pg);

    if (hrecs->pg_end)
        free(hrecs->pg_end);

    if (hrecs->type_pool)
        pool_destroy(hrecs->type_pool);

    if (hrecs->tag_pool)
        pool_destroy(hrecs->tag_pool);

    if (hrecs->str_pool)
        string_pool_destroy(hrecs->str_pool);

    if (hrecs->type_order)
        free(hrecs->type_order);

    if (hrecs->ID_buf)
        free(hrecs->ID_buf);

    free(hrecs);
}

/*
 * Internal method already used by the CRAM code
 * Returns the first header item matching 'type'. If ID is non-NULL it checks
 * for the tag ID: and compares against the specified ID.
 *
 * Returns NULL if no type/ID is found
 */
sam_hrec_type_t *sam_hrecs_find_type_id(sam_hrecs_t *hrecs, const char *type,
                                     const char *ID_key, const char *ID_value) {
    if (!hrecs || !type)
        return NULL;
    sam_hrec_type_t *t1, *t2;
    khint_t k;

    /* Special case for types we have prebuilt hashes on */
    if (ID_key) {
        if (!ID_value)
            return NULL;

        if (type[0]   == 'S' && type[1]   == 'Q' &&
            ID_key[0] == 'S' && ID_key[1] == 'N') {
            k = kh_get(m_s2i, hrecs->ref_hash, ID_value);
            return k != kh_end(hrecs->ref_hash)
                ? hrecs->ref[kh_val(hrecs->ref_hash, k)].ty
                : NULL;
        }

        if (type[0]   == 'R' && type[1]   == 'G' &&
            ID_key[0] == 'I' && ID_key[1] == 'D') {
            k = kh_get(m_s2i, hrecs->rg_hash, ID_value);
            return k != kh_end(hrecs->rg_hash)
                ? hrecs->rg[kh_val(hrecs->rg_hash, k)].ty
                : NULL;
        }

        if (type[0]   == 'P' && type[1]   == 'G' &&
            ID_key[0] == 'I' && ID_key[1] == 'D') {
            k = kh_get(m_s2i, hrecs->pg_hash, ID_value);
            return k != kh_end(hrecs->pg_hash)
                ? hrecs->pg[kh_val(hrecs->pg_hash, k)].ty
                : NULL;
        }
    }

    k = kh_get(sam_hrecs_t, hrecs->h, TYPEKEY(type));
    if (k == kh_end(hrecs->h))
        return NULL;

    if (!ID_key)
        return kh_val(hrecs->h, k);

    t1 = t2 = kh_val(hrecs->h, k);
    do {
        sam_hrec_tag_t *tag;
        for (tag = t1->tag; tag; tag = tag->next) {
            if (tag->str[0] == ID_key[0] && tag->str[1] == ID_key[1]) {
                const char *cp1 = tag->str+3;
                const char *cp2 = ID_value;
                while (*cp1 && *cp1 == *cp2)
                    cp1++, cp2++;
                if (*cp2 || *cp1)
                    continue;
                return t1;
            }
        }
        t1 = t1->next;
    } while (t1 != t2);

    return NULL;
}

/*
 * Adds or updates tag key,value pairs in a header line.
 * Eg for adding M5 tags to @SQ lines or updating sort order for the
 * @HD line.
 *
 * va_list contains multiple key,value pairs ending in NULL.
 *
 * Returns 0 on success
 *        -1 on failure
 */
int sam_hrecs_vupdate(sam_hrecs_t *hrecs, sam_hrec_type_t *type, va_list ap) {
    if (!hrecs)
        return -1;

    for (;;) {
        char *k, *v, *str;
        sam_hrec_tag_t *tag, *prev = NULL;

        if (!(k = (char *)va_arg(ap, char *)))
            break;
        if (!(v = va_arg(ap, char *)))
            v = "";

        tag = sam_hrecs_find_key(type, k, &prev);
        if (!tag) {
            if (!(tag = pool_alloc(hrecs->tag_pool)))
                return -1;
            if (prev)
                prev->next = tag;
            else
                type->tag = tag;

            tag->next = NULL;
        }

        tag->len = 3 + strlen(v);
        str = string_alloc(hrecs->str_pool, hts_add_sat2(tag->len, 1));
        if (!str)
            return -1;

        if (snprintf(str, tag->len+1, "%2.2s:%s", k, v) < 0)
            return -1;

        tag->str = str;
    }

    hrecs->dirty = 1; //mark text as dirty and force a rebuild

    return 0;
}

/*
 * Adds or updates tag key,value pairs in a header line.
 * Eg for adding M5 tags to @SQ lines or updating sort order for the
 * @HD line.
 *
 * Specify multiple key,value pairs ending in NULL.
 *
 * Returns 0 on success
 *        -1 on failure
 */
static int sam_hrecs_update(sam_hrecs_t *hrecs, sam_hrec_type_t *type, ...) {
    va_list args;
    int res;
    va_start(args, type);
    res = sam_hrecs_vupdate(hrecs, type, args);
    va_end(args);
    return res;
}

/*
 * Looks for a specific key in a single sam header line identified by *type.
 * If prev is non-NULL it also fills this out with the previous tag, to
 * permit use in key removal. *prev is set to NULL when the tag is the first
 * key in the list. When a tag isn't found, prev (if non NULL) will be the last
 * tag in the existing list.
 *
 * Returns the tag pointer on success
 *         NULL on failure
 */
sam_hrec_tag_t *sam_hrecs_find_key(sam_hrec_type_t *type,
                                   const char *key,
                                   sam_hrec_tag_t **prev) {
    sam_hrec_tag_t *tag, *p = NULL;
    if (!type)
        return NULL;

    for (tag = type->tag; tag; p = tag, tag = tag->next) {
        if (tag->str[0] == key[0] && tag->str[1] == key[1]) {
            if (prev)
                *prev = p;
            return tag;
        }
    }

    if (prev)
        *prev = p;

    return NULL;
}

int sam_hrecs_remove_key(sam_hrecs_t *hrecs,
                         sam_hrec_type_t *type,
                         const char *key) {
    sam_hrec_tag_t *tag, *prev;
    if (!hrecs)
        return -1;
    tag = sam_hrecs_find_key(type, key, &prev);
    if (!tag)
        return 0; // Not there anyway

    if (type->type == TYPEKEY("SQ") && tag->str[0] == 'A' && tag->str[1] == 'N') {
        assert(tag->len >= 3);
        sam_hrec_tag_t *sn_tag = sam_hrecs_find_key(type, "SN", NULL);
        if (sn_tag) {
            assert(sn_tag->len >= 3);
            khint_t k = kh_get(m_s2i, hrecs->ref_hash, sn_tag->str + 3);
            if (k != kh_end(hrecs->ref_hash))
                sam_hrecs_remove_ref_altnames(hrecs, kh_val(hrecs->ref_hash, k), tag->str + 3);
        }
    }

    if (!prev) { //first tag
        type->tag = tag->next;
    } else {
        prev->next = tag->next;
    }
    pool_free(hrecs->tag_pool, tag);
    hrecs->dirty = 1; //mark text as dirty and force a rebuild

    return 1;
}

/*
 * Looks up a read-group by name and returns a pointer to the start of the
 * associated tag list.
 *
 * Returns NULL on failure
 */
sam_hrec_rg_t *sam_hrecs_find_rg(sam_hrecs_t *hrecs, const char *rg) {
    khint_t k = kh_get(m_s2i, hrecs->rg_hash, rg);
    return k == kh_end(hrecs->rg_hash)
        ? NULL
        : &hrecs->rg[kh_val(hrecs->rg_hash, k)];
}

#if DEBUG_HEADER
void sam_hrecs_dump(sam_hrecs_t *hrecs) {
    khint_t k;
    int i;

    printf("===DUMP===\n");
    for (k = kh_begin(hrecs->h); k != kh_end(hrecs->h); k++) {
        sam_hrec_type_t *t1, *t2;
        char c[2];
        int idx = 0;

        if (!kh_exist(hrecs->h, k))
            continue;

        t1 = t2 = kh_val(hrecs->h, k);
        c[0] = kh_key(hrecs->h, k)>>8;
        c[1] = kh_key(hrecs->h, k)&0xff;
        printf("Type %.2s\n", c);

        do {
            sam_hrec_tag_t *tag;
            printf(">>>%d ", idx++);
            for (tag = t1->tag; tag; tag=tag->next) {
                if (strncmp(c, "CO", 2))
                    printf("\"%.2s\":\"%.*s\"\t", tag->str, tag->len-3, tag->str+3);
                else
                    printf("%s", tag->str);
            }
            putchar('\n');
            t1 = t1->next;
        } while (t1 != t2);
    }

    /* Dump out PG chains */
    printf("\n@PG chains:\n");
    for (i = 0; i < hrecs->npg_end; i++) {
        int j;
        printf("  %d:", i);
        for (j = hrecs->pg_end[i]; j != -1; j = hrecs->pg[j].prev_id) {
            printf("%s%d(%.*s)",
                   j == hrecs->pg_end[i] ? " " : "->",
                   j, hrecs->pg[j].name_len, hrecs->pg[j].name);
        }
        printf("\n");
    }

    puts("===END DUMP===");
}
#endif

/*
 * Returns the sort order:
 */
enum sam_sort_order sam_hrecs_sort_order(sam_hrecs_t *hrecs) {
    khint_t k;
    enum sam_sort_order so;

    so = ORDER_UNKNOWN;
    k = kh_get(sam_hrecs_t, hrecs->h, TYPEKEY("HD"));
    if (k != kh_end(hrecs->h)) {
        sam_hrec_type_t *ty = kh_val(hrecs->h, k);
        sam_hrec_tag_t *tag;
        for (tag = ty->tag; tag; tag = tag->next) {
            if (tag->str[0] == 'S' && tag->str[1] == 'O') {
                if (strcmp(tag->str+3, "unsorted") == 0)
                    so = ORDER_UNSORTED;
                else if (strcmp(tag->str+3, "queryname") == 0)
                    so = ORDER_NAME;
                else if (strcmp(tag->str+3, "coordinate") == 0)
                    so = ORDER_COORD;
                else if (strcmp(tag->str+3, "unknown") != 0)
                    hts_log_error("Unknown sort order field: %s", tag->str+3);
            }
        }
    }

    return so;
}

enum sam_group_order sam_hrecs_group_order(sam_hrecs_t *hrecs) {
    khint_t k;
    enum sam_group_order go;

    go = ORDER_NONE;
    k = kh_get(sam_hrecs_t, hrecs->h, TYPEKEY("HD"));
    if (k != kh_end(hrecs->h)) {
        sam_hrec_type_t *ty = kh_val(hrecs->h, k);
        sam_hrec_tag_t *tag;
        for (tag = ty->tag; tag; tag = tag->next) {
            if (tag->str[0] == 'G' && tag->str[1] == 'O') {
                if (strcmp(tag->str+3, "query") == 0)
                    go = ORDER_QUERY;
                else if (strcmp(tag->str+3, "reference") == 0)
                    go = ORDER_REFERENCE;
            }
        }
    }

    return go;
}


/*
 *
 * Removes all the instances of a tag.
 * Returns number removed (can be 0) or -1 on a fail.
 *
 */
int sam_hdr_remove_tag_all(sam_hdr_t *bh, const char *type, const char *tag_id) {
    sam_hrecs_t *hrecs;
    sam_hrec_type_t *t1, *t2;
    khint_t k;
    int ret = 0;
    int deleted = 0;

    if (!bh || !type || !tag_id)
        return -1;

    if (!(hrecs = bh->hrecs)) {
        if (sam_hdr_fill_hrecs(bh) != 0)
            return -1;
        hrecs = bh->hrecs;
    }

    k = kh_get(sam_hrecs_t, hrecs->h, TYPEKEY(type));

    if (k == kh_end(hrecs->h))
        return 0;

    t1 = t2 = kh_val(hrecs->h, k);

    do {
        if ((ret = sam_hrecs_remove_key(hrecs, t1, tag_id)) == 1) {
            deleted++;
        } else if (ret == -1) {
            return -1;
        }

        t1 = t1->next;
    } while (t1 != t2);

    return deleted;
}