#define HTS_BUILDING_LIBRARY // Enables HTSLIB_EXPORT, see htslib/hts_defs.h

#include <config.h>

#include <stdlib.h>

#include <string.h>

#include <stdio.h>

#include <assert.h>

#include <errno.h>

#include "htslib/tbx.h"

#include "htslib/bgzf.h"

#include "htslib/hts_alloc.h"

#include "htslib/hts_endian.h"

#include "hts_internal.h"

#include "bgzf_internal.h"

#include "htslib/khash.h"

KHASH_DECLARE(s2i, kh_cstr_t, int64_t)

const tbx_conf_t tbx_conf_gff = { 0, 1, 4, 5, '#', 0 };

const tbx_conf_t tbx_conf_bed = { TBX_UCSC, 1, 2, 3, '#', 0 };

const tbx_conf_t tbx_conf_psltbl = { TBX_UCSC, 15, 17, 18, '#', 0 };

const tbx_conf_t tbx_conf_sam = { TBX_SAM, 3, 4, 0, '@', 0 };

const tbx_conf_t tbx_conf_vcf = { TBX_VCF, 1, 2, 0, '#', 0 };

const tbx_conf_t tbx_conf_gaf = { TBX_GAF, 1, 6, 0, '#', 0 };

typedef struct {

int64_t beg, end;

char *ss, *se;

int tid;

} tbx_intv_t;

static inline int get_tid(tbx_t *tbx, const char *ss, int is_add)

{

khint_t k;

khash_t(s2i) *d;

if ((tbx->conf.preset&0xffff) == TBX_GAF) return(0);

if (tbx->dict == 0) tbx->dict = kh_init(s2i);

if (!tbx->dict) return -1; // Out of memory

d = (khash_t(s2i)*)tbx->dict;

if (is_add) {

int absent;

k = kh_put(s2i, d, ss, &absent);

if (absent < 0) {

return -1; // Out of memory

} else if (absent) {

char *ss_dup = strdup(ss);

if (ss_dup) {

kh_key(d, k) = ss_dup;

kh_val(d, k) = kh_size(d) - 1;

} else {

kh_del(s2i, d, k);

return -1; // Out of memory

}

}

} else k = kh_get(s2i, d, ss);

return k == kh_end(d)? -1 : kh_val(d, k);

}

int tbx_name2id(tbx_t *tbx, const char *ss)

{

return get_tid(tbx, ss, 0);

}

int tbx_parse1(const tbx_conf_t *conf, size_t len, char *line, tbx_intv_t *intv)

{

size_t i, b = 0;

int id = 1, getlen = 0, alcnt = 0, use_svlen = 0, lenpos = -1;

// The last column that can affect the result. Scanning stops once it

// has been seen, so trailing columns (e.g. the samples of a large VCF)

// are not examined at all.

int last_col = conf->sc > conf->bc ? conf->sc : conf->bc;

char *s, *t;

uint8_t svlenals[8192];

int64_t reflen = 0, svlen = 0, fmtlen = 0, tmp = 0;

if (conf->ec > last_col) last_col = conf->ec;

if ((conf->preset&0xffff) == TBX_VCF && last_col < 8)

last_col = 8; // INFO

else if ((conf->preset&0xffff) == TBX_SAM && last_col < 6)

last_col = 6; // CIGAR

intv->ss = intv->se = 0; intv->beg = intv->end = -1;

for (i = 0; i <= len; ++i) {

if (line[i] == '\t' || line[i] == 0) {

if (id == conf->sc) {

intv->ss = line + b; intv->se = line + i;

} else if (id == conf->bc) {

// here ->beg is 0-based.

if ((conf->preset&0xffff) == TBX_GAF){

// if gaf find the smallest and largest node id

char *t;

int64_t nodeid = -1;

for (s = line + b + 1; s < line + i;) {

nodeid = strtoll(s, &t, 0);

if(intv->beg == -1){

intv->beg = intv->end = nodeid;

} else {

if(nodeid < intv->beg){

intv->beg = nodeid;

}

if(nodeid > intv->end){

intv->end = nodeid;

}

}

s = t + 1;

}

} else {

intv->beg = strtoll(line + b, &s, 0);

if (conf->bc <= conf->ec) // don't overwrite an already set end point

intv->end = intv->beg;

if ( s==line+b ) return -1; // expected int

if (!(conf->preset&TBX_UCSC))

--intv->beg;

else if (conf->bc <= conf->ec)

++intv->end;

if (intv->beg < 0) {

hts_log_warning("Coordinate <= 0 detected. "

"Did you forget to use the -0 option?");

intv->beg = 0;

}

if (intv->end < 1) intv->end = 1;

}

} else {

if ((conf->preset&0xffff) == TBX_GENERIC) {

if (id == conf->ec)

{

intv->end = strtoll(line + b, &s, 0);

if ( s==line+b ) return -1; // expected int

}

} else if ((conf->preset&0xffff) == TBX_SAM) {

if (id == 6) { // CIGAR

int l = 0;

char *t;

for (s = line + b; s < line + i;) {

long x = strtol(s, &t, 10);

char op = toupper_c(*t);

if (op == 'M' || op == 'D' || op == 'N') l += x;

s = t + 1;

}

if (l == 0) l = 1;

intv->end = intv->beg + l;

}

} else if ((conf->preset&0xffff) == TBX_VCF) {

if (id == 4) { //ref allele

if (b < i) intv->end = intv->beg + (i - b);

++alcnt;

reflen = i - b;

} if (id == 5) { //alt allele

int lastbyte = 0, c = line[i];

svlenals[lastbyte] = 0;

line[i] = 0;

s = line + b;

do {

t = strchr(s, ',');

if (alcnt >> 3 != lastbyte) { //initialize insals

lastbyte = alcnt >> 3;

svlenals[lastbyte] = 0;

}

++alcnt;

if (t) {

*t = 0;

}

if (svlen_on_ref_for_vcf_alt(s, -1)) {

// Need to check SVLEN for this ALT

svlenals[lastbyte] |= 1 << ((alcnt - 1) & 7);

use_svlen = 1;

} else if (!strcmp("<*>", s) ||

!strcmp("<NON_REF>", s)) { //note gvcf

getlen = 1;

}

if (t) {

*t = ',';

s = t + 1;

}

} while (t && alcnt < 65536); //max allcnt is 65535

line[i] = c;

} else if (id == 8) { //INFO, look for "END=" / "SVLEN"

int c = line[i], d = 1;

line[i] = 0;

s = strstr(line + b, "END=");

if (s == line + b) s += 4;

else if (s) {

s = strstr(line + b, ";END=");

if (s) s += 5;

}

if (s && *s != '.') {

long long end = strtoll(s, &s, 0);

if (end <= intv->beg) {

static int reported = 0;

if (!reported) {

int l = intv->ss ? (int) (intv->se - intv->ss) : 0;

hts_log_warning("VCF INFO/END=%lld is smaller than POS at %.*s:%"PRIhts_pos"\n"

"This tag will be ignored. "

"Note: only one invalid END tag will be reported.",

end, l >= 0 ? l : 0,

intv->ss ? intv->ss : "",

intv->beg);

reported = 1;

}

} else {

intv->end = end;

}

}

// SVLEN only matters when a <DEL>-like symbolic

// allele refers to it, as in vcf.c:get_rlen

s = NULL;

if (use_svlen) {

s = strstr(line + b, "SVLEN=");

if (s == line + b) s += 6; //at start of info

else if (s) { //not at the start

s = strstr(line + b, ";SVLEN=");

if (s) s += 7;

}

}

while (s && d < alcnt) {

t = strchr(s, ',');

if (svlenals[d >> 3] & (1 << (d & 7))) {

// <DEL> symbolic allele

tmp = atoll(s);

tmp = tmp < 0 ? llabs(tmp) : tmp;

} else {

tmp = 1;

}

svlen = svlen < tmp ? tmp : svlen;

s = t ? t + 1 : NULL;

++d;

}

line[i] = c;

} else if (getlen && id == 9 ) { //FORMAT

int c = line[i], pos = -1;

line[i] = 0;

s = line + b;

while (s) {

++pos;

if (!(t = strchr(s, ':'))) { //no further fields

if (!strcmp(s, "LEN")) {

lenpos = pos;

}

break; //not present at all!

} else {

*t = '\0';

if (!strcmp(s, "LEN")) {

lenpos = pos;

*t = ':';

break;

}

*t = ':';

s = t + 1; //check next one

}

}

line[i] = c;

if (lenpos == -1) { //not present

break;

}

} else if (id > 9 && getlen && lenpos != -1) {

//get LEN from sample

int c = line[i], d = 0;

line[i] = 0; tmp = 0;

s = line + b;

for (d = 0; d <= lenpos; ++d) {

if (d == lenpos) {

tmp = atoll(s);

break;

}

if ((t = strchr(s, ':'))) {

s = t + 1;

} else {

break; //not in sync with fmt def!

}

}

fmtlen = fmtlen < tmp ? tmp : fmtlen;

line[i] = c;

}

}

}

b = i + 1; //beginning if current field

++id;

// Nothing past last_col can change the result, unless a gVCF

// <*>/<NON_REF> allele means the FORMAT LEN field must be checked

if (id > last_col

&& !((conf->preset&0xffff) == TBX_VCF && getlen)) break;

}

}

if ((conf->preset&0xffff) == TBX_VCF) {

tmp = reflen < svlen ?

svlen < fmtlen ? fmtlen : svlen :

reflen < fmtlen ? fmtlen : reflen ;

tmp += intv->beg;

intv->end = intv->end < tmp ? tmp : intv->end;

//NOTE: 'end' calculation be in sync with end/rlen in vcf.c:get_rlen

}

if (intv->ss == 0 || intv->se == 0 || intv->beg < 0 || intv->end < 0) return -1;

return 0;

}

static inline int get_intv(tbx_t *tbx, kstring_t *str, tbx_intv_t *intv, int is_add)

{

if (tbx_parse1(&tbx->conf, str->l, str->s, intv) == 0) {

int c = *intv->se;

*intv->se = '\0';

if ((tbx->conf.preset&0xffff) == TBX_GAF){

intv->tid = 0;

} else {

intv->tid = get_tid(tbx, intv->ss, is_add);

}

*intv->se = c;

if (intv->tid < 0) return -2; // get_tid out of memory

return (intv->beg >= 0 && intv->end >= 0)? 0 : -1;

} else {

char *type = NULL;

switch (tbx->conf.preset&0xffff)

{

case TBX_SAM: type = "TBX_SAM"; break;

case TBX_VCF: type = "TBX_VCF"; break;

case TBX_GAF: type = "TBX_GAF"; break;

case TBX_UCSC: type = "TBX_UCSC"; break;

default: type = "TBX_GENERIC"; break;

}

if (hts_is_utf16_text(str))

hts_log_error("Failed to parse %s: offending line appears to be encoded as UTF-16", type);

else

hts_log_error("Failed to parse %s: was wrong -p [type] used?\nThe offending line was: \"%s\"",

type, str->s);

return -1;

}

}

int tbx_readrec(BGZF *fp, void *tbxv, void *sv, int *tid, hts_pos_t *beg, hts_pos_t *end)

{

tbx_t *tbx = (tbx_t *) tbxv;

kstring_t *s = (kstring_t *) sv;

int ret;

// Get a line until either EOF or a non-meta character

do {

ret = bgzf_getline(fp, '\n', s);

} while (ret >= 0 && s->l && *s->s == tbx->conf.meta_char);

// Parse line

if (ret >= 0) {

tbx_intv_t intv;

if (get_intv(tbx, s, &intv, 0) < 0)

return -2;

*tid = intv.tid; *beg = intv.beg; *end = intv.end;

}

return ret;

}

typedef struct tbx_wrapper {

void *result_ptr;

tbx_t *tbx;

} tbx_wrapper;

int tbx_itr_next1(htsFile *htsfp, tbx_t *tbx, hts_itr_t *iter, void *r)

{

if (!htsfp->is_bgzf) {

hts_log_error("Only bgzf compressed files can be used with iterators");

errno = EINVAL;

return -2;

}

if (iter->multi) {

tbx_wrapper tmp = { r, tbx };

return hts_itr_multi_next(htsfp, iter, &tmp);

} else {

return hts_itr_next(htsfp->fp.bgzf, iter, r, tbx);

}

}

static int tbx_multi_readrec(BGZF *fp, void *fpv, void *r,

int *tid, hts_pos_t *beg, hts_pos_t *end)

{

tbx_wrapper *tmp = (tbx_wrapper *) r;

return tbx_readrec(fp, tmp->tbx, tmp->result_ptr, tid, beg, end);

}

static int tbx_set_meta(tbx_t *tbx)

{

int i, l = 0, l_nm;

uint32_t x[7];

char **name;

uint8_t *meta;

khint_t k;

khash_t(s2i) *d = (khash_t(s2i)*)tbx->dict;

memcpy(x, &tbx->conf, 24);

name = hts_malloc_p(sizeof(char*), kh_size(d));

if (!name) return -1;

for (k = kh_begin(d), l = 0; k != kh_end(d); ++k) {

if (!kh_exist(d, k)) continue;

name[kh_val(d, k)] = (char*)kh_key(d, k);

l += strlen(kh_key(d, k)) + 1; // +1 to include '\0'

}

l_nm = x[6] = l;

meta = hts_malloc_ps(sizeof(*meta), l_nm, 28);

if (!meta) { free(name); return -1; }

if (ed_is_big())

for (i = 0; i < 7; ++i)

x[i] = ed_swap_4(x[i]);

memcpy(meta, x, 28);

for (l = 28, i = 0; i < (int)kh_size(d); ++i) {

int x = strlen(name[i]) + 1;

memcpy(meta + l, name[i], x);

l += x;

}

free(name);

hts_idx_set_meta(tbx->idx, l, meta, 0);

return 0;

}

static void adjust_max_ref_len_vcf(const char *str, int64_t *max_ref_len)

{

const char *ptr;

int64_t len;

if (strncmp(str, "##contig", 8) != 0) return;

ptr = strstr(str + 8, "length");

if (!ptr) return;

for (ptr += 6; *ptr == ' ' || *ptr == '='; ptr++) {}

len = strtoll(ptr, NULL, 10);

if (*max_ref_len < len) *max_ref_len = len;

}

static void adjust_max_ref_len_sam(const char *str, int64_t *max_ref_len)

{

const char *ptr;

int64_t len;

if (strncmp(str, "@SQ", 3) != 0) return;

ptr = strstr(str + 3, "\tLN:");

if (!ptr) return;

ptr += 4;

len = strtoll(ptr, NULL, 10);

if (*max_ref_len < len) *max_ref_len = len;

}

tbx_t *tbx_index(BGZF *fp, int min_shift, const tbx_conf_t *conf)

{

tbx_t *tbx;

kstring_t str;

int ret, first = 0, n_lvls, fmt;

int64_t lineno = 0;

uint64_t last_off = 0;

tbx_intv_t intv;

int64_t max_ref_len = 0;

str.s = 0; str.l = str.m = 0;

tbx = (tbx_t*)calloc(1, sizeof(tbx_t));

if (!tbx) return NULL;

tbx->conf = *conf;

if (min_shift > 0) n_lvls = (TBX_MAX_SHIFT - min_shift + 2) / 3, fmt = HTS_FMT_CSI;

else min_shift = 14, n_lvls = 5, fmt = HTS_FMT_TBI;

while ((ret = bgzf_getline(fp, '\n', &str)) >= 0) {

++lineno;

if (str.s[0] == tbx->conf.meta_char && fmt == HTS_FMT_CSI) {

switch (tbx->conf.preset) {

case TBX_SAM:

adjust_max_ref_len_sam(str.s, &max_ref_len); break;

case TBX_VCF:

adjust_max_ref_len_vcf(str.s, &max_ref_len); break;

default:

break;

}

}

if (lineno <= tbx->conf.line_skip || str.s[0] == tbx->conf.meta_char) {

last_off = bgzf_tell(fp);

continue;

}

if (first == 0) {

if (fmt == HTS_FMT_CSI) {

if (max_ref_len) {

hts_adjust_csi_settings(max_ref_len, &min_shift, &n_lvls);

} else {

// This will give a maximum reference length of at

// least 100Gbases for min_shift >= 10, and the

// maximum possible for min_shift < 10.

const int max_n_lvls = 9; // To prevent bin number overflow

n_lvls = (min_shift < 10

? max_n_lvls

: (min_shift < 25

? max_n_lvls - (min_shift - 10) / 3

: 4));

}

}

tbx->idx = hts_idx_init(0, fmt, last_off, min_shift, n_lvls);

if (!tbx->idx) goto fail;

first = 1;

}

ret = get_intv(tbx, &str, &intv, 1);

if (ret < 0) goto fail; // Out of memory or unparsable lines

if (hts_idx_push(tbx->idx, intv.tid, intv.beg, intv.end,

bgzf_tell(fp), 1) < 0) {

goto fail;

}

}

if (ret < -1) goto fail;

if ( !tbx->idx ) tbx->idx = hts_idx_init(0, fmt, last_off, min_shift, n_lvls); // empty file

if (!tbx->idx) goto fail;

if ( !tbx->dict ) tbx->dict = kh_init(s2i);

if (!tbx->dict) goto fail;

if (hts_idx_finish(tbx->idx, bgzf_tell(fp)) != 0) goto fail;

if (tbx_set_meta(tbx) != 0) goto fail;

free(str.s);

return tbx;

fail:

free(str.s);

tbx_destroy(tbx);

return NULL;

}

void tbx_destroy(tbx_t *tbx)

{

khash_t(s2i) *d = (khash_t(s2i)*)tbx->dict;

if (d != NULL)

{

khint_t k;

for (k = kh_begin(d); k != kh_end(d); ++k)

if (kh_exist(d, k)) free((char*)kh_key(d, k));

}

hts_idx_destroy(tbx->idx);

kh_destroy(s2i, d);

free(tbx);

}

int tbx_index_build3(const char *fn, const char *fnidx, int min_shift, int n_threads, const tbx_conf_t *conf)

{

tbx_t *tbx;

BGZF *fp;

int ret;

if ((fp = bgzf_open(fn, "r")) == 0) return -1;

if ( n_threads ) bgzf_mt(fp, n_threads, 256);

if ( bgzf_compression(fp) != bgzf ) { bgzf_close(fp); return -2; }

tbx = tbx_index(fp, min_shift, conf);

bgzf_close(fp);

if ( !tbx ) return -1;

ret = hts_idx_save_as(tbx->idx, fn, fnidx, min_shift > 0? HTS_FMT_CSI : HTS_FMT_TBI);

tbx_destroy(tbx);

return ret;

}

int tbx_index_build2(const char *fn, const char *fnidx, int min_shift, const tbx_conf_t *conf)

{

return tbx_index_build3(fn, fnidx, min_shift, 0, conf);

}

int tbx_index_build(const char *fn, int min_shift, const tbx_conf_t *conf)

{

return tbx_index_build3(fn, NULL, min_shift, 0, conf);

}

static tbx_t *index_load(const char *fn, const char *fnidx, int flags)

{

tbx_t *tbx;

uint8_t *meta;

char *nm, *p;

uint32_t l_meta, l_nm;

tbx = (tbx_t*)calloc(1, sizeof(tbx_t));

if (!tbx)

return NULL;

tbx->idx = hts_idx_load3(fn, fnidx, HTS_FMT_TBI, flags);

if ( !tbx->idx )

{

free(tbx);

return NULL;

}

meta = hts_idx_get_meta(tbx->idx, &l_meta);

if ( !meta || l_meta < 28) goto invalid;

tbx->conf.preset = le_to_i32(&meta[0]);

tbx->conf.sc = le_to_i32(&meta[4]);

tbx->conf.bc = le_to_i32(&meta[8]);

tbx->conf.ec = le_to_i32(&meta[12]);

tbx->conf.meta_char = le_to_i32(&meta[16]);

tbx->conf.line_skip = le_to_i32(&meta[20]);

l_nm = le_to_u32(&meta[24]);

if (l_nm > l_meta - 28) goto invalid;

p = nm = (char*)meta + 28;

// This assumes meta is NUL-terminated, so we can merrily strlen away.

// hts_idx_load_local() assures this for us by adding a NUL on the end

// of whatever it reads.

for (; p - nm < l_nm; p += strlen(p) + 1) {

if (get_tid(tbx, p, 1) < 0) {

hts_log_error("%s", strerror(errno));

goto fail;

}

}

return tbx;

invalid:

hts_log_error("Invalid index header for %s", fnidx ? fnidx : fn);

fail:

tbx_destroy(tbx);

return NULL;

}

tbx_t *tbx_index_load3(const char *fn, const char *fnidx, int flags)

{

return index_load(fn, fnidx, flags);

}

tbx_t *tbx_index_load2(const char *fn, const char *fnidx)

{

return index_load(fn, fnidx, 1);

}

tbx_t *tbx_index_load(const char *fn)

{

return index_load(fn, NULL, 1);

}

const char **tbx_seqnames(tbx_t *tbx, int *n)

{

khash_t(s2i) *d = (khash_t(s2i)*)tbx->dict;

if (d == NULL)

{

*n = 0;

return calloc(1, sizeof(char *));

}

int tid, m = kh_size(d);

const char **names = (const char**) calloc(m,sizeof(const char*));

khint_t k;

if (!names) {

*n = 0;

return NULL;

}

for (k=kh_begin(d); k<kh_end(d); k++)

{

if ( !kh_exist(d,k) ) continue;

tid = kh_val(d,k);

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;

}

static int tbx_name2id_wrapper(void *vhdr, const char *ref)

{

return tbx_name2id((tbx_t *) vhdr, ref);

}

hts_itr_t *tbx_itr_querys1(tbx_t *tbx, const char *region)

{

return hts_itr_querys(tbx->idx, region, tbx_name2id_wrapper, tbx,

hts_itr_query, tbx_readrec);

}

hts_itr_t *tbx_itr_regarray(tbx_t *tbx, char **regarray, unsigned int regcount)

{

hts_itr_t *itr = NULL;

hts_reglist_t *r_list = NULL;

int r_count = 0;

r_list = hts_reglist_create(regarray, regcount, &r_count, tbx,

tbx_name2id_wrapper);

if (!r_list)

return NULL;

itr = hts_itr_regions(tbx->idx, r_list, r_count, tbx_name2id_wrapper, tbx,

hts_itr_multi_bam, tbx_multi_readrec,

bgzf_pseek, bgzf_ptell);

return itr;

}

hts_itr_t *tbx_itr_regions(const tbx_t *tbx, hts_reglist_t *reglist,

unsigned int regcount)

{

return hts_itr_regions(tbx->idx, reglist, regcount, tbx_name2id_wrapper,

(void *) tbx, hts_itr_multi_bam, tbx_multi_readrec,

bgzf_pseek, bgzf_ptell);

}