Skip to content

Commit 6d2002b

Browse files
committed
Resructure of main vcf_parse function.
This can speed things up by 6-7% on basic single-sample files. The previous loop caused lots of branch prediction misses due to the counter 'i' being used to do 8 different parts of code depending on token number. Additionally it's got better error checking now as previously running out of tokens early just did a return 0 rather than complaining about missing columns.
1 parent bcb25ab commit 6d2002b

1 file changed

Lines changed: 142 additions & 71 deletions

File tree

vcf.c

Lines changed: 142 additions & 71 deletions
Original file line numberDiff line numberDiff line change
@@ -46,9 +46,29 @@ DEALINGS IN THE SOFTWARE. */
4646
#include "htslib/khash_str2int.h"
4747
#include "htslib/kstring.h"
4848
#include "htslib/sam.h"
49-
5049
#include "htslib/khash.h"
50+
51+
#if 0
52+
// This helps on Intel a bit, often 6-7% faster VCF parsing.
53+
// Conversely sometimes harms AMD Zen4 as ~9% slower.
54+
// Possibly related to IPC differences. However for now it's just a
55+
// curiousity we ignore and stick with the simpler code.
56+
//
57+
// Left here as a hint for future explorers.
58+
static inline int xstreq(const char *a, const char *b) {
59+
while (*a && *a == *b)
60+
a++, b++;
61+
return *a == *b;
62+
}
63+
64+
#define KHASH_MAP_INIT_XSTR(name, khval_t) \
65+
KHASH_INIT(name, kh_cstr_t, khval_t, 1, kh_str_hash_func, xstreq)
66+
67+
KHASH_MAP_INIT_XSTR(vdict, bcf_idinfo_t)
68+
#else
5169
KHASH_MAP_INIT_STR(vdict, bcf_idinfo_t)
70+
#endif
71+
5272
typedef khash_t(vdict) vdict_t;
5373

5474
KHASH_MAP_INIT_STR(hdict, bcf_hrec_t*)
@@ -1373,7 +1393,8 @@ bcf_hdr_t *bcf_hdr_init(const char *mode)
13731393
for (i = 0; i < 3; ++i) {
13741394
if ((h->dict[i] = kh_init(vdict)) == NULL) goto fail;
13751395
// Supersize the hash to make collisions very unlikely
1376-
if (kh_resize(vdict, h->dict[i], 32768) < 0) goto fail;
1396+
static int dsize[3] = {16384,16384,2048}; // info, contig, format
1397+
if (kh_resize(vdict, h->dict[i], dsize[i]) < 0) goto fail;
13771398
}
13781399

13791400
bcf_hdr_aux_t *aux = (bcf_hdr_aux_t*)calloc(1,sizeof(bcf_hdr_aux_t));
@@ -3328,6 +3349,10 @@ static int vcf_parse_info(kstring_t *str, const bcf_hdr_t *h, bcf1_t *v, char *p
33283349
// key=value;
33293350
*next_equals = 0;
33303351
from = val = next_equals+1;
3352+
// Prefetching d->keys[hash] helps here provided we avoid
3353+
// computing hash twice (needs API change), but not universally.
3354+
// It may be significant for other applications though so it's
3355+
// something to consider for the future.
33313356
} else {
33323357
// key;
33333358
val = NULL;
@@ -3491,7 +3516,7 @@ static int vcf_parse_info(kstring_t *str, const bcf_hdr_t *h, bcf1_t *v, char *p
34913516

34923517
int vcf_parse(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v)
34933518
{
3494-
int i = 0, ret = -2, overflow = 0;
3519+
int ret = -2, overflow = 0;
34953520
char *p, *q, *r, *t;
34963521
kstring_t *str;
34973522
khint_t k;
@@ -3528,79 +3553,125 @@ int vcf_parse(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v)
35283553
bcf_clear1(v);
35293554
str = &v->shared;
35303555
memset(&aux, 0, sizeof(ks_tokaux_t));
3531-
for (p = kstrtok(s->s, "\t", &aux), i = 0; p; p = kstrtok(0, 0, &aux), ++i) {
3532-
q = (char*)aux.p;
3533-
*q = 0;
3534-
if (i == 0) { // CHROM
3535-
vdict_t *d = (vdict_t*)h->dict[BCF_DT_CTG];
3536-
k = kh_get(vdict, d, p);
3537-
if (k == kh_end(d))
3538-
{
3539-
hts_log_warning("Contig '%s' is not defined in the header. (Quick workaround: index the file with tabix.)", p);
3540-
v->errcode = BCF_ERR_CTG_UNDEF;
3541-
if ((k = fix_chromosome(h, d, p)) == kh_end(d)) {
3542-
hts_log_error("Could not add dummy header for contig '%s'", p);
3543-
v->errcode |= BCF_ERR_CTG_INVALID;
3556+
3557+
// CHROM
3558+
if (!(p = kstrtok(s->s, "\t", &aux)))
3559+
goto err;
3560+
*(q = (char*)aux.p) = 0;
3561+
3562+
vdict_t *d = (vdict_t*)h->dict[BCF_DT_CTG];
3563+
k = kh_get(vdict, d, p);
3564+
if (k == kh_end(d)) {
3565+
hts_log_warning("Contig '%s' is not defined in the header. (Quick workaround: index the file with tabix.)", p);
3566+
v->errcode = BCF_ERR_CTG_UNDEF;
3567+
if ((k = fix_chromosome(h, d, p)) == kh_end(d)) {
3568+
hts_log_error("Could not add dummy header for contig '%s'", p);
3569+
v->errcode |= BCF_ERR_CTG_INVALID;
3570+
goto err;
3571+
}
3572+
}
3573+
v->rid = kh_val(d, k).id;
3574+
3575+
// POS
3576+
if (!(p = kstrtok(0, 0, &aux)))
3577+
goto err;
3578+
*(q = (char*)aux.p) = 0;
3579+
3580+
overflow = 0;
3581+
char *tmp = p;
3582+
v->pos = hts_str2uint(p, &p, 63, &overflow);
3583+
if (overflow) {
3584+
hts_log_error("Position value '%s' is too large", tmp);
3585+
goto err;
3586+
} else if ( *p ) {
3587+
hts_log_error("Could not parse the position '%s'", tmp);
3588+
goto err;
3589+
} else {
3590+
v->pos -= 1;
3591+
}
3592+
if (v->pos >= INT32_MAX)
3593+
v->unpacked |= BCF_IS_64BIT;
3594+
3595+
// ID
3596+
if (!(p = kstrtok(0, 0, &aux)))
3597+
goto err;
3598+
*(q = (char*)aux.p) = 0;
3599+
3600+
if (NOT_DOT(p)) bcf_enc_vchar(str, q - p, p);
3601+
else bcf_enc_size(str, 0, BCF_BT_CHAR);
3602+
3603+
// REF
3604+
if (!(p = kstrtok(0, 0, &aux)))
3605+
goto err;
3606+
*(q = (char*)aux.p) = 0;
3607+
3608+
bcf_enc_vchar(str, q - p, p);
3609+
v->n_allele = 1, v->rlen = q - p;
3610+
3611+
// ALT
3612+
if (!(p = kstrtok(0, 0, &aux)))
3613+
goto err;
3614+
*(q = (char*)aux.p) = 0;
3615+
3616+
if (NOT_DOT(p)) {
3617+
for (r = t = p;; ++r) {
3618+
if (*r == ',' || *r == 0) {
3619+
if (v->n_allele == UINT16_MAX) {
3620+
hts_log_error("Too many ALT alleles at %s:%"PRIhts_pos,
3621+
bcf_seqname_safe(h,v), v->pos+1);
3622+
v->errcode |= BCF_ERR_LIMITS;
35443623
goto err;
35453624
}
3625+
bcf_enc_vchar(str, r - t, t);
3626+
t = r + 1;
3627+
++v->n_allele;
35463628
}
3547-
v->rid = kh_val(d, k).id;
3548-
} else if (i == 1) { // POS
3549-
overflow = 0;
3550-
char *tmp = p;
3551-
v->pos = hts_str2uint(p, &p, 63, &overflow);
3552-
if (overflow) {
3553-
hts_log_error("Position value '%s' is too large", tmp);
3554-
goto err;
3555-
} else if ( *p ) {
3556-
hts_log_error("Could not parse the position '%s'", tmp);
3557-
goto err;
3558-
} else {
3559-
v->pos -= 1;
3560-
}
3561-
if (v->pos >= INT32_MAX)
3562-
v->unpacked |= BCF_IS_64BIT;
3563-
} else if (i == 2) { // ID
3564-
if (NOT_DOT(p)) bcf_enc_vchar(str, q - p, p);
3565-
else bcf_enc_size(str, 0, BCF_BT_CHAR);
3566-
} else if (i == 3) { // REF
3567-
bcf_enc_vchar(str, q - p, p);
3568-
v->n_allele = 1, v->rlen = q - p;
3569-
} else if (i == 4) { // ALT
3570-
if (NOT_DOT(p)) {
3571-
for (r = t = p;; ++r) {
3572-
if (*r == ',' || *r == 0) {
3573-
if (v->n_allele == UINT16_MAX) {
3574-
hts_log_error("Too many ALT alleles at %s:%"PRIhts_pos,
3575-
bcf_seqname_safe(h,v), v->pos+1);
3576-
v->errcode |= BCF_ERR_LIMITS;
3577-
goto err;
3578-
}
3579-
bcf_enc_vchar(str, r - t, t);
3580-
t = r + 1;
3581-
++v->n_allele;
3582-
}
3583-
if (r == q) break;
3584-
}
3585-
}
3586-
} else if (i == 5) { // QUAL
3587-
if (NOT_DOT(p)) v->qual = atof(p);
3588-
else bcf_float_set_missing(v->qual);
3589-
if ( v->max_unpack && !(v->max_unpack>>1) ) goto end; // BCF_UN_STR
3590-
} else if (i == 6) { // FILTER
3591-
if (NOT_DOT(p)) {
3592-
if (vcf_parse_filter(str, h, v, p, q)) goto err;
3593-
} else bcf_enc_vint(str, 0, 0, -1);
3594-
if ( v->max_unpack && !(v->max_unpack>>2) ) goto end; // BCF_UN_FLT
3595-
} else if (i == 7) { // INFO
3596-
if (NOT_DOT(p)) {
3597-
if (vcf_parse_info(str, h, v, p, q)) goto err;
3598-
}
3599-
if ( v->max_unpack && !(v->max_unpack>>3) ) goto end;
3600-
} else if (i == 8) {// FORMAT
3601-
return vcf_parse_format(s, h, v, p, q) == 0 ? 0 : -2;
3629+
if (r == q) break;
3630+
}
3631+
}
3632+
3633+
// QUAL
3634+
if (!(p = kstrtok(0, 0, &aux)))
3635+
goto err;
3636+
*(q = (char*)aux.p) = 0;
3637+
3638+
if (NOT_DOT(p)) v->qual = atof(p);
3639+
else bcf_float_set_missing(v->qual);
3640+
if ( v->max_unpack && !(v->max_unpack>>1) ) goto end; // BCF_UN_STR
3641+
3642+
// FILTER
3643+
if (!(p = kstrtok(0, 0, &aux)))
3644+
goto err;
3645+
*(q = (char*)aux.p) = 0;
3646+
3647+
if (NOT_DOT(p)) {
3648+
if (vcf_parse_filter(str, h, v, p, q)) {
3649+
goto err;
3650+
}
3651+
} else bcf_enc_vint(str, 0, 0, -1);
3652+
if ( v->max_unpack && !(v->max_unpack>>2) ) goto end; // BCF_UN_FLT
3653+
3654+
// INFO
3655+
if (!(p = kstrtok(0, 0, &aux)))
3656+
goto err;
3657+
*(q = (char*)aux.p) = 0;
3658+
3659+
if (NOT_DOT(p)) {
3660+
if (vcf_parse_info(str, h, v, p, q)) {
3661+
goto err;
36023662
}
36033663
}
3664+
if ( v->max_unpack && !(v->max_unpack>>3) ) goto end;
3665+
3666+
// FORMAT; optional
3667+
p = kstrtok(0, 0, &aux);
3668+
if (p) {
3669+
*(q = (char*)aux.p) = 0;
3670+
3671+
return vcf_parse_format(s, h, v, p, q) == 0 ? 0 : -2;
3672+
} else {
3673+
return 0;
3674+
}
36043675

36053676
end:
36063677
ret = 0;

0 commit comments

Comments
 (0)