Mercurial > repos > rliterman > csp2

annotate CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/lib/python3.8/site-packages/Bio/bgzf.py @ 68:5028fdace37b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 2e9511a184a1ca667c7be0c6321a36dc4e3d116d
author jpayne
date Tue, 18 Mar 2025 16:23:26 -0400
parents
children
rev   line source
jpayne@68 1 #!/usr/bin/env python
jpayne@68 2 # Copyright 2010-2018 by Peter Cock.
jpayne@68 3 # All rights reserved.
jpayne@68 4 #
jpayne@68 5 # This file is part of the Biopython distribution and governed by your
jpayne@68 6 # choice of the "Biopython License Agreement" or the "BSD 3-Clause License".
jpayne@68 7 # Please see the LICENSE file that should have been included as part of this
jpayne@68 8 # package.
jpayne@68 9 r"""Read and write BGZF compressed files (the GZIP variant used in BAM).
jpayne@68 10
jpayne@68 11 The SAM/BAM file format (Sequence Alignment/Map) comes in a plain text
jpayne@68 12 format (SAM), and a compressed binary format (BAM). The latter uses a
jpayne@68 13 modified form of gzip compression called BGZF (Blocked GNU Zip Format),
jpayne@68 14 which can be applied to any file format to provide compression with
jpayne@68 15 efficient random access. BGZF is described together with the SAM/BAM
jpayne@68 16 file format at https://samtools.sourceforge.net/SAM1.pdf
jpayne@68 17
jpayne@68 18 Please read the text below about 'virtual offsets' before using BGZF
jpayne@68 19 files for random access.
jpayne@68 20
jpayne@68 21 Aim of this module
jpayne@68 22 ------------------
jpayne@68 23 The Python gzip library can be used to read BGZF files, since for
jpayne@68 24 decompression they are just (specialised) gzip files. What this
jpayne@68 25 module aims to facilitate is random access to BGZF files (using the
jpayne@68 26 'virtual offset' idea), and writing BGZF files (which means using
jpayne@68 27 suitably sized gzip blocks and writing the extra 'BC' field in the
jpayne@68 28 gzip headers). As in the gzip library, the zlib library is used
jpayne@68 29 internally.
jpayne@68 30
jpayne@68 31 In addition to being required for random access to and writing of
jpayne@68 32 BAM files, the BGZF format can also be used on other sequential
jpayne@68 33 data (in the sense of one record after another), such as most of
jpayne@68 34 the sequence data formats supported in Bio.SeqIO (like FASTA,
jpayne@68 35 FASTQ, GenBank, etc) or large MAF alignments.
jpayne@68 36
jpayne@68 37 The Bio.SeqIO indexing functions use this module to support BGZF files.
jpayne@68 38
jpayne@68 39 Technical Introduction to BGZF
jpayne@68 40 ------------------------------
jpayne@68 41 The gzip file format allows multiple compressed blocks, each of which
jpayne@68 42 could be a stand alone gzip file. As an interesting bonus, this means
jpayne@68 43 you can use Unix ``cat`` to combine two or more gzip files into one by
jpayne@68 44 concatenating them. Also, each block can have one of several compression
jpayne@68 45 levels (including uncompressed, which actually takes up a little bit
jpayne@68 46 more space due to the gzip header).
jpayne@68 47
jpayne@68 48 What the BAM designers realised was that while random access to data
jpayne@68 49 stored in traditional gzip files was slow, breaking the file into
jpayne@68 50 gzip blocks would allow fast random access to each block. To access
jpayne@68 51 a particular piece of the decompressed data, you just need to know
jpayne@68 52 which block it starts in (the offset of the gzip block start), and
jpayne@68 53 how far into the (decompressed) contents of the block you need to
jpayne@68 54 read.
jpayne@68 55
jpayne@68 56 One problem with this is finding the gzip block sizes efficiently.
jpayne@68 57 You can do it with a standard gzip file, but it requires every block
jpayne@68 58 to be decompressed -- and that would be rather slow. Additionally
jpayne@68 59 typical gzip files may use very large blocks.
jpayne@68 60
jpayne@68 61 All that differs in BGZF is that compressed size of each gzip block
jpayne@68 62 is limited to 2^16 bytes, and an extra 'BC' field in the gzip header
jpayne@68 63 records this size. Traditional decompression tools can ignore this,
jpayne@68 64 and unzip the file just like any other gzip file.
jpayne@68 65
jpayne@68 66 The point of this is you can look at the first BGZF block, find out
jpayne@68 67 how big it is from this 'BC' header, and thus seek immediately to
jpayne@68 68 the second block, and so on.
jpayne@68 69
jpayne@68 70 The BAM indexing scheme records read positions using a 64 bit
jpayne@68 71 'virtual offset', comprising ``coffset << 16 | uoffset``, where ``coffset``
jpayne@68 72 is the file offset of the BGZF block containing the start of the read
jpayne@68 73 (unsigned integer using up to 64-16 = 48 bits), and ``uoffset`` is the
jpayne@68 74 offset within the (decompressed) block (unsigned 16 bit integer).
jpayne@68 75
jpayne@68 76 This limits you to BAM files where the last block starts by 2^48
jpayne@68 77 bytes, or 256 petabytes, and the decompressed size of each block
jpayne@68 78 is at most 2^16 bytes, or 64kb. Note that this matches the BGZF
jpayne@68 79 'BC' field size which limits the compressed size of each block to
jpayne@68 80 2^16 bytes, allowing for BAM files to use BGZF with no gzip
jpayne@68 81 compression (useful for intermediate files in memory to reduce
jpayne@68 82 CPU load).
jpayne@68 83
jpayne@68 84 Warning about namespaces
jpayne@68 85 ------------------------
jpayne@68 86 It is considered a bad idea to use "from XXX import ``*``" in Python, because
jpayne@68 87 it pollutes the namespace. This is a real issue with Bio.bgzf (and the
jpayne@68 88 standard Python library gzip) because they contain a function called open
jpayne@68 89 i.e. Suppose you do this:
jpayne@68 90
jpayne@68 91 >>> from Bio.bgzf import *
jpayne@68 92 >>> print(open.__module__)
jpayne@68 93 Bio.bgzf
jpayne@68 94
jpayne@68 95 Or,
jpayne@68 96
jpayne@68 97 >>> from gzip import *
jpayne@68 98 >>> print(open.__module__)
jpayne@68 99 gzip
jpayne@68 100
jpayne@68 101 Notice that the open function has been replaced. You can "fix" this if you
jpayne@68 102 need to by importing the built-in open function:
jpayne@68 103
jpayne@68 104 >>> from builtins import open
jpayne@68 105
jpayne@68 106 However, what we recommend instead is to use the explicit namespace, e.g.
jpayne@68 107
jpayne@68 108 >>> from Bio import bgzf
jpayne@68 109 >>> print(bgzf.open.__module__)
jpayne@68 110 Bio.bgzf
jpayne@68 111
jpayne@68 112
jpayne@68 113 Examples
jpayne@68 114 --------
jpayne@68 115 This is an ordinary GenBank file compressed using BGZF, so it can
jpayne@68 116 be decompressed using gzip,
jpayne@68 117
jpayne@68 118 >>> import gzip
jpayne@68 119 >>> handle = gzip.open("GenBank/NC_000932.gb.bgz", "r")
jpayne@68 120 >>> assert 0 == handle.tell()
jpayne@68 121 >>> line = handle.readline()
jpayne@68 122 >>> assert 80 == handle.tell()
jpayne@68 123 >>> line = handle.readline()
jpayne@68 124 >>> assert 143 == handle.tell()
jpayne@68 125 >>> data = handle.read(70000)
jpayne@68 126 >>> assert 70143 == handle.tell()
jpayne@68 127 >>> handle.close()
jpayne@68 128
jpayne@68 129 We can also access the file using the BGZF reader - but pay
jpayne@68 130 attention to the file offsets which will be explained below:
jpayne@68 131
jpayne@68 132 >>> handle = BgzfReader("GenBank/NC_000932.gb.bgz", "r")
jpayne@68 133 >>> assert 0 == handle.tell()
jpayne@68 134 >>> print(handle.readline().rstrip())
jpayne@68 135 LOCUS NC_000932 154478 bp DNA circular PLN 15-APR-2009
jpayne@68 136 >>> assert 80 == handle.tell()
jpayne@68 137 >>> print(handle.readline().rstrip())
jpayne@68 138 DEFINITION Arabidopsis thaliana chloroplast, complete genome.
jpayne@68 139 >>> assert 143 == handle.tell()
jpayne@68 140 >>> data = handle.read(70000)
jpayne@68 141 >>> assert 987828735 == handle.tell()
jpayne@68 142 >>> print(handle.readline().rstrip())
jpayne@68 143 f="GeneID:844718"
jpayne@68 144 >>> print(handle.readline().rstrip())
jpayne@68 145 CDS complement(join(84337..84771,85454..85843))
jpayne@68 146 >>> offset = handle.seek(make_virtual_offset(55074, 126))
jpayne@68 147 >>> print(handle.readline().rstrip())
jpayne@68 148 68521 tatgtcattc gaaattgtat aaagacaact cctatttaat agagctattt gtgcaagtat
jpayne@68 149 >>> handle.close()
jpayne@68 150
jpayne@68 151 Notice the handle's offset looks different as a BGZF file. This
jpayne@68 152 brings us to the key point about BGZF, which is the block structure:
jpayne@68 153
jpayne@68 154 >>> handle = open("GenBank/NC_000932.gb.bgz", "rb")
jpayne@68 155 >>> for values in BgzfBlocks(handle):
jpayne@68 156 ... print("Raw start %i, raw length %i; data start %i, data length %i" % values)
jpayne@68 157 Raw start 0, raw length 15073; data start 0, data length 65536
jpayne@68 158 Raw start 15073, raw length 17857; data start 65536, data length 65536
jpayne@68 159 Raw start 32930, raw length 22144; data start 131072, data length 65536
jpayne@68 160 Raw start 55074, raw length 22230; data start 196608, data length 65536
jpayne@68 161 Raw start 77304, raw length 14939; data start 262144, data length 43478
jpayne@68 162 Raw start 92243, raw length 28; data start 305622, data length 0
jpayne@68 163 >>> handle.close()
jpayne@68 164
jpayne@68 165 In this example the first three blocks are 'full' and hold 65536 bytes
jpayne@68 166 of uncompressed data. The fourth block isn't full and holds 43478 bytes.
jpayne@68 167 Finally there is a special empty fifth block which takes 28 bytes on
jpayne@68 168 disk and serves as an 'end of file' (EOF) marker. If this is missing,
jpayne@68 169 it is possible your BGZF file is incomplete.
jpayne@68 170
jpayne@68 171 By reading ahead 70,000 bytes we moved into the second BGZF block,
jpayne@68 172 and at that point the BGZF virtual offsets start to look different
jpayne@68 173 to a simple offset into the decompressed data as exposed by the gzip
jpayne@68 174 library.
jpayne@68 175
jpayne@68 176 As an example, consider seeking to the decompressed position 196734.
jpayne@68 177 Since 196734 = 65536 + 65536 + 65536 + 126 = 65536*3 + 126, this
jpayne@68 178 is equivalent to jumping the first three blocks (which in this
jpayne@68 179 specific example are all size 65536 after decompression - which
jpayne@68 180 does not always hold) and starting at byte 126 of the fourth block
jpayne@68 181 (after decompression). For BGZF, we need to know the fourth block's
jpayne@68 182 offset of 55074 and the offset within the block of 126 to get the
jpayne@68 183 BGZF virtual offset.
jpayne@68 184
jpayne@68 185 >>> print(55074 << 16 | 126)
jpayne@68 186 3609329790
jpayne@68 187 >>> print(bgzf.make_virtual_offset(55074, 126))
jpayne@68 188 3609329790
jpayne@68 189
jpayne@68 190 Thus for this BGZF file, decompressed position 196734 corresponds
jpayne@68 191 to the virtual offset 3609329790. However, another BGZF file with
jpayne@68 192 different contents would have compressed more or less efficiently,
jpayne@68 193 so the compressed blocks would be different sizes. What this means
jpayne@68 194 is the mapping between the uncompressed offset and the compressed
jpayne@68 195 virtual offset depends on the BGZF file you are using.
jpayne@68 196
jpayne@68 197 If you are accessing a BGZF file via this module, just use the
jpayne@68 198 handle.tell() method to note the virtual offset of a position you
jpayne@68 199 may later want to return to using handle.seek().
jpayne@68 200
jpayne@68 201 The catch with BGZF virtual offsets is while they can be compared
jpayne@68 202 (which offset comes first in the file), you cannot safely subtract
jpayne@68 203 them to get the size of the data between them, nor add/subtract
jpayne@68 204 a relative offset.
jpayne@68 205
jpayne@68 206 Of course you can parse this file with Bio.SeqIO using BgzfReader,
jpayne@68 207 although there isn't any benefit over using gzip.open(...), unless
jpayne@68 208 you want to index BGZF compressed sequence files:
jpayne@68 209
jpayne@68 210 >>> from Bio import SeqIO
jpayne@68 211 >>> handle = BgzfReader("GenBank/NC_000932.gb.bgz")
jpayne@68 212 >>> record = SeqIO.read(handle, "genbank")
jpayne@68 213 >>> handle.close()
jpayne@68 214 >>> print(record.id)
jpayne@68 215 NC_000932.1
jpayne@68 216
jpayne@68 217 Text Mode
jpayne@68 218 ---------
jpayne@68 219
jpayne@68 220 Like the standard library gzip.open(...), the BGZF code defaults to opening
jpayne@68 221 files in binary mode.
jpayne@68 222
jpayne@68 223 You can request the file be opened in text mode, but beware that this is hard
jpayne@68 224 coded to the simple "latin1" (aka "iso-8859-1") encoding (which includes all
jpayne@68 225 the ASCII characters), which works well with most Western European languages.
jpayne@68 226 However, it is not fully compatible with the more widely used UTF-8 encoding.
jpayne@68 227
jpayne@68 228 In variable width encodings like UTF-8, some single characters in the unicode
jpayne@68 229 text output are represented by multiple bytes in the raw binary form. This is
jpayne@68 230 problematic with BGZF, as we cannot always decode each block in isolation - a
jpayne@68 231 single unicode character could be split over two blocks. This can even happen
jpayne@68 232 with fixed width unicode encodings, as the BGZF block size is not fixed.
jpayne@68 233
jpayne@68 234 Therefore, this module is currently restricted to only support single byte
jpayne@68 235 unicode encodings, such as ASCII, "latin1" (which is a superset of ASCII), or
jpayne@68 236 potentially other character maps (not implemented).
jpayne@68 237
jpayne@68 238 Furthermore, unlike the default text mode on Python 3, we do not attempt to
jpayne@68 239 implement universal new line mode. This transforms the various operating system
jpayne@68 240 new line conventions like Windows (CR LF or "\r\n"), Unix (just LF, "\n"), or
jpayne@68 241 old Macs (just CR, "\r"), into just LF ("\n"). Here we have the same problem -
jpayne@68 242 is "\r" at the end of a block an incomplete Windows style new line?
jpayne@68 243
jpayne@68 244 Instead, you will get the CR ("\r") and LF ("\n") characters as is.
jpayne@68 245
jpayne@68 246 If your data is in UTF-8 or any other incompatible encoding, you must use
jpayne@68 247 binary mode, and decode the appropriate fragments yourself.
jpayne@68 248 """
jpayne@68 249
jpayne@68 250 import struct
jpayne@68 251 import sys
jpayne@68 252 import zlib
jpayne@68 253
jpayne@68 254 from builtins import open as _open
jpayne@68 255
jpayne@68 256 _bgzf_magic = b"\x1f\x8b\x08\x04"
jpayne@68 257 _bgzf_header = b"\x1f\x8b\x08\x04\x00\x00\x00\x00\x00\xff\x06\x00\x42\x43\x02\x00"
jpayne@68 258 _bgzf_eof = b"\x1f\x8b\x08\x04\x00\x00\x00\x00\x00\xff\x06\x00BC\x02\x00\x1b\x00\x03\x00\x00\x00\x00\x00\x00\x00\x00\x00"
jpayne@68 259 _bytes_BC = b"BC"
jpayne@68 260
jpayne@68 261
jpayne@68 262 def open(filename, mode="rb"):
jpayne@68 263 r"""Open a BGZF file for reading, writing or appending.
jpayne@68 264
jpayne@68 265 If text mode is requested, in order to avoid multi-byte characters, this is
jpayne@68 266 hard coded to use the "latin1" encoding, and "\r" and "\n" are passed as is
jpayne@68 267 (without implementing universal new line mode).
jpayne@68 268
jpayne@68 269 If your data is in UTF-8 or any other incompatible encoding, you must use
jpayne@68 270 binary mode, and decode the appropriate fragments yourself.
jpayne@68 271 """
jpayne@68 272 if "r" in mode.lower():
jpayne@68 273 return BgzfReader(filename, mode)
jpayne@68 274 elif "w" in mode.lower() or "a" in mode.lower():
jpayne@68 275 return BgzfWriter(filename, mode)
jpayne@68 276 else:
jpayne@68 277 raise ValueError(f"Bad mode {mode!r}")
jpayne@68 278
jpayne@68 279
jpayne@68 280 def make_virtual_offset(block_start_offset, within_block_offset):
jpayne@68 281 """Compute a BGZF virtual offset from block start and within block offsets.
jpayne@68 282
jpayne@68 283 The BAM indexing scheme records read positions using a 64 bit
jpayne@68 284 'virtual offset', comprising in C terms:
jpayne@68 285
jpayne@68 286 block_start_offset << 16 | within_block_offset
jpayne@68 287
jpayne@68 288 Here block_start_offset is the file offset of the BGZF block
jpayne@68 289 start (unsigned integer using up to 64-16 = 48 bits), and
jpayne@68 290 within_block_offset within the (decompressed) block (unsigned
jpayne@68 291 16 bit integer).
jpayne@68 292
jpayne@68 293 >>> make_virtual_offset(0, 0)
jpayne@68 294 0
jpayne@68 295 >>> make_virtual_offset(0, 1)
jpayne@68 296 1
jpayne@68 297 >>> make_virtual_offset(0, 2**16 - 1)
jpayne@68 298 65535
jpayne@68 299 >>> make_virtual_offset(0, 2**16)
jpayne@68 300 Traceback (most recent call last):
jpayne@68 301 ...
jpayne@68 302 ValueError: Require 0 <= within_block_offset < 2**16, got 65536
jpayne@68 303
jpayne@68 304 >>> 65536 == make_virtual_offset(1, 0)
jpayne@68 305 True
jpayne@68 306 >>> 65537 == make_virtual_offset(1, 1)
jpayne@68 307 True
jpayne@68 308 >>> 131071 == make_virtual_offset(1, 2**16 - 1)
jpayne@68 309 True
jpayne@68 310
jpayne@68 311 >>> 6553600000 == make_virtual_offset(100000, 0)
jpayne@68 312 True
jpayne@68 313 >>> 6553600001 == make_virtual_offset(100000, 1)
jpayne@68 314 True
jpayne@68 315 >>> 6553600010 == make_virtual_offset(100000, 10)
jpayne@68 316 True
jpayne@68 317
jpayne@68 318 >>> make_virtual_offset(2**48, 0)
jpayne@68 319 Traceback (most recent call last):
jpayne@68 320 ...
jpayne@68 321 ValueError: Require 0 <= block_start_offset < 2**48, got 281474976710656
jpayne@68 322
jpayne@68 323 """
jpayne@68 324 if within_block_offset < 0 or within_block_offset >= 65536:
jpayne@68 325 raise ValueError(
jpayne@68 326 "Require 0 <= within_block_offset < 2**16, got %i" % within_block_offset
jpayne@68 327 )
jpayne@68 328 if block_start_offset < 0 or block_start_offset >= 281474976710656:
jpayne@68 329 raise ValueError(
jpayne@68 330 "Require 0 <= block_start_offset < 2**48, got %i" % block_start_offset
jpayne@68 331 )
jpayne@68 332 return (block_start_offset << 16) | within_block_offset
jpayne@68 333
jpayne@68 334
jpayne@68 335 def split_virtual_offset(virtual_offset):
jpayne@68 336 """Divides a 64-bit BGZF virtual offset into block start & within block offsets.
jpayne@68 337
jpayne@68 338 >>> (100000, 0) == split_virtual_offset(6553600000)
jpayne@68 339 True
jpayne@68 340 >>> (100000, 10) == split_virtual_offset(6553600010)
jpayne@68 341 True
jpayne@68 342
jpayne@68 343 """
jpayne@68 344 start = virtual_offset >> 16
jpayne@68 345 return start, virtual_offset ^ (start << 16)
jpayne@68 346
jpayne@68 347
jpayne@68 348 def BgzfBlocks(handle):
jpayne@68 349 """Low level debugging function to inspect BGZF blocks.
jpayne@68 350
jpayne@68 351 Expects a BGZF compressed file opened in binary read mode using
jpayne@68 352 the builtin open function. Do not use a handle from this bgzf
jpayne@68 353 module or the gzip module's open function which will decompress
jpayne@68 354 the file.
jpayne@68 355
jpayne@68 356 Returns the block start offset (see virtual offsets), the block
jpayne@68 357 length (add these for the start of the next block), and the
jpayne@68 358 decompressed length of the blocks contents (limited to 65536 in
jpayne@68 359 BGZF), as an iterator - one tuple per BGZF block.
jpayne@68 360
jpayne@68 361 >>> from builtins import open
jpayne@68 362 >>> handle = open("SamBam/ex1.bam", "rb")
jpayne@68 363 >>> for values in BgzfBlocks(handle):
jpayne@68 364 ... print("Raw start %i, raw length %i; data start %i, data length %i" % values)
jpayne@68 365 Raw start 0, raw length 18239; data start 0, data length 65536
jpayne@68 366 Raw start 18239, raw length 18223; data start 65536, data length 65536
jpayne@68 367 Raw start 36462, raw length 18017; data start 131072, data length 65536
jpayne@68 368 Raw start 54479, raw length 17342; data start 196608, data length 65536
jpayne@68 369 Raw start 71821, raw length 17715; data start 262144, data length 65536
jpayne@68 370 Raw start 89536, raw length 17728; data start 327680, data length 65536
jpayne@68 371 Raw start 107264, raw length 17292; data start 393216, data length 63398
jpayne@68 372 Raw start 124556, raw length 28; data start 456614, data length 0
jpayne@68 373 >>> handle.close()
jpayne@68 374
jpayne@68 375 Indirectly we can tell this file came from an old version of
jpayne@68 376 samtools because all the blocks (except the final one and the
jpayne@68 377 dummy empty EOF marker block) are 65536 bytes. Later versions
jpayne@68 378 avoid splitting a read between two blocks, and give the header
jpayne@68 379 its own block (useful to speed up replacing the header). You
jpayne@68 380 can see this in ex1_refresh.bam created using samtools 0.1.18:
jpayne@68 381
jpayne@68 382 samtools view -b ex1.bam > ex1_refresh.bam
jpayne@68 383
jpayne@68 384 >>> handle = open("SamBam/ex1_refresh.bam", "rb")
jpayne@68 385 >>> for values in BgzfBlocks(handle):
jpayne@68 386 ... print("Raw start %i, raw length %i; data start %i, data length %i" % values)
jpayne@68 387 Raw start 0, raw length 53; data start 0, data length 38
jpayne@68 388 Raw start 53, raw length 18195; data start 38, data length 65434
jpayne@68 389 Raw start 18248, raw length 18190; data start 65472, data length 65409
jpayne@68 390 Raw start 36438, raw length 18004; data start 130881, data length 65483
jpayne@68 391 Raw start 54442, raw length 17353; data start 196364, data length 65519
jpayne@68 392 Raw start 71795, raw length 17708; data start 261883, data length 65411
jpayne@68 393 Raw start 89503, raw length 17709; data start 327294, data length 65466
jpayne@68 394 Raw start 107212, raw length 17390; data start 392760, data length 63854
jpayne@68 395 Raw start 124602, raw length 28; data start 456614, data length 0
jpayne@68 396 >>> handle.close()
jpayne@68 397
jpayne@68 398 The above example has no embedded SAM header (thus the first block
jpayne@68 399 is very small at just 38 bytes of decompressed data), while the next
jpayne@68 400 example does (a larger block of 103 bytes). Notice that the rest of
jpayne@68 401 the blocks show the same sizes (they contain the same read data):
jpayne@68 402
jpayne@68 403 >>> handle = open("SamBam/ex1_header.bam", "rb")
jpayne@68 404 >>> for values in BgzfBlocks(handle):
jpayne@68 405 ... print("Raw start %i, raw length %i; data start %i, data length %i" % values)
jpayne@68 406 Raw start 0, raw length 104; data start 0, data length 103
jpayne@68 407 Raw start 104, raw length 18195; data start 103, data length 65434
jpayne@68 408 Raw start 18299, raw length 18190; data start 65537, data length 65409
jpayne@68 409 Raw start 36489, raw length 18004; data start 130946, data length 65483
jpayne@68 410 Raw start 54493, raw length 17353; data start 196429, data length 65519
jpayne@68 411 Raw start 71846, raw length 17708; data start 261948, data length 65411
jpayne@68 412 Raw start 89554, raw length 17709; data start 327359, data length 65466
jpayne@68 413 Raw start 107263, raw length 17390; data start 392825, data length 63854
jpayne@68 414 Raw start 124653, raw length 28; data start 456679, data length 0
jpayne@68 415 >>> handle.close()
jpayne@68 416
jpayne@68 417 """
jpayne@68 418 if isinstance(handle, BgzfReader):
jpayne@68 419 raise TypeError("Function BgzfBlocks expects a binary handle")
jpayne@68 420 data_start = 0
jpayne@68 421 while True:
jpayne@68 422 start_offset = handle.tell()
jpayne@68 423 try:
jpayne@68 424 block_length, data = _load_bgzf_block(handle)
jpayne@68 425 except StopIteration:
jpayne@68 426 break
jpayne@68 427 data_len = len(data)
jpayne@68 428 yield start_offset, block_length, data_start, data_len
jpayne@68 429 data_start += data_len
jpayne@68 430
jpayne@68 431
jpayne@68 432 def _load_bgzf_block(handle, text_mode=False):
jpayne@68 433 """Load the next BGZF block of compressed data (PRIVATE).
jpayne@68 434
jpayne@68 435 Returns a tuple (block size and data), or at end of file
jpayne@68 436 will raise StopIteration.
jpayne@68 437 """
jpayne@68 438 magic = handle.read(4)
jpayne@68 439 if not magic:
jpayne@68 440 # End of file - should we signal this differently now?
jpayne@68 441 # See https://www.python.org/dev/peps/pep-0479/
jpayne@68 442 raise StopIteration
jpayne@68 443 if magic != _bgzf_magic:
jpayne@68 444 raise ValueError(
jpayne@68 445 r"A BGZF (e.g. a BAM file) block should start with "
jpayne@68 446 r"%r, not %r; handle.tell() now says %r"
jpayne@68 447 % (_bgzf_magic, magic, handle.tell())
jpayne@68 448 )
jpayne@68 449 gzip_mod_time, gzip_extra_flags, gzip_os, extra_len = struct.unpack(
jpayne@68 450 "<LBBH", handle.read(8)
jpayne@68 451 )
jpayne@68 452
jpayne@68 453 block_size = None
jpayne@68 454 x_len = 0
jpayne@68 455 while x_len < extra_len:
jpayne@68 456 subfield_id = handle.read(2)
jpayne@68 457 subfield_len = struct.unpack("<H", handle.read(2))[0] # uint16_t
jpayne@68 458 subfield_data = handle.read(subfield_len)
jpayne@68 459 x_len += subfield_len + 4
jpayne@68 460 if subfield_id == _bytes_BC:
jpayne@68 461 if subfield_len != 2:
jpayne@68 462 raise ValueError("Wrong BC payload length")
jpayne@68 463 if block_size is not None:
jpayne@68 464 raise ValueError("Two BC subfields?")
jpayne@68 465 block_size = struct.unpack("<H", subfield_data)[0] + 1 # uint16_t
jpayne@68 466 if x_len != extra_len:
jpayne@68 467 raise ValueError(f"x_len and extra_len differ {x_len}, {extra_len}")
jpayne@68 468 if block_size is None:
jpayne@68 469 raise ValueError("Missing BC, this isn't a BGZF file!")
jpayne@68 470 # Now comes the compressed data, CRC, and length of uncompressed data.
jpayne@68 471 deflate_size = block_size - 1 - extra_len - 19
jpayne@68 472 d = zlib.decompressobj(-15) # Negative window size means no headers
jpayne@68 473 data = d.decompress(handle.read(deflate_size)) + d.flush()
jpayne@68 474 expected_crc = handle.read(4)
jpayne@68 475 expected_size = struct.unpack("<I", handle.read(4))[0]
jpayne@68 476 if expected_size != len(data):
jpayne@68 477 raise RuntimeError("Decompressed to %i, not %i" % (len(data), expected_size))
jpayne@68 478 # Should cope with a mix of Python platforms...
jpayne@68 479 crc = zlib.crc32(data)
jpayne@68 480 if crc < 0:
jpayne@68 481 crc = struct.pack("<i", crc)
jpayne@68 482 else:
jpayne@68 483 crc = struct.pack("<I", crc)
jpayne@68 484 if expected_crc != crc:
jpayne@68 485 raise RuntimeError(f"CRC is {crc}, not {expected_crc}")
jpayne@68 486 if text_mode:
jpayne@68 487 # Note ISO-8859-1 aka Latin-1 preserves first 256 chars
jpayne@68 488 # (i.e. ASCII), but critically is a single byte encoding
jpayne@68 489 return block_size, data.decode("latin-1")
jpayne@68 490 else:
jpayne@68 491 return block_size, data
jpayne@68 492
jpayne@68 493
jpayne@68 494 class BgzfReader:
jpayne@68 495 r"""BGZF reader, acts like a read only handle but seek/tell differ.
jpayne@68 496
jpayne@68 497 Let's use the BgzfBlocks function to have a peek at the BGZF blocks
jpayne@68 498 in an example BAM file,
jpayne@68 499
jpayne@68 500 >>> from builtins import open
jpayne@68 501 >>> handle = open("SamBam/ex1.bam", "rb")
jpayne@68 502 >>> for values in BgzfBlocks(handle):
jpayne@68 503 ... print("Raw start %i, raw length %i; data start %i, data length %i" % values)
jpayne@68 504 Raw start 0, raw length 18239; data start 0, data length 65536
jpayne@68 505 Raw start 18239, raw length 18223; data start 65536, data length 65536
jpayne@68 506 Raw start 36462, raw length 18017; data start 131072, data length 65536
jpayne@68 507 Raw start 54479, raw length 17342; data start 196608, data length 65536
jpayne@68 508 Raw start 71821, raw length 17715; data start 262144, data length 65536
jpayne@68 509 Raw start 89536, raw length 17728; data start 327680, data length 65536
jpayne@68 510 Raw start 107264, raw length 17292; data start 393216, data length 63398
jpayne@68 511 Raw start 124556, raw length 28; data start 456614, data length 0
jpayne@68 512 >>> handle.close()
jpayne@68 513
jpayne@68 514 Now let's see how to use this block information to jump to
jpayne@68 515 specific parts of the decompressed BAM file:
jpayne@68 516
jpayne@68 517 >>> handle = BgzfReader("SamBam/ex1.bam", "rb")
jpayne@68 518 >>> assert 0 == handle.tell()
jpayne@68 519 >>> magic = handle.read(4)
jpayne@68 520 >>> assert 4 == handle.tell()
jpayne@68 521
jpayne@68 522 So far nothing so strange, we got the magic marker used at the
jpayne@68 523 start of a decompressed BAM file, and the handle position makes
jpayne@68 524 sense. Now however, let's jump to the end of this block and 4
jpayne@68 525 bytes into the next block by reading 65536 bytes,
jpayne@68 526
jpayne@68 527 >>> data = handle.read(65536)
jpayne@68 528 >>> len(data)
jpayne@68 529 65536
jpayne@68 530 >>> assert 1195311108 == handle.tell()
jpayne@68 531
jpayne@68 532 Expecting 4 + 65536 = 65540 were you? Well this is a BGZF 64-bit
jpayne@68 533 virtual offset, which means:
jpayne@68 534
jpayne@68 535 >>> split_virtual_offset(1195311108)
jpayne@68 536 (18239, 4)
jpayne@68 537
jpayne@68 538 You should spot 18239 as the start of the second BGZF block, while
jpayne@68 539 the 4 is the offset into this block. See also make_virtual_offset,
jpayne@68 540
jpayne@68 541 >>> make_virtual_offset(18239, 4)
jpayne@68 542 1195311108
jpayne@68 543
jpayne@68 544 Let's jump back to almost the start of the file,
jpayne@68 545
jpayne@68 546 >>> make_virtual_offset(0, 2)
jpayne@68 547 2
jpayne@68 548 >>> handle.seek(2)
jpayne@68 549 2
jpayne@68 550 >>> handle.close()
jpayne@68 551
jpayne@68 552 Note that you can use the max_cache argument to limit the number of
jpayne@68 553 BGZF blocks cached in memory. The default is 100, and since each
jpayne@68 554 block can be up to 64kb, the default cache could take up to 6MB of
jpayne@68 555 RAM. The cache is not important for reading through the file in one
jpayne@68 556 pass, but is important for improving performance of random access.
jpayne@68 557 """
jpayne@68 558
jpayne@68 559 def __init__(self, filename=None, mode="r", fileobj=None, max_cache=100):
jpayne@68 560 r"""Initialize the class for reading a BGZF file.
jpayne@68 561
jpayne@68 562 You would typically use the top level ``bgzf.open(...)`` function
jpayne@68 563 which will call this class internally. Direct use is discouraged.
jpayne@68 564
jpayne@68 565 Either the ``filename`` (string) or ``fileobj`` (input file object in
jpayne@68 566 binary mode) arguments must be supplied, but not both.
jpayne@68 567
jpayne@68 568 Argument ``mode`` controls if the data will be returned as strings in
jpayne@68 569 text mode ("rt", "tr", or default "r"), or bytes binary mode ("rb"
jpayne@68 570 or "br"). The argument name matches the built-in ``open(...)`` and
jpayne@68 571 standard library ``gzip.open(...)`` function.
jpayne@68 572
jpayne@68 573 If text mode is requested, in order to avoid multi-byte characters,
jpayne@68 574 this is hard coded to use the "latin1" encoding, and "\r" and "\n"
jpayne@68 575 are passed as is (without implementing universal new line mode). There
jpayne@68 576 is no ``encoding`` argument.
jpayne@68 577
jpayne@68 578 If your data is in UTF-8 or any other incompatible encoding, you must
jpayne@68 579 use binary mode, and decode the appropriate fragments yourself.
jpayne@68 580
jpayne@68 581 Argument ``max_cache`` controls the maximum number of BGZF blocks to
jpayne@68 582 cache in memory. Each can be up to 64kb thus the default of 100 blocks
jpayne@68 583 could take up to 6MB of RAM. This is important for efficient random
jpayne@68 584 access, a small value is fine for reading the file in one pass.
jpayne@68 585 """
jpayne@68 586 # TODO - Assuming we can seek, check for 28 bytes EOF empty block
jpayne@68 587 # and if missing warn about possible truncation (as in samtools)?
jpayne@68 588 if max_cache < 1:
jpayne@68 589 raise ValueError("Use max_cache with a minimum of 1")
jpayne@68 590 # Must open the BGZF file in binary mode, but we may want to
jpayne@68 591 # treat the contents as either text or binary (unicode or
jpayne@68 592 # bytes under Python 3)
jpayne@68 593 if filename and fileobj:
jpayne@68 594 raise ValueError("Supply either filename or fileobj, not both")
jpayne@68 595 # Want to reject output modes like w, a, x, +
jpayne@68 596 if mode.lower() not in ("r", "tr", "rt", "rb", "br"):
jpayne@68 597 raise ValueError(
jpayne@68 598 "Must use a read mode like 'r' (default), 'rt', or 'rb' for binary"
jpayne@68 599 )
jpayne@68 600 # If an open file was passed, make sure it was opened in binary mode.
jpayne@68 601 if fileobj:
jpayne@68 602 if fileobj.read(0) != b"":
jpayne@68 603 raise ValueError("fileobj not opened in binary mode")
jpayne@68 604 handle = fileobj
jpayne@68 605 else:
jpayne@68 606 handle = _open(filename, "rb")
jpayne@68 607 self._text = "b" not in mode.lower()
jpayne@68 608 if self._text:
jpayne@68 609 self._newline = "\n"
jpayne@68 610 else:
jpayne@68 611 self._newline = b"\n"
jpayne@68 612 self._handle = handle
jpayne@68 613 self.max_cache = max_cache
jpayne@68 614 self._buffers = {}
jpayne@68 615 self._block_start_offset = None
jpayne@68 616 self._block_raw_length = None
jpayne@68 617 self._load_block(handle.tell())
jpayne@68 618
jpayne@68 619 def _load_block(self, start_offset=None):
jpayne@68 620 if start_offset is None:
jpayne@68 621 # If the file is being read sequentially, then _handle.tell()
jpayne@68 622 # should be pointing at the start of the next block.
jpayne@68 623 # However, if seek has been used, we can't assume that.
jpayne@68 624 start_offset = self._block_start_offset + self._block_raw_length
jpayne@68 625 if start_offset == self._block_start_offset:
jpayne@68 626 self._within_block_offset = 0
jpayne@68 627 return
jpayne@68 628 elif start_offset in self._buffers:
jpayne@68 629 # Already in cache
jpayne@68 630 self._buffer, self._block_raw_length = self._buffers[start_offset]
jpayne@68 631 self._within_block_offset = 0
jpayne@68 632 self._block_start_offset = start_offset
jpayne@68 633 return
jpayne@68 634 # Must hit the disk... first check cache limits,
jpayne@68 635 while len(self._buffers) >= self.max_cache:
jpayne@68 636 # TODO - Implement LRU cache removal?
jpayne@68 637 self._buffers.popitem()
jpayne@68 638 # Now load the block
jpayne@68 639 handle = self._handle
jpayne@68 640 if start_offset is not None:
jpayne@68 641 handle.seek(start_offset)
jpayne@68 642 self._block_start_offset = handle.tell()
jpayne@68 643 try:
jpayne@68 644 block_size, self._buffer = _load_bgzf_block(handle, self._text)
jpayne@68 645 except StopIteration:
jpayne@68 646 # EOF
jpayne@68 647 block_size = 0
jpayne@68 648 if self._text:
jpayne@68 649 self._buffer = ""
jpayne@68 650 else:
jpayne@68 651 self._buffer = b""
jpayne@68 652 self._within_block_offset = 0
jpayne@68 653 self._block_raw_length = block_size
jpayne@68 654 # Finally save the block in our cache,
jpayne@68 655 self._buffers[self._block_start_offset] = self._buffer, block_size
jpayne@68 656
jpayne@68 657 def tell(self):
jpayne@68 658 """Return a 64-bit unsigned BGZF virtual offset."""
jpayne@68 659 if 0 < self._within_block_offset and self._within_block_offset == len(
jpayne@68 660 self._buffer
jpayne@68 661 ):
jpayne@68 662 # Special case where we're right at the end of a (non empty) block.
jpayne@68 663 # For non-maximal blocks could give two possible virtual offsets,
jpayne@68 664 # but for a maximal block can't use 65536 as the within block
jpayne@68 665 # offset. Therefore for consistency, use the next block and a
jpayne@68 666 # within block offset of zero.
jpayne@68 667 return (self._block_start_offset + self._block_raw_length) << 16
jpayne@68 668 else:
jpayne@68 669 # return make_virtual_offset(self._block_start_offset,
jpayne@68 670 # self._within_block_offset)
jpayne@68 671 # TODO - Include bounds checking as in make_virtual_offset?
jpayne@68 672 return (self._block_start_offset << 16) | self._within_block_offset
jpayne@68 673
jpayne@68 674 def seek(self, virtual_offset):
jpayne@68 675 """Seek to a 64-bit unsigned BGZF virtual offset."""
jpayne@68 676 # Do this inline to avoid a function call,
jpayne@68 677 # start_offset, within_block = split_virtual_offset(virtual_offset)
jpayne@68 678 start_offset = virtual_offset >> 16
jpayne@68 679 within_block = virtual_offset ^ (start_offset << 16)
jpayne@68 680 if start_offset != self._block_start_offset:
jpayne@68 681 # Don't need to load the block if already there
jpayne@68 682 # (this avoids a function call since _load_block would do nothing)
jpayne@68 683 self._load_block(start_offset)
jpayne@68 684 if start_offset != self._block_start_offset:
jpayne@68 685 raise ValueError("start_offset not loaded correctly")
jpayne@68 686 if within_block > len(self._buffer):
jpayne@68 687 if not (within_block == 0 and len(self._buffer) == 0):
jpayne@68 688 raise ValueError(
jpayne@68 689 "Within offset %i but block size only %i"
jpayne@68 690 % (within_block, len(self._buffer))
jpayne@68 691 )
jpayne@68 692 self._within_block_offset = within_block
jpayne@68 693 # assert virtual_offset == self.tell(), \
jpayne@68 694 # "Did seek to %i (%i, %i), but tell says %i (%i, %i)" \
jpayne@68 695 # % (virtual_offset, start_offset, within_block,
jpayne@68 696 # self.tell(), self._block_start_offset,
jpayne@68 697 # self._within_block_offset)
jpayne@68 698 return virtual_offset
jpayne@68 699
jpayne@68 700 def read(self, size=-1):
jpayne@68 701 """Read method for the BGZF module."""
jpayne@68 702 if size < 0:
jpayne@68 703 raise NotImplementedError("Don't be greedy, that could be massive!")
jpayne@68 704
jpayne@68 705 result = "" if self._text else b""
jpayne@68 706 while size and self._block_raw_length:
jpayne@68 707 if self._within_block_offset + size <= len(self._buffer):
jpayne@68 708 # This may leave us right at the end of a block
jpayne@68 709 # (lazy loading, don't load the next block unless we have too)
jpayne@68 710 data = self._buffer[
jpayne@68 711 self._within_block_offset : self._within_block_offset + size
jpayne@68 712 ]
jpayne@68 713 self._within_block_offset += size
jpayne@68 714 if not data:
jpayne@68 715 raise ValueError("Must be at least 1 byte")
jpayne@68 716 result += data
jpayne@68 717 break
jpayne@68 718 else:
jpayne@68 719 data = self._buffer[self._within_block_offset :]
jpayne@68 720 size -= len(data)
jpayne@68 721 self._load_block() # will reset offsets
jpayne@68 722 result += data
jpayne@68 723
jpayne@68 724 return result
jpayne@68 725
jpayne@68 726 def readline(self):
jpayne@68 727 """Read a single line for the BGZF file."""
jpayne@68 728 result = "" if self._text else b""
jpayne@68 729 while self._block_raw_length:
jpayne@68 730 i = self._buffer.find(self._newline, self._within_block_offset)
jpayne@68 731 # Three cases to consider,
jpayne@68 732 if i == -1:
jpayne@68 733 # No newline, need to read in more data
jpayne@68 734 data = self._buffer[self._within_block_offset :]
jpayne@68 735 self._load_block() # will reset offsets
jpayne@68 736 result += data
jpayne@68 737 elif i + 1 == len(self._buffer):
jpayne@68 738 # Found new line, but right at end of block (SPECIAL)
jpayne@68 739 data = self._buffer[self._within_block_offset :]
jpayne@68 740 # Must now load the next block to ensure tell() works
jpayne@68 741 self._load_block() # will reset offsets
jpayne@68 742 if not data:
jpayne@68 743 raise ValueError("Must be at least 1 byte")
jpayne@68 744 result += data
jpayne@68 745 break
jpayne@68 746 else:
jpayne@68 747 # Found new line, not at end of block (easy case, no IO)
jpayne@68 748 data = self._buffer[self._within_block_offset : i + 1]
jpayne@68 749 self._within_block_offset = i + 1
jpayne@68 750 # assert data.endswith(self._newline)
jpayne@68 751 result += data
jpayne@68 752 break
jpayne@68 753
jpayne@68 754 return result
jpayne@68 755
jpayne@68 756 def __next__(self):
jpayne@68 757 """Return the next line."""
jpayne@68 758 line = self.readline()
jpayne@68 759 if not line:
jpayne@68 760 raise StopIteration
jpayne@68 761 return line
jpayne@68 762
jpayne@68 763 def __iter__(self):
jpayne@68 764 """Iterate over the lines in the BGZF file."""
jpayne@68 765 return self
jpayne@68 766
jpayne@68 767 def close(self):
jpayne@68 768 """Close BGZF file."""
jpayne@68 769 self._handle.close()
jpayne@68 770 self._buffer = None
jpayne@68 771 self._block_start_offset = None
jpayne@68 772 self._buffers = None
jpayne@68 773
jpayne@68 774 def seekable(self):
jpayne@68 775 """Return True indicating the BGZF supports random access."""
jpayne@68 776 return True
jpayne@68 777
jpayne@68 778 def isatty(self):
jpayne@68 779 """Return True if connected to a TTY device."""
jpayne@68 780 return False
jpayne@68 781
jpayne@68 782 def fileno(self):
jpayne@68 783 """Return integer file descriptor."""
jpayne@68 784 return self._handle.fileno()
jpayne@68 785
jpayne@68 786 def __enter__(self):
jpayne@68 787 """Open a file operable with WITH statement."""
jpayne@68 788 return self
jpayne@68 789
jpayne@68 790 def __exit__(self, type, value, traceback):
jpayne@68 791 """Close a file with WITH statement."""
jpayne@68 792 self.close()
jpayne@68 793
jpayne@68 794
jpayne@68 795 class BgzfWriter:
jpayne@68 796 """Define a BGZFWriter object."""
jpayne@68 797
jpayne@68 798 def __init__(self, filename=None, mode="w", fileobj=None, compresslevel=6):
jpayne@68 799 """Initilize the class."""
jpayne@68 800 if filename and fileobj:
jpayne@68 801 raise ValueError("Supply either filename or fileobj, not both")
jpayne@68 802 # If an open file was passed, make sure it was opened in binary mode.
jpayne@68 803 if fileobj:
jpayne@68 804 if fileobj.read(0) != b"":
jpayne@68 805 raise ValueError("fileobj not opened in binary mode")
jpayne@68 806 handle = fileobj
jpayne@68 807 else:
jpayne@68 808 if "w" not in mode.lower() and "a" not in mode.lower():
jpayne@68 809 raise ValueError(f"Must use write or append mode, not {mode!r}")
jpayne@68 810 if "a" in mode.lower():
jpayne@68 811 handle = _open(filename, "ab")
jpayne@68 812 else:
jpayne@68 813 handle = _open(filename, "wb")
jpayne@68 814 self._text = "b" not in mode.lower()
jpayne@68 815 self._handle = handle
jpayne@68 816 self._buffer = b""
jpayne@68 817 self.compresslevel = compresslevel
jpayne@68 818
jpayne@68 819 def _write_block(self, block):
jpayne@68 820 """Write provided data to file as a single BGZF compressed block (PRIVATE)."""
jpayne@68 821 # print("Saving %i bytes" % len(block))
jpayne@68 822 if len(block) > 65536:
jpayne@68 823 raise ValueError(f"{len(block)} Block length > 65536")
jpayne@68 824 # Giving a negative window bits means no gzip/zlib headers,
jpayne@68 825 # -15 used in samtools
jpayne@68 826 c = zlib.compressobj(
jpayne@68 827 self.compresslevel, zlib.DEFLATED, -15, zlib.DEF_MEM_LEVEL, 0
jpayne@68 828 )
jpayne@68 829 compressed = c.compress(block) + c.flush()
jpayne@68 830 del c
jpayne@68 831 if len(compressed) > 65536:
jpayne@68 832 raise RuntimeError(
jpayne@68 833 "TODO - Didn't compress enough, try less data in this block"
jpayne@68 834 )
jpayne@68 835 crc = zlib.crc32(block)
jpayne@68 836 # Should cope with a mix of Python platforms...
jpayne@68 837 if crc < 0:
jpayne@68 838 crc = struct.pack("<i", crc)
jpayne@68 839 else:
jpayne@68 840 crc = struct.pack("<I", crc)
jpayne@68 841 bsize = struct.pack("<H", len(compressed) + 25) # includes -1
jpayne@68 842 crc = struct.pack("<I", zlib.crc32(block) & 0xFFFFFFFF)
jpayne@68 843 uncompressed_length = struct.pack("<I", len(block))
jpayne@68 844 # Fixed 16 bytes,
jpayne@68 845 # gzip magic bytes (4) mod time (4),
jpayne@68 846 # gzip flag (1), os (1), extra length which is six (2),
jpayne@68 847 # sub field which is BC (2), sub field length of two (2),
jpayne@68 848 # Variable data,
jpayne@68 849 # 2 bytes: block length as BC sub field (2)
jpayne@68 850 # X bytes: the data
jpayne@68 851 # 8 bytes: crc (4), uncompressed data length (4)
jpayne@68 852 data = _bgzf_header + bsize + compressed + crc + uncompressed_length
jpayne@68 853 self._handle.write(data)
jpayne@68 854
jpayne@68 855 def write(self, data):
jpayne@68 856 """Write method for the class."""
jpayne@68 857 # TODO - Check bytes vs unicode
jpayne@68 858 if isinstance(data, str):
jpayne@68 859 # When reading we can't cope with multi-byte characters
jpayne@68 860 # being split between BGZF blocks, so we restrict to a
jpayne@68 861 # single byte encoding - like ASCII or latin-1.
jpayne@68 862 # On output we could probably allow any encoding, as we
jpayne@68 863 # don't care about splitting unicode characters between blocks
jpayne@68 864 data = data.encode("latin-1")
jpayne@68 865 # block_size = 2**16 = 65536
jpayne@68 866 data_len = len(data)
jpayne@68 867 if len(self._buffer) + data_len < 65536:
jpayne@68 868 # print("Cached %r" % data)
jpayne@68 869 self._buffer += data
jpayne@68 870 else:
jpayne@68 871 # print("Got %r, writing out some data..." % data)
jpayne@68 872 self._buffer += data
jpayne@68 873 while len(self._buffer) >= 65536:
jpayne@68 874 self._write_block(self._buffer[:65536])
jpayne@68 875 self._buffer = self._buffer[65536:]
jpayne@68 876
jpayne@68 877 def flush(self):
jpayne@68 878 """Flush data explicitally."""
jpayne@68 879 while len(self._buffer) >= 65536:
jpayne@68 880 self._write_block(self._buffer[:65535])
jpayne@68 881 self._buffer = self._buffer[65535:]
jpayne@68 882 self._write_block(self._buffer)
jpayne@68 883 self._buffer = b""
jpayne@68 884 self._handle.flush()
jpayne@68 885
jpayne@68 886 def close(self):
jpayne@68 887 """Flush data, write 28 bytes BGZF EOF marker, and close BGZF file.
jpayne@68 888
jpayne@68 889 samtools will look for a magic EOF marker, just a 28 byte empty BGZF
jpayne@68 890 block, and if it is missing warns the BAM file may be truncated. In
jpayne@68 891 addition to samtools writing this block, so too does bgzip - so this
jpayne@68 892 implementation does too.
jpayne@68 893 """
jpayne@68 894 if self._buffer:
jpayne@68 895 self.flush()
jpayne@68 896 self._handle.write(_bgzf_eof)
jpayne@68 897 self._handle.flush()
jpayne@68 898 self._handle.close()
jpayne@68 899
jpayne@68 900 def tell(self):
jpayne@68 901 """Return a BGZF 64-bit virtual offset."""
jpayne@68 902 return make_virtual_offset(self._handle.tell(), len(self._buffer))
jpayne@68 903
jpayne@68 904 def seekable(self):
jpayne@68 905 """Return True indicating the BGZF supports random access."""
jpayne@68 906 # Not seekable, but we do support tell...
jpayne@68 907 return False
jpayne@68 908
jpayne@68 909 def isatty(self):
jpayne@68 910 """Return True if connected to a TTY device."""
jpayne@68 911 return False
jpayne@68 912
jpayne@68 913 def fileno(self):
jpayne@68 914 """Return integer file descriptor."""
jpayne@68 915 return self._handle.fileno()
jpayne@68 916
jpayne@68 917 def __enter__(self):
jpayne@68 918 """Open a file operable with WITH statement."""
jpayne@68 919 return self
jpayne@68 920
jpayne@68 921 def __exit__(self, type, value, traceback):
jpayne@68 922 """Close a file with WITH statement."""
jpayne@68 923 self.close()
jpayne@68 924
jpayne@68 925
jpayne@68 926 if __name__ == "__main__":
jpayne@68 927 if len(sys.argv) > 1:
jpayne@68 928 print("Call this with no arguments and pipe uncompressed data in on stdin")
jpayne@68 929 print("and it will produce BGZF compressed data on stdout. e.g.")
jpayne@68 930 print("")
jpayne@68 931 print("./bgzf.py < example.fastq > example.fastq.bgz")
jpayne@68 932 print("")
jpayne@68 933 print("The extension convention of *.bgz is to distinugish these from *.gz")
jpayne@68 934 print("used for standard gzipped files without the block structure of BGZF.")
jpayne@68 935 print("You can use the standard gunzip command to decompress BGZF files,")
jpayne@68 936 print("if it complains about the extension try something like this:")
jpayne@68 937 print("")
jpayne@68 938 print("cat example.fastq.bgz | gunzip > example.fastq")
jpayne@68 939 print("")
jpayne@68 940 print("See also the tool bgzip that comes with samtools")
jpayne@68 941 sys.exit(0)
jpayne@68 942
jpayne@68 943 # Ensure we have binary mode handles
jpayne@68 944 # (leave stderr as default text mode)
jpayne@68 945 stdin = sys.stdin.buffer
jpayne@68 946 stdout = sys.stdout.buffer
jpayne@68 947
jpayne@68 948 sys.stderr.write("Producing BGZF output from stdin...\n")
jpayne@68 949 w = BgzfWriter(fileobj=stdout)
jpayne@68 950 while True:
jpayne@68 951 data = stdin.read(65536)
jpayne@68 952 w.write(data)
jpayne@68 953 if not data:
jpayne@68 954 break
jpayne@68 955 # Doing close will write an empty BGZF block as EOF marker:
jpayne@68 956 w.close()
jpayne@68 957 sys.stderr.write("BGZF data produced\n")