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)