ezbamqc: ezBAMQC/src/htslib/vcf.c comparison

comparison ezBAMQC/src/htslib/vcf.c @ 0:dfa3745e5fd8

Uploaded

author	youngkim
date	Thu, 24 Mar 2016 17:12:52 -0400
parents
children

comparison

equal deleted inserted replaced

--1:000000000000
+:dfa3745e5fd8
+/*  vcf.c -- VCF/BCF API functions.
+Copyright (C) 2012, 2013 Broad Institute.
+Copyright (C) 2012-2014 Genome Research Ltd.
+Portions copyright (C) 2014 Intel Corporation.
+Author: Heng Li <lh3@sanger.ac.uk>
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+DEALINGS IN THE SOFTWARE.  */
+#include <zlib.h>
+#include <stdio.h>
+#include <ctype.h>
+#include <assert.h>
+#include <string.h>
+#include <stdlib.h>
+#include <limits.h>
+#include "htslib/kstring.h"
+#include "htslib/bgzf.h"
+#include "htslib/vcf.h"
+#include "htslib/tbx.h"
+#include "htslib/hfile.h"
+#include "htslib/khash_str2int.h"
+#include "htslib/khash.h"
+KHASH_MAP_INIT_STR(vdict, bcf_idinfo_t)
+typedef khash_t(vdict) vdict_t;
+#include "htslib/kseq.h"
+KSTREAM_DECLARE(gzFile, gzread)
+uint32_t bcf_float_missing    = 0x7F800001;
+uint32_t bcf_float_vector_end = 0x7F800002;
+uint8_t bcf_type_shift[] = { 0, 0, 1, 2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
+static bcf_idinfo_t bcf_idinfo_def = { .info = { 15, 15, 15 }, .hrec = { NULL, NULL, NULL}, .id = -1 };
+/*************************
+*** VCF header parser ***
+*************************/
+int bcf_hdr_sync(bcf_hdr_t *h);
+int bcf_hdr_add_sample(bcf_hdr_t *h, const char *s)
+{
+if ( !s ) return 0;
+const char *ss = s;
+while ( !*ss && isspace(*ss) ) ss++;
+if ( !*ss )
+{
+fprintf(stderr,"[E::%s] Empty sample name: trailing spaces/tabs in the header line?\n", __func__);
+abort();
+}
+vdict_t *d = (vdict_t*)h->dict[BCF_DT_SAMPLE];
+int ret;
+char *sdup = strdup(s);
+int k = kh_put(vdict, d, sdup, &ret);
+if (ret) { // absent
+kh_val(d, k) = bcf_idinfo_def;
+kh_val(d, k).id = kh_size(d) - 1;
+} else {
+if (hts_verbose >= 2)
+{
+fprintf(stderr, "[E::%s] Duplicated sample name '%s'\n", __func__, s);
+abort();
+}
+free(sdup);
+return -1;
+}
+int n = kh_size(d);
+h->samples = (char**) realloc(h->samples,sizeof(char*)*n);
+h->samples[n-1] = sdup;
+h->dirty = 1;
+return 0;
+}
+int bcf_hdr_parse_sample_line(bcf_hdr_t *h, const char *str)
+{
+int ret = 0;
+int i = 0;
+const char *p, *q;
+// add samples
+for (p = q = str;; ++q) {
+if (*q != '\t' && *q != 0 && *q != '\n') continue;
+if (++i > 9) {
+char *s = (char*)malloc(q - p + 1);
+strncpy(s, p, q - p);
+s[q - p] = 0;
+if ( bcf_hdr_add_sample(h,s) < 0 ) ret = -1;
+free(s);
+}
+if (*q == 0 || *q == '\n') break;
+p = q + 1;
+}
+bcf_hdr_add_sample(h,NULL);
+return ret;
+}
+int bcf_hdr_sync(bcf_hdr_t *h)
+{
+int i;
+for (i = 0; i < 3; i++)
+{
+vdict_t *d = (vdict_t*)h->dict[i];
+khint_t k;
+// find out the largest id, there may be holes because of IDX
+int max_id = -1;
+for (k=kh_begin(d); k<kh_end(d); k++)
+{
+if (!kh_exist(d,k)) continue;
+if ( max_id < kh_val(d,k).id ) max_id = kh_val(d,k).id;
+}
+if ( max_id >= h->n[i] )
+{
+h->id[i] = (bcf_idpair_t*)realloc(h->id[i], (max_id+1)*sizeof(bcf_idpair_t));
+for (k=h->n[i]; k<=max_id; k++)
+{
+h->id[i][k].key = NULL;
+h->id[i][k].val = NULL;
+}
+h->n[i] = max_id+1;
+}
+for (k=kh_begin(d); k<kh_end(d); k++)
+{
+if (!kh_exist(d,k)) continue;
+h->id[i][kh_val(d,k).id].key = kh_key(d,k);
+h->id[i][kh_val(d,k).id].val = &kh_val(d,k);
+}
+}
+h->dirty = 0;
+return 0;
+}
+void bcf_hrec_destroy(bcf_hrec_t *hrec)
+{
+free(hrec->key);
+if ( hrec->value ) free(hrec->value);
+int i;
+for (i=0; i<hrec->nkeys; i++)
+{
+free(hrec->keys[i]);
+free(hrec->vals[i]);
+}
+free(hrec->keys);
+free(hrec->vals);
+free(hrec);
+}
+// Copies all fields except IDX.
+bcf_hrec_t *bcf_hrec_dup(bcf_hrec_t *hrec)
+{
+bcf_hrec_t *out = (bcf_hrec_t*) calloc(1,sizeof(bcf_hrec_t));
+out->type = hrec->type;
+if ( hrec->key ) out->key = strdup(hrec->key);
+if ( hrec->value ) out->value = strdup(hrec->value);
+out->nkeys = hrec->nkeys;
+out->keys = (char**) malloc(sizeof(char*)*hrec->nkeys);
+out->vals = (char**) malloc(sizeof(char*)*hrec->nkeys);
+int i, j = 0;
+for (i=0; i<hrec->nkeys; i++)
+{
+if ( hrec->keys[i] && !strcmp("IDX",hrec->keys[i]) ) continue;
+if ( hrec->keys[i] ) out->keys[j] = strdup(hrec->keys[i]);
+if ( hrec->vals[i] ) out->vals[j] = strdup(hrec->vals[i]);
+j++;
+}
+if ( i!=j ) out->nkeys -= i-j;   // IDX was omitted
+return out;
+}
+void bcf_hrec_debug(FILE *fp, bcf_hrec_t *hrec)
+{
+fprintf(fp, "key=[%s] value=[%s]", hrec->key, hrec->value?hrec->value:"");
+int i;
+for (i=0; i<hrec->nkeys; i++)
+fprintf(fp, "\t[%s]=[%s]", hrec->keys[i],hrec->vals[i]);
+fprintf(fp, "\n");
+}
+void bcf_header_debug(bcf_hdr_t *hdr)
+{
+int i, j;
+for (i=0; i<hdr->nhrec; i++)
+{
+if ( !hdr->hrec[i]->value )
+{
+fprintf(stderr, "##%s=<", hdr->hrec[i]->key);
+fprintf(stderr,"%s=%s", hdr->hrec[i]->keys[0], hdr->hrec[i]->vals[0]);
+for (j=1; j<hdr->hrec[i]->nkeys; j++)
+fprintf(stderr,",%s=%s", hdr->hrec[i]->keys[j], hdr->hrec[i]->vals[j]);
+fprintf(stderr,">\n");
+}
+else
+fprintf(stderr,"##%s=%s\n", hdr->hrec[i]->key,hdr->hrec[i]->value);
+}
+}
+void bcf_hrec_add_key(bcf_hrec_t *hrec, const char *str, int len)
+{
+int n = ++hrec->nkeys;
+hrec->keys = (char**) realloc(hrec->keys, sizeof(char*)*n);
+hrec->vals = (char**) realloc(hrec->vals, sizeof(char*)*n);
+assert( len );
+hrec->keys[n-1] = (char*) malloc((len+1)*sizeof(char));
+memcpy(hrec->keys[n-1],str,len);
+hrec->keys[n-1][len] = 0;
+hrec->vals[n-1] = NULL;
+}
+void bcf_hrec_set_val(bcf_hrec_t *hrec, int i, const char *str, int len, int is_quoted)
+{
+if ( !str ) { hrec->vals[i] = NULL; return; }
+if ( hrec->vals[i] ) free(hrec->vals[i]);
+if ( is_quoted )
+{
+hrec->vals[i] = (char*) malloc((len+3)*sizeof(char));
+hrec->vals[i][0] = '"';
+memcpy(&hrec->vals[i][1],str,len);
+hrec->vals[i][len+1] = '"';
+hrec->vals[i][len+2] = 0;
+}
+else
+{
+hrec->vals[i] = (char*) malloc((len+1)*sizeof(char));
+memcpy(hrec->vals[i],str,len);
+hrec->vals[i][len] = 0;
+}
+}
+void hrec_add_idx(bcf_hrec_t *hrec, int idx)
+{
+int n = ++hrec->nkeys;
+hrec->keys = (char**) realloc(hrec->keys, sizeof(char*)*n);
+hrec->vals = (char**) realloc(hrec->vals, sizeof(char*)*n);
+hrec->keys[n-1] = strdup("IDX");
+kstring_t str = {0,0,0};
+kputw(idx, &str);
+hrec->vals[n-1] = str.s;
+}
+int bcf_hrec_find_key(bcf_hrec_t *hrec, const char *key)
+{
+int i;
+for (i=0; i<hrec->nkeys; i++)
+if ( !strcasecmp(key,hrec->keys[i]) ) return i;
+return -1;
+}
+static inline int is_escaped(const char *min, const char *str)
+{
+int n = 0;
+while ( --str>=min && *str=='\\' ) n++;
+return n%2;
+}
+bcf_hrec_t *bcf_hdr_parse_line(const bcf_hdr_t *h, const char *line, int *len)
+{
+const char *p = line;
+if (p[0] != '#' || p[1] != '#') { *len = 0; return NULL; }
+p += 2;
+const char *q = p;
+while ( *q && *q!='=' ) q++;
+int n = q-p;
+if ( *q!='=' || !n ) { *len = q-line+1; return NULL; } // wrong format
+bcf_hrec_t *hrec = (bcf_hrec_t*) calloc(1,sizeof(bcf_hrec_t));
+hrec->key = (char*) malloc(sizeof(char)*(n+1));
+memcpy(hrec->key,p,n);
+hrec->key[n] = 0;
+p = ++q;
+if ( *p!='<' ) // generic field, e.g. ##samtoolsVersion=0.1.18-r579
+{
+while ( *q && *q!='\n' ) q++;
+hrec->value = (char*) malloc((q-p+1)*sizeof(char));
+memcpy(hrec->value, p, q-p);
+hrec->value[q-p] = 0;
+*len = q-line+1;
+return hrec;
+}
+// structured line, e.g. ##INFO=<ID=PV1,Number=1,Type=Float,Description="P-value for baseQ bias">
+int nopen = 1;
+while ( *q && *q!='\n' && nopen )
+{
+p = ++q;
+while ( *q && isalnum(*q) ) q++;
+n = q-p;
+if ( *q!='=' || !n )
+{
+// wrong format
+while ( *q && *q!='\n' ) q++;
+kstring_t tmp = {0,0,0};
+kputsn(line,q-line,&tmp);
+fprintf(stderr,"Could not parse the header line: \"%s\"\n", tmp.s);
+free(tmp.s);
+*len = q-line+1;
+bcf_hrec_destroy(hrec);
+return NULL;
+}
+bcf_hrec_add_key(hrec, p, q-p);
+p = ++q;
+int quoted = *p=='"' ? 1 : 0;
+if ( quoted ) p++, q++;
+while (1)
+{
+if ( !*q ) break;
+if ( quoted ) { if ( *q=='"' && !is_escaped(p,q) ) break; }
+else
+{
+if ( *q=='<' ) nopen++;
+if ( *q=='>' ) nopen--;
+if ( !nopen ) break;
+if ( *q==',' && nopen==1 ) break;
+}
+q++;
+}
+bcf_hrec_set_val(hrec, hrec->nkeys-1, p, q-p, quoted);
+if ( quoted ) q++;
+if ( *q=='>' ) { nopen--; q++; }
+}
+*len = q-line+1;
+return hrec;
+}
+// returns: 1 when hdr needs to be synced, 0 otherwise
+int bcf_hdr_register_hrec(bcf_hdr_t *hdr, bcf_hrec_t *hrec)
+{
+// contig
+int i,j,k, ret;
+char *str;
+if ( !strcmp(hrec->key, "contig") )
+{
+hrec->type = BCF_HL_CTG;
+// Get the contig ID ($str) and length ($j)
+i = bcf_hrec_find_key(hrec,"length");
+if ( i<0 ) j = 0;
+else if ( sscanf(hrec->vals[i],"%d",&j)!=1 ) return 0;
+i = bcf_hrec_find_key(hrec,"ID");
+if ( i<0 ) return 0;
+str = strdup(hrec->vals[i]);
+// Register in the dictionary
+vdict_t *d = (vdict_t*)hdr->dict[BCF_DT_CTG];
+k = kh_put(vdict, d, str, &ret);
+if ( !ret ) { free(str); return 0; }    // already present
+int idx = bcf_hrec_find_key(hrec,"IDX");
+if ( idx!=-1 )
+{
+char *tmp = hrec->vals[idx];
+idx = strtol(hrec->vals[idx], &tmp, 10);
+if ( *tmp )
+{
+fprintf(stderr,"[%s:%d %s] Error parsing the IDX tag, skipping.\n", __FILE__,__LINE__,__FUNCTION__);
+return 0;
+}
+}
+else
+{
+idx = kh_size(d) - 1;
+hrec_add_idx(hrec, idx);
+}
+kh_val(d, k) = bcf_idinfo_def;
+kh_val(d, k).id = idx;
+kh_val(d, k).info[0] = j;
+kh_val(d, k).hrec[0] = hrec;
+return 1;
+}
+if ( !strcmp(hrec->key, "INFO") ) hrec->type = BCF_HL_INFO;
+else if ( !strcmp(hrec->key, "FILTER") ) hrec->type = BCF_HL_FLT;
+else if ( !strcmp(hrec->key, "FORMAT") ) hrec->type = BCF_HL_FMT;
+else if ( hrec->nkeys>0 ) { hrec->type = BCF_HL_STR; return 1; }
+else return 0;
+// INFO/FILTER/FORMAT
+char *id = NULL;
+int type = -1, num = -1, var = -1, idx = -1;
+for (i=0; i<hrec->nkeys; i++)
+{
+if ( !strcmp(hrec->keys[i], "ID") ) id = hrec->vals[i];
+else if ( !strcmp(hrec->keys[i], "IDX") )
+{
+char *tmp = hrec->vals[i];
+idx = strtol(hrec->vals[i], &tmp, 10);
+if ( *tmp )
+{
+fprintf(stderr,"[%s:%d %s] Error parsing the IDX tag, skipping.\n", __FILE__,__LINE__,__FUNCTION__);
+return 0;
+}
+}
+else if ( !strcmp(hrec->keys[i], "Type") )
+{
+if ( !strcmp(hrec->vals[i], "Integer") ) type = BCF_HT_INT;
+else if ( !strcmp(hrec->vals[i], "Float") ) type = BCF_HT_REAL;
+else if ( !strcmp(hrec->vals[i], "String") ) type = BCF_HT_STR;
+else if ( !strcmp(hrec->vals[i], "Character") ) type = BCF_HT_STR;
+else if ( !strcmp(hrec->vals[i], "Flag") ) type = BCF_HT_FLAG;
+else
+{
+fprintf(stderr, "[E::%s] The type \"%s\" not supported, assuming \"String\"\n", __func__, hrec->vals[i]);
+type = BCF_HT_STR;
+}
+}
+else if ( !strcmp(hrec->keys[i], "Number") )
+{
+if ( !strcmp(hrec->vals[i],"A") ) var = BCF_VL_A;
+else if ( !strcmp(hrec->vals[i],"R") ) var = BCF_VL_R;
+else if ( !strcmp(hrec->vals[i],"G") ) var = BCF_VL_G;
+else if ( !strcmp(hrec->vals[i],".") ) var = BCF_VL_VAR;
+else
+{
+sscanf(hrec->vals[i],"%d",&num);
+var = BCF_VL_FIXED;
+}
+if (var != BCF_VL_FIXED) num = 0xfffff;
+}
+}
+uint32_t info = (uint32_t)num<<12 | var<<8 | type<<4 | hrec->type;
+if ( !id ) return 0;
+str = strdup(id);
+vdict_t *d = (vdict_t*)hdr->dict[BCF_DT_ID];
+k = kh_put(vdict, d, str, &ret);
+if ( !ret )
+{
+// already present
+free(str);
+if ( kh_val(d, k).hrec[info&0xf] ) return 0;
+kh_val(d, k).info[info&0xf] = info;
+kh_val(d, k).hrec[info&0xf] = hrec;
+if ( idx==-1 ) hrec_add_idx(hrec, kh_val(d, k).id);
+return 1;
+}
+kh_val(d, k) = bcf_idinfo_def;
+kh_val(d, k).info[info&0xf] = info;
+kh_val(d, k).hrec[info&0xf] = hrec;
+kh_val(d, k).id = idx==-1 ? kh_size(d) - 1 : idx;
+if ( idx==-1 ) hrec_add_idx(hrec, kh_val(d, k).id);
+return 1;
+}
+int bcf_hdr_add_hrec(bcf_hdr_t *hdr, bcf_hrec_t *hrec)
+{
+if ( !hrec ) return 0;
+hrec->type = BCF_HL_GEN;
+if ( !bcf_hdr_register_hrec(hdr,hrec) )
+{
+// If one of the hashed field, then it is already present
+if ( hrec->type != BCF_HL_GEN )
+{
+bcf_hrec_destroy(hrec);
+return 0;
+}
+// Is one of the generic fields and already present?
+int i;
+for (i=0; i<hdr->nhrec; i++)
+{
+if ( hdr->hrec[i]->type!=BCF_HL_GEN ) continue;
+if ( !strcmp(hdr->hrec[i]->key,hrec->key) && !strcmp(hrec->key,"fileformat") ) break;
+if ( !strcmp(hdr->hrec[i]->key,hrec->key) && !strcmp(hdr->hrec[i]->value,hrec->value) ) break;
+}
+if ( i<hdr->nhrec )
+{
+bcf_hrec_destroy(hrec);
+return 0;
+}
+}
+// New record, needs to be added
+int n = ++hdr->nhrec;
+hdr->hrec = (bcf_hrec_t**) realloc(hdr->hrec, n*sizeof(bcf_hrec_t*));
+hdr->hrec[n-1] = hrec;
+hdr->dirty = 1;
+return hrec->type==BCF_HL_GEN ? 0 : 1;
+}
+/*
+*  Note that while querying of FLT,INFO,FMT,CTG lines is fast (the keys are hashed),
+*  the STR,GEN lines are searched for linearly in a linked list of all header lines.
+*  This may become a problem for VCFs with huge headers, we might need to build a
+*  dictionary for these lines as well.
+*/
+bcf_hrec_t *bcf_hdr_get_hrec(const bcf_hdr_t *hdr, int type, const char *key, const char *value, const char *str_class)
+{
+int i;
+if ( type==BCF_HL_GEN )
+{
+for (i=0; i<hdr->nhrec; i++)
+{
+if ( hdr->hrec[i]->type!=type ) continue;
+if ( strcmp(hdr->hrec[i]->key,key) ) continue;
+if ( !value || !strcmp(hdr->hrec[i]->value,value) ) return hdr->hrec[i];
+}
+return NULL;
+}
+else if ( type==BCF_HL_STR )
+{
+for (i=0; i<hdr->nhrec; i++)
+{
+if ( hdr->hrec[i]->type!=type ) continue;
+if ( strcmp(hdr->hrec[i]->key,str_class) ) continue;
+int j = bcf_hrec_find_key(hdr->hrec[i],key);
+if ( j>=0 && !strcmp(hdr->hrec[i]->vals[j],value) ) return hdr->hrec[i];
+}
+return NULL;
+}
+vdict_t *d = type==BCF_HL_CTG ? (vdict_t*)hdr->dict[BCF_DT_CTG] : (vdict_t*)hdr->dict[BCF_DT_ID];
+khint_t k = kh_get(vdict, d, value);
+if ( k == kh_end(d) ) return NULL;
+return kh_val(d, k).hrec[type==BCF_HL_CTG?0:type];
+}
+void bcf_hdr_check_sanity(bcf_hdr_t *hdr)
+{
+static int PL_warned = 0, GL_warned = 0;
+if ( !PL_warned )
+{
+int id = bcf_hdr_id2int(hdr, BCF_DT_ID, "PL");
+if ( bcf_hdr_idinfo_exists(hdr,BCF_HL_FMT,id) && bcf_hdr_id2length(hdr,BCF_HL_FMT,id)!=BCF_VL_G )
+{
+fprintf(stderr,"[W::%s] PL should be declared as Number=G\n", __func__);
+PL_warned = 1;
+}
+}
+if ( !GL_warned )
+{
+int id = bcf_hdr_id2int(hdr, BCF_HL_FMT, "GL");
+if ( bcf_hdr_idinfo_exists(hdr,BCF_HL_FMT,id) && bcf_hdr_id2length(hdr,BCF_HL_FMT,id)!=BCF_VL_G )
+{
+fprintf(stderr,"[W::%s] GL should be declared as Number=G\n", __func__);
+PL_warned = 1;
+}
+}
+}
+int bcf_hdr_parse(bcf_hdr_t *hdr, char *htxt)
+{
+int len, needs_sync = 0;
+char *p = htxt;
+// Check sanity: "fileformat" string must come as first
+bcf_hrec_t *hrec = bcf_hdr_parse_line(hdr,p,&len);
+if ( !hrec || !hrec->key || strcasecmp(hrec->key,"fileformat") )
+fprintf(stderr, "[W::%s] The first line should be ##fileformat; is the VCF/BCF header broken?\n", __func__);
+needs_sync += bcf_hdr_add_hrec(hdr, hrec);
+// The filter PASS must appear first in the dictionary
+hrec = bcf_hdr_parse_line(hdr,"##FILTER=<ID=PASS,Description=\"All filters passed\">",&len);
+needs_sync += bcf_hdr_add_hrec(hdr, hrec);
+// Parse the whole header
+while ( (hrec=bcf_hdr_parse_line(hdr,p,&len)) )
+{
+needs_sync += bcf_hdr_add_hrec(hdr, hrec);
+p += len;
+}
+int ret = bcf_hdr_parse_sample_line(hdr,p);
+bcf_hdr_sync(hdr);
+bcf_hdr_check_sanity(hdr);
+return ret;
+}
+int bcf_hdr_append(bcf_hdr_t *hdr, const char *line)
+{
+int len;
+bcf_hrec_t *hrec = bcf_hdr_parse_line(hdr, (char*) line, &len);
+if ( !hrec ) return -1;
+bcf_hdr_add_hrec(hdr, hrec);
+return 0;
+}
+void bcf_hdr_remove(bcf_hdr_t *hdr, int type, const char *key)
+{
+int i;
+bcf_hrec_t *hrec;
+while (1)
+{
+if ( type==BCF_HL_FLT || type==BCF_HL_INFO || type==BCF_HL_FMT || type== BCF_HL_CTG )
+{
+hrec = bcf_hdr_get_hrec(hdr, type, "ID", key, NULL);
+if ( !hrec ) return;
+for (i=0; i<hdr->nhrec; i++)
+if ( hdr->hrec[i]==hrec ) break;
+assert( i<hdr->nhrec );
+vdict_t *d = type==BCF_HL_CTG ? (vdict_t*)hdr->dict[BCF_DT_CTG] : (vdict_t*)hdr->dict[BCF_DT_ID];
+khint_t k = kh_get(vdict, d, key);
+kh_val(d, k).hrec[type==BCF_HL_CTG?0:type] = NULL;
+}
+else
+{
+for (i=0; i<hdr->nhrec; i++)
+{
+if ( hdr->hrec[i]->type!=type ) continue;
+if ( type==BCF_HL_GEN )
+{
+if ( !strcmp(hdr->hrec[i]->key,key) ) break;
+}
+else
+{
+// not all structured lines have ID, we could be more sophisticated as in bcf_hdr_get_hrec()
+int j = bcf_hrec_find_key(hdr->hrec[i], "ID");
+if ( j>=0 && !strcmp(hdr->hrec[i]->vals[j],key) ) break;
+}
+}
+if ( i==hdr->nhrec ) return;
+hrec = hdr->hrec[i];
+}
+hdr->nhrec--;
+if ( i < hdr->nhrec )
+memmove(&hdr->hrec[i],&hdr->hrec[i+1],(hdr->nhrec-i)*sizeof(bcf_hrec_t*));
+bcf_hrec_destroy(hrec);
+hdr->dirty = 1;
+}
+}
+int bcf_hdr_printf(bcf_hdr_t *hdr, const char *fmt, ...)
+{
+va_list ap;
+va_start(ap, fmt);
+int n = vsnprintf(NULL, 0, fmt, ap) + 2;
+va_end(ap);
+char *line = (char*)malloc(n);
+va_start(ap, fmt);
+vsnprintf(line, n, fmt, ap);
+va_end(ap);
+int ret = bcf_hdr_append(hdr, line);
+free(line);
+return ret;
+}
+/**********************
+*** BCF header I/O ***
+**********************/
+const char *bcf_hdr_get_version(const bcf_hdr_t *hdr)
+{
+bcf_hrec_t *hrec = bcf_hdr_get_hrec(hdr, BCF_HL_GEN, "fileformat", NULL, NULL);
+if ( !hrec )
+{
+fprintf(stderr,"No version string found, assuming VCFv4.2\n");
+return "VCFv4.2";
+}
+return hrec->value;
+}
+void bcf_hdr_set_version(bcf_hdr_t *hdr, const char *version)
+{
+bcf_hrec_t *hrec = bcf_hdr_get_hrec(hdr, BCF_HL_GEN, "fileformat", NULL, NULL);
+if ( !hrec )
+{
+int len;
+kstring_t str = {0,0,0};
+ksprintf(&str,"##fileformat=%s", version);
+hrec = bcf_hdr_parse_line(hdr, str.s, &len);
+free(str.s);
+}
+else
+{
+free(hrec->value);
+hrec->value = strdup(version);
+}
+hdr->dirty = 1;
+}
+bcf_hdr_t *bcf_hdr_init(const char *mode)
+{
+int i;
+bcf_hdr_t *h;
+h = (bcf_hdr_t*)calloc(1, sizeof(bcf_hdr_t));
+for (i = 0; i < 3; ++i)
+h->dict[i] = kh_init(vdict);
+if ( strchr(mode,'w') )
+{
+bcf_hdr_append(h, "##fileformat=VCFv4.2");
+// The filter PASS must appear first in the dictionary
+bcf_hdr_append(h, "##FILTER=<ID=PASS,Description=\"All filters passed\">");
+}
+return h;
+}
+void bcf_hdr_destroy(bcf_hdr_t *h)
+{
+int i;
+khint_t k;
+for (i = 0; i < 3; ++i) {
+vdict_t *d = (vdict_t*)h->dict[i];
+if (d == 0) continue;
+for (k = kh_begin(d); k != kh_end(d); ++k)
+if (kh_exist(d, k)) free((char*)kh_key(d, k));
+kh_destroy(vdict, d);
+free(h->id[i]);
+}
+for (i=0; i<h->nhrec; i++)
+bcf_hrec_destroy(h->hrec[i]);
+if (h->nhrec) free(h->hrec);
+if (h->samples) free(h->samples);
+free(h->keep_samples);
+free(h->transl[0]); free(h->transl[1]);
+free(h->mem.s);
+free(h);
+}
+bcf_hdr_t *bcf_hdr_read(htsFile *hfp)
+{
+if (hfp->format.format == vcf)
+return vcf_hdr_read(hfp);
+BGZF *fp = hfp->fp.bgzf;
+uint8_t magic[5];
+bcf_hdr_t *h;
+h = bcf_hdr_init("r");
+if ( bgzf_read(fp, magic, 5)<0 )
+{
+fprintf(stderr,"[%s:%d %s] Failed to read the header (reading BCF in text mode?)\n", __FILE__,__LINE__,__FUNCTION__);
+return NULL;
+}
+if (strncmp((char*)magic, "BCF\2\2", 5) != 0)
+{
+if (!strncmp((char*)magic, "BCF", 3))
+fprintf(stderr,"[%s:%d %s] invalid BCF2 magic string: only BCFv2.2 is supported.\n", __FILE__,__LINE__,__FUNCTION__);
+else if (hts_verbose >= 2)
+fprintf(stderr, "[E::%s] invalid BCF2 magic string\n", __func__);
+bcf_hdr_destroy(h);
+return 0;
+}
+int hlen;
+char *htxt;
+bgzf_read(fp, &hlen, 4);
+htxt = (char*)malloc(hlen);
+bgzf_read(fp, htxt, hlen);
+bcf_hdr_parse(h, htxt);
+free(htxt);
+return h;
+}
+int bcf_hdr_write(htsFile *hfp, bcf_hdr_t *h)
+{
+if ( h->dirty ) bcf_hdr_sync(h);
+if (hfp->format.format == vcf || hfp->format.format == text_format)
+return vcf_hdr_write(hfp, h);
+int hlen;
+char *htxt = bcf_hdr_fmt_text(h, 1, &hlen);
+hlen++; // include the \0 byte
+BGZF *fp = hfp->fp.bgzf;
+if ( bgzf_write(fp, "BCF\2\2", 5) !=5 ) return -1;
+if ( bgzf_write(fp, &hlen, 4) !=4 ) return -1;
+if ( bgzf_write(fp, htxt, hlen) != hlen ) return -1;
+free(htxt);
+return 0;
+}
+/********************
+*** BCF site I/O ***
+********************/
+bcf1_t *bcf_init1()
+{
+bcf1_t *v;
+v = (bcf1_t*)calloc(1, sizeof(bcf1_t));
+return v;
+}
+void bcf_clear(bcf1_t *v)
+{
+int i;
+for (i=0; i<v->d.m_info; i++)
+{
+if ( v->d.info[i].vptr_free )
+{
+free(v->d.info[i].vptr - v->d.info[i].vptr_off);
+v->d.info[i].vptr_free = 0;
+}
+}
+for (i=0; i<v->d.m_fmt; i++)
+{
+if ( v->d.fmt[i].p_free )
+{
+free(v->d.fmt[i].p - v->d.fmt[i].p_off);
+v->d.fmt[i].p_free = 0;
+}
+}
+v->rid = v->pos = v->rlen = v->unpacked = 0;
+bcf_float_set_missing(v->qual);
+v->n_info = v->n_allele = v->n_fmt = v->n_sample = 0;
+v->shared.l = v->indiv.l = 0;
+v->d.var_type = -1;
+v->d.shared_dirty = 0;
+v->d.indiv_dirty  = 0;
+v->d.n_flt = 0;
+v->errcode = 0;
+if (v->d.m_als) v->d.als[0] = 0;
+if (v->d.m_id) v->d.id[0] = 0;
+}
+void bcf_empty1(bcf1_t *v)
+{
+bcf_clear1(v);
+free(v->d.id);
+free(v->d.als);
+free(v->d.allele); free(v->d.flt); free(v->d.info); free(v->d.fmt);
+if (v->d.var ) free(v->d.var);
+free(v->shared.s); free(v->indiv.s);
+}
+void bcf_destroy1(bcf1_t *v)
+{
+bcf_empty1(v);
+free(v);
+}
+static inline int bcf_read1_core(BGZF *fp, bcf1_t *v)
+{
+uint32_t x[8];
+int ret;
+if ((ret = bgzf_read(fp, x, 32)) != 32) {
+if (ret == 0) return -1;
+return -2;
+}
+bcf_clear1(v);
+x[0] -= 24; // to exclude six 32-bit integers
+ks_resize(&v->shared, x[0]);
+ks_resize(&v->indiv, x[1]);
+memcpy(v, x + 2, 16);
+v->n_allele = x[6]>>16; v->n_info = x[6]&0xffff;
+v->n_fmt = x[7]>>24; v->n_sample = x[7]&0xffffff;
+v->shared.l = x[0], v->indiv.l = x[1];
+// silent fix of broken BCFs produced by earlier versions of bcf_subset, prior to and including bd6ed8b4
+if ( (!v->indiv.l || !v->n_sample) && v->n_fmt ) v->n_fmt = 0;
+bgzf_read(fp, v->shared.s, v->shared.l);
+bgzf_read(fp, v->indiv.s, v->indiv.l);
+return 0;
+}
+#define bit_array_size(n) ((n)/8+1)
+#define bit_array_set(a,i)   ((a)[(i)/8] |=   1 << ((i)%8))
+#define bit_array_clear(a,i) ((a)[(i)/8] &= ~(1 << ((i)%8)))
+#define bit_array_test(a,i)  ((a)[(i)/8] &   (1 << ((i)%8)))
+static inline uint8_t *bcf_unpack_fmt_core1(uint8_t *ptr, int n_sample, bcf_fmt_t *fmt);
+int bcf_subset_format(const bcf_hdr_t *hdr, bcf1_t *rec)
+{
+if ( !hdr->keep_samples ) return 0;
+if ( !bcf_hdr_nsamples(hdr) )
+{
+rec->indiv.l = rec->n_sample = 0;
+return 0;
+}
+int i, j;
+uint8_t *ptr = (uint8_t*)rec->indiv.s, *dst = NULL, *src;
+bcf_dec_t *dec = &rec->d;
+hts_expand(bcf_fmt_t, rec->n_fmt, dec->m_fmt, dec->fmt);
+for (i=0; i<dec->m_fmt; ++i) dec->fmt[i].p_free = 0;
+for (i=0; i<rec->n_fmt; i++)
+{
+ptr = bcf_unpack_fmt_core1(ptr, rec->n_sample, &dec->fmt[i]);
+src = dec->fmt[i].p - dec->fmt[i].size;
+if ( dst )
+{
+memmove(dec->fmt[i-1].p + dec->fmt[i-1].p_len, dec->fmt[i].p - dec->fmt[i].p_off, dec->fmt[i].p_off);
+dec->fmt[i].p = dec->fmt[i-1].p + dec->fmt[i-1].p_len + dec->fmt[i].p_off;
+}
+dst = dec->fmt[i].p;
+for (j=0; j<hdr->nsamples_ori; j++)
+{
+src += dec->fmt[i].size;
+if ( !bit_array_test(hdr->keep_samples,j) ) continue;
+memmove(dst, src, dec->fmt[i].size);
+dst += dec->fmt[i].size;
+}
+rec->indiv.l -= dec->fmt[i].p_len - (dst - dec->fmt[i].p);
+dec->fmt[i].p_len = dst - dec->fmt[i].p;
+}
+rec->unpacked |= BCF_UN_FMT;
+rec->n_sample = bcf_hdr_nsamples(hdr);
+return 0;
+}
+int bcf_read(htsFile *fp, const bcf_hdr_t *h, bcf1_t *v)
+{
+if (fp->format.format == vcf) return vcf_read(fp,h,v);
+int ret = bcf_read1_core(fp->fp.bgzf, v);
+if ( ret!=0 || !h->keep_samples ) return ret;
+return bcf_subset_format(h,v);
+}
+int bcf_readrec(BGZF *fp, void *null, void *vv, int *tid, int *beg, int *end)
+{
+bcf1_t *v = (bcf1_t *) vv;
+int ret;
+if ((ret = bcf_read1_core(fp, v)) >= 0)
+*tid = v->rid, *beg = v->pos, *end = v->pos + v->rlen;
+return ret;
+}
+static inline void bcf1_sync_id(bcf1_t *line, kstring_t *str)
+{
+// single typed string
+if ( line->d.id && strcmp(line->d.id, ".") ) bcf_enc_vchar(str, strlen(line->d.id), line->d.id);
+else bcf_enc_size(str, 0, BCF_BT_CHAR);
+}
+static inline void bcf1_sync_alleles(bcf1_t *line, kstring_t *str)
+{
+// list of typed strings
+int i;
+for (i=0; i<line->n_allele; i++)
+bcf_enc_vchar(str, strlen(line->d.allele[i]), line->d.allele[i]);
+if ( !line->rlen && line->n_allele ) line->rlen = strlen(line->d.allele[0]);
+}
+static inline void bcf1_sync_filter(bcf1_t *line, kstring_t *str)
+{
+// typed vector of integers
+if ( line->d.n_flt ) bcf_enc_vint(str, line->d.n_flt, line->d.flt, -1);
+else bcf_enc_vint(str, 0, 0, -1);
+}
+static inline void bcf1_sync_info(bcf1_t *line, kstring_t *str)
+{
+// pairs of typed vectors
+int i, irm = -1;
+for (i=0; i<line->n_info; i++)
+{
+bcf_info_t *info = &line->d.info[i];
+if ( !info->vptr )
+{
+// marked for removal
+if ( irm < 0 ) irm = i;
+continue;
+}
+kputsn_(info->vptr - info->vptr_off, info->vptr_len + info->vptr_off, str);
+if ( irm >=0 )
+{
+bcf_info_t tmp = line->d.info[irm]; line->d.info[irm] = line->d.info[i]; line->d.info[i] = tmp;
+while ( irm<=i && line->d.info[irm].vptr ) irm++;
+}
+}
+if ( irm>=0 ) line->n_info = irm;
+}
+static int bcf1_sync(bcf1_t *line)
+{
+char *shared_ori = line->shared.s;
+size_t prev_len;
+kstring_t tmp = {0,0,0};
+if ( !line->shared.l )
+{
+// New line created via API, BCF data blocks do not exist. Get it ready for BCF output
+tmp = line->shared;
+bcf1_sync_id(line, &tmp);
+line->unpack_size[0] = tmp.l; prev_len = tmp.l;
+bcf1_sync_alleles(line, &tmp);
+line->unpack_size[1] = tmp.l - prev_len; prev_len = tmp.l;
+bcf1_sync_filter(line, &tmp);
+line->unpack_size[2] = tmp.l - prev_len;
+bcf1_sync_info(line, &tmp);
+line->shared = tmp;
+}
+else if ( line->d.shared_dirty )
+{
+// The line was edited, update the BCF data block, ptr_ori points
+// to the original unchanged BCF data.
+uint8_t *ptr_ori = (uint8_t *) line->shared.s;
+assert( line->unpacked & BCF_UN_STR );
+// ID: single typed string
+if ( line->d.shared_dirty & BCF1_DIRTY_ID )
+bcf1_sync_id(line, &tmp);
+else
+kputsn_(ptr_ori, line->unpack_size[0], &tmp);
+ptr_ori += line->unpack_size[0];
+line->unpack_size[0] = tmp.l; prev_len = tmp.l;
+// REF+ALT: list of typed strings
+if ( line->d.shared_dirty & BCF1_DIRTY_ALS )
+bcf1_sync_alleles(line, &tmp);
+else
+{
+kputsn_(ptr_ori, line->unpack_size[1], &tmp);
+if ( !line->rlen && line->n_allele ) line->rlen = strlen(line->d.allele[0]);
+}
+ptr_ori += line->unpack_size[1];
+line->unpack_size[1] = tmp.l - prev_len; prev_len = tmp.l;
+if ( line->unpacked & BCF_UN_FLT )
+{
+// FILTER: typed vector of integers
+if ( line->d.shared_dirty & BCF1_DIRTY_FLT )
+bcf1_sync_filter(line, &tmp);
+else if ( line->d.n_flt )
+kputsn_(ptr_ori, line->unpack_size[2], &tmp);
+else
+bcf_enc_vint(&tmp, 0, 0, -1);
+ptr_ori += line->unpack_size[2];
+line->unpack_size[2] = tmp.l - prev_len;
+if ( line->unpacked & BCF_UN_INFO )
+{
+// INFO: pairs of typed vectors
+if ( line->d.shared_dirty & BCF1_DIRTY_INF )
+{
+bcf1_sync_info(line, &tmp);
+ptr_ori = (uint8_t*)line->shared.s + line->shared.l;
+}
+}
+}
+int size = line->shared.l - (size_t)ptr_ori + (size_t)line->shared.s;
+if ( size ) kputsn_(ptr_ori, size, &tmp);
+free(line->shared.s);
+line->shared = tmp;
+}
+if ( line->shared.s != shared_ori && line->unpacked & BCF_UN_INFO )
+{
+// Reallocated line->shared.s block invalidated line->d.info[].vptr pointers
+size_t off_new = line->unpack_size[0] + line->unpack_size[1] + line->unpack_size[2];
+int i;
+for (i=0; i<line->n_info; i++)
+{
+uint8_t *vptr_free = line->d.info[i].vptr_free ? line->d.info[i].vptr - line->d.info[i].vptr_off : NULL;
+line->d.info[i].vptr = (uint8_t*) line->shared.s + off_new + line->d.info[i].vptr_off;
+off_new += line->d.info[i].vptr_len + line->d.info[i].vptr_off;
+if ( vptr_free )
+{
+free(vptr_free);
+line->d.info[i].vptr_free = 0;
+}
+}
+}
+if ( line->n_sample && line->n_fmt && (!line->indiv.l || line->d.indiv_dirty) )
+{
+// The genotype fields changed or are not present
+tmp.l = tmp.m = 0; tmp.s = NULL;
+int i, irm = -1;
+for (i=0; i<line->n_fmt; i++)
+{
+bcf_fmt_t *fmt = &line->d.fmt[i];
+if ( !fmt->p )
+{
+// marked for removal
+if ( irm < 0 ) irm = i;
+continue;
+}
+kputsn_(fmt->p - fmt->p_off, fmt->p_len + fmt->p_off, &tmp);
+if ( irm >=0 )
+{
+bcf_fmt_t tfmt = line->d.fmt[irm]; line->d.fmt[irm] = line->d.fmt[i]; line->d.fmt[i] = tfmt;
+while ( irm<=i && line->d.fmt[irm].p ) irm++;
+}
+}
+if ( irm>=0 ) line->n_fmt = irm;
+free(line->indiv.s);
+line->indiv = tmp;
+// Reallocated line->indiv.s block invalidated line->d.fmt[].p pointers
+size_t off_new = 0;
+for (i=0; i<line->n_fmt; i++)
+{
+uint8_t *p_free = line->d.fmt[i].p_free ? line->d.fmt[i].p - line->d.fmt[i].p_off : NULL;
+line->d.fmt[i].p = (uint8_t*) line->indiv.s + off_new + line->d.fmt[i].p_off;
+off_new += line->d.fmt[i].p_len + line->d.fmt[i].p_off;
+if ( p_free )
+{
+free(p_free);
+line->d.fmt[i].p_free = 0;
+}
+}
+}
+if ( !line->n_sample ) line->n_fmt = 0;
+line->d.shared_dirty = line->d.indiv_dirty = 0;
+return 0;
+}
+bcf1_t *bcf_copy(bcf1_t *dst, bcf1_t *src)
+{
+bcf1_sync(src);
+bcf_clear(dst);
+dst->rid  = src->rid;
+dst->pos  = src->pos;
+dst->rlen = src->rlen;
+dst->qual = src->qual;
+dst->n_info = src->n_info; dst->n_allele = src->n_allele;
+dst->n_fmt = src->n_fmt; dst->n_sample = src->n_sample;
+dst->shared.m = dst->shared.l = src->shared.l;
+dst->shared.s = (char*) malloc(dst->shared.l);
+memcpy(dst->shared.s,src->shared.s,dst->shared.l);
+dst->indiv.m = dst->indiv.l = src->indiv.l;
+dst->indiv.s = (char*) malloc(dst->indiv.l);
+memcpy(dst->indiv.s,src->indiv.s,dst->indiv.l);
+return dst;
+}
+bcf1_t *bcf_dup(bcf1_t *src)
+{
+bcf1_t *out = bcf_init1();
+return bcf_copy(out, src);
+}
+int bcf_write(htsFile *hfp, const bcf_hdr_t *h, bcf1_t *v)
+{
+if ( h->dirty )
+{
+// we could as well call bcf_hdr_sync here, not sure
+fprintf(stderr,"FIXME: dirty header not synced\n");
+exit(1);
+}
+if ( bcf_hdr_nsamples(h)!=v->n_sample )
+{
+fprintf(stderr,"[%s:%d %s] Broken VCF record, the number of columns at %s:%d does not match the number of samples (%d vs %d).\n",
+__FILE__,__LINE__,__FUNCTION__,bcf_seqname(h,v),v->pos+1, v->n_sample,bcf_hdr_nsamples(h));
+return -1;
+}
+if ( hfp->format.format == vcf || hfp->format.format == text_format )
+return vcf_write(hfp,h,v);
+if ( v->errcode )
+{
+// vcf_parse1() encountered a new contig or tag, undeclared in the
+// header.  At this point, the header must have been printed,
+// proceeding would lead to a broken BCF file. Errors must be checked
+// and cleared by the caller before we can proceed.
+fprintf(stderr,"[%s:%d %s] Unchecked error (%d), exiting.\n", __FILE__,__LINE__,__FUNCTION__,v->errcode);
+exit(1);
+}
+bcf1_sync(v);   // check if the BCF record was modified
+BGZF *fp = hfp->fp.bgzf;
+uint32_t x[8];
+x[0] = v->shared.l + 24; // to include six 32-bit integers
+x[1] = v->indiv.l;
+memcpy(x + 2, v, 16);
+x[6] = (uint32_t)v->n_allele<<16 | v->n_info;
+x[7] = (uint32_t)v->n_fmt<<24 | v->n_sample;
+if ( bgzf_write(fp, x, 32) != 32 ) return -1;
+if ( bgzf_write(fp, v->shared.s, v->shared.l) != v->shared.l ) return -1;
+if ( bgzf_write(fp, v->indiv.s, v->indiv.l) != v->indiv.l ) return -1;
+return 0;
+}
+/**********************
+*** VCF header I/O ***
+**********************/
+bcf_hdr_t *vcf_hdr_read(htsFile *fp)
+{
+kstring_t txt, *s = &fp->line;
+bcf_hdr_t *h;
+h = bcf_hdr_init("r");
+txt.l = txt.m = 0; txt.s = 0;
+while (hts_getline(fp, KS_SEP_LINE, s) >= 0) {
+if (s->l == 0) continue;
+if (s->s[0] != '#') {
+if (hts_verbose >= 2)
+fprintf(stderr, "[E::%s] no sample line\n", __func__);
+free(txt.s);
+bcf_hdr_destroy(h);
+return 0;
+}
+if (s->s[1] != '#' && fp->fn_aux) { // insert contigs here
+int dret;
+gzFile f;
+kstream_t *ks;
+kstring_t tmp;
+tmp.l = tmp.m = 0; tmp.s = 0;
+f = gzopen(fp->fn_aux, "r");
+ks = ks_init(f);
+while (ks_getuntil(ks, 0, &tmp, &dret) >= 0) {
+int c;
+kputs("##contig=<ID=", &txt); kputs(tmp.s, &txt);
+ks_getuntil(ks, 0, &tmp, &dret);
+kputs(",length=", &txt); kputw(atol(tmp.s), &txt);
+kputsn(">\n", 2, &txt);
+if (dret != '\n')
+while ((c = ks_getc(ks)) != '\n' && c != -1); // skip the rest of the line
+}
+free(tmp.s);
+ks_destroy(ks);
+gzclose(f);
+}
+kputsn(s->s, s->l, &txt);
+kputc('\n', &txt);
+if (s->s[1] != '#') break;
+}
+if ( !txt.s )
+{
+fprintf(stderr,"[%s:%d %s] Could not read the header\n", __FILE__,__LINE__,__FUNCTION__);
+return NULL;
+}
+bcf_hdr_parse(h, txt.s);
+// check tabix index, are all contigs listed in the header? add the missing ones
+tbx_t *idx = tbx_index_load(fp->fn);
+if ( idx )
+{
+int i, n, need_sync = 0;
+const char **names = tbx_seqnames(idx, &n);
+for (i=0; i<n; i++)
+{
+bcf_hrec_t *hrec = bcf_hdr_get_hrec(h, BCF_HL_CTG, "ID", (char*) names[i], NULL);
+if ( hrec ) continue;
+hrec = (bcf_hrec_t*) calloc(1,sizeof(bcf_hrec_t));
+hrec->key = strdup("contig");
+bcf_hrec_add_key(hrec, "ID", strlen("ID"));
+bcf_hrec_set_val(hrec, hrec->nkeys-1, (char*) names[i], strlen(names[i]), 0);
+bcf_hdr_add_hrec(h, hrec);
+need_sync = 1;
+}
+free(names);
+tbx_destroy(idx);
+if ( need_sync )
+bcf_hdr_sync(h);
+}
+free(txt.s);
+return h;
+}
+int bcf_hdr_set(bcf_hdr_t *hdr, const char *fname)
+{
+int i, n;
+char **lines = hts_readlines(fname, &n);
+if ( !lines ) return 1;
+for (i=0; i<n-1; i++)
+{
+int k;
+bcf_hrec_t *hrec = bcf_hdr_parse_line(hdr,lines[i],&k);
+if ( hrec ) bcf_hdr_add_hrec(hdr, hrec);
+free(lines[i]);
+}
+bcf_hdr_parse_sample_line(hdr,lines[n-1]);
+free(lines[n-1]);
+free(lines);
+bcf_hdr_sync(hdr);
+return 0;
+}
+static void _bcf_hrec_format(const bcf_hrec_t *hrec, int is_bcf, kstring_t *str)
+{
+if ( !hrec->value )
+{
+int j, nout = 0;
+ksprintf(str, "##%s=<", hrec->key);
+for (j=0; j<hrec->nkeys; j++)
+{
+// do not output IDX if output is VCF
+if ( !is_bcf && !strcmp("IDX",hrec->keys[j]) ) continue;
+if ( nout ) kputc(',',str);
+ksprintf(str,"%s=%s", hrec->keys[j], hrec->vals[j]);
+nout++;
+}
+ksprintf(str,">\n");
+}
+else
+ksprintf(str,"##%s=%s\n", hrec->key,hrec->value);
+}
+void bcf_hrec_format(const bcf_hrec_t *hrec, kstring_t *str)
+{
+_bcf_hrec_format(hrec,0,str);
+}
+char *bcf_hdr_fmt_text(const bcf_hdr_t *hdr, int is_bcf, int *len)
+{
+int i;
+kstring_t txt = {0,0,0};
+for (i=0; i<hdr->nhrec; i++)
+_bcf_hrec_format(hdr->hrec[i], is_bcf, &txt);
+ksprintf(&txt,"#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO");
+if ( bcf_hdr_nsamples(hdr) )
+{
+ksprintf(&txt,"\tFORMAT");
+for (i=0; i<bcf_hdr_nsamples(hdr); i++)
+ksprintf(&txt,"\t%s", hdr->samples[i]);
+}
+ksprintf(&txt,"\n");
+if ( len ) *len = txt.l;
+return txt.s;
+}
+const char **bcf_hdr_seqnames(const bcf_hdr_t *h, int *n)
+{
+vdict_t *d = (vdict_t*)h->dict[BCF_DT_CTG];
+int tid, m = kh_size(d);
+const char **names = (const char**) calloc(m,sizeof(const char*));
+khint_t k;
+for (k=kh_begin(d); k<kh_end(d); k++)
+{
+if ( !kh_exist(d,k) ) continue;
+tid = kh_val(d,k).id;
+assert( tid<m );
+names[tid] = kh_key(d,k);
+}
+// sanity check: there should be no gaps
+for (tid=0; tid<m; tid++)
+assert(names[tid]);
+*n = m;
+return names;
+}
+int vcf_hdr_write(htsFile *fp, const bcf_hdr_t *h)
+{
+int hlen;
+char *htxt = bcf_hdr_fmt_text(h, 0, &hlen);
+while (hlen && htxt[hlen-1] == 0) --hlen; // kill trailing zeros
+int ret;
+if ( fp->format.compression!=no_compression )
+ret = bgzf_write(fp->fp.bgzf, htxt, hlen);
+else
+ret = hwrite(fp->fp.hfile, htxt, hlen);
+free(htxt);
+return ret<0 ? -1 : 0;
+}
+/***********************
+*** Typed value I/O ***
+***********************/
+void bcf_enc_vint(kstring_t *s, int n, int32_t *a, int wsize)
+{
+int32_t max = INT32_MIN + 1, min = INT32_MAX;
+int i;
+if (n == 0) bcf_enc_size(s, 0, BCF_BT_NULL);
+else if (n == 1) bcf_enc_int1(s, a[0]);
+else {
+if (wsize <= 0) wsize = n;
+for (i = 0; i < n; ++i) {
+if (a[i] == bcf_int32_missing || a[i] == bcf_int32_vector_end ) continue;
+if (max < a[i]) max = a[i];
+if (min > a[i]) min = a[i];
+}
+if (max <= INT8_MAX && min > bcf_int8_vector_end) {
+bcf_enc_size(s, wsize, BCF_BT_INT8);
+for (i = 0; i < n; ++i)
+if ( a[i]==bcf_int32_vector_end ) kputc(bcf_int8_vector_end, s);
+else if ( a[i]==bcf_int32_missing ) kputc(bcf_int8_missing, s);
+else kputc(a[i], s);
+} else if (max <= INT16_MAX && min > bcf_int16_vector_end) {
+bcf_enc_size(s, wsize, BCF_BT_INT16);
+for (i = 0; i < n; ++i)
+{
+int16_t x;
+if ( a[i]==bcf_int32_vector_end ) x = bcf_int16_vector_end;
+else if ( a[i]==bcf_int32_missing ) x = bcf_int16_missing;
+else x = a[i];
+kputsn((char*)&x, 2, s);
+}
+} else {
+bcf_enc_size(s, wsize, BCF_BT_INT32);
+for (i = 0; i < n; ++i) {
+int32_t x = a[i];
+kputsn((char*)&x, 4, s);
+}
+}
+}
+}
+void bcf_enc_vfloat(kstring_t *s, int n, float *a)
+{
+bcf_enc_size(s, n, BCF_BT_FLOAT);
+kputsn((char*)a, n << 2, s);
+}
+void bcf_enc_vchar(kstring_t *s, int l, const char *a)
+{
+bcf_enc_size(s, l, BCF_BT_CHAR);
+kputsn(a, l, s);
+}
+void bcf_fmt_array(kstring_t *s, int n, int type, void *data)
+{
+int j = 0;
+if (n == 0) {
+kputc('.', s);
+return;
+}
+if (type == BCF_BT_CHAR)
+{
+char *p = (char*)data;
+for (j = 0; j < n && *p; ++j, ++p)
+{
+if ( *p==bcf_str_missing ) kputc('.', s);
+else kputc(*p, s);
+}
+}
+else
+{
+#define BRANCH(type_t, is_missing, is_vector_end, kprint) { \
+type_t *p = (type_t *) data; \
+for (j=0; j<n; j++) \
+{ \
+if ( is_vector_end ) break; \
+if ( j ) kputc(',', s); \
+if ( is_missing ) kputc('.', s); \
+else kprint; \
+} \
+}
+switch (type) {
+case BCF_BT_INT8:  BRANCH(int8_t,  p[j]==bcf_int8_missing,  p[j]==bcf_int8_vector_end,  kputw(p[j], s)); break;
+case BCF_BT_INT16: BRANCH(int16_t, p[j]==bcf_int16_missing, p[j]==bcf_int16_vector_end, kputw(p[j], s)); break;
+case BCF_BT_INT32: BRANCH(int32_t, p[j]==bcf_int32_missing, p[j]==bcf_int32_vector_end, kputw(p[j], s)); break;
+case BCF_BT_FLOAT: BRANCH(float,   bcf_float_is_missing(p[j]), bcf_float_is_vector_end(p[j]), ksprintf(s, "%g", p[j])); break;
+default: fprintf(stderr,"todo: type %d\n", type); exit(1); break;
+}
+#undef BRANCH
+}
+}
+uint8_t *bcf_fmt_sized_array(kstring_t *s, uint8_t *ptr)
+{
+int x, type;
+x = bcf_dec_size(ptr, &ptr, &type);
+bcf_fmt_array(s, x, type, ptr);
+return ptr + (x << bcf_type_shift[type]);
+}
+/********************
+*** VCF site I/O ***
+********************/
+typedef struct {
+int key, max_m, size, offset;
+uint32_t is_gt:1, max_g:15, max_l:16;
+uint32_t y;
+uint8_t *buf;
+} fmt_aux_t;
+static inline void align_mem(kstring_t *s)
+{
+if (s->l&7) {
+uint64_t zero = 0;
+int l = ((s->l + 7)>>3<<3) - s->l;
+kputsn((char*)&zero, l, s);
+}
+}
+// p,q is the start and the end of the FORMAT field
+int _vcf_parse_format(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v, char *p, char *q)
+{
+if ( !bcf_hdr_nsamples(h) ) return 0;
+char *r, *t;
+int j, l, m, g;
+khint_t k;
+ks_tokaux_t aux1;
+vdict_t *d = (vdict_t*)h->dict[BCF_DT_ID];
+kstring_t *mem = (kstring_t*)&h->mem;
+mem->l = 0;
+// count the number of format fields
+for (r = p, v->n_fmt = 1; *r; ++r)
+if (*r == ':') ++v->n_fmt;
+char *end = s->s + s->l;
+if ( q>=end )
+{
+fprintf(stderr,"[%s:%d %s] Error: FORMAT column with no sample columns starting at %s:%d\n", __FILE__,__LINE__,__FUNCTION__,s->s,v->pos+1);
+return -1;
+}
+fmt_aux_t *fmt = (fmt_aux_t*)alloca(v->n_fmt * sizeof(fmt_aux_t));
+// get format information from the dictionary
+for (j = 0, t = kstrtok(p, ":", &aux1); t; t = kstrtok(0, 0, &aux1), ++j) {
+*(char*)aux1.p = 0;
+k = kh_get(vdict, d, t);
+if (k == kh_end(d) || kh_val(d, k).info[BCF_HL_FMT] == 15) {
+fprintf(stderr, "[W::%s] FORMAT '%s' is not defined in the header, assuming Type=String\n", __func__, t);
+kstring_t tmp = {0,0,0};
+int l;
+ksprintf(&tmp, "##FORMAT=<ID=%s,Number=1,Type=String,Description=\"Dummy\">", t);
+bcf_hrec_t *hrec = bcf_hdr_parse_line(h,tmp.s,&l);
+free(tmp.s);
+if ( bcf_hdr_add_hrec((bcf_hdr_t*)h, hrec) ) bcf_hdr_sync((bcf_hdr_t*)h);
+k = kh_get(vdict, d, t);
+v->errcode = BCF_ERR_TAG_UNDEF;
+}
+fmt[j].max_l = fmt[j].max_m = fmt[j].max_g = 0;
+fmt[j].key = kh_val(d, k).id;
+fmt[j].is_gt = !strcmp(t, "GT");
+fmt[j].y = h->id[0][fmt[j].key].val->info[BCF_HL_FMT];
+}
+// compute max
+int n_sample_ori = -1;
+r = q + 1;  // r: position in the format string
+m = l = g = 1, v->n_sample = 0;  // m: max vector size, l: max field len, g: max number of alleles
+while ( r<end )
+{
+// can we skip some samples?
+if ( h->keep_samples )
+{
+n_sample_ori++;
+if ( !bit_array_test(h->keep_samples,n_sample_ori) )
+{
+while ( *r!='\t' && r<end ) r++;
+if ( *r=='\t' ) { *r = 0; r++; }
+continue;
+}
+}
+// collect fmt stats: max vector size, length, number of alleles
+j = 0;  // j-th format field
+for (;;)
+{
+if ( *r == '\t' ) *r = 0;
+if ( *r == ':' || !*r )  // end of field or end of sample
+{
+if (fmt[j].max_m < m) fmt[j].max_m = m;
+if (fmt[j].max_l < l - 1) fmt[j].max_l = l - 1;
+if (fmt[j].is_gt && fmt[j].max_g < g) fmt[j].max_g = g;
+l = 0, m = g = 1;
+if ( *r==':' )
+{
+j++;
+if ( j>=v->n_fmt )
+{
+fprintf(stderr,"Incorrect number of FORMAT fields at %s:%d\n", h->id[BCF_DT_CTG][v->rid].key,v->pos+1);
+exit(1);
+}
+}
+else break;
+}
+else if ( *r== ',' ) m++;
+else if ( fmt[j].is_gt && (*r == '|' || *r == '/') ) g++;
+if ( r>=end ) break;
+r++; l++;
+}
+v->n_sample++;
+if ( v->n_sample == bcf_hdr_nsamples(h) ) break;
+r++;
+}
+// allocate memory for arrays
+for (j = 0; j < v->n_fmt; ++j) {
+fmt_aux_t *f = &fmt[j];
+if ( !f->max_m ) f->max_m = 1;  // omitted trailing format field
+if ((f->y>>4&0xf) == BCF_HT_STR) {
+f->size = f->is_gt? f->max_g << 2 : f->max_l;
+} else if ((f->y>>4&0xf) == BCF_HT_REAL || (f->y>>4&0xf) == BCF_HT_INT) {
+f->size = f->max_m << 2;
+} else
+{
+fprintf(stderr, "[E::%s] the format type %d currently not supported\n", __func__, f->y>>4&0xf);
+abort(); // I do not know how to do with Flag in the genotype fields
+}
+align_mem(mem);
+f->offset = mem->l;
+ks_resize(mem, mem->l + v->n_sample * f->size);
+mem->l += v->n_sample * f->size;
+}
+for (j = 0; j < v->n_fmt; ++j)
+fmt[j].buf = (uint8_t*)mem->s + fmt[j].offset;
+// fill the sample fields; at beginning of the loop, t points to the first char of a format
+n_sample_ori = -1;
+t = q + 1; m = 0;   // m: sample id
+while ( t<end )
+{
+// can we skip some samples?
+if ( h->keep_samples )
+{
+n_sample_ori++;
+if ( !bit_array_test(h->keep_samples,n_sample_ori) )
+{
+while ( *t && t<end ) t++;
+t++;
+continue;
+}
+}
+if ( m == bcf_hdr_nsamples(h) ) break;
+j = 0; // j-th format field, m-th sample
+while ( *t )
+{
+fmt_aux_t *z = &fmt[j];
+if ((z->y>>4&0xf) == BCF_HT_STR) {
+if (z->is_gt) { // genotypes
+int32_t is_phased = 0, *x = (int32_t*)(z->buf + z->size * m);
+for (l = 0;; ++t) {
+if (*t == '.') ++t, x[l++] = is_phased;
+else x[l++] = (strtol(t, &t, 10) + 1) << 1 | is_phased;
+#if THOROUGH_SANITY_CHECKS
+assert( 0 );    // success of strtol,strtod not checked
+#endif
+is_phased = (*t == '|');
+if (*t == ':' || *t == 0) break;
+}
+if ( !l ) x[l++] = 0;   // An empty field, insert missing value
+for (; l < z->size>>2; ++l) x[l] = bcf_int32_vector_end;
+} else {
+char *x = (char*)z->buf + z->size * m;
+for (r = t, l = 0; *t != ':' && *t; ++t) x[l++] = *t;
+for (; l < z->size; ++l) x[l] = 0;
+}
+} else if ((z->y>>4&0xf) == BCF_HT_INT) {
+int32_t *x = (int32_t*)(z->buf + z->size * m);
+for (l = 0;; ++t) {
+if (*t == '.') x[l++] = bcf_int32_missing, ++t; // ++t to skip "."
+else x[l++] = strtol(t, &t, 10);
+if (*t == ':' || *t == 0) break;
+}
+if ( !l ) x[l++] = bcf_int32_missing;
+for (; l < z->size>>2; ++l) x[l] = bcf_int32_vector_end;
+} else if ((z->y>>4&0xf) == BCF_HT_REAL) {
+float *x = (float*)(z->buf + z->size * m);
+for (l = 0;; ++t) {
+if (*t == '.' && !isdigit(t[1])) bcf_float_set_missing(x[l++]), ++t; // ++t to skip "."
+else x[l++] = strtod(t, &t);
+if (*t == ':' || *t == 0) break;
+}
+if ( !l ) bcf_float_set_missing(x[l++]);    // An empty field, insert missing value
+for (; l < z->size>>2; ++l) bcf_float_set_vector_end(x[l]);
+} else abort();
+if (*t == 0) {
+for (++j; j < v->n_fmt; ++j) { // fill end-of-vector values
+z = &fmt[j];
+if ((z->y>>4&0xf) == BCF_HT_STR) {
+if (z->is_gt) {
+int32_t *x = (int32_t*)(z->buf + z->size * m);
+x[0] = bcf_int32_missing;
+for (l = 1; l < z->size>>2; ++l) x[l] = bcf_int32_vector_end;
+} else {
+char *x = (char*)z->buf + z->size * m;
+if ( z->size ) x[0] = '.';
+for (l = 1; l < z->size; ++l) x[l] = 0;
+}
+} else if ((z->y>>4&0xf) == BCF_HT_INT) {
+int32_t *x = (int32_t*)(z->buf + z->size * m);
+x[0] = bcf_int32_missing;
+for (l = 1; l < z->size>>2; ++l) x[l] = bcf_int32_vector_end;
+} else if ((z->y>>4&0xf) == BCF_HT_REAL) {
+float *x = (float*)(z->buf + z->size * m);
+bcf_float_set_missing(x[0]);
+for (l = 1; l < z->size>>2; ++l) bcf_float_set_vector_end(x[l]);
+}
+}
+break;
+}
+else
+{
+if (*t == ':') ++j;
+t++;
+}
+}
+m++; t++;
+}
+// write individual genotype information
+kstring_t *str = &v->indiv;
+int i;
+if (v->n_sample > 0) {
+for (i = 0; i < v->n_fmt; ++i) {
+fmt_aux_t *z = &fmt[i];
+bcf_enc_int1(str, z->key);
+if ((z->y>>4&0xf) == BCF_HT_STR && !z->is_gt) {
+bcf_enc_size(str, z->size, BCF_BT_CHAR);
+kputsn((char*)z->buf, z->size * v->n_sample, str);
+} else if ((z->y>>4&0xf) == BCF_HT_INT || z->is_gt) {
+bcf_enc_vint(str, (z->size>>2) * v->n_sample, (int32_t*)z->buf, z->size>>2);
+} else {
+bcf_enc_size(str, z->size>>2, BCF_BT_FLOAT);
+kputsn((char*)z->buf, z->size * v->n_sample, str);
+}
+}
+}
+if ( v->n_sample!=bcf_hdr_nsamples(h) )
+{
+fprintf(stderr,"[%s:%d %s] Number of columns at %s:%d does not match the number of samples (%d vs %d).\n",
+__FILE__,__LINE__,__FUNCTION__,bcf_seqname(h,v),v->pos+1, v->n_sample,bcf_hdr_nsamples(h));
+v->errcode |= BCF_ERR_NCOLS;
+return -1;
+}
+return 0;
+}
+int vcf_parse(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v)
+{
+int i = 0;
+char *p, *q, *r, *t;
+kstring_t *str;
+khint_t k;
+ks_tokaux_t aux;
+bcf_clear1(v);
+str = &v->shared;
+memset(&aux, 0, sizeof(ks_tokaux_t));
+for (p = kstrtok(s->s, "\t", &aux), i = 0; p; p = kstrtok(0, 0, &aux), ++i) {
+q = (char*)aux.p;
+*q = 0;
+if (i == 0) { // CHROM
+vdict_t *d = (vdict_t*)h->dict[BCF_DT_CTG];
+k = kh_get(vdict, d, p);
+if (k == kh_end(d))
+{
+// Simple error recovery for chromosomes not defined in the header. It will not help when VCF header has
+// been already printed, but will enable tools like vcfcheck to proceed.
+fprintf(stderr, "[W::%s] contig '%s' is not defined in the header. (Quick workaround: index the file with tabix.)\n", __func__, p);
+kstring_t tmp = {0,0,0};
+int l;
+ksprintf(&tmp, "##contig=<ID=%s>", p);
+bcf_hrec_t *hrec = bcf_hdr_parse_line(h,tmp.s,&l);
+free(tmp.s);
+if ( bcf_hdr_add_hrec((bcf_hdr_t*)h, hrec) ) bcf_hdr_sync((bcf_hdr_t*)h);
+k = kh_get(vdict, d, p);
+v->errcode = BCF_ERR_CTG_UNDEF;
+}
+v->rid = kh_val(d, k).id;
+} else if (i == 1) { // POS
+v->pos = atoi(p) - 1;
+} else if (i == 2) { // ID
+if (strcmp(p, ".")) bcf_enc_vchar(str, q - p, p);
+else bcf_enc_size(str, 0, BCF_BT_CHAR);
+} else if (i == 3) { // REF
+bcf_enc_vchar(str, q - p, p);
+v->n_allele = 1, v->rlen = q - p;
+} else if (i == 4) { // ALT
+if (strcmp(p, ".")) {
+for (r = t = p;; ++r) {
+if (*r == ',' || *r == 0) {
+bcf_enc_vchar(str, r - t, t);
+t = r + 1;
+++v->n_allele;
+}
+if (r == q) break;
+}
+}
+} else if (i == 5) { // QUAL
+if (strcmp(p, ".")) v->qual = atof(p);
+else memcpy(&v->qual, &bcf_float_missing, 4);
+if ( v->max_unpack && !(v->max_unpack>>1) ) return 0; // BCF_UN_STR
+} else if (i == 6) { // FILTER
+if (strcmp(p, ".")) {
+int32_t *a;
+int n_flt = 1, i;
+ks_tokaux_t aux1;
+vdict_t *d = (vdict_t*)h->dict[BCF_DT_ID];
+// count the number of filters
+if (*(q-1) == ';') *(q-1) = 0;
+for (r = p; *r; ++r)
+if (*r == ';') ++n_flt;
+a = (int32_t*)alloca(n_flt * sizeof(int32_t));
+// add filters
+for (t = kstrtok(p, ";", &aux1), i = 0; t; t = kstrtok(0, 0, &aux1)) {
+*(char*)aux1.p = 0;
+k = kh_get(vdict, d, t);
+if (k == kh_end(d))
+{
+// Simple error recovery for FILTERs not defined in the header. It will not help when VCF header has
+// been already printed, but will enable tools like vcfcheck to proceed.
+fprintf(stderr, "[W::%s] FILTER '%s' is not defined in the header\n", __func__, t);
+kstring_t tmp = {0,0,0};
+int l;
+ksprintf(&tmp, "##FILTER=<ID=%s,Description=\"Dummy\">", t);
+bcf_hrec_t *hrec = bcf_hdr_parse_line(h,tmp.s,&l);
+free(tmp.s);
+if ( bcf_hdr_add_hrec((bcf_hdr_t*)h, hrec) ) bcf_hdr_sync((bcf_hdr_t*)h);
+k = kh_get(vdict, d, t);
+v->errcode = BCF_ERR_TAG_UNDEF;
+}
+a[i++] = kh_val(d, k).id;
+}
+n_flt = i;
+bcf_enc_vint(str, n_flt, a, -1);
+} else bcf_enc_vint(str, 0, 0, -1);
+if ( v->max_unpack && !(v->max_unpack>>2) ) return 0;    // BCF_UN_FLT
+} else if (i == 7) { // INFO
+char *key;
+vdict_t *d = (vdict_t*)h->dict[BCF_DT_ID];
+v->n_info = 0;
+if (strcmp(p, ".")) {
+if (*(q-1) == ';') *(q-1) = 0;
+for (r = key = p;; ++r) {
+int c;
+char *val, *end;
+if (*r != ';' && *r != '=' && *r != 0) continue;
+val = end = 0;
+c = *r; *r = 0;
+if (c == '=') {
+val = r + 1;
+for (end = val; *end != ';' && *end != 0; ++end);
+c = *end; *end = 0;
+} else end = r;
+if ( !*key ) { if (c==0) break; r = end; key = r + 1; continue; }  // faulty VCF, ";;" in the INFO
+k = kh_get(vdict, d, key);
+if (k == kh_end(d) || kh_val(d, k).info[BCF_HL_INFO] == 15)
+{
+fprintf(stderr, "[W::%s] INFO '%s' is not defined in the header, assuming Type=String\n", __func__, key);
+kstring_t tmp = {0,0,0};
+int l;
+ksprintf(&tmp, "##INFO=<ID=%s,Number=1,Type=String,Description=\"Dummy\">", key);
+bcf_hrec_t *hrec = bcf_hdr_parse_line(h,tmp.s,&l);
+free(tmp.s);
+if ( bcf_hdr_add_hrec((bcf_hdr_t*)h, hrec) ) bcf_hdr_sync((bcf_hdr_t*)h);
+k = kh_get(vdict, d, key);
+v->errcode = BCF_ERR_TAG_UNDEF;
+}
+uint32_t y = kh_val(d, k).info[BCF_HL_INFO];
+++v->n_info;
+bcf_enc_int1(str, kh_val(d, k).id);
+if (val == 0) {
+bcf_enc_size(str, 0, BCF_BT_NULL);
+} else if ((y>>4&0xf) == BCF_HT_FLAG || (y>>4&0xf) == BCF_HT_STR) { // if Flag has a value, treat it as a string
+bcf_enc_vchar(str, end - val, val);
+} else { // int/float value/array
+int i, n_val;
+char *t, *te;
+for (t = val, n_val = 1; *t; ++t) // count the number of values
+if (*t == ',') ++n_val;
+if ((y>>4&0xf) == BCF_HT_INT) {
+int32_t *z;
+z = (int32_t*)alloca(n_val * sizeof(int32_t));
+for (i = 0, t = val; i < n_val; ++i, ++t)
+{
+z[i] = strtol(t, &te, 10);
+if ( te==t ) // conversion failed
+{
+z[i] = bcf_int32_missing;
+while ( *te && *te!=',' ) te++;
+}
+t = te;
+}
+bcf_enc_vint(str, n_val, z, -1);
+if (strcmp(key, "END") == 0) v->rlen = z[0] - v->pos;
+} else if ((y>>4&0xf) == BCF_HT_REAL) {
+float *z;
+z = (float*)alloca(n_val * sizeof(float));
+for (i = 0, t = val; i < n_val; ++i, ++t)
+{
+z[i] = strtod(t, &te);
+if ( te==t ) // conversion failed
+{
+bcf_float_set_missing(z[i]);
+while ( *te && *te!=',' ) te++;
+}
+t = te;
+}
+bcf_enc_vfloat(str, n_val, z);
+}
+}
+if (c == 0) break;
+r = end;
+key = r + 1;
+}
+}
+if ( v->max_unpack && !(v->max_unpack>>3) ) return 0;
+} else if (i == 8) // FORMAT
+return _vcf_parse_format(s, h, v, p, q);
+}
+return 0;
+}
+int vcf_read(htsFile *fp, const bcf_hdr_t *h, bcf1_t *v)
+{
+int ret;
+ret = hts_getline(fp, KS_SEP_LINE, &fp->line);
+if (ret < 0) return -1;
+return vcf_parse1(&fp->line, h, v);
+}
+static inline uint8_t *bcf_unpack_fmt_core1(uint8_t *ptr, int n_sample, bcf_fmt_t *fmt)
+{
+uint8_t *ptr_start = ptr;
+fmt->id = bcf_dec_typed_int1(ptr, &ptr);
+fmt->n = bcf_dec_size(ptr, &ptr, &fmt->type);
+fmt->size = fmt->n << bcf_type_shift[fmt->type];
+fmt->p = ptr;
+fmt->p_off  = ptr - ptr_start;
+fmt->p_free = 0;
+ptr += n_sample * fmt->size;
+fmt->p_len = ptr - fmt->p;
+return ptr;
+}
+static inline uint8_t *bcf_unpack_info_core1(uint8_t *ptr, bcf_info_t *info)
+{
+uint8_t *ptr_start = ptr;
+info->key = bcf_dec_typed_int1(ptr, &ptr);
+info->len = bcf_dec_size(ptr, &ptr, &info->type);
+info->vptr = ptr;
+info->vptr_off  = ptr - ptr_start;
+info->vptr_free = 0;
+info->v1.i = 0;
+if (info->len == 1) {
+if (info->type == BCF_BT_INT8 || info->type == BCF_BT_CHAR) info->v1.i = *(int8_t*)ptr;
+else if (info->type == BCF_BT_INT32) info->v1.i = *(int32_t*)ptr;
+else if (info->type == BCF_BT_FLOAT) info->v1.f = *(float*)ptr;
+else if (info->type == BCF_BT_INT16) info->v1.i = *(int16_t*)ptr;
+}
+ptr += info->len << bcf_type_shift[info->type];
+info->vptr_len = ptr - info->vptr;
+return ptr;
+}
+int bcf_unpack(bcf1_t *b, int which)
+{
+if ( !b->shared.l ) return 0; // Building a new BCF record from scratch
+uint8_t *ptr = (uint8_t*)b->shared.s, *ptr_ori;
+int *offset, i;
+bcf_dec_t *d = &b->d;
+if (which & BCF_UN_FLT) which |= BCF_UN_STR;
+if (which & BCF_UN_INFO) which |= BCF_UN_SHR;
+if ((which&BCF_UN_STR) && !(b->unpacked&BCF_UN_STR))
+{
+kstring_t tmp;
+// ID
+tmp.l = 0; tmp.s = d->id; tmp.m = d->m_id;
+ptr_ori = ptr;
+ptr = bcf_fmt_sized_array(&tmp, ptr);
+b->unpack_size[0] = ptr - ptr_ori;
+kputc('\0', &tmp);
+d->id = tmp.s; d->m_id = tmp.m;
+// REF and ALT are in a single block (d->als) and d->alleles are pointers into this block
+tmp.l = 0; tmp.s = d->als; tmp.m = d->m_als;
+offset = (int*)alloca(b->n_allele * sizeof(int));
+ptr_ori = ptr;
+for (i = 0; i < b->n_allele; ++i) {
+offset[i] = tmp.l;
+ptr = bcf_fmt_sized_array(&tmp, ptr);
+kputc('\0', &tmp);
+}
+b->unpack_size[1] = ptr - ptr_ori;
+d->als = tmp.s; d->m_als = tmp.m;
+hts_expand(char*, b->n_allele, d->m_allele, d->allele); // NM: hts_expand() is a macro
+for (i = 0; i < b->n_allele; ++i)
+d->allele[i] = d->als + offset[i];
+b->unpacked |= BCF_UN_STR;
+}
+if ((which&BCF_UN_FLT) && !(b->unpacked&BCF_UN_FLT)) { // FILTER
+ptr = (uint8_t*)b->shared.s + b->unpack_size[0] + b->unpack_size[1];
+ptr_ori = ptr;
+if (*ptr>>4) {
+int type;
+d->n_flt = bcf_dec_size(ptr, &ptr, &type);
+hts_expand(int, d->n_flt, d->m_flt, d->flt);
+for (i = 0; i < d->n_flt; ++i)
+d->flt[i] = bcf_dec_int1(ptr, type, &ptr);
+} else ++ptr, d->n_flt = 0;
+b->unpack_size[2] = ptr - ptr_ori;
+b->unpacked |= BCF_UN_FLT;
+}
+if ((which&BCF_UN_INFO) && !(b->unpacked&BCF_UN_INFO)) { // INFO
+ptr = (uint8_t*)b->shared.s + b->unpack_size[0] + b->unpack_size[1] + b->unpack_size[2];
+hts_expand(bcf_info_t, b->n_info, d->m_info, d->info);
+for (i = 0; i < d->m_info; ++i) d->info[i].vptr_free = 0;
+for (i = 0; i < b->n_info; ++i)
+ptr = bcf_unpack_info_core1(ptr, &d->info[i]);
+b->unpacked |= BCF_UN_INFO;
+}
+if ((which&BCF_UN_FMT) && b->n_sample && !(b->unpacked&BCF_UN_FMT)) { // FORMAT
+ptr = (uint8_t*)b->indiv.s;
+hts_expand(bcf_fmt_t, b->n_fmt, d->m_fmt, d->fmt);
+for (i = 0; i < d->m_fmt; ++i) d->fmt[i].p_free = 0;
+for (i = 0; i < b->n_fmt; ++i)
+ptr = bcf_unpack_fmt_core1(ptr, b->n_sample, &d->fmt[i]);
+b->unpacked |= BCF_UN_FMT;
+}
+return 0;
+}
+int vcf_format(const bcf_hdr_t *h, const bcf1_t *v, kstring_t *s)
+{
+int i;
+bcf_unpack((bcf1_t*)v, BCF_UN_ALL);
+kputs(h->id[BCF_DT_CTG][v->rid].key, s); // CHROM
+kputc('\t', s); kputw(v->pos + 1, s); // POS
+kputc('\t', s); kputs(v->d.id ? v->d.id : ".", s); // ID
+kputc('\t', s); // REF
+if (v->n_allele > 0) kputs(v->d.allele[0], s);
+else kputc('.', s);
+kputc('\t', s); // ALT
+if (v->n_allele > 1) {
+for (i = 1; i < v->n_allele; ++i) {
+if (i > 1) kputc(',', s);
+kputs(v->d.allele[i], s);
+}
+} else kputc('.', s);
+kputc('\t', s); // QUAL
+if ( bcf_float_is_missing(v->qual) ) kputc('.', s); // QUAL
+else ksprintf(s, "%g", v->qual);
+kputc('\t', s); // FILTER
+if (v->d.n_flt) {
+for (i = 0; i < v->d.n_flt; ++i) {
+if (i) kputc(';', s);
+kputs(h->id[BCF_DT_ID][v->d.flt[i]].key, s);
+}
+} else kputc('.', s);
+kputc('\t', s); // INFO
+if (v->n_info) {
+int first = 1;
+for (i = 0; i < v->n_info; ++i) {
+bcf_info_t *z = &v->d.info[i];
+if ( !z->vptr ) continue;
+if ( !first ) kputc(';', s); first = 0;
+kputs(h->id[BCF_DT_ID][z->key].key, s);
+if (z->len <= 0) continue;
+kputc('=', s);
+if (z->len == 1)
+{
+switch (z->type)
+{
+case BCF_BT_INT8:  if ( z->v1.i==bcf_int8_missing ) kputc('.', s); else kputw(z->v1.i, s); break;
+case BCF_BT_INT16: if ( z->v1.i==bcf_int16_missing ) kputc('.', s); else kputw(z->v1.i, s); break;
+case BCF_BT_INT32: if ( z->v1.i==bcf_int32_missing ) kputc('.', s); else kputw(z->v1.i, s); break;
+case BCF_BT_FLOAT: if ( bcf_float_is_missing(z->v1.f) ) kputc('.', s); else ksprintf(s, "%g", z->v1.f); break;
+case BCF_BT_CHAR:  kputc(z->v1.i, s); break;
+default: fprintf(stderr,"todo: type %d\n", z->type); exit(1); break;
+}
+}
+else bcf_fmt_array(s, z->len, z->type, z->vptr);
+}
+if ( first ) kputc('.', s);
+} else kputc('.', s);
+// FORMAT and individual information
+if (v->n_sample)
+{
+int i,j;
+if ( v->n_fmt)
+{
+int gt_i = -1;
+bcf_fmt_t *fmt = v->d.fmt;
+int first = 1;
+for (i = 0; i < (int)v->n_fmt; ++i) {
+if ( !fmt[i].p ) continue;
+kputc(!first ? ':' : '\t', s); first = 0;
+if ( fmt[i].id<0 ) //!bcf_hdr_idinfo_exists(h,BCF_HL_FMT,fmt[i].id) )
+{
+fprintf(stderr, "[E::%s] invalid BCF, the FORMAT tag id=%d not present in the header.\n", __func__, fmt[i].id);
+abort();
+}
+kputs(h->id[BCF_DT_ID][fmt[i].id].key, s);
+if (strcmp(h->id[BCF_DT_ID][fmt[i].id].key, "GT") == 0) gt_i = i;
+}
+if ( first ) kputs("\t.", s);
+for (j = 0; j < v->n_sample; ++j) {
+kputc('\t', s);
+first = 1;
+for (i = 0; i < (int)v->n_fmt; ++i) {
+bcf_fmt_t *f = &fmt[i];
+if ( !f->p ) continue;
+if (!first) kputc(':', s); first = 0;
+if (gt_i == i)
+bcf_format_gt(f,j,s);
+else
+bcf_fmt_array(s, f->n, f->type, f->p + j * f->size);
+}
+if ( first ) kputc('.', s);
+}
+}
+else
+for (j=0; j<=v->n_sample; j++)
+kputs("\t.", s);
+}
+kputc('\n', s);
+return 0;
+}
+int vcf_write_line(htsFile *fp, kstring_t *line)
+{
+int ret;
+if ( line->s[line->l-1]!='\n' ) kputc('\n',line);
+if ( fp->format.compression!=no_compression )
+ret = bgzf_write(fp->fp.bgzf, line->s, line->l);
+else
+ret = hwrite(fp->fp.hfile, line->s, line->l);
+return ret==line->l ? 0 : -1;
+}
+int vcf_write(htsFile *fp, const bcf_hdr_t *h, bcf1_t *v)
+{
+int ret;
+fp->line.l = 0;
+vcf_format1(h, v, &fp->line);
+if ( fp->format.compression!=no_compression )
+ret = bgzf_write(fp->fp.bgzf, fp->line.s, fp->line.l);
+else
+ret = hwrite(fp->fp.hfile, fp->line.s, fp->line.l);
+return ret==fp->line.l ? 0 : -1;
+}
+/************************
+* Data access routines *
+************************/
+int bcf_hdr_id2int(const bcf_hdr_t *h, int which, const char *id)
+{
+khint_t k;
+vdict_t *d = (vdict_t*)h->dict[which];
+k = kh_get(vdict, d, id);
+return k == kh_end(d)? -1 : kh_val(d, k).id;
+}
+/********************
+*** BCF indexing ***
+********************/
+hts_idx_t *bcf_index(htsFile *fp, int min_shift)
+{
+int n_lvls, i;
+bcf1_t *b;
+hts_idx_t *idx;
+bcf_hdr_t *h;
+int64_t max_len = 0, s;
+h = bcf_hdr_read(fp);
+if ( !h ) return NULL;
+int nids = 0;
+for (i = 0; i < h->n[BCF_DT_CTG]; ++i)
+{
+if ( !h->id[BCF_DT_CTG][i].val ) continue;
+if ( max_len < h->id[BCF_DT_CTG][i].val->info[0] ) max_len = h->id[BCF_DT_CTG][i].val->info[0];
+nids++;
+}
+if ( !max_len ) max_len = ((int64_t)1<<31) - 1;  // In case contig line is broken.
+max_len += 256;
+for (n_lvls = 0, s = 1<<min_shift; max_len > s; ++n_lvls, s <<= 3);
+idx = hts_idx_init(nids, HTS_FMT_CSI, bgzf_tell(fp->fp.bgzf), min_shift, n_lvls);
+b = bcf_init1();
+while (bcf_read1(fp,h, b) >= 0) {
+int ret;
+ret = hts_idx_push(idx, b->rid, b->pos, b->pos + b->rlen, bgzf_tell(fp->fp.bgzf), 1);
+if (ret < 0)
+{
+bcf_destroy1(b);
+hts_idx_destroy(idx);
+return NULL;
+}
+}
+hts_idx_finish(idx, bgzf_tell(fp->fp.bgzf));
+bcf_destroy1(b);
+bcf_hdr_destroy(h);
+return idx;
+}
+int bcf_index_build(const char *fn, int min_shift)
+{
+htsFile *fp;
+hts_idx_t *idx;
+if ((fp = hts_open(fn, "rb")) == 0) return -1;
+if ( fp->format.compression!=bgzf ) { hts_close(fp); return -1; }
+idx = bcf_index(fp, min_shift);
+hts_close(fp);
+if ( !idx ) return -1;
+hts_idx_save(idx, fn, HTS_FMT_CSI);
+hts_idx_destroy(idx);
+return 0;
+}
+/*****************
+*** Utilities ***
+*****************/
+int bcf_hdr_combine(bcf_hdr_t *dst, const bcf_hdr_t *src)
+{
+int i, ndst_ori = dst->nhrec, need_sync = 0, ret = 0;
+for (i=0; i<src->nhrec; i++)
+{
+if ( src->hrec[i]->type==BCF_HL_GEN && src->hrec[i]->value )
+{
+int j;
+for (j=0; j<ndst_ori; j++)
+{
+if ( dst->hrec[j]->type!=BCF_HL_GEN ) continue;
+// Checking only the key part of generic lines, otherwise
+// the VCFs are too verbose. Should we perhaps add a flag
+// to bcf_hdr_combine() and make this optional?
+if ( !strcmp(src->hrec[i]->key,dst->hrec[j]->key) ) break;
+}
+if ( j>=ndst_ori )
+need_sync += bcf_hdr_add_hrec(dst, bcf_hrec_dup(src->hrec[i]));
+}
+else if ( src->hrec[i]->type==BCF_HL_STR )
+{
+// NB: we are ignoring fields without ID
+int j = bcf_hrec_find_key(src->hrec[i],"ID");
+if ( j>=0 )
+{
+bcf_hrec_t *rec = bcf_hdr_get_hrec(dst, src->hrec[i]->type, "ID", src->hrec[i]->vals[j], src->hrec[i]->key);
+if ( !rec )
+need_sync += bcf_hdr_add_hrec(dst, bcf_hrec_dup(src->hrec[i]));
+}
+}
+else
+{
+int j = bcf_hrec_find_key(src->hrec[i],"ID");
+assert( j>=0 ); // this should always be true for valid VCFs
+bcf_hrec_t *rec = bcf_hdr_get_hrec(dst, src->hrec[i]->type, "ID", src->hrec[i]->vals[j], NULL);
+if ( !rec )
+need_sync += bcf_hdr_add_hrec(dst, bcf_hrec_dup(src->hrec[i]));
+else if ( src->hrec[i]->type==BCF_HL_INFO || src->hrec[i]->type==BCF_HL_FMT )
+{
+// Check that both records are of the same type. The bcf_hdr_id2length
+// macro cannot be used here because dst header is not synced yet.
+vdict_t *d_src = (vdict_t*)src->dict[BCF_DT_ID];
+vdict_t *d_dst = (vdict_t*)dst->dict[BCF_DT_ID];
+khint_t k_src  = kh_get(vdict, d_src, src->hrec[i]->vals[0]);
+khint_t k_dst  = kh_get(vdict, d_dst, src->hrec[i]->vals[0]);
+if ( (kh_val(d_src,k_src).info[rec->type]>>8 & 0xf) != (kh_val(d_dst,k_dst).info[rec->type]>>8 & 0xf) )
+{
+fprintf(stderr,"Warning: trying to combine \"%s\" tag definitions of different lengths\n", src->hrec[i]->vals[0]);
+ret |= 1;
+}
+if ( (kh_val(d_src,k_src).info[rec->type]>>4 & 0xf) != (kh_val(d_dst,k_dst).info[rec->type]>>4 & 0xf) )
+{
+fprintf(stderr,"Warning: trying to combine \"%s\" tag definitions of different types\n", src->hrec[i]->vals[0]);
+ret |= 1;
+}
+}
+}
+}
+if ( need_sync ) bcf_hdr_sync(dst);
+return ret;
+}
+int bcf_translate(const bcf_hdr_t *dst_hdr, bcf_hdr_t *src_hdr, bcf1_t *line)
+{
+int i;
+if ( line->errcode )
+{
+fprintf(stderr,"[%s:%d %s] Unchecked error (%d), exiting.\n", __FILE__,__LINE__,__FUNCTION__,line->errcode);
+exit(1);
+}
+if ( src_hdr->ntransl==-1 ) return 0;    // no need to translate, all tags have the same id
+if ( !src_hdr->ntransl )  // called for the first time, see what needs translating
+{
+int dict;
+for (dict=0; dict<2; dict++)    // BCF_DT_ID and BCF_DT_CTG
+{
+src_hdr->transl[dict] = (int*) malloc(src_hdr->n[dict]*sizeof(int));
+for (i=0; i<src_hdr->n[dict]; i++)
+{
+if ( !src_hdr->id[dict][i].key || !dst_hdr->id[dict][i].key )    // gap left after removed BCF header lines
+src_hdr->transl[dict][i] = -1;
+else if ( i>=dst_hdr->n[dict] || strcmp(src_hdr->id[dict][i].key,dst_hdr->id[dict][i].key) )
+{
+src_hdr->transl[dict][i] = bcf_hdr_id2int(dst_hdr,dict,src_hdr->id[dict][i].key);
+src_hdr->ntransl++;
+}
+else
+src_hdr->transl[dict][i] = -1;
+}
+}
+if ( !src_hdr->ntransl )
+{
+free(src_hdr->transl[0]); src_hdr->transl[0] = NULL;
+free(src_hdr->transl[1]); src_hdr->transl[1] = NULL;
+src_hdr->ntransl = -1;
+}
+if ( src_hdr->ntransl==-1 ) return 0;
+}
+bcf_unpack(line,BCF_UN_ALL);
+// CHROM
+if ( src_hdr->transl[BCF_DT_CTG][line->rid] >=0 ) line->rid = src_hdr->transl[BCF_DT_CTG][line->rid];
+// FILTER
+for (i=0; i<line->d.n_flt; i++)
+{
+int src_id = line->d.flt[i];
+if ( src_hdr->transl[BCF_DT_ID][src_id] >=0 )
+line->d.flt[i] = src_hdr->transl[BCF_DT_ID][src_id];
+line->d.shared_dirty |= BCF1_DIRTY_FLT;
+}
+// INFO
+for (i=0; i<line->n_info; i++)
+{
+int src_id = line->d.info[i].key;
+int dst_id = src_hdr->transl[BCF_DT_ID][src_id];
+if ( dst_id<0 ) continue;
+int src_size = src_id>>7 ? ( src_id>>15 ? BCF_BT_INT32 : BCF_BT_INT16) : BCF_BT_INT8;
+int dst_size = dst_id>>7 ? ( dst_id>>15 ? BCF_BT_INT32 : BCF_BT_INT16) : BCF_BT_INT8;
+if ( src_size==dst_size )   // can overwrite
+{
+line->d.info[i].key = dst_id;
+uint8_t *vptr = line->d.info[i].vptr - line->d.info[i].vptr_off;
+if ( dst_size==BCF_BT_INT8 ) { vptr[1] = (uint8_t)dst_id; }
+else if ( dst_size==BCF_BT_INT16 ) { *(uint16_t*)vptr = (uint16_t)dst_id; }
+else { *(uint32_t*)vptr = (uint32_t)dst_id; }
+}
+else    // must realloc
+{
+bcf_info_t *info = &line->d.info[i];
+assert( !info->vptr_free );
+kstring_t str = {0,0,0};
+bcf_enc_int1(&str, dst_id);
+bcf_enc_size(&str, info->len,info->type);
+info->vptr_off = str.l;
+kputsn((char*)info->vptr, info->vptr_len, &str);
+info->vptr = (uint8_t*)str.s + info->vptr_off;
+info->vptr_free = 1;
+info->key = dst_id;
+line->d.shared_dirty |= BCF1_DIRTY_INF;
+}
+}
+// FORMAT
+for (i=0; i<line->n_fmt; i++)
+{
+int src_id = line->d.fmt[i].id;
+int dst_id = src_hdr->transl[BCF_DT_ID][src_id];
+if ( dst_id<0 ) continue;
+int src_size = src_id>>7 ? ( src_id>>15 ? BCF_BT_INT32 : BCF_BT_INT16) : BCF_BT_INT8;
+int dst_size = dst_id>>7 ? ( dst_id>>15 ? BCF_BT_INT32 : BCF_BT_INT16) : BCF_BT_INT8;
+if ( src_size==dst_size )   // can overwrite
+{
+line->d.fmt[i].id = dst_id;
+uint8_t *p = line->d.fmt[i].p - line->d.fmt[i].p_off;    // pointer to the vector size (4bits) and BT type (4bits)
+if ( dst_size==BCF_BT_INT8 ) { p[1] = dst_id; }
+else if ( dst_size==BCF_BT_INT16 ) { uint8_t *x = (uint8_t*) &dst_id; p[1] = x[0]; p[2] = x[1]; }
+else { uint8_t *x = (uint8_t*) &dst_id; p[1] = x[0]; p[2] = x[1]; p[3] = x[2]; p[4] = x[3]; }
+}
+else    // must realloc
+{
+bcf_fmt_t *fmt = &line->d.fmt[i];
+assert( !fmt->p_free );
+kstring_t str = {0,0,0};
+bcf_enc_int1(&str, dst_id);
+bcf_enc_size(&str, fmt->n, fmt->type);
+fmt->p_off = str.l;
+kputsn((char*)fmt->p, fmt->p_len, &str);
+fmt->p = (uint8_t*)str.s + fmt->p_off;
+fmt->p_free = 1;
+fmt->id = dst_id;
+line->d.indiv_dirty = 1;
+}
+}
+return 0;
+}
+bcf_hdr_t *bcf_hdr_dup(const bcf_hdr_t *hdr)
+{
+bcf_hdr_t *hout = bcf_hdr_init("r");
+char *htxt = bcf_hdr_fmt_text(hdr, 1, NULL);
+bcf_hdr_parse(hout, htxt);
+free(htxt);
+return hout;
+}
+bcf_hdr_t *bcf_hdr_subset(const bcf_hdr_t *h0, int n, char *const* samples, int *imap)
+{
+int hlen;
+void *names_hash = khash_str2int_init();
+char *htxt = bcf_hdr_fmt_text(h0, 1, &hlen);
+kstring_t str;
+bcf_hdr_t *h;
+str.l = str.m = 0; str.s = 0;
+h = bcf_hdr_init("w");
+bcf_hdr_set_version(h,bcf_hdr_get_version(h0));
+int j;
+for (j=0; j<n; j++) imap[j] = -1;
+if ( bcf_hdr_nsamples(h0) > 0) {
+char *p;
+int i = 0, end = n? 8 : 7;
+while ((p = strstr(htxt, "#CHROM\t")) != 0)
+if (p > htxt && *(p-1) == '\n') break;
+while ((p = strchr(p, '\t')) != 0 && i < end) ++i, ++p;
+if (i != end) {
+free(h); free(str.s);
+return 0; // malformated header
+}
+kputsn(htxt, p - htxt, &str);
+for (i = 0; i < n; ++i) {
+if ( khash_str2int_has_key(names_hash,samples[i]) )
+{
+fprintf(stderr,"[E::bcf_hdr_subset] Duplicate sample name \"%s\".\n", samples[i]);
+free(str.s);
+free(htxt);
+khash_str2int_destroy(names_hash);
+bcf_hdr_destroy(h);
+return NULL;
+}
+imap[i] = bcf_hdr_id2int(h0, BCF_DT_SAMPLE, samples[i]);
+if (imap[i] < 0) continue;
+kputc('\t', &str);
+kputs(samples[i], &str);
+khash_str2int_inc(names_hash,samples[i]);
+}
+} else kputsn(htxt, hlen, &str);
+while (str.l && (!str.s[str.l-1] || str.s[str.l-1]=='\n') ) str.l--; // kill trailing zeros and newlines
+kputc('\n',&str);
+bcf_hdr_parse(h, str.s);
+free(str.s);
+free(htxt);
+khash_str2int_destroy(names_hash);
+return h;
+}
+int bcf_hdr_set_samples(bcf_hdr_t *hdr, const char *samples, int is_file)
+{
+if ( samples && !strcmp("-",samples) ) return 0;            // keep all samples
+hdr->nsamples_ori = bcf_hdr_nsamples(hdr);
+if ( !samples ) { bcf_hdr_nsamples(hdr) = 0; return 0; }    // exclude all samples
+int i, narr = bit_array_size(bcf_hdr_nsamples(hdr));
+hdr->keep_samples = (uint8_t*) calloc(narr,1);
+if ( samples[0]=='^' )
+for (i=0; i<bcf_hdr_nsamples(hdr); i++) bit_array_set(hdr->keep_samples,i);
+int idx, n, ret = 0;
+char **smpls = hts_readlist(samples[0]=='^'?samples+1:samples, is_file, &n);
+if ( !smpls ) return -1;
+for (i=0; i<n; i++)
+{
+idx = bcf_hdr_id2int(hdr,BCF_DT_SAMPLE,smpls[i]);
+if ( idx<0 )
+{
+if ( !ret ) ret = i+1;
+continue;
+}
+assert( idx<bcf_hdr_nsamples(hdr) );
+if (  samples[0]=='^' )
+bit_array_clear(hdr->keep_samples, idx);
+else
+bit_array_set(hdr->keep_samples, idx);
+}
+for (i=0; i<n; i++) free(smpls[i]);
+free(smpls);
+bcf_hdr_nsamples(hdr) = 0;
+for (i=0; i<hdr->nsamples_ori; i++)
+if ( bit_array_test(hdr->keep_samples,i) ) bcf_hdr_nsamples(hdr)++;
+if ( !bcf_hdr_nsamples(hdr) ) { free(hdr->keep_samples); hdr->keep_samples=NULL; }
+else
+{
+char **samples = (char**) malloc(sizeof(char*)*bcf_hdr_nsamples(hdr));
+idx = 0;
+for (i=0; i<hdr->nsamples_ori; i++)
+if ( bit_array_test(hdr->keep_samples,i) ) samples[idx++] = strdup(hdr->samples[i]);
+free(hdr->samples);
+hdr->samples = samples;
+// delete original samples from the dictionary
+vdict_t *d = (vdict_t*)hdr->dict[BCF_DT_SAMPLE];
+int k;
+for (k = kh_begin(d); k != kh_end(d); ++k)
+if (kh_exist(d, k)) free((char*)kh_key(d, k));
+kh_destroy(vdict, d);
+// add the subset back
+hdr->dict[BCF_DT_SAMPLE] = d = kh_init(vdict);
+for (i=0; i<bcf_hdr_nsamples(hdr); i++)
+{
+int ignore, k = kh_put(vdict, d, hdr->samples[i], &ignore);
+kh_val(d, k) = bcf_idinfo_def;
+kh_val(d, k).id = kh_size(d) - 1;
+}
+bcf_hdr_sync(hdr);
+}
+return ret;
+}
+int bcf_subset(const bcf_hdr_t *h, bcf1_t *v, int n, int *imap)
+{
+kstring_t ind;
+ind.s = 0; ind.l = ind.m = 0;
+if (n) {
+bcf_fmt_t *fmt;
+int i, j;
+fmt = (bcf_fmt_t*)alloca(v->n_fmt * sizeof(bcf_fmt_t));
+uint8_t *ptr = (uint8_t*)v->indiv.s;
+for (i = 0; i < v->n_fmt; ++i)
+ptr = bcf_unpack_fmt_core1(ptr, v->n_sample, &fmt[i]);
+for (i = 0; i < (int)v->n_fmt; ++i) {
+bcf_fmt_t *f = &fmt[i];
+bcf_enc_int1(&ind, f->id);
+bcf_enc_size(&ind, f->n, f->type);
+for (j = 0; j < n; ++j)
+if (imap[j] >= 0) kputsn((char*)(f->p + imap[j] * f->size), f->size, &ind);
+}
+for (i = j = 0; j < n; ++j) if (imap[j] >= 0) ++i;
+v->n_sample = i;
+} else v->n_sample = 0;
+if ( !v->n_sample ) v->n_fmt = 0;
+free(v->indiv.s);
+v->indiv = ind;
+v->unpacked &= ~BCF_UN_FMT;    // only BCF is ready for output, VCF will need to unpack again
+return 0;
+}
+int bcf_is_snp(bcf1_t *v)
+{
+int i;
+bcf_unpack(v, BCF_UN_STR);
+for (i = 0; i < v->n_allele; ++i)
+{
+if ( v->d.allele[i][1]==0 ) continue;
+// mpileup's <X> allele, see also below. This is not completely satisfactory,
+// a general library is here narrowly tailored to fit samtools.
+if ( v->d.allele[i][0]=='<' && v->d.allele[i][1]=='X' && v->d.allele[i][2]=='>' ) continue;
+break;
+}
+return i == v->n_allele;
+}
+static void bcf_set_variant_type(const char *ref, const char *alt, variant_t *var)
+{
+// The most frequent case
+if ( !ref[1] && !alt[1] )
+{
+if ( *alt == '.' || *ref==*alt ) { var->n = 0; var->type = VCF_REF; return; }
+if ( *alt == 'X' ) { var->n = 0; var->type = VCF_REF; return; }  // mpileup's X allele shouldn't be treated as variant
+var->n = 1; var->type = VCF_SNP; return;
+}
+if ( alt[0]=='<' )
+{
+if ( alt[1]=='X' && alt[2]=='>' ) { var->n = 0; var->type = VCF_REF; return; }  // mpileup's X allele shouldn't be treated as variant
+var->type = VCF_OTHER;
+return;
+}
+const char *r = ref, *a = alt;
+while (*r && *a && *r==*a ) { r++; a++; }
+if ( *a && !*r )
+{
+while ( *a ) a++;
+var->n = (a-alt)-(r-ref); var->type = VCF_INDEL; return;
+}
+else if ( *r && !*a )
+{
+while ( *r ) r++;
+var->n = (a-alt)-(r-ref); var->type = VCF_INDEL; return;
+}
+else if ( !*r && !*a )
+{
+var->n = 0; var->type = VCF_REF; return;
+}
+const char *re = r, *ae = a;
+while ( re[1] ) re++;
+while ( ae[1] ) ae++;
+while ( *re==*ae && re>r && ae>a ) { re--; ae--; }
+if ( ae==a )
+{
+if ( re==r ) { var->n = 1; var->type = VCF_SNP; return; }
+var->n = -(re-r);
+if ( *re==*ae ) { var->type = VCF_INDEL; return; }
+var->type = VCF_OTHER; return;
+}
+else if ( re==r )
+{
+var->n = ae-a;
+if ( *re==*ae ) { var->type = VCF_INDEL; return; }
+var->type = VCF_OTHER; return;
+}
+var->type = ( re-r == ae-a ) ? VCF_MNP : VCF_OTHER;
+var->n = ( re-r > ae-a ) ? -(re-r+1) : ae-a+1;
+// should do also complex events, SVs, etc...
+}
+static void bcf_set_variant_types(bcf1_t *b)
+{
+if ( !(b->unpacked & BCF_UN_STR) ) bcf_unpack(b, BCF_UN_STR);
+bcf_dec_t *d = &b->d;
+if ( d->n_var < b->n_allele )
+{
+d->var = (variant_t *) realloc(d->var, sizeof(variant_t)*b->n_allele);
+d->n_var = b->n_allele;
+}
+int i;
+b->d.var_type = 0;
+for (i=1; i<b->n_allele; i++)
+{
+bcf_set_variant_type(d->allele[0],d->allele[i], &d->var[i]);
+b->d.var_type |= d->var[i].type;
+//fprintf(stderr,"[set_variant_type] %d   %s %s -> %d %d .. %d\n", b->pos+1,d->allele[0],d->allele[i],d->var[i].type,d->var[i].n, b->d.var_type);
+}
+}
+int bcf_get_variant_types(bcf1_t *rec)
+{
+if ( rec->d.var_type==-1 ) bcf_set_variant_types(rec);
+return rec->d.var_type;
+}
+int bcf_get_variant_type(bcf1_t *rec, int ith_allele)
+{
+if ( rec->d.var_type==-1 ) bcf_set_variant_types(rec);
+return rec->d.var[ith_allele].type;
+}
+int bcf_update_info(const bcf_hdr_t *hdr, bcf1_t *line, const char *key, const void *values, int n, int type)
+{
+// Is the field already present?
+int i, inf_id = bcf_hdr_id2int(hdr,BCF_DT_ID,key);
+if ( !bcf_hdr_idinfo_exists(hdr,BCF_HL_INFO,inf_id) ) return -1;    // No such INFO field in the header
+if ( !(line->unpacked & BCF_UN_INFO) ) bcf_unpack(line, BCF_UN_INFO);
+for (i=0; i<line->n_info; i++)
+if ( inf_id==line->d.info[i].key ) break;
+bcf_info_t *inf = i==line->n_info ? NULL : &line->d.info[i];
+if ( !n || (type==BCF_HT_STR && !values) )
+{
+if ( inf )
+{
+// Mark the tag for removal, free existing memory if necessary
+if ( inf->vptr_free )
+{
+free(inf->vptr - inf->vptr_off);
+inf->vptr_free = 0;
+}
+line->d.shared_dirty |= BCF1_DIRTY_INF;
+inf->vptr = NULL;
+}
+return 0;
+}
+// Encode the values and determine the size required to accommodate the values
+kstring_t str = {0,0,0};
+bcf_enc_int1(&str, inf_id);
+if ( type==BCF_HT_INT )
+bcf_enc_vint(&str, n, (int32_t*)values, -1);
+else if ( type==BCF_HT_REAL )
+bcf_enc_vfloat(&str, n, (float*)values);
+else if ( type==BCF_HT_FLAG || type==BCF_HT_STR )
+{
+if ( values==NULL )
+bcf_enc_size(&str, 0, BCF_BT_NULL);
+else
+bcf_enc_vchar(&str, strlen((char*)values), (char*)values);
+}
+else
+{
+fprintf(stderr, "[E::%s] the type %d not implemented yet\n", __func__, type);
+abort();
+}
+// Is the INFO tag already present
+if ( inf )
+{
+// Is it big enough to accommodate new block?
+if ( str.l <= inf->vptr_len + inf->vptr_off )
+{
+if ( str.l != inf->vptr_len + inf->vptr_off ) line->d.shared_dirty |= BCF1_DIRTY_INF;
+uint8_t *ptr = inf->vptr - inf->vptr_off;
+memcpy(ptr, str.s, str.l);
+free(str.s);
+int vptr_free = inf->vptr_free;
+bcf_unpack_info_core1(ptr, inf);
+inf->vptr_free = vptr_free;
+}
+else
+{
+assert( !inf->vptr_free );  // fix the caller or improve here: this has been modified before
+bcf_unpack_info_core1((uint8_t*)str.s, inf);
+inf->vptr_free = 1;
+line->d.shared_dirty |= BCF1_DIRTY_INF;
+}
+}
+else
+{
+// The tag is not present, create new one
+line->n_info++;
+hts_expand0(bcf_info_t, line->n_info, line->d.m_info , line->d.info);
+inf = &line->d.info[line->n_info-1];
+bcf_unpack_info_core1((uint8_t*)str.s, inf);
+inf->vptr_free = 1;
+line->d.shared_dirty |= BCF1_DIRTY_INF;
+}
+line->unpacked |= BCF_UN_INFO;
+return 0;
+}
+int bcf_update_format_string(const bcf_hdr_t *hdr, bcf1_t *line, const char *key, const char **values, int n)
+{
+if ( !n )
+return bcf_update_format(hdr,line,key,NULL,0,BCF_HT_STR);
+int i, max_len = 0;
+for (i=0; i<n; i++)
+{
+int len = strlen(values[i]);
+if ( len > max_len ) max_len = len;
+}
+char *out = (char*) malloc(max_len*n);
+if ( !out ) return -2;
+for (i=0; i<n; i++)
+{
+char *dst = out+i*max_len;
+const char *src = values[i];
+int j = 0;
+while ( src[j] ) { dst[j] = src[j]; j++; }
+for (; j<max_len; j++) dst[j] = 0;
+}
+int ret = bcf_update_format(hdr,line,key,out,max_len*n,BCF_HT_STR);
+free(out);
+return ret;
+}
+int bcf_update_format(const bcf_hdr_t *hdr, bcf1_t *line, const char *key, const void *values, int n, int type)
+{
+// Is the field already present?
+int i, fmt_id = bcf_hdr_id2int(hdr,BCF_DT_ID,key);
+if ( !bcf_hdr_idinfo_exists(hdr,BCF_HL_FMT,fmt_id) )
+{
+if ( !n ) return 0;
+return -1;  // the key not present in the header
+}
+if ( !(line->unpacked & BCF_UN_FMT) ) bcf_unpack(line, BCF_UN_FMT);
+for (i=0; i<line->n_fmt; i++)
+if ( line->d.fmt[i].id==fmt_id ) break;
+bcf_fmt_t *fmt = i==line->n_fmt ? NULL : &line->d.fmt[i];
+if ( !n )
+{
+if ( fmt )
+{
+// Mark the tag for removal, free existing memory if necessary
+if ( fmt->p_free )
+{
+free(fmt->p - fmt->p_off);
+fmt->p_free = 0;
+}
+line->d.indiv_dirty = 1;
+fmt->p = NULL;
+}
+return 0;
+}
+line->n_sample = bcf_hdr_nsamples(hdr);
+int nps = n / line->n_sample;  // number of values per sample
+assert( nps && nps*line->n_sample==n );     // must be divisible by n_sample
+// Encode the values and determine the size required to accommodate the values
+kstring_t str = {0,0,0};
+bcf_enc_int1(&str, fmt_id);
+if ( type==BCF_HT_INT )
+bcf_enc_vint(&str, n, (int32_t*)values, nps);
+else if ( type==BCF_HT_REAL )
+{
+bcf_enc_size(&str, nps, BCF_BT_FLOAT);
+kputsn((char*)values, nps*line->n_sample*sizeof(float), &str);
+}
+else if ( type==BCF_HT_STR )
+{
+bcf_enc_size(&str, nps, BCF_BT_CHAR);
+kputsn((char*)values, nps*line->n_sample, &str);
+}
+else
+{
+fprintf(stderr, "[E::%s] the type %d not implemented yet\n", __func__, type);
+abort();
+}
+if ( !fmt )
+{
+// Not present, new format field
+line->n_fmt++;
+hts_expand0(bcf_fmt_t, line->n_fmt, line->d.m_fmt, line->d.fmt);
+// Special case: VCF specification requires that GT is always first
+if ( line->n_fmt > 1 && key[0]=='G' && key[1]=='T' && !key[2] )
+{
+for (i=line->n_fmt-1; i>0; i--)
+line->d.fmt[i] = line->d.fmt[i-1];
+fmt = &line->d.fmt[0];
+}
+else
+fmt = &line->d.fmt[line->n_fmt-1];
+bcf_unpack_fmt_core1((uint8_t*)str.s, line->n_sample, fmt);
+line->d.indiv_dirty = 1;
+fmt->p_free = 1;
+}
+else
+{
+// The tag is already present, check if it is big enough to accomodate the new block
+if ( str.l <= fmt->p_len + fmt->p_off )
+{
+// good, the block is big enough
+if ( str.l != fmt->p_len + fmt->p_off ) line->d.indiv_dirty = 1;
+uint8_t *ptr = fmt->p - fmt->p_off;
+memcpy(ptr, str.s, str.l);
+free(str.s);
+int p_free = fmt->p_free;
+bcf_unpack_fmt_core1(ptr, line->n_sample, fmt);
+fmt->p_free = p_free;
+}
+else
+{
+assert( !fmt->p_free );  // fix the caller or improve here: this has been modified before
+bcf_unpack_fmt_core1((uint8_t*)str.s, line->n_sample, fmt);
+fmt->p_free = 1;
+line->d.indiv_dirty = 1;
+}
+}
+line->unpacked |= BCF_UN_FMT;
+return 0;
+}
+int bcf_update_filter(const bcf_hdr_t *hdr, bcf1_t *line, int *flt_ids, int n)
+{
+if ( !(line->unpacked & BCF_UN_FLT) ) bcf_unpack(line, BCF_UN_FLT);
+line->d.shared_dirty |= BCF1_DIRTY_FLT;
+line->d.n_flt = n;
+if ( !n ) return 0;
+hts_expand(int, line->d.n_flt, line->d.m_flt, line->d.flt);
+int i;
+for (i=0; i<n; i++)
+line->d.flt[i] = flt_ids[i];
+return 0;
+}
+int bcf_add_filter(const bcf_hdr_t *hdr, bcf1_t *line, int flt_id)
+{
+if ( !(line->unpacked & BCF_UN_FLT) ) bcf_unpack(line, BCF_UN_FLT);
+int i;
+for (i=0; i<line->d.n_flt; i++)
+if ( flt_id==line->d.flt[i] ) break;
+if ( i<line->d.n_flt ) return 0;    // this filter is already set
+line->d.shared_dirty |= BCF1_DIRTY_FLT;
+if ( flt_id==0 )    // set to PASS
+line->d.n_flt = 1;
+else if ( line->d.n_flt==1 && line->d.flt[0]==0 )
+line->d.n_flt = 1;
+else
+line->d.n_flt++;
+hts_expand(int, line->d.n_flt, line->d.m_flt, line->d.flt);
+line->d.flt[line->d.n_flt-1] = flt_id;
+return 1;
+}
+int bcf_remove_filter(const bcf_hdr_t *hdr, bcf1_t *line, int flt_id, int pass)
+{
+if ( !(line->unpacked & BCF_UN_FLT) ) bcf_unpack(line, BCF_UN_FLT);
+int i;
+for (i=0; i<line->d.n_flt; i++)
+if ( flt_id==line->d.flt[i] ) break;
+if ( i==line->d.n_flt ) return 0;   // the filter is not present
+line->d.shared_dirty |= BCF1_DIRTY_FLT;
+if ( i!=line->d.n_flt-1 ) memmove(line->d.flt+i,line->d.flt+i+1,(line->d.n_flt-i-1)*sizeof(*line->d.flt));
+line->d.n_flt--;
+if ( !line->d.n_flt && pass ) bcf_add_filter(hdr,line,0);
+return 0;
+}
+int bcf_has_filter(const bcf_hdr_t *hdr, bcf1_t *line, char *filter)
+{
+if ( filter[0]=='.' && !filter[1] ) filter = "PASS";
+int id = bcf_hdr_id2int(hdr, BCF_DT_ID, filter);
+if ( !bcf_hdr_idinfo_exists(hdr,BCF_HL_FLT,id) ) return -1;  // not defined in the header
+if ( !(line->unpacked & BCF_UN_FLT) ) bcf_unpack(line, BCF_UN_FLT);
+if ( id==0 && !line->d.n_flt) return 1; // PASS
+int i;
+for (i=0; i<line->d.n_flt; i++)
+if ( line->d.flt[i]==id ) return 1;
+return 0;
+}
+static inline int _bcf1_sync_alleles(const bcf_hdr_t *hdr, bcf1_t *line, int nals)
+{
+line->d.shared_dirty |= BCF1_DIRTY_ALS;
+line->n_allele = nals;
+hts_expand(char*, line->n_allele, line->d.m_allele, line->d.allele);
+char *als = line->d.als;
+int n = 0;
+while (n<nals)
+{
+line->d.allele[n] = als;
+while ( *als ) als++;
+als++;
+n++;
+}
+return 0;
+}
+int bcf_update_alleles(const bcf_hdr_t *hdr, bcf1_t *line, const char **alleles, int nals)
+{
+kstring_t tmp = {0,0,0};
+char *free_old = NULL;
+// If the supplied alleles are not pointers to line->d.als, the existing block can be reused.
+int i;
+for (i=0; i<nals; i++)
+if ( alleles[i]>=line->d.als && alleles[i]<line->d.als+line->d.m_als ) break;
+if ( i==nals )
+{
+// all alleles point elsewhere, reuse the existing block
+tmp.l = 0; tmp.s = line->d.als; tmp.m = line->d.m_als;
+}
+else
+free_old = line->d.als;
+for (i=0; i<nals; i++)
+{
+kputs(alleles[i], &tmp);
+kputc(0, &tmp);
+}
+line->d.als = tmp.s; line->d.m_als = tmp.m;
+free(free_old);
+return _bcf1_sync_alleles(hdr,line,nals);
+}
+int bcf_update_alleles_str(const bcf_hdr_t *hdr, bcf1_t *line, const char *alleles_string)
+{
+kstring_t tmp;
+tmp.l = 0; tmp.s = line->d.als; tmp.m = line->d.m_als;
+kputs(alleles_string, &tmp);
+line->d.als = tmp.s; line->d.m_als = tmp.m;
+int nals = 1;
+char *t = line->d.als;
+while (*t)
+{
+if ( *t==',' ) { *t = 0; nals++; }
+t++;
+}
+return _bcf1_sync_alleles(hdr, line, nals);
+}
+int bcf_update_id(const bcf_hdr_t *hdr, bcf1_t *line, const char *id)
+{
+kstring_t tmp;
+tmp.l = 0; tmp.s = line->d.id; tmp.m = line->d.m_id;
+if ( id )
+kputs(id, &tmp);
+else
+kputs(".", &tmp);
+line->d.id = tmp.s; line->d.m_id = tmp.m;
+line->d.shared_dirty |= BCF1_DIRTY_ID;
+return 0;
+}
+bcf_fmt_t *bcf_get_fmt(const bcf_hdr_t *hdr, bcf1_t *line, const char *key)
+{
+int id = bcf_hdr_id2int(hdr, BCF_DT_ID, key);
+if ( !bcf_hdr_idinfo_exists(hdr,BCF_HL_FMT,id) ) return NULL;   // no such FMT field in the header
+return bcf_get_fmt_id(line, id);
+}
+bcf_info_t *bcf_get_info(const bcf_hdr_t *hdr, bcf1_t *line, const char *key)
+{
+int id = bcf_hdr_id2int(hdr, BCF_DT_ID, key);
+if ( !bcf_hdr_idinfo_exists(hdr,BCF_HL_INFO,id) ) return NULL;   // no such INFO field in the header
+return bcf_get_info_id(line, id);
+}
+bcf_fmt_t *bcf_get_fmt_id(bcf1_t *line, const int id)
+{
+int i;
+if ( !(line->unpacked & BCF_UN_FMT) ) bcf_unpack(line, BCF_UN_FMT);
+for (i=0; i<line->n_fmt; i++)
+{
+if ( line->d.fmt[i].id==id ) return &line->d.fmt[i];
+}
+return NULL;
+}
+bcf_info_t *bcf_get_info_id(bcf1_t *line, const int id)
+{
+int i;
+if ( !(line->unpacked & BCF_UN_INFO) ) bcf_unpack(line, BCF_UN_INFO);
+for (i=0; i<line->n_info; i++)
+{
+if ( line->d.info[i].key==id ) return &line->d.info[i];
+}
+return NULL;
+}
+int bcf_get_info_values(const bcf_hdr_t *hdr, bcf1_t *line, const char *tag, void **dst, int *ndst, int type)
+{
+int i,j, tag_id = bcf_hdr_id2int(hdr, BCF_DT_ID, tag);
+if ( !bcf_hdr_idinfo_exists(hdr,BCF_HL_INFO,tag_id) ) return -1;    // no such INFO field in the header
+if ( bcf_hdr_id2type(hdr,BCF_HL_INFO,tag_id)!=type ) return -2;     // expected different type
+if ( !(line->unpacked & BCF_UN_INFO) ) bcf_unpack(line, BCF_UN_INFO);
+for (i=0; i<line->n_info; i++)
+if ( line->d.info[i].key==tag_id ) break;
+if ( i==line->n_info ) return ( type==BCF_HT_FLAG ) ? 0 : -3;       // the tag is not present in this record
+if ( type==BCF_HT_FLAG ) return 1;
+bcf_info_t *info = &line->d.info[i];
+if ( type==BCF_HT_STR )
+{
+if ( *ndst < info->len+1 )
+{
+*ndst = info->len + 1;
+*dst  = realloc(*dst, *ndst);
+}
+memcpy(*dst,info->vptr,info->len);
+((uint8_t*)*dst)[info->len] = 0;
+return info->len;
+}
+// Make sure the buffer is big enough
+int size1 = type==BCF_HT_INT ? sizeof(int32_t) : sizeof(float);
+if ( *ndst < info->len )
+{
+*ndst = info->len;
+*dst  = realloc(*dst, *ndst * size1);
+}
+if ( info->len == 1 )
+{
+if ( info->type==BCF_BT_FLOAT ) *((float*)*dst) = info->v1.f;
+else
+{
+#define BRANCH(type_t, missing) { \
+if ( info->v1.i==missing ) *((int32_t*)*dst) = bcf_int32_missing; \
+else *((int32_t*)*dst) = info->v1.i; \
+}
+switch (info->type)
+{
+case BCF_BT_INT8:  BRANCH(int8_t,  bcf_int8_missing ); break;
+case BCF_BT_INT16: BRANCH(int16_t, bcf_int16_missing); break;
+case BCF_BT_INT32: BRANCH(int32_t, bcf_int32_missing); break;
+}
+#undef BRANCH
+}
+return 1;
+}
+#define BRANCH(type_t, is_missing, is_vector_end, set_missing, out_type_t) { \
+out_type_t *tmp = (out_type_t *) *dst; \
+type_t *p = (type_t *) info->vptr; \
+for (j=0; j<info->len; j++) \
+{ \
+if ( is_vector_end ) return j; \
+if ( is_missing ) set_missing; \
+else *tmp = p[j]; \
+tmp++; \
+} \
+return j; \
+}
+switch (info->type) {
+case BCF_BT_INT8:  BRANCH(int8_t,  p[j]==bcf_int8_missing,  p[j]==bcf_int8_vector_end,  *tmp=bcf_int32_missing, int32_t); break;
+case BCF_BT_INT16: BRANCH(int16_t, p[j]==bcf_int16_missing, p[j]==bcf_int16_vector_end, *tmp=bcf_int32_missing, int32_t); break;
+case BCF_BT_INT32: BRANCH(int32_t, p[j]==bcf_int32_missing, p[j]==bcf_int32_vector_end, *tmp=bcf_int32_missing, int32_t); break;
+case BCF_BT_FLOAT: BRANCH(float,   bcf_float_is_missing(p[j]), bcf_float_is_vector_end(p[j]), bcf_float_set_missing(*tmp), float); break;
+default: fprintf(stderr,"TODO: %s:%d .. info->type=%d\n", __FILE__,__LINE__, info->type); exit(1);
+}
+#undef BRANCH
+return -4;  // this can never happen
+}
+int bcf_get_format_string(const bcf_hdr_t *hdr, bcf1_t *line, const char *tag, char ***dst, int *ndst)
+{
+int i,tag_id = bcf_hdr_id2int(hdr, BCF_DT_ID, tag);
+if ( !bcf_hdr_idinfo_exists(hdr,BCF_HL_FMT,tag_id) ) return -1;    // no such FORMAT field in the header
+if ( bcf_hdr_id2type(hdr,BCF_HL_FMT,tag_id)!=BCF_HT_STR ) return -2;     // expected different type
+if ( !(line->unpacked & BCF_UN_FMT) ) bcf_unpack(line, BCF_UN_FMT);
+for (i=0; i<line->n_fmt; i++)
+if ( line->d.fmt[i].id==tag_id ) break;
+if ( i==line->n_fmt ) return -3;                               // the tag is not present in this record
+bcf_fmt_t *fmt = &line->d.fmt[i];
+int nsmpl = bcf_hdr_nsamples(hdr);
+if ( !*dst )
+{
+*dst = (char**) malloc(sizeof(char*)*nsmpl);
+if ( !*dst ) return -4;     // could not alloc
+(*dst)[0] = NULL;
+}
+int n = (fmt->n+1)*nsmpl;
+if ( *ndst < n )
+{
+(*dst)[0] = realloc((*dst)[0], n);
+if ( !(*dst)[0] ) return -4;    // could not alloc
+*ndst = n;
+}
+for (i=0; i<nsmpl; i++)
+{
+uint8_t *src = fmt->p + i*fmt->n;
+uint8_t *tmp = (uint8_t*)(*dst)[0] + i*(fmt->n+1);
+memcpy(tmp,src,fmt->n);
+tmp[fmt->n] = 0;
+(*dst)[i] = (char*) tmp;
+}
+return n;
+}
+int bcf_get_format_values(const bcf_hdr_t *hdr, bcf1_t *line, const char *tag, void **dst, int *ndst, int type)
+{
+int i,j, tag_id = bcf_hdr_id2int(hdr, BCF_DT_ID, tag);
+if ( !bcf_hdr_idinfo_exists(hdr,BCF_HL_FMT,tag_id) ) return -1;    // no such FORMAT field in the header
+if ( tag[0]=='G' && tag[1]=='T' && tag[2]==0 )
+{
+// Ugly: GT field is considered to be a string by the VCF header but BCF represents it as INT.
+if ( bcf_hdr_id2type(hdr,BCF_HL_FMT,tag_id)!=BCF_HT_STR ) return -2;
+}
+else if ( bcf_hdr_id2type(hdr,BCF_HL_FMT,tag_id)!=type ) return -2;     // expected different type
+if ( !(line->unpacked & BCF_UN_FMT) ) bcf_unpack(line, BCF_UN_FMT);
+for (i=0; i<line->n_fmt; i++)
+if ( line->d.fmt[i].id==tag_id ) break;
+if ( i==line->n_fmt ) return -3;                               // the tag is not present in this record
+bcf_fmt_t *fmt = &line->d.fmt[i];
+if ( type==BCF_HT_STR )
+{
+int n = fmt->n*bcf_hdr_nsamples(hdr);
+if ( *ndst < n )
+{
+*dst  = realloc(*dst, n);
+if ( !*dst ) return -4;     // could not alloc
+*ndst = n;
+}
+memcpy(*dst,fmt->p,n);
+return n;
+}
+// Make sure the buffer is big enough
+int nsmpl = bcf_hdr_nsamples(hdr);
+int size1 = type==BCF_HT_INT ? sizeof(int32_t) : sizeof(float);
+if ( *ndst < fmt->n*nsmpl )
+{
+*ndst = fmt->n*nsmpl;
+*dst  = realloc(*dst, *ndst*size1);
+if ( !dst ) return -4;     // could not alloc
+}
+#define BRANCH(type_t, is_missing, is_vector_end, set_missing, set_vector_end, out_type_t) { \
+out_type_t *tmp = (out_type_t *) *dst; \
+type_t *p = (type_t*) fmt->p; \
+for (i=0; i<nsmpl; i++) \
+{ \
+for (j=0; j<fmt->n; j++) \
+{ \
+if ( is_missing ) set_missing; \
+else if ( is_vector_end ) { set_vector_end; break; } \
+else *tmp = p[j]; \
+tmp++; \
+} \
+for (; j<fmt->n; j++) { set_vector_end; tmp++; } \
+p = (type_t *)((char *)p + fmt->size); \
+} \
+}
+switch (fmt->type) {
+case BCF_BT_INT8:  BRANCH(int8_t,  p[j]==bcf_int8_missing,  p[j]==bcf_int8_vector_end,  *tmp=bcf_int32_missing, *tmp=bcf_int32_vector_end, int32_t); break;
+case BCF_BT_INT16: BRANCH(int16_t, p[j]==bcf_int16_missing, p[j]==bcf_int16_vector_end, *tmp=bcf_int32_missing, *tmp=bcf_int32_vector_end, int32_t); break;
+case BCF_BT_INT32: BRANCH(int32_t, p[j]==bcf_int32_missing, p[j]==bcf_int32_vector_end, *tmp=bcf_int32_missing, *tmp=bcf_int32_vector_end, int32_t); break;
+case BCF_BT_FLOAT: BRANCH(float,   bcf_float_is_missing(p[j]), bcf_float_is_vector_end(p[j]), bcf_float_set_missing(*tmp), bcf_float_set_vector_end(*tmp), float); break;
+default: fprintf(stderr,"TODO: %s:%d .. fmt->type=%d\n", __FILE__,__LINE__, fmt->type); exit(1);
+}
+#undef BRANCH
+return nsmpl*fmt->n;
+}

Mercurial > repos > youngkim > ezbamqc

comparison ezBAMQC/src/htslib/vcf.c @ 0:dfa3745e5fd8