Tabix reader written 100% in Python
Project description
tabixpy
Tabix parser writtern in Python3.
Website
https://pypi.org/project/tabixpy/
Install
pip install tabixpy
Usage
import tabixpy
ingz = "example.vcf.gz"
data = tabixpy.read_tabix(ingz)
tabixpy.save(data, ingz, compress=True)
Tabix
https://samtools.github.io/hts-specs/tabix.pdf
Field Description Type Value
---------------------------------------------------------------------------------------
magic Magic string char[4] TBI 1
n_ref # sequences int32_t
format Format (0: generic; 1: SAM; 2: VCF) int32_t
col_seq Column for the sequence name int32_t
col_beg Column for the start of a region int32_t
col_end Column for the end of a region int32_t
meta Leading character for comment lines int32_t
skip # lines to skip at the beginning int32_t
l_nm Length of concatenated sequence names int32_t
names Concatenated names, each zero terminated char[l_nm]
======================= List of indices (n=n_ref ) =======================
n_bin # distinct bins (for the binning index) int32_t
======================= List of distinct bins (n=n_bin) =======================
bin Distinct bin number uint32_t
n_chunk # chunks int32_t
======================= List of chunks (n=n_chunk) =======================
cnk_beg Virtual file offset of the start of the chunk uint64_t
cnk_end Virtual file offset of the end of the chunk uint64_t
n_intv # 16kb intervals (for the linear index) int32_t
======================= List of distinct intervals (n=n_intv) =======================
ioff File offset of the first record in the interval uint64_t
n_no_coor (optional) # unmapped reads without coordinates set uint64_t
Notes
-
The index file is BGZF compressed.
-
All integers are little-endian.
-
When (format&0x10000) is true, the coordinate follows the BED rule (i.e. half-closed-half-open and zero based); otherwise, the coordinate follows the GFF rule (closed and one based).
-
For the SAM format, the end of a region equals POS plus the reference length in the alignment, inferred from CIGAR. For the VCF format, the end of a region equals POS plus the size of the deletion.
-
Field col beg may equal col end, and in this case, the end of a region is end=beg+1.
-
Example:
- For GFF, format=0 , col seq=1, col beg=4, col end=5, meta=‘#’ and skip=0.
- For BED, format=0x10000, col seq=1, col beg=2, col end=3, meta=‘#’ and skip=0.
-
Given a zero-based, half-closed and half-open region [beg, end), the bin number is calculated with the following C function:
int reg2bin(int beg, int end) {
--end;
if (beg>>14 == end>>14) return ((1<<15)-1)/7 + (beg>>14);
if (beg>>17 == end>>17) return ((1<<12)-1)/7 + (beg>>17);
if (beg>>20 == end>>20) return ((1<< 9)-1)/7 + (beg>>20);
if (beg>>23 == end>>23) return ((1<< 6)-1)/7 + (beg>>23);
if (beg>>26 == end>>26) return ((1<< 3)-1)/7 + (beg>>26);
return 0;
}
- The list of bins that may overlap a region [beg, end) can be obtained with the following C function:
#define MAX_BIN (((1<<18)-1)/7)
int reg2bins(int rbeg, int rend, uint16_t list[MAX_BIN]) {
int i = 0, k;
--rend;
list[i++] = 0;
for (k = 1 + (rbeg>>26); k <= 1 + (rend>>26); ++k) list[i++] = k;
for (k = 9 + (rbeg>>23); k <= 9 + (rend>>23); ++k) list[i++] = k;
for (k = 73 + (rbeg>>20); k <= 73 + (rend>>20); ++k) list[i++] = k;
for (k = 585 + (rbeg>>17); k <= 585 + (rend>>17); ++k) list[i++] = k;
for (k = 4681 + (rbeg>>14); k <= 4681 + (rend>>14); ++k) list[i++] = k;
return i; // #elements in list[]
}
Offset calculation
block_bytes_mask = 0x0000000000000FFFF
real_file_offset = chunk >> 16
block_bytes_offset = chunk & block_bytes_mask
begin
block 0000000000000000000000000110100011010101100111111001011010001111 450,260,604,559
mask 0000000000000000000000000000000000000000000000001111111111111111
real offset 0000000000000000000000000000000000000000011010001101010110011111 6,870,431
block byte offset 0000000000000000000000000000000000000000000000001001011010001111 38,543
end
block 0000000000000000000000000110101001010101110101110011001010010100 456,706,699,924
mask 0000000000000000000000000000000000000000000000001111111111111111
real offset 0000000000000000000000000000000000000000011010100101010111010111 6,968,791
block byte offset 0000000000000000000000000000000000000000000000000011001010010100 12,948
BGZIP
http://samtools.github.io/hts-specs/SAMv1.pdf
The random access method to be described next limits the uncompressed contents
of each BGZF block to a maximum of 216 bytes of data. Thus while ISIZE is
stored as a uint32 t as per the gzip format, in BGZF it is limited to the range
[0, 65536]. BSIZE can represent BGZF block sizes in the range [1, 65536],
though typically BSIZE will be rather less than ISIZE due to compression.
4.1.1 Random access
BGZF files support random access through the BAM file index. To achieve this,
the BAM file index uses virtual file offsets into the BGZF file. Each virtual
file offset is an unsigned 64-bit integer, defined as: coffset<<16|uoffset,
where coffset is an unsigned byte offset into the BGZF file to the beginning of
a BGZF block, and uoffset is an unsigned byte offset into the uncompressed data
stream represented by that BGZF block. Virtual file offsets can be compared,
but subtraction between virtual file offsets and addition between a virtual
offset and an integer are both disallowed.
BGZIP TABIX
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3042176/
https://samtools.github.io/hts-specs/tabix.pdf
2.1 Sorting and BGZF compression
Bgzip compresses the data file in the BGZF format, which is the concatenation
of a series of gzip blocks with each block holding at most 216 bytes of
uncompressed data.
In the compressed file, each uncompressed byte in the text data file is
assigned a unique 64-bit virtual file offset where the higher 48 bits keep the
real file offset of the start of the gzip block the byte falls in, and the
lower 16 bits store the offset of the byte inside the gzip block. Given a
virtual file offset, one can directly seek to the start of the gzip block using
the higher 48 bits, decompress the block and retrieve the byte pointed by the
lower 16 bits of the virtual offset. Random access can thus be achieved without
the help of additional index structures. As gzip works with concatenated gzip
files, it can also seamlessly decompress a BGZF file. The detailed description
of the BGZF format is described in the SAM specification.
2.2 Coupled binning and linear indices
Tabix builds two types of indices for a data file: a binning index and a linear
index. We can actually achieve fast retrieval with only one of them. However,
using the binning index alone may incur many unnecessary seek calls, while
using the linear index alone has bad worst-case performance (when some records
span very long distances). Using them together avoids their weakness.
2.2.1 The binning index
The basic idea of binning is to cluster records into large intervals, called
bins. A record is assigned to bin k if k is the bin of the smallest size that
fully contains the record. For each bin, we keep in the index file the virtual
file offsets of all records assigned to the bin. When we search for records
overlapping a query interval, we first collect bins overlapping the interval
and then test each record in the collected bins for overlaps.
In principle, bins can be selected freely as long as each record can be
assigned to a bin. In the Tabix binning index, we adopt a multilevel binning
scheme where bins at same level are non-overlapping and of the same size. In
Tabix, each bin k, 0<=k<=37,449, represents a half-close-half-open interval
[(k-ol)sl, (k-ol+1)sl), where l = [log2(7k+1)/3] is the level of the bin,
sl = 2^(29-3l) is the size of the bin at level l and ol = (2^3l - 1)/7 is the
offset at l. In this scheme, bin 0 spans 512 Mb, 1-8 span 64 Mb, 9-72 8 Mb,
73-584 1 Mb, 585-4680 128 kb and 4681-37449 span 16 kb intervals. The scheme is
very similar to the UCSC binning (Kent et al., 2002) except that in UCSC,
0<=k<=4681 and therefore the smallest bin size is 128 kb.
2.2.2 The linear index
With the binning index alone, we frequently need to seek to records assigned to
bins at lower levels, in particular bin 0, the bin spanning the entire
sequence. However, frequently records at lower levels do not overlap the query
interval, which leads to unsuccessful seeks and hurts performance especially
when data are transferred over network. The linear index is to reduce
unnecessary seek calls in this case (Table 1).
In the linear index, we keep for each tiling 16 kb window the virtual file
offset of the leftmost record (i.e. having the smallest start coordinate) that
overlaps the window. When we search for records overlapping a query interval,
we will know from the index the leftmost record that possibly overlaps the
query interval. Records having smaller coordinates than this leftmost record
can be skipped and unsuccessful seek calls can be saved.
When the data files are sorted by coordinates, records assigned to the same bin
tend to be adjacent. Thus in the index file, we do not need to keep the virtual
file offset of each record, but only to keep the start offset of a chunk of
records assigned to the same bin.
Schema
https://www.liquid-technologies.com/online-json-schema-validator
https://www.jsonschemavalidator.net/
Example output
{
"n_ref": 1,
"format": 2,
"col_seq": 1,
"col_beg": 2,
"col_end": 0,
"meta": "#",
"skip": 0,
"l_nm": 11,
"names": [
"SL2.50ch00"
],
"refs": [{
"ref_n": 0,
"ref_name": "SL2.50ch00",
"n_bin": 86,
"bins": [{
"bin_n": 0,
"bin": 4681,
"n_chunk": 1,
"chunks": {
"chunk_begin": [{
"bin_n": 0,
"chunk_n": 0,
"real": 0,
"bytes": 29542,
"block_len": 8031,
"bin_pos": 280,
"first_pos": 280,
"last_pos": 1506
}],
"chunk_end": [{
"bin_n": 0,
"chunk_n": 0,
"real": 124525,
"bytes": 19630,
"block_len": 9015,
"bin_pos": 16388,
"first_pos": 16141,
"last_pos": 17830
}]
},
"chunks_begin": {
"bin_n": 0,
"chunk_n": 0,
"real": 0,
"bytes": 29542,
"block_len": 8031,
"bin_pos": 280,
"first_pos": 280,
"last_pos": 1506
},
"chunks_end": {
"bin_n": 0,
"chunk_n": 0,
"real": 124525,
"bytes": 19630,
"block_len": 9015,
"bin_pos": 16388,
"first_pos": 16141,
"last_pos": 17830
}
},
{
"bin_n": 85,
"bin": 4766,
"n_chunk": 1,
"chunks": {
"chunk_begin": [{
"bin_n": 84,
"chunk_n": 0,
"real": 7021611,
"bytes": 4631,
"block_len": 6621,
"bin_pos": 1392700,
"first_pos": 1392519,
"last_pos": 1393974
}],
"chunk_end": [{
"bin_n": 85,
"chunk_n": 0,
"real": 7039684,
"bytes": 0,
"block_len": -1,
"bin_pos": -1,
"first_pos": -1,
"last_pos": -1
}]
},
"chunks_begin": {
"bin_n": 84,
"chunk_n": 0,
"real": 7021611,
"bytes": 4631,
"block_len": 6621,
"bin_pos": 1392700,
"first_pos": 1392519,
"last_pos": 1393974
},
"chunks_end": {
"bin_n": 85,
"chunk_n": 0,
"real": 7039684,
"bytes": 0,
"block_len": -1,
"bin_pos": -1,
"first_pos": -1,
"last_pos": -1
}
}
],
"bins_begin": {
"bin_n": 0,
"chunk_n": 0,
"real": 0,
"bytes": 29542,
"block_len": 8031,
"bin_pos": 280,
"first_pos": 280,
"last_pos": 1506
},
"bins_end": {
"bin_n": 84,
"chunk_n": 0,
"real": 7021611,
"bytes": 4631,
"block_len": 6621,
"bin_pos": 1392700,
"first_pos": 1392519,
"last_pos": 1393974
},
"first_block": {
"bin_n": 0,
"chunk_n": 0,
"real": 0,
"bytes": 29542,
"block_len": 8031,
"bin_pos": 280,
"first_pos": 280,
"last_pos": 1506
},
"last_block": {
"bin_n": -1,
"chunk_n": -1,
"real": 7035849,
"bytes": 0,
"block_len": 3835,
"bin_pos": -1,
"first_pos": 1395124,
"last_pos": 1395638
},
"n_intv": 86,
"intvs": [{
"bin_n": 0,
"chunk_n": 0,
"real": 0,
"bytes": 29542,
"block_len": 8031,
"bin_pos": 280,
"first_pos": 280,
"last_pos": 1506
},
{
"bin_n": 84,
"chunk_n": 0,
"real": 7021611,
"bytes": 4631,
"block_len": 6621,
"bin_pos": 1392700,
"first_pos": 1392519,
"last_pos": 1393974
}
]
}],
"n_no_coor": null,
"__format_name__": "TBJ",
"__format_ver__": 5
}
Timming
2020-06-12 11:30:34,716 - tabixpy - INFO - reading annotated_tomato_150.100000.vcf.gz
2020-06-12 11:30:34,738 - tabixpy - INFO - saving annotated_tomato_150.100000.vcf.gz.tbj
,024
2020-06-12 11:31:16,506 - tabixpy - INFO - reading annotated_tomato_150.vcf.bgz
2020-06-12 11:31:24,152 - tabixpy - INFO - saving annotated_tomato_150.vcf.bgz.tbj
8,646
File Sizes
TBI Tabix Index
TBK Binary TabixPy 'all chunks' index
TBJ Compressed JSON Tabix index
TBJ.json Uncompressed JSON Tabix Index
6.8M annotated_tomato_150.100000.vcf.gz
1.1K annotated_tomato_150.100000.vcf.gz.tbi
5.9K annotated_tomato_150.100000.vcf.gz.tbk 6.5X
8.9K annotated_tomato_150.100000.vcf.gz.tbj 8.1X
104K annotated_tomato_150.100000.vcf.gz.tbj.json 94.5X
44M annotated_tomato_150.SL2.50ch00-01-02.vcf.gz
4.3K annotated_tomato_150.SL2.50ch00-01-02.vcf.gz.tbi
40K annotated_tomato_150.SL2.50ch00-01-02.vcf.gz.tbk 9.3X
39K annotated_tomato_150.SL2.50ch00-01-02.vcf.gz.tbj 9.1X
468K annotated_tomato_150.SL2.50ch00-01-02.vcf.gz.tbj.json 108.3X
5.6G annotated_tomato_150.vcf.bgz
727K annotated_tomato_150.vcf.bgz.tbi
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.
Source Distribution
tabixpy-1.5.2.tar.gz
(53.9 MB
view hashes)
Built Distribution
tabixpy-1.5.2-py3-none-any.whl
(53.9 MB
view hashes)
