-
Notifications
You must be signed in to change notification settings - Fork 10
/
vcf.h
339 lines (283 loc) · 9.27 KB
/
vcf.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
#ifndef BCF_H
#define BCF_H
#include <stdint.h>
#include <limits.h>
#include <assert.h>
#include "bgzf.h"
#include "hts.h"
/*****************
* Header struct *
*****************/
#define BCF_HL_FLT 0 // header line
#define BCF_HL_INFO 1
#define BCF_HL_FMT 2
#define BCF_HL_CTG 3
#define BCF_HT_FLAG 0 // header type
#define BCF_HT_INT 1
#define BCF_HT_REAL 2
#define BCF_HT_STR 3
#define BCF_VL_FIXED 0 // variable length
#define BCF_VL_VAR 1
#define BCF_VL_A 2
#define BCF_VL_G 3
#define BCF_VL_R 4
/* === Dictionary ===
The header keeps three dictonaries. The first keeps IDs in the
"FILTER/INFO/FORMAT" lines, the second keeps the sequence names and lengths
in the "contig" lines and the last keeps the sample names. bcf_hdr_t::dict[]
is the actual hash table, which is opaque to the end users. In the hash
table, the key is the ID or sample name as a C string and the value is a
bcf_idinfo_t struct. bcf_hdr_t::id[] points to key-value pairs in the hash
table in the order that they appear in the VCF header. bcf_hdr_t::n[] is the
size of the hash table or, equivalently, the length of the id[] arrays.
*/
#define BCF_DT_ID 0 // dictionary type
#define BCF_DT_CTG 1
#define BCF_DT_SAMPLE 2
typedef struct {
uint32_t info[3]; // for each number => Number:20, var:4, Type:4, ColType:4
int id;
} bcf_idinfo_t;
typedef struct {
const char *key;
const bcf_idinfo_t *val;
} bcf_idpair_t;
typedef struct {
int32_t l_text, m_text, n[3];
bcf_idpair_t *id[3];
void *dict[3]; // ID dictionary, contig dict and sample dict
char *text;
kstring_t mem;
} bcf_hdr_t;
extern uint8_t bcf_type_shift[];
/**************
* VCF record *
**************/
#define BCF_BT_NULL 0
#define BCF_BT_INT8 1
#define BCF_BT_INT16 2
#define BCF_BT_INT32 3
#define BCF_BT_FLOAT 5
#define BCF_BT_CHAR 7
typedef struct {
int id, n, type, size; // bcf_hdr_t::id[BCF_DT_ID][$id].key is the key in string; $size is the per-sample size
uint8_t *p; // point to the data array
} bcf_fmt_t;
typedef struct {
int key, type, len; // bcf_hdr_t::id[BCF_DT_ID][$key].key is the key in string; $len: the length of the vector
union {
int32_t i; // integer value
float f; // float value
} v1; // only set if $len==1; for easier access
uint8_t *vptr; // point to data array, excluding sized bytes
} bcf_info_t;
typedef struct {
int m_fmt, m_info, m_str, m_allele, m_flt; // allocated size (high-water mark); do not change
int n_flt; // # FILTER fields
char *id, **allele; // ID, REF and ALT; allele[0] is the REF; all null terminated
int *flt; // filter keys in the dictionary
bcf_info_t *info; // INFO
bcf_fmt_t *fmt; // FORMAT and individual sample
} bcf_dec_t;
typedef struct {
int32_t rid; // CHROM
int32_t pos; // POS
int32_t rlen; // length of REF
float qual; // QUAL
uint32_t n_info:16, n_allele:16;
uint32_t n_fmt:8, n_sample:24;
kstring_t shared, indiv;
bcf_dec_t d; // lazy evaluation: $d is not generated by bcf_read1(), but by explicitly calling bcf_unpack()
int unpacked; // allow calling bcf_unpack() repeatedly without redoing the work
uint8_t *unpack_ptr; // keep the place before which unpack is done
} bcf1_t;
#define bcf_int8_missing INT8_MIN
#define bcf_int16_missing INT16_MIN
#define bcf_int32_missing INT32_MIN
#define bcf_int8_end (INT8_MIN +1)
#define bcf_int16_end (INT16_MIN+1)
#define bcf_int32_end (INT32_MIN+1)
/*******
* API *
*******/
#ifdef __cplusplus
extern "C" {
#endif
/***************
*** BCF I/O ***
***************/
bcf_hdr_t *bcf_hdr_init(void);
int bcf_hdr_parse(bcf_hdr_t *h);
/**
* Read BCF header
*
* @param fp BGZF file pointer; file offset must be placed at the beginning
*
* @return BCF header struct
*/
bcf_hdr_t *bcf_hdr_read(BGZF *fp);
/**
* Write BCF header to BCF
*
* @param fp BGZF file pointer; file offset placed at the beginning
* @param h BCF header
*/
void bcf_hdr_write(BGZF *fp, const bcf_hdr_t *h);
/** Destroy a BCF header struct */
void bcf_hdr_destroy(bcf_hdr_t *h);
/** Initialize a bcf1_t object; equivalent to calloc(1, sizeof(bcf1_t)) */
bcf1_t *bcf_init1();
/** Deallocate a bcf1_t object */
void bcf_destroy1(bcf1_t *v);
void bcf_clear1(bcf1_t *v);
/**
* Read one BCF record
*
* @param fp BGZF file pointer
* @param v BCF record read from $fp
*
* @return 0 on success; -1 on normal file end; <-1 on error
*/
int bcf_read1(BGZF *fp, bcf1_t *v);
/**
* Write one BCF record
*
* @param fp BGZF file pointer
* @param v BCF record to write
*
* @return
*/
int bcf_write1(BGZF *fp, const bcf1_t *v);
/** Helper function for the bcf_iter_next() macro; ignore it */
int bcf_readrec(BGZF *fp, void *null, bcf1_t *v, int *tid, int *beg, int *end);
#define BCF_UN_STR 1 // up to ALT inclusive
#define BCF_UN_FLT 2 // up to FILTER
#define BCF_UN_INFO 4 // up to INFO
#define BCF_UN_SHR (BCF_UN_STR|BCF_UN_FLT|BCF_UN_INFO) // all shared information
#define BCF_UN_FMT 8 // unpack format and each sample
#define BCF_UN_IND BCF_UN_FMT // a synonymous of BCF_UN_FMT
#define BCF_UN_ALL (BCF_UN_SHR|BCF_UN_FMT) // everything
/**
* Unpack/decode a BCF record (fill the bcf1_t::d field
*/
int bcf_unpack(bcf1_t *b, int which); // to unpack everything, set $which to BCF_UN_ALL
int bcf_id2int(const bcf_hdr_t *h, int which, const char *id);
int bcf_name2id(const bcf_hdr_t *h, const char *id);
void bcf_fmt_array(kstring_t *s, int n, int type, void *data);
uint8_t *bcf_fmt_sized_array(kstring_t *s, uint8_t *ptr);
void bcf_enc_vchar(kstring_t *s, int l, const char *a);
void bcf_enc_vint(kstring_t *s, int n, const int32_t *a, int wsize);
void bcf_enc_vfloat(kstring_t *s, int n, float *a);
/*****************
*** BCF index ***
*****************/
#define bcf_itr_destroy(iter) hts_itr_destroy(iter)
#define bcf_itr_queryi(idx, tid, beg, end) hts_itr_query((idx), (tid), (beg), (end))
#define bcf_itr_querys(idx, hdr, s) hts_itr_querys((idx), (s), (hts_name2id_f)(bcf_name2id), (hdr))
#define bcf_itr_next(fp, itr, r) hts_itr_next((fp), (itr), (r), (hts_readrec_f)(bcf_readrec), 0)
#define bcf_index_load(fn) hts_idx_load(fn, HTS_FMT_CSI)
int bcf_index_build(const char *fn, int min_shift);
int bcf_seekn(BGZF *fp, const hts_idx_t *idx, int64_t n);
/***************
*** VCF I/O ***
***************/
typedef htsFile vcfFile;
#define vcf_open(fn, mode, fn_ref) hts_open((fn), (mode), (fn_ref)) // strchr(mode, 'b')!=0 for BCF; otherwise VCF
#define vcf_close(fp) hts_close(fp)
bcf_hdr_t *vcf_hdr_read(htsFile *fp);
void vcf_hdr_write(htsFile *fp, const bcf_hdr_t *h);
int vcf_parse1(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v);
int vcf_format1(const bcf_hdr_t *h, const bcf1_t *v, kstring_t *s);
int vcf_read1(htsFile *fp, const bcf_hdr_t *h, bcf1_t *v);
int vcf_write1(htsFile *fp, const bcf_hdr_t *h, const bcf1_t *v);
/*************************
*** VCF/BCF utilities ***
*************************/
int bcf_hdr_append(bcf_hdr_t *h, const char *line);
bcf_hdr_t *bcf_hdr_subset(const bcf_hdr_t *h0, int n, char *const* samples, int *imap);
int bcf_subset(const bcf_hdr_t *h, bcf1_t *v, int n, int *imap);
int bcf_is_snp(bcf1_t *v);
int bcf_is_filtered(bcf1_t *b);
void bcf_get_ref_alt1(const bcf1_t *b, int *l_ref, char **ref, int *l_alt, char **alt);
char *bcf_get_alt1(const bcf1_t *b, int *len);
void bcfcpy(bcf1_t *dst, const bcf1_t *src);
int bcfcmp(const bcf1_t *a, const bcf1_t *b);
int bcfcpy_min(bcf1_t *b, const bcf1_t *b0, const char *alt2);
int bcf_append_info_ints(const bcf_hdr_t *h, bcf1_t *b, const char *key, int n, const int32_t *val);
#ifdef __cplusplus
}
#endif
/*******************
* Typed value I/O *
*******************/
#include "kstring.h"
static inline void bcf_enc_size(kstring_t *s, int size, int type)
{
if (size >= 15) {
kputc(15<<4|type, s);
if (size >= 128) {
if (size >= 32768) {
int32_t x = size;
kputc(1<<4|BCF_BT_INT32, s);
kputsn((char*)&x, 4, s);
} else {
int16_t x = size;
kputc(1<<4|BCF_BT_INT16, s);
kputsn((char*)&x, 2, s);
}
} else {
kputc(1<<4|BCF_BT_INT8, s);
kputc(size, s);
}
} else kputc(size<<4|type, s);
}
static inline int bcf_enc_inttype(long x)
{
if (x <= INT8_MAX && x > INT8_MIN) return BCF_BT_INT8;
if (x <= INT16_MAX && x > INT16_MIN) return BCF_BT_INT16;
return BCF_BT_INT32;
}
static inline void bcf_enc_int1(kstring_t *s, int32_t x)
{
if (x == INT32_MIN) {
bcf_enc_size(s, 1, BCF_BT_INT8);
kputc(INT8_MIN, s);
} else if (x <= INT8_MAX && x > INT8_MIN) {
bcf_enc_size(s, 1, BCF_BT_INT8);
kputc(x, s);
} else if (x <= INT16_MAX && x > INT16_MIN) {
int16_t z = x;
bcf_enc_size(s, 1, BCF_BT_INT16);
kputsn((char*)&z, 2, s);
} else {
int32_t z = x;
bcf_enc_size(s, 1, BCF_BT_INT32);
kputsn((char*)&z, 4, s);
}
}
static inline int32_t bcf_dec_int1(const uint8_t *p, int type, uint8_t **q)
{
if (type == BCF_BT_INT8) {
*q = (uint8_t*)p + 1;
return *(int8_t*)p;
} else if (type == BCF_BT_INT16) {
*q = (uint8_t*)p + 2;
return *(int16_t*)p;
} else {
*q = (uint8_t*)p + 4;
return *(int32_t*)p;
}
}
static inline int32_t bcf_dec_typed_int1(const uint8_t *p, uint8_t **q)
{
return bcf_dec_int1(p + 1, *p&0xf, q);
}
static inline int32_t bcf_dec_size(const uint8_t *p, uint8_t **q, int *type)
{
*type = *p & 0xf;
if (*p>>4 != 15) {
*q = (uint8_t*)p + 1;
return *p>>4;
} else return bcf_dec_typed_int1(p + 1, q);
}
#endif