6 BCF specification
下面详细讲解 VCF 规范第 6 节 BCF specification。这一节讲的是 VCF 的二进制压缩版本 BCF2。你可以把它理解成:VCF 是给人看的文本格式,BCF2 是给机器高效读写的二进制格式。两者表达的信息基本等价,但 BCF2 更小、更快,并且可以被 tabix 索引。
6 BCF specification:为什么需要 BCF2
VCF 很灵活,能容纳多样本、多种字段,社区广泛使用。但它有两个大缺点:
- 文件大:文本格式占磁盘空间。100 个外显子组就是几 GB,上千个样本的外显子组 VCF 轻松上百 GB。
- 解析慢:文本解析非常耗时。
BCF2 的目标就是解决这两个问题:
- 它是 二进制、压缩的 VCF 等价物。
- 可以用 tabix 索引。
- 可以从磁盘或流中高效解码。
- 为了效率,BCF2 只支持 VCF 的一个子集:所有 INFO 和 GENOTYPE 字段必须在头部有完整类型定义。
- 例如,如果 INFO 里有
AC字段,头部必须声明AC是等位基因索引数组,类型为整数。 - 遇到未知或未完全指定的字段,转换 VCF 到 BCF2 时应该报错。
- 例如,如果 INFO 里有
6.1 Overall file organization:整体文件组织
一个 BCF2 文件由三部分组成:
- 强制头部:二进制头部 + VCF 文本头部。
- 一系列 BGZF 压缩块:每个块包含若干二进制 BCF2 记录。
- 一个空 BGZF 块:表示记录结束。
BGZF 是一种块压缩格式,允许 tabix 索引。BCF2 文件通常以 .bcf 为扩展名。
BCF2 头部遵循 VCF 规范,但有几个扩展/限制:
- 所有 contig 定义必须完整。记录不能引用头部没有的 contig。
- 所有 INFO 和 GENOTYPE 字段必须在头部完全类型化。
- 遇到未知或未完全指定的字段,转换时应报错。
6.2 Header:BCF2 头部
BCF2 头部包含以下字段:
| 字段 | 类型 | 说明 |
|---|---|---|
magic | char[3] | 字符 "BCF" |
major_version | uint8_t | 主版本号,BCF2 为 2 |
minor_version | uint8_t | 次版本号,BCF2.1 为 1 |
l_text | uint32_t | text 字段长度,包括结尾的 NULL 字符 |
text | char[l_text] | VCF 格式头部文本,以 NULL 结尾 |
magic和版本号用于快速判断文件是不是 BCF2.1。text字段包含标准 VCF 头部行,从##fileformat=VCFv4.2到#CHROM ...,最后以 NULL 字符结尾。- BCF2 支持字符串字典,用
##dictionary=SO,S1,...,SN定义,逗号分隔的有序字符串列表。
6.2.1 Dictionary of strings:字符串字典
BCF2 中大多数字符串值不直接存储,而是用整数引用字典中的字符串。这样可以大幅减少重复字符串的存储空间。
例如,VCF 头部:
##INFO=<ID=ASP,Number=0,Type=Flag,Description="X">
##INFO=<ID=RSPOS,Number=1,Type=Integer,Description="Y">
##INFO=<ID=dbSNPBuildID,Number=1,Type=Integer,Description="Z">
##contig=<ID=20,length=62435964,assembly=B36,md5=...,species="Homo sapiens">
#CHROM POS ID REF ALT QUAL FILTER INFO
20 10144 rs144773400 TA T . PASS ASP;RSPOS=10145,dbSNPBuildID=134
20 10228 rs143255646 TA T . PASS ASP;RSPOS=10229;dbSNPBuildID=134
在 BCF2 中,ASP、RSPOS、dbSNPBuildID 这些键名会被字典中的整数偏移代替。
6.3.1 Site encoding:位点编码
BCF2 记录的核心字段如下:
| 字段 | 类型 | 说明 |
|---|---|---|
l_shared | uint32_t | 从 CHROM 到 INFO 结束的数据长度 |
l_indiv | uint32_t | FORMAT 和个体基因型字段的数据长度 |
CHROM | int32_t | 强制 contig 字典中的偏移 |
POS | int32_t | 0 起始的最左坐标 |
rlen | int32_t | 记录投影到参考序列的长度。必须是 REF 等位基因长度,或符号等位基因声明的长度(尊重 END 属性) |
QUAL | float | 变异质量;缺失值为 0x7F800001 |
n_info | uint16_t | 该记录中 INFO 字段数量 |
n_allele | uint16_t | REF+ALT 等位基因数量 |
n_sample | uint24_t | 样本数量,3 字节小端。必须等于头部中的样本数 |
n_fmt | uint8_t | FORMAT 键数量 |
ID | typed string | 变异标识符;缺失值为 0x07 |
REF+ALT | 类型化字符串列表 | 第一个是 REF(强制),后面是 n_allele - 1 个 ALT |
FILTER | 类型化整数向量 | 整数偏移向量,每个 FILTER 值一个;. 编码为 MISSING |
INFO | 键/值对 | n_info 对类型化向量。第一个值是 INFO 键在字典中的偏移,第二个值是字段值 |
Genotype values | 见下文 | 见下文 |
6.3.2 Genotype encoding:基因型编码
VCF 中基因型按样本组织,每个样本一列。BCF2 中则 按字段组织:每个 FORMAT 字段有一个向量,向量中每个元素对应一个样本。
例如 VCF:
FORMAT NA00001 NA00002 NA00003
GT:GQ:DP 0/0:48:1 0/1:9:8 1/1:43:5
在 BCF2 中等价于:
GT = 0/0, 0/1, 1/1
GQ = 48, 9, 43
DP = 1, 8, 5
每个基因型字段编码为一个三元组:
| 字段 | 类型 | 说明 |
|---|---|---|
fnt_key | typed int | FORMAT 键在字典中的偏移 |
fnt_type | uint8_t+ | 每个个体值的类型字节,可能后跟类型化整数表示向量长度。对基因型值,它只出现一次,位于基因型字段值数组之前 |
fnt_values | 按类型 | 所有个体的值连接成向量。每个值类型相同。变长向量用缺失值填充;字符串存储为 char 向量 |
值始终隐含是 N 个值的向量,N 是样本数。类型字节表示 N 长度向量中每个值的类型。对于原子值,大小 = 1。如果类型字段表示值本身是向量(如 PL 字段),则外层向量的每个值都是向量。这种编码在基因型字段向量中每个值长度和类型相同时很高效。
字段顺序没有严格定义,但最好遵循输入 VCF/BCF2 的顺序。
如果没有样本记录,基因型块大小为 0。
6.3.3 Type encoding:类型编码
BCF2 中所有值都强类型。类型信息编码在值前面的前缀字节中。这个字节包含:
- 低层类型(如 int32、float);
- 元素数量。
编码方式:
| 位 | 含义 |
|---|---|
| 5,6,7,8 位 | 即将到来的类型的元素数量。原子值大小必须为 1。如果大小为 15,表示向量有 15 个或更多元素,后续 BCF2 字节流包含一个类型化整数表示真实大小。如果大小在 2–14 之间,则省略这个整数,后续流立即开始第一个值。大小为 0 表示 MISSING |
| 1,2,3,4 位 | 类型 |
低 4 位编码类型:
| 低 4 位 | 十六进制 | 原子类型 |
|---|---|---|
| 0 | 0x00 | MISSING 值,无显式类型 |
| 1 | 0x?1 | Integer [8 bit] |
| 2 | 0x?2 | Integer [16 bit] |
| 3 | 0x?3 | Integer [32 bit] |
| 5 | 0x?5 | Float [32 bit] |
| 7 | 0x?7 | Character,ASCII 编码 8 位 |
字符串类型 = 7,带适当长度。单字符是长度为 1 的字符串。类型码 4、6、8–15 保留。
- 整数可编码为 8、16、32 位,小端。缺失值:
0x80、0x8000、0x8000000。 - 浮点数单精度 32 位,IEEE-754-1985。NaN 表示缺失:
- qNaN:
0x7FC00000 - MISSING:
0x7F800001
- qNaN:
- 字符值用单字符字符串编码。
- Flag 值只能出现在 INFO 字段,编码为 1 元素 INT8,值 1。
字符串编码有两种:
- INFO、FORMAT、FILTER 键:用字典偏移。
- 字符串值(ID、REF、ALT、INFO、FORMAT 字段中的字符串):类型化 ASCII 字节数组,不以 NUL 结尾,长度由类型描述符给出。
例子:编码 "ACAC":
- 类型字节:字符串类型
0x07+ 内联大小 4 →0x40 | 0x07 = 0x47 - 后跟 ASCII:
0x41 0x43 0x41 0x43 - 最终:
0x47 0x41 0x43 0x41 0x43
例子:编码 "VariantCallFormatSampleText"(长度 27):
- 类型字节:
0xF0 | 0x07 = 0xF7(溢出大小) - 后跟类型化大小 27:
0x11 0x1B - 后跟 ASCII 字节。
缺失值 . 就是大小 0 的字符串:0x07。
字符串列表折叠:VCF 中有时有字符串列表(如氨基酸变化),BCF2 不支持字符串向量(字符向量已是字符串),所以折叠成逗号分隔字符串,编码为普通 BCF2 字符向量,读取时再展开。因为 VCF 字符串不能包含逗号,所以安全。折叠字符串开头加逗号,便于判断。
例子:X=[A,B,C,D] 编码为字符串 ",A,B,C,D",大小 8,类型字节 0x87,后跟 ASCII 0x2C 0x41 0x2C 0x42 0x2C 0x43 0x2C 0x44。
向量:
- 固定长度向量:类型字节指定长度,后跟值。
- 缺失值:可以部分缺失,每个值设为对应 MISSING;或者整个向量缺失,类型字节大小为 0,类型为 MISSING。
- 混合长度向量:取最大长度,每个样本向量左对齐,不足处填 MISSING。编码器/解码器自动添加/删除。
- 例子:样本 A
[1],样本 B[2,3],编码为[1, MISSING]和[2, 3]。 - 字节:
0x11 0x03(X 字典偏移 3),0x21(每个值是 2 元素 INT8),0x01 0x80(A),0x02 0x03(B)。
- 例子:样本 A
GT 字段编码:
- 类型化整数向量(8、16 或 32 位),元素数量等于该位点所有样本最大倍性。
- 每个整数:
(allele + 1) << 1 | phasedallele是等位基因编号;如果是.,设为 -1,高位全 0。phased0 或 1。
- 向量按最大倍性填充 MISSING。
- 例子:
0/1:(0+1)<<1 | 0 = 0x02,(1+1)<<1 | 0 = 0x04→0x02 0x04- 三个样本
0/1, 1/1, 0/0:0x02 04 04 04 02 02 0|1:(1+1)<<1 | 1 = 0x05,前面标准首字节0x02→0x02 05./.:0x00 00- 单倍体
0:0x02 - 单倍体
1:0x04 - 四倍体
0/1/2:0x02 04 06 - 四倍体
0/1|2:0x02 04 07 - 单倍体
0和二倍体0/1:0x02 80 02 04(单倍体补 MISSING)
6.4 Encoding a VCF record example:完整编码示例
用一个具体 VCF 记录:
CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001 NA00002 NA00003
chr1 101 rs123 A C 30.1 PASS HM3;AC=3;AN=6;AA=C GT:GQ:DP:AD:PL 0/0:10:32:32,0:0,10,100 0/1:10:48:32,16:10,0,100 1/1:10:64:0,64:100,10,0
假设 chr1 是 contig 列表中第二个(偏移 1)。
6.4.1 编码 CHROM 和 POS
CHROM偏移:0x01000000(32 位小端)POS:0 起始,所以 101 →0x64000000rlen:REF 长度 1 →0x01000000QUAL:30.1 的 32 位浮点 →0x41 0xF0 0xCC 0xCDn_info:4 →0x0400n_allele:2 →0x0200n_sample:3 →0x030000(24 位小端)n_fmt:5 →0x05
6.4.7 编码 Genotypes
假设 FORMAT 键偏移:GT=1,GQ=2,DP=3,AD=4,PL=5。
- GT:
0x1101 0x21 0x0202020404040x1101:键 GT 偏移 1,类型 INT80x21:2 元素 INT8 向量- 值:
0x02 0x02(0/0),0x02 0x04(0/1),0x04 0x04(1/1)
- GQ:
0x1102 0x11 0x0A0A0A(值 10) - DP:
0x1103 0x11 0x203040(值 32,48,64) - AD:
0x1104 0x21 0x300030200040 - PL:
0x1105 0x31 0x000A640A0064640A00
Putting it all together
最终字节块:
0x33000000 l_shared = 51
0x2A000000 l_indiv = 42
0x01000000 CHROM 偏移 1
0x64000000 POS 101
0x01000000 rlen = 1
0x41 0xF0 0xCC 0xCD QUAL = 30.1
0x0400 n_info = 4
0x0200 n_allele = 2
0x030000 n_sample = 3
0x05 n_fmt = 5
0x57 0x72 0x73 0x31 0x32 0x33 ID = rs123
0x17 0x41 REF = A
0x17 0x43 ALT = C
0x11 0x00 FILTER = PASS
0x11 0x50 0x00 HM3 flag
0x11 0x51 0x11 0x03 AC=3
0x11 0x52 0x11 0x06 AN=6
0x11 0x53 0x17 0x43 AA=C
0x1101 0x21 0x020202040404 GT
0x1102 0x11 0x0A0A0A GQ
0x1103 0x11 0x203040 DP
0x1104 0x21 0x300030200040 AD
0x1105 0x31 0x000A640A0064640A00 PL
总共只有 96 字节!这就是 BCF2 高效的原因。
6.5 BCF2 block gzip and indexing:BGZF 压缩与索引
原始二进制记录随后被编码进 BGZF 块,遵循 SAM 规范第 3 节的 BGZF 压缩格式。BCF2 记录也可以是未压缩的,但通常仍用 BGZF 压缩,只是压缩级别为 0。例如,通过管道用 samtools/bcftools 流式处理 BCF2 时,解码/编码成本可能使不压缩更合理。
- GATK 目前生成原始 BCF2 文件(完全不是 BGZF 压缩),但未来会改变。
- BCF2 文件预期用与 BAM 相同的索引方案(SAM 规范第 4 节)进行索引,即 tabix 索引。
总结
BCF2 是 VCF 的二进制压缩版本,核心特点:
- 文件结构:头部 + BGZF 压缩块 + 空块结束。
- 头部:magic
"BCF"、版本、VCF 文本头部、字符串字典。 - 字符串字典:字符串用整数偏移引用,减少重复。
- 记录编码:按字段组织,类型强指定。核心字段:
l_shared、l_indiv、CHROM、POS、rlen、QUAL、n_info、n_allele、n_sample、n_fmt、ID、REF+ALT、FILTER、INFO、Genotype。 - 类型系统:类型描述符字节,高 4 位长度,低 4 位类型。支持 INT8/16/32、FLOAT、CHAR、MISSING。
- 基因型编码:按字段编码,每个字段是 N 个样本的向量。GT 编码为
(allele+1)<<1 | phased。 - 向量:固定长度和混合长度,缺失值用 MISSING 填充。
- 压缩与索引:BGZF 压缩,tabix 索引。
理解 BCF2 的关键是:它和 VCF 表达相同信息,但用二进制、字典、类型化向量来节省空间和加速解析。