Mercurial > repos > rliterman > csp2

annotate CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/lib/python3.8/site-packages/pysam/libchtslib.pxd @ 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 # cython: language_level=3
jpayne@68 2 from libc.stdint cimport int8_t, int16_t, int32_t, int64_t
jpayne@68 3 from libc.stdint cimport uint8_t, uint16_t, uint32_t, uint64_t
jpayne@68 4 from libc.stdlib cimport malloc, calloc, realloc, free
jpayne@68 5 from libc.string cimport memcpy, memcmp, strncpy, strlen, strdup
jpayne@68 6 from libc.stdio cimport FILE, printf
jpayne@68 7 from posix.types cimport off_t
jpayne@68 8
jpayne@68 9 cdef extern from "Python.h":
jpayne@68 10 FILE* PyFile_AsFile(object)
jpayne@68 11
jpayne@68 12
jpayne@68 13 # cython does not wrap stdarg
jpayne@68 14 cdef extern from "stdarg.h":
jpayne@68 15 ctypedef struct va_list:
jpayne@68 16 pass
jpayne@68 17
jpayne@68 18
jpayne@68 19 cdef extern from "htslib/kstring.h" nogil:
jpayne@68 20 ctypedef struct kstring_t:
jpayne@68 21 size_t l, m
jpayne@68 22 char *s
jpayne@68 23
jpayne@68 24 int kputc(int c, kstring_t *s)
jpayne@68 25 int kputw(int c, kstring_t *s)
jpayne@68 26 int kputl(long c, kstring_t *s)
jpayne@68 27 int ksprintf(kstring_t *s, const char *fmt, ...)
jpayne@68 28
jpayne@68 29
jpayne@68 30 cdef extern from "htslib_util.h" nogil:
jpayne@68 31 int hts_set_verbosity(int verbosity)
jpayne@68 32 int hts_get_verbosity()
jpayne@68 33
jpayne@68 34 ctypedef uint32_t khint32_t
jpayne@68 35 ctypedef uint32_t khint_t
jpayne@68 36 ctypedef khint_t khiter_t
jpayne@68 37
jpayne@68 38 # Used to manage BCF Header info
jpayne@68 39 ctypedef struct vdict_t:
jpayne@68 40 khint_t n_buckets, size, n_occupied, upper_bound
jpayne@68 41 khint32_t *flags
jpayne@68 42 const char *keys
jpayne@68 43 bcf_idinfo_t *vals
jpayne@68 44
jpayne@68 45 # Used to manage indexed contigs in Tabix
jpayne@68 46 ctypedef struct s2i_t:
jpayne@68 47 khint_t n_buckets, size, n_occupied, upper_bound
jpayne@68 48 khint32_t *flags
jpayne@68 49 const char *keys
jpayne@68 50 int64_t *vals
jpayne@68 51
jpayne@68 52 # Generic khash methods
jpayne@68 53 khint_t kh_size(void *d)
jpayne@68 54 khint_t kh_begin(void *d)
jpayne@68 55 khint_t kh_end(void *d)
jpayne@68 56 int kh_exist(void *d, khiter_t i)
jpayne@68 57
jpayne@68 58 # Specialized khash methods for vdict
jpayne@68 59 khint_t kh_get_vdict(vdict_t *d, const char *key)
jpayne@68 60 const char *kh_key_vdict "kh_key" (vdict_t *d, khint_t i)
jpayne@68 61 bcf_idinfo_t kh_val_vdict "kh_val" (vdict_t *d, khint_t i)
jpayne@68 62
jpayne@68 63
jpayne@68 64 cdef extern from "htslib/hfile.h" nogil:
jpayne@68 65 ctypedef struct hFILE
jpayne@68 66
jpayne@68 67 # @abstract Open the named file or URL as a stream
jpayne@68 68 # @return An hFILE pointer, or NULL (with errno set) if an error occurred.
jpayne@68 69 hFILE *hopen(const char *filename, const char *mode, ...)
jpayne@68 70
jpayne@68 71 # @abstract Associate a stream with an existing open file descriptor
jpayne@68 72 # @return An hFILE pointer, or NULL (with errno set) if an error occurred.
jpayne@68 73 # @notes For socket descriptors (on Windows), mode should contain 's'.
jpayne@68 74 hFILE *hdopen(int fd, const char *mode)
jpayne@68 75
jpayne@68 76 # @abstract Report whether the file name or URL denotes remote storage
jpayne@68 77 # @return 0 if local, 1 if remote.
jpayne@68 78 # @notes "Remote" means involving e.g. explicit network access, with the
jpayne@68 79 # implication that callers may wish to cache such files' contents locally.
jpayne@68 80 int hisremote(const char *filename)
jpayne@68 81
jpayne@68 82 # @abstract Flush (for output streams) and close the stream
jpayne@68 83 # @return 0 if successful, or EOF (with errno set) if an error occurred.
jpayne@68 84 int hclose(hFILE *fp)
jpayne@68 85
jpayne@68 86 # @abstract Close the stream, without flushing or propagating errors
jpayne@68 87 # @notes For use while cleaning up after an error only. Preserves errno.
jpayne@68 88 void hclose_abruptly(hFILE *fp)
jpayne@68 89
jpayne@68 90 # @abstract Return the stream's error indicator
jpayne@68 91 # @return Non-zero (in fact, an errno value) if an error has occurred.
jpayne@68 92 # @notes This would be called herror() and return true/false to parallel
jpayne@68 93 # ferror(3), but a networking-related herror(3) function already exists. */
jpayne@68 94 int herrno(hFILE *fp)
jpayne@68 95
jpayne@68 96 # @abstract Clear the stream's error indicator
jpayne@68 97 void hclearerr(hFILE *fp)
jpayne@68 98
jpayne@68 99 # @abstract Reposition the read/write stream offset
jpayne@68 100 # @return The resulting offset within the stream (as per lseek(2)),
jpayne@68 101 # or negative if an error occurred.
jpayne@68 102 off_t hseek(hFILE *fp, off_t offset, int whence)
jpayne@68 103
jpayne@68 104 # @abstract Report the current stream offset
jpayne@68 105 # @return The offset within the stream, starting from zero.
jpayne@68 106 off_t htell(hFILE *fp)
jpayne@68 107
jpayne@68 108 # @abstract Read one character from the stream
jpayne@68 109 # @return The character read, or EOF on end-of-file or error
jpayne@68 110 int hgetc(hFILE *fp)
jpayne@68 111
jpayne@68 112 # Read from the stream until the delimiter, up to a maximum length
jpayne@68 113 # @param buffer The buffer into which bytes will be written
jpayne@68 114 # @param size The size of the buffer
jpayne@68 115 # @param delim The delimiter (interpreted as an `unsigned char`)
jpayne@68 116 # @param fp The file stream
jpayne@68 117 # @return The number of bytes read, or negative on error.
jpayne@68 118 # @since 1.4
jpayne@68 119 #
jpayne@68 120 # Bytes will be read into the buffer up to and including a delimiter, until
jpayne@68 121 # EOF is reached, or _size-1_ bytes have been written, whichever comes first.
jpayne@68 122 # The string will then be terminated with a NUL byte (`\0`).
jpayne@68 123 ssize_t hgetdelim(char *buffer, size_t size, int delim, hFILE *fp)
jpayne@68 124
jpayne@68 125 # Read a line from the stream, up to a maximum length
jpayne@68 126 # @param buffer The buffer into which bytes will be written
jpayne@68 127 # @param size The size of the buffer
jpayne@68 128 # @param fp The file stream
jpayne@68 129 # @return The number of bytes read, or negative on error.
jpayne@68 130 # @since 1.4
jpayne@68 131 #
jpayne@68 132 # Specialization of hgetdelim() for a `\n` delimiter.
jpayne@68 133 ssize_t hgetln(char *buffer, size_t size, hFILE *fp)
jpayne@68 134
jpayne@68 135 # Read a line from the stream, up to a maximum length
jpayne@68 136 # @param buffer The buffer into which bytes will be written
jpayne@68 137 # @param size The size of the buffer (must be > 1 to be useful)
jpayne@68 138 # @param fp The file stream
jpayne@68 139 # @return _buffer_ on success, or `NULL` if an error occurred.
jpayne@68 140 # @since 1.4
jpayne@68 141 #
jpayne@68 142 # This function can be used as a replacement for `fgets(3)`, or together with
jpayne@68 143 # kstring's `kgetline()` to read arbitrarily-long lines into a _kstring_t_.
jpayne@68 144 char *hgets(char *buffer, int size, hFILE *fp)
jpayne@68 145
jpayne@68 146 # @abstract Peek at characters to be read without removing them from buffers
jpayne@68 147 # @param fp The file stream
jpayne@68 148 # @param buffer The buffer to which the peeked bytes will be written
jpayne@68 149 # @param nbytes The number of bytes to peek at; limited by the size of the
jpayne@68 150 # internal buffer, which could be as small as 4K.
jpayne@68 151 # @return The number of bytes peeked, which may be less than nbytes if EOF
jpayne@68 152 # is encountered; or negative, if there was an I/O error.
jpayne@68 153 # @notes The characters peeked at remain in the stream's internal buffer,
jpayne@68 154 # and will be returned by later hread() etc calls.
jpayne@68 155 ssize_t hpeek(hFILE *fp, void *buffer, size_t nbytes)
jpayne@68 156
jpayne@68 157 # @abstract Read a block of characters from the file
jpayne@68 158 # @return The number of bytes read, or negative if an error occurred.
jpayne@68 159 # @notes The full nbytes requested will be returned, except as limited
jpayne@68 160 # by EOF or I/O errors.
jpayne@68 161 ssize_t hread(hFILE *fp, void *buffer, size_t nbytes)
jpayne@68 162
jpayne@68 163 # @abstract Write a character to the stream
jpayne@68 164 # @return The character written, or EOF if an error occurred.
jpayne@68 165 int hputc(int c, hFILE *fp)
jpayne@68 166
jpayne@68 167 # @abstract Write a string to the stream
jpayne@68 168 # @return 0 if successful, or EOF if an error occurred.
jpayne@68 169 int hputs(const char *text, hFILE *fp)
jpayne@68 170
jpayne@68 171 # @abstract Write a block of characters to the file
jpayne@68 172 # @return Either nbytes, or negative if an error occurred.
jpayne@68 173 # @notes In the absence of I/O errors, the full nbytes will be written.
jpayne@68 174 ssize_t hwrite(hFILE *fp, const void *buffer, size_t nbytes)
jpayne@68 175
jpayne@68 176 # @abstract For writing streams, flush buffered output to the underlying stream
jpayne@68 177 # @return 0 if successful, or EOF if an error occurred.
jpayne@68 178 int hflush(hFILE *fp)
jpayne@68 179
jpayne@68 180
jpayne@68 181 cdef extern from "htslib/bgzf.h" nogil:
jpayne@68 182 ctypedef struct bgzf_mtaux_t
jpayne@68 183 ctypedef struct bgzidx_t
jpayne@68 184 ctypedef struct z_stream
jpayne@68 185
jpayne@68 186 ctypedef struct BGZF:
jpayne@68 187 unsigned errcode
jpayne@68 188 unsigned is_write
jpayne@68 189 int is_be
jpayne@68 190 int compress_level
jpayne@68 191 int is_compressed
jpayne@68 192 int is_gzip
jpayne@68 193 int cache_size
jpayne@68 194 int64_t block_address
jpayne@68 195 int64_t uncompressed_address
jpayne@68 196 void *uncompressed_block
jpayne@68 197 void *compressed_block
jpayne@68 198 void *cache
jpayne@68 199 hFILE *fp
jpayne@68 200 bgzf_mtaux_t *mt
jpayne@68 201 bgzidx_t *idx
jpayne@68 202 int idx_build_otf
jpayne@68 203 z_stream *gz_stream
jpayne@68 204
jpayne@68 205 #*****************
jpayne@68 206 # Basic routines *
jpayne@68 207 # *****************/
jpayne@68 208
jpayne@68 209 # Open an existing file descriptor for reading or writing.
jpayne@68 210 #
jpayne@68 211 # @param fd file descriptor
jpayne@68 212 # @param mode mode matching /[rwag][u0-9]+/: 'r' for reading, 'w' for
jpayne@68 213 # writing, 'a' for appending, 'g' for gzip rather than BGZF
jpayne@68 214 # compression (with 'w' only), and digit specifies the zlib
jpayne@68 215 # compression level.
jpayne@68 216 # Note that there is a distinction between 'u' and '0': the
jpayne@68 217 # first yields plain uncompressed output whereas the latter
jpayne@68 218 # outputs uncompressed data wrapped in the zlib format.
jpayne@68 219 # @return BGZF file handler; 0 on error
jpayne@68 220
jpayne@68 221 BGZF* bgzf_dopen(int fd, const char *mode)
jpayne@68 222 BGZF* bgzf_fdopen(int fd, const char *mode) # for backward compatibility
jpayne@68 223
jpayne@68 224 # Open the specified file for reading or writing.
jpayne@68 225 BGZF* bgzf_open(const char* path, const char *mode)
jpayne@68 226
jpayne@68 227 # Open an existing hFILE stream for reading or writing.
jpayne@68 228 BGZF* bgzf_hopen(hFILE *fp, const char *mode)
jpayne@68 229
jpayne@68 230 # Close the BGZF and free all associated resources.
jpayne@68 231 #
jpayne@68 232 # @param fp BGZF file handler
jpayne@68 233 # @return 0 on success and -1 on error
jpayne@68 234 int bgzf_close(BGZF *fp)
jpayne@68 235
jpayne@68 236 # Read up to _length_ bytes from the file storing into _data_.
jpayne@68 237 #
jpayne@68 238 # @param fp BGZF file handler
jpayne@68 239 # @param data data array to read into
jpayne@68 240 # @param length size of data to read
jpayne@68 241 # @return number of bytes actually read; 0 on end-of-file and -1 on error
jpayne@68 242 ssize_t bgzf_read(BGZF *fp, void *data, size_t length)
jpayne@68 243
jpayne@68 244 # Write _length_ bytes from _data_ to the file. If no I/O errors occur,
jpayne@68 245 # the complete _length_ bytes will be written (or queued for writing).
jpayne@68 246 #
jpayne@68 247 # @param fp BGZF file handler
jpayne@68 248 # @param data data array to write
jpayne@68 249 # @param length size of data to write
jpayne@68 250 # @return number of bytes written (i.e., _length_); negative on error
jpayne@68 251 ssize_t bgzf_write(BGZF *fp, const void *data, size_t length)
jpayne@68 252
jpayne@68 253 # Read up to _length_ bytes directly from the underlying stream without
jpayne@68 254 # decompressing. Bypasses BGZF blocking, so must be used with care in
jpayne@68 255 # specialised circumstances only.
jpayne@68 256 #
jpayne@68 257 # @param fp BGZF file handler
jpayne@68 258 # @param data data array to read into
jpayne@68 259 # @param length number of raw bytes to read
jpayne@68 260 # @return number of bytes actually read; 0 on end-of-file and -1 on error
jpayne@68 261 ssize_t bgzf_raw_read(BGZF *fp, void *data, size_t length)
jpayne@68 262
jpayne@68 263 # Write _length_ bytes directly to the underlying stream without
jpayne@68 264 # compressing. Bypasses BGZF blocking, so must be used with care
jpayne@68 265 # in specialised circumstances only.
jpayne@68 266 #
jpayne@68 267 # @param fp BGZF file handler
jpayne@68 268 # @param data data array to write
jpayne@68 269 # @param length number of raw bytes to write
jpayne@68 270 # @return number of bytes actually written; -1 on error
jpayne@68 271 ssize_t bgzf_raw_write(BGZF *fp, const void *data, size_t length)
jpayne@68 272
jpayne@68 273 # Write the data in the buffer to the file.
jpayne@68 274 int bgzf_flush(BGZF *fp)
jpayne@68 275
jpayne@68 276 # Return a virtual file pointer to the current location in the file.
jpayne@68 277 # No interpretation of the value should be made, other than a subsequent
jpayne@68 278 # call to bgzf_seek can be used to position the file at the same point.
jpayne@68 279 # Return value is non-negative on success.
jpayne@68 280 int64_t bgzf_tell(BGZF *fp)
jpayne@68 281
jpayne@68 282 # Set the file to read from the location specified by _pos_.
jpayne@68 283 #
jpayne@68 284 # @param fp BGZF file handler
jpayne@68 285 # @param pos virtual file offset returned by bgzf_tell()
jpayne@68 286 # @param whence must be SEEK_SET (cimported from libc.stdio / posix.unistd)
jpayne@68 287 # @return 0 on success and -1 on error
jpayne@68 288 # /
jpayne@68 289 int64_t bgzf_seek(BGZF *fp, int64_t pos, int whence)
jpayne@68 290
jpayne@68 291 # Check if the BGZF end-of-file (EOF) marker is present
jpayne@68 292 #
jpayne@68 293 # @param fp BGZF file handler opened for reading
jpayne@68 294 # @return 1 if the EOF marker is present and correct
jpayne@68 295 # 2 if it can't be checked, e.g., because fp isn't seekable
jpayne@68 296 # 0 if the EOF marker is absent
jpayne@68 297 # -1 (with errno set) on error
jpayne@68 298 int bgzf_check_EOF(BGZF *fp)
jpayne@68 299
jpayne@68 300 # Check if a file is in the BGZF format
jpayne@68 301 #
jpayne@68 302 # @param fn file name
jpayne@68 303 # @return 1 if _fn_ is BGZF; 0 if not or on I/O error
jpayne@68 304 int bgzf_is_bgzf(const char *fn)
jpayne@68 305
jpayne@68 306 #*********************
jpayne@68 307 # Advanced routines *
jpayne@68 308 #*********************
jpayne@68 309
jpayne@68 310 # Set the cache size. Only effective when compiled with -DBGZF_CACHE.
jpayne@68 311 #
jpayne@68 312 # @param fp BGZF file handler
jpayne@68 313 # @param size size of cache in bytes; 0 to disable caching (default)
jpayne@68 314 void bgzf_set_cache_size(BGZF *fp, int size)
jpayne@68 315
jpayne@68 316 # Flush the file if the remaining buffer size is smaller than _size_
jpayne@68 317 # @return 0 if flushing succeeded or was not needed; negative on error
jpayne@68 318 int bgzf_flush_try(BGZF *fp, ssize_t size)
jpayne@68 319
jpayne@68 320 # Read one byte from a BGZF file. It is faster than bgzf_read()
jpayne@68 321 # @param fp BGZF file handler
jpayne@68 322 # @return byte read; -1 on end-of-file or error
jpayne@68 323 int bgzf_getc(BGZF *fp)
jpayne@68 324
jpayne@68 325 # Read one line from a BGZF file. It is faster than bgzf_getc()
jpayne@68 326 #
jpayne@68 327 # @param fp BGZF file handler
jpayne@68 328 # @param delim delimiter
jpayne@68 329 # @param str string to write to; must be initialized
jpayne@68 330 # @return length of the string; 0 on end-of-file; negative on error
jpayne@68 331 int bgzf_getline(BGZF *fp, int delim, kstring_t *str)
jpayne@68 332
jpayne@68 333 # Read the next BGZF block.
jpayne@68 334 int bgzf_read_block(BGZF *fp)
jpayne@68 335
jpayne@68 336 # Enable multi-threading (only effective on writing and when the
jpayne@68 337 # library was compiled with -DBGZF_MT)
jpayne@68 338 #
jpayne@68 339 # @param fp BGZF file handler; must be opened for writing
jpayne@68 340 # @param n_threads #threads used for writing
jpayne@68 341 # @param n_sub_blks #blocks processed by each thread; a value 64-256 is recommended
jpayne@68 342 int bgzf_mt(BGZF *fp, int n_threads, int n_sub_blks)
jpayne@68 343
jpayne@68 344
jpayne@68 345 # Compress a single BGZF block.
jpayne@68 346 #
jpayne@68 347 # @param dst output buffer (must have size >= BGZF_MAX_BLOCK_SIZE)
jpayne@68 348 # @param dlen size of output buffer; updated on return to the number
jpayne@68 349 # of bytes actually written to dst
jpayne@68 350 # @param src buffer to be compressed
jpayne@68 351 # @param slen size of data to compress (must be <= BGZF_BLOCK_SIZE)
jpayne@68 352 # @param level compression level
jpayne@68 353 # @return 0 on success and negative on error
jpayne@68 354 #
jpayne@68 355 int bgzf_compress(void *dst, size_t *dlen, const void *src, size_t slen, int level)
jpayne@68 356
jpayne@68 357 #*******************
jpayne@68 358 # bgzidx routines *
jpayne@68 359 # BGZF at the uncompressed offset
jpayne@68 360 #
jpayne@68 361 # @param fp BGZF file handler; must be opened for reading
jpayne@68 362 # @param uoffset file offset in the uncompressed data
jpayne@68 363 # @param where SEEK_SET (cimported from libc.stdio) supported atm
jpayne@68 364 #
jpayne@68 365 # Returns 0 on success and -1 on error.
jpayne@68 366 int bgzf_useek(BGZF *fp, long uoffset, int where)
jpayne@68 367
jpayne@68 368 # Position in uncompressed BGZF
jpayne@68 369 #
jpayne@68 370 # @param fp BGZF file handler; must be opened for reading
jpayne@68 371 #
jpayne@68 372 # Returns the current offset on success and -1 on error.
jpayne@68 373 long bgzf_utell(BGZF *fp)
jpayne@68 374
jpayne@68 375 # Tell BGZF to build index while compressing.
jpayne@68 376 #
jpayne@68 377 # @param fp BGZF file handler; can be opened for reading or writing.
jpayne@68 378 #
jpayne@68 379 # Returns 0 on success and -1 on error.
jpayne@68 380 int bgzf_index_build_init(BGZF *fp)
jpayne@68 381
jpayne@68 382 # Load BGZF index
jpayne@68 383 #
jpayne@68 384 # @param fp BGZF file handler
jpayne@68 385 # @param bname base name
jpayne@68 386 # @param suffix suffix to add to bname (can be NULL)
jpayne@68 387 #
jpayne@68 388 # Returns 0 on success and -1 on error.
jpayne@68 389 int bgzf_index_load(BGZF *fp, const char *bname, const char *suffix)
jpayne@68 390
jpayne@68 391 # Save BGZF index
jpayne@68 392 #
jpayne@68 393 # @param fp BGZF file handler
jpayne@68 394 # @param bname base name
jpayne@68 395 # @param suffix suffix to add to bname (can be NULL)
jpayne@68 396 #
jpayne@68 397 # Returns 0 on success and -1 on error.
jpayne@68 398 int bgzf_index_dump(BGZF *fp, const char *bname, const char *suffix)
jpayne@68 399
jpayne@68 400
jpayne@68 401 cdef extern from "htslib/hts.h" nogil:
jpayne@68 402 uint32_t kroundup32(uint32_t x)
jpayne@68 403
jpayne@68 404 ctypedef struct cram_fd
jpayne@68 405
jpayne@68 406 union FilePointerUnion:
jpayne@68 407 BGZF *bgzf
jpayne@68 408 cram_fd *cram
jpayne@68 409 hFILE *hfile
jpayne@68 410 void *voidp
jpayne@68 411
jpayne@68 412 enum htsFormatCategory:
jpayne@68 413 unknown_category
jpayne@68 414 sequence_data # Sequence data -- SAM, BAM, CRAM, etc
jpayne@68 415 variant_data # Variant calling data -- VCF, BCF, etc
jpayne@68 416 index_file # Index file associated with some data file
jpayne@68 417 region_list # Coordinate intervals or regions -- BED, etc
jpayne@68 418 category_maximum
jpayne@68 419
jpayne@68 420 enum htsExactFormat:
jpayne@68 421 unknown_format
jpayne@68 422 binary_format
jpayne@68 423 text_format
jpayne@68 424 sam, bam, bai, cram, crai, vcf, bcf, csi, gzi, tbi, bed
jpayne@68 425 format_maximum
jpayne@68 426
jpayne@68 427 enum htsCompression:
jpayne@68 428 no_compression, gzip, bgzf, custom
jpayne@68 429 compression_maximum
jpayne@68 430
jpayne@68 431 cdef enum hts_fmt_option:
jpayne@68 432 CRAM_OPT_DECODE_MD,
jpayne@68 433 CRAM_OPT_PREFIX,
jpayne@68 434 CRAM_OPT_VERBOSITY,
jpayne@68 435 CRAM_OPT_SEQS_PER_SLICE,
jpayne@68 436 CRAM_OPT_SLICES_PER_CONTAINER,
jpayne@68 437 CRAM_OPT_RANGE,
jpayne@68 438 CRAM_OPT_VERSION,
jpayne@68 439 CRAM_OPT_EMBED_REF,
jpayne@68 440 CRAM_OPT_IGNORE_MD5,
jpayne@68 441 CRAM_OPT_REFERENCE,
jpayne@68 442 CRAM_OPT_MULTI_SEQ_PER_SLICE,
jpayne@68 443 CRAM_OPT_NO_REF,
jpayne@68 444 CRAM_OPT_USE_BZIP2,
jpayne@68 445 CRAM_OPT_SHARED_REF,
jpayne@68 446 CRAM_OPT_NTHREADS,
jpayne@68 447 CRAM_OPT_THREAD_POOL,
jpayne@68 448 CRAM_OPT_USE_LZMA,
jpayne@68 449 CRAM_OPT_USE_RANS,
jpayne@68 450 CRAM_OPT_REQUIRED_FIELDS,
jpayne@68 451 HTS_OPT_COMPRESSION_LEVEL,
jpayne@68 452 HTS_OPT_NTHREADS,
jpayne@68 453
jpayne@68 454 ctypedef struct htsVersion:
jpayne@68 455 short major, minor
jpayne@68 456
jpayne@68 457 ctypedef struct htsFormat:
jpayne@68 458 htsFormatCategory category
jpayne@68 459 htsExactFormat format
jpayne@68 460 htsVersion version
jpayne@68 461 htsCompression compression
jpayne@68 462 short compression_level
jpayne@68 463 void *specific
jpayne@68 464
jpayne@68 465 ctypedef struct htsFile:
jpayne@68 466 uint8_t is_bin
jpayne@68 467 uint8_t is_write
jpayne@68 468 uint8_t is_be
jpayne@68 469 uint8_t is_cram
jpayne@68 470 int64_t lineno
jpayne@68 471 kstring_t line
jpayne@68 472 char *fn
jpayne@68 473 char *fn_aux
jpayne@68 474 FilePointerUnion fp
jpayne@68 475 htsFormat format
jpayne@68 476
jpayne@68 477 int hts_verbose
jpayne@68 478
jpayne@68 479 cdef union hts_opt_val_union:
jpayne@68 480 int i
jpayne@68 481 char *s
jpayne@68 482
jpayne@68 483 ctypedef struct hts_opt:
jpayne@68 484 char *arg
jpayne@68 485 hts_fmt_option opt
jpayne@68 486 hts_opt_val_union val
jpayne@68 487 void *next
jpayne@68 488
jpayne@68 489 # @abstract Parses arg and appends it to the option list.
jpayne@68 490 # @return 0 on success and -1 on failure
jpayne@68 491 int hts_opt_add(hts_opt **opts, const char *c_arg)
jpayne@68 492
jpayne@68 493 # @abstract Applies an hts_opt option list to a given htsFile.
jpayne@68 494 # @return 0 on success and -1 on failure
jpayne@68 495 int hts_opt_apply(htsFile *fp, hts_opt *opts)
jpayne@68 496
jpayne@68 497 # @abstract Frees an hts_opt list.
jpayne@68 498 void hts_opt_free(hts_opt *opts)
jpayne@68 499
jpayne@68 500 # @abstract Table for converting a nucleotide character to 4-bit encoding.
jpayne@68 501 # The input character may be either an IUPAC ambiguity code, '=' for 0, or
jpayne@68 502 # '0'/'1'/'2'/'3' for a result of 1/2/4/8. The result is encoded as 1/2/4/8
jpayne@68 503 # for A/C/G/T or combinations of these bits for ambiguous bases.
jpayne@68 504 const unsigned char *seq_nt16_table
jpayne@68 505
jpayne@68 506 # @abstract Table for converting a 4-bit encoded nucleotide to an IUPAC
jpayne@68 507 # ambiguity code letter (or '=' when given 0).
jpayne@68 508 const char *seq_nt16_str
jpayne@68 509
jpayne@68 510 # @abstract Table for converting a 4-bit encoded nucleotide to about 2 bits.
jpayne@68 511 # Returns 0/1/2/3 for 1/2/4/8 (i.e., A/C/G/T), or 4 otherwise (0 or ambiguous).
jpayne@68 512 const int *seq_nt16_int
jpayne@68 513
jpayne@68 514 # @abstract Get the htslib version number
jpayne@68 515 # @return For released versions, a string like "N.N[.N]"; or git describe
jpayne@68 516 # output if using a library built within a Git repository.
jpayne@68 517 const char *hts_version()
jpayne@68 518
jpayne@68 519 # @abstract Determine format by peeking at the start of a file
jpayne@68 520 # @param fp File opened for reading, positioned at the beginning
jpayne@68 521 # @param fmt Format structure that will be filled out on return
jpayne@68 522 # @return 0 for success, or negative if an error occurred.
jpayne@68 523 int hts_detect_format(hFILE *fp, htsFormat *fmt)
jpayne@68 524
jpayne@68 525 # @abstract Get a human-readable description of the file format
jpayne@68 526 # @return Description string, to be freed by the caller after use.
jpayne@68 527 char *hts_format_description(const htsFormat *format)
jpayne@68 528
jpayne@68 529 # @abstract Open a SAM/BAM/CRAM/VCF/BCF/etc file
jpayne@68 530 # @param fn The file name or "-" for stdin/stdout
jpayne@68 531 # @param mode Mode matching / [rwa][bceguxz0-9]* /
jpayne@68 532 # @discussion
jpayne@68 533 # With 'r' opens for reading; any further format mode letters are ignored
jpayne@68 534 # as the format is detected by checking the first few bytes or BGZF blocks
jpayne@68 535 # of the file. With 'w' or 'a' opens for writing or appending, with format
jpayne@68 536 # specifier letters:
jpayne@68 537 # b binary format (BAM, BCF, etc) rather than text (SAM, VCF, etc)
jpayne@68 538 # c CRAM format
jpayne@68 539 # g gzip compressed
jpayne@68 540 # u uncompressed
jpayne@68 541 # z bgzf compressed
jpayne@68 542 # [0-9] zlib compression level
jpayne@68 543 # and with non-format option letters (for any of 'r'/'w'/'a'):
jpayne@68 544 # e close the file on exec(2) (opens with O_CLOEXEC, where supported)
jpayne@68 545 # x create the file exclusively (opens with O_EXCL, where supported)
jpayne@68 546 # Note that there is a distinction between 'u' and '0': the first yields
jpayne@68 547 # plain uncompressed output whereas the latter outputs uncompressed data
jpayne@68 548 # wrapped in the zlib format.
jpayne@68 549 # @example
jpayne@68 550 # [rw]b .. compressed BCF, BAM, FAI
jpayne@68 551 # [rw]bu .. uncompressed BCF
jpayne@68 552 # [rw]z .. compressed VCF
jpayne@68 553 # [rw] .. uncompressed VCF
jpayne@68 554 htsFile *hts_open(const char *fn, const char *mode)
jpayne@68 555
jpayne@68 556 # @abstract Open a SAM/BAM/CRAM/VCF/BCF/etc file
jpayne@68 557 # @param fn The file name or "-" for stdin/stdout
jpayne@68 558 # @param mode Open mode, as per hts_open()
jpayne@68 559 # @param fmt Optional format specific parameters
jpayne@68 560 # @discussion
jpayne@68 561 # See hts_open() for description of fn and mode.
jpayne@68 562 # // TODO Update documentation for s/opts/fmt/
jpayne@68 563 # Opts contains a format string (sam, bam, cram, vcf, bcf) which will,
jpayne@68 564 # if defined, override mode. Opts also contains a linked list of hts_opt
jpayne@68 565 # structures to apply to the open file handle. These can contain things
jpayne@68 566 # like pointers to the reference or information on compression levels,
jpayne@68 567 # block sizes, etc.
jpayne@68 568 htsFile *hts_open_format(const char *fn, const char *mode, const htsFormat *fmt)
jpayne@68 569
jpayne@68 570 # @abstract Open an existing stream as a SAM/BAM/CRAM/VCF/BCF/etc file
jpayne@68 571 # @param fp The already-open file handle
jpayne@68 572 # @param fn The file name or "-" for stdin/stdout
jpayne@68 573 # @param mode Open mode, as per hts_open()
jpayne@68 574 htsFile *hts_hopen(hFILE *fp, const char *fn, const char *mode)
jpayne@68 575
jpayne@68 576 # @abstract For output streams, flush any buffered data
jpayne@68 577 # @param fp The file handle to be flushed
jpayne@68 578 # @return 0 for success, or negative if an error occurred.
jpayne@68 579 # @since 1.14
jpayne@68 580 int hts_flush(htsFile *fp)
jpayne@68 581
jpayne@68 582 # @abstract Close a file handle, flushing buffered data for output streams
jpayne@68 583 # @param fp The file handle to be closed
jpayne@68 584 # @return 0 for success, or negative if an error occurred.
jpayne@68 585 int hts_close(htsFile *fp)
jpayne@68 586
jpayne@68 587 # @abstract Returns the file's format information
jpayne@68 588 # @param fp The file handle
jpayne@68 589 # @return Read-only pointer to the file's htsFormat.
jpayne@68 590 const htsFormat *hts_get_format(htsFile *fp)
jpayne@68 591
jpayne@68 592 # @ abstract Returns a string containing the file format extension.
jpayne@68 593 # @ param format Format structure containing the file type.
jpayne@68 594 # @ return A string ("sam", "bam", etc) or "?" for unknown formats.
jpayne@68 595 const char *hts_format_file_extension(const htsFormat *format)
jpayne@68 596
jpayne@68 597 # @abstract Sets a specified CRAM option on the open file handle.
jpayne@68 598 # @param fp The file handle open the open file.
jpayne@68 599 # @param opt The CRAM_OPT_* option.
jpayne@68 600 # @param ... Optional arguments, dependent on the option used.
jpayne@68 601 # @return 0 for success, or negative if an error occurred.
jpayne@68 602 int hts_set_opt(htsFile *fp, hts_fmt_option opt, ...)
jpayne@68 603
jpayne@68 604 int hts_getline(htsFile *fp, int delimiter, kstring_t *str)
jpayne@68 605 char **hts_readlines(const char *fn, int *_n)
jpayne@68 606
jpayne@68 607 # @abstract Parse comma-separated list or read list from a file
jpayne@68 608 # @param list File name or comma-separated list
jpayne@68 609 # @param is_file
jpayne@68 610 # @param _n Size of the output array (number of items read)
jpayne@68 611 # @return NULL on failure or pointer to newly allocated array of
jpayne@68 612 # strings
jpayne@68 613 char **hts_readlist(const char *fn, int is_file, int *_n)
jpayne@68 614
jpayne@68 615 # @abstract Create extra threads to aid compress/decompression for this file
jpayne@68 616 # @param fp The file handle
jpayne@68 617 # @param n The number of worker threads to create
jpayne@68 618 # @return 0 for success, or negative if an error occurred.
jpayne@68 619 # @notes THIS THREADING API IS LIKELY TO CHANGE IN FUTURE.
jpayne@68 620 int hts_set_threads(htsFile *fp, int n)
jpayne@68 621
jpayne@68 622 # @abstract Set .fai filename for a file opened for reading
jpayne@68 623 # @return 0 for success, negative on failure
jpayne@68 624 # @discussion
jpayne@68 625 # Called before *_hdr_read(), this provides the name of a .fai file
jpayne@68 626 # used to provide a reference list if the htsFile contains no @SQ headers.
jpayne@68 627 int hts_set_fai_filename(htsFile *fp, const char *fn_aux)
jpayne@68 628
jpayne@68 629 int8_t HTS_IDX_NOCOOR
jpayne@68 630 int8_t HTS_IDX_START
jpayne@68 631 int8_t HTS_IDX_REST
jpayne@68 632 int8_t HTS_IDX_NONE
jpayne@68 633
jpayne@68 634 int8_t HTS_FMT_CSI
jpayne@68 635 int8_t HTS_FMT_BAI
jpayne@68 636 int8_t HTS_FMT_TBI
jpayne@68 637 int8_t HTS_FMT_CRAI
jpayne@68 638
jpayne@68 639 BGZF *hts_get_bgzfp(htsFile *fp)
jpayne@68 640
jpayne@68 641 ctypedef struct hts_idx_t
jpayne@68 642
jpayne@68 643 ctypedef struct hts_pair64_t:
jpayne@68 644 uint64_t u, v
jpayne@68 645
jpayne@68 646 ctypedef int hts_readrec_func(BGZF *fp, void *data, void *r, int *tid, int *beg, int *end)
jpayne@68 647
jpayne@68 648 ctypedef struct hts_bins_t:
jpayne@68 649 int n, m
jpayne@68 650 int *a
jpayne@68 651
jpayne@68 652 ctypedef struct hts_itr_t:
jpayne@68 653 uint32_t read_rest
jpayne@68 654 uint32_t finished
jpayne@68 655 int tid, bed, end, n_off, i
jpayne@68 656 int curr_tid, curr_beg, curr_end
jpayne@68 657 uint64_t curr_off
jpayne@68 658 hts_pair64_t *off
jpayne@68 659 hts_readrec_func *readfunc
jpayne@68 660 hts_bins_t bins
jpayne@68 661
jpayne@68 662 hts_idx_t *hts_idx_init(int n, int fmt, uint64_t offset0, int min_shift, int n_lvls)
jpayne@68 663 void hts_idx_destroy(hts_idx_t *idx)
jpayne@68 664 int hts_idx_push(hts_idx_t *idx, int tid, int beg, int end, uint64_t offset, int is_mapped)
jpayne@68 665 void hts_idx_finish(hts_idx_t *idx, uint64_t final_offset)
jpayne@68 666
jpayne@68 667 #### Save an index to a file
jpayne@68 668 # @param idx Index to be written
jpayne@68 669 # @param fn Input BAM/BCF/etc filename, to which .bai/.csi/etc will be added
jpayne@68 670 # @param fmt One of the HTS_FMT_* index formats
jpayne@68 671 # @return 0 if successful, or negative if an error occurred.
jpayne@68 672 int hts_idx_save(const hts_idx_t *idx, const char *fn, int fmt)
jpayne@68 673
jpayne@68 674 #### Save an index to a specific file
jpayne@68 675 # @param idx Index to be written
jpayne@68 676 # @param fn Input BAM/BCF/etc filename
jpayne@68 677 # @param fnidx Output filename, or NULL to add .bai/.csi/etc to @a fn
jpayne@68 678 # @param fmt One of the HTS_FMT_* index formats
jpayne@68 679 # @return 0 if successful, or negative if an error occurred.
jpayne@68 680 int hts_idx_save_as(const hts_idx_t *idx, const char *fn, const char *fnidx, int fmt)
jpayne@68 681
jpayne@68 682 #### Load an index file
jpayne@68 683 # @param fn BAM/BCF/etc filename, to which .bai/.csi/etc will be added or
jpayne@68 684 # the extension substituted, to search for an existing index file
jpayne@68 685 # @param fmt One of the HTS_FMT_* index formats
jpayne@68 686 # @return The index, or NULL if an error occurred.
jpayne@68 687 hts_idx_t *hts_idx_load(const char *fn, int fmt)
jpayne@68 688
jpayne@68 689 #### Load a specific index file
jpayne@68 690 # @param fn Input BAM/BCF/etc filename
jpayne@68 691 # @param fnidx The input index filename
jpayne@68 692 # @return The index, or NULL if an error occurred.
jpayne@68 693 hts_idx_t *hts_idx_load2(const char *fn, const char *fnidx)
jpayne@68 694
jpayne@68 695 #### Load a specific index file
jpayne@68 696 # @param fn Input BAM/BCF/etc filename
jpayne@68 697 # @param fnidx The input index filename
jpayne@68 698 # @param fmt One of the HTS_FMT_* index formats
jpayne@68 699 # @param flags Flags to alter behaviour (see description)
jpayne@68 700 # @return The index, or NULL if an error occurred.
jpayne@68 701 hts_idx_t *hts_idx_load3(const char *fn, const char *fnidx, int fmt, int flags)
jpayne@68 702
jpayne@68 703 int HTS_IDX_SAVE_REMOTE
jpayne@68 704 int HTS_IDX_SILENT_FAIL
jpayne@68 705
jpayne@68 706 uint8_t *hts_idx_get_meta(hts_idx_t *idx, uint32_t *l_meta)
jpayne@68 707 void hts_idx_set_meta(hts_idx_t *idx, int l_meta, uint8_t *meta, int is_copy)
jpayne@68 708
jpayne@68 709 int hts_idx_get_stat(const hts_idx_t* idx, int tid,
jpayne@68 710 uint64_t* mapped, uint64_t* unmapped)
jpayne@68 711
jpayne@68 712 uint64_t hts_idx_get_n_no_coor(const hts_idx_t* idx)
jpayne@68 713
jpayne@68 714 int HTS_PARSE_THOUSANDS_SEP # Ignore ',' separators within numbers
jpayne@68 715
jpayne@68 716 # Parse a numeric string
jpayne@68 717 # The number may be expressed in scientific notation, and optionally may
jpayne@68 718 # contain commas in the integer part (before any decimal point or E notation).
jpayne@68 719 # @param str String to be parsed
jpayne@68 720 # @param strend If non-NULL, set on return to point to the first character
jpayne@68 721 # in @a str after those forming the parsed number
jpayne@68 722 # @param flags Or'ed-together combination of HTS_PARSE_* flags
jpayne@68 723 # @return Converted value of the parsed number.
jpayne@68 724 #
jpayne@68 725 # When @a strend is NULL, a warning will be printed (if hts_verbose is 2
jpayne@68 726 # or more) if there are any trailing characters after the number.
jpayne@68 727 long long hts_parse_decimal(const char *str, char **strend, int flags)
jpayne@68 728
jpayne@68 729 # Parse a "CHR:START-END"-style region string
jpayne@68 730 # @param str String to be parsed
jpayne@68 731 # @param beg Set on return to the 0-based start of the region
jpayne@68 732 # @param end Set on return to the 1-based end of the region
jpayne@68 733 # @return Pointer to the colon or '\0' after the reference sequence name,
jpayne@68 734 # or NULL if @a str could not be parsed.
jpayne@68 735 const char *hts_parse_reg(const char *str, int *beg, int *end)
jpayne@68 736
jpayne@68 737 hts_itr_t *hts_itr_query(const hts_idx_t *idx, int tid, int beg, int end, hts_readrec_func *readrec)
jpayne@68 738 void hts_itr_destroy(hts_itr_t *iter)
jpayne@68 739
jpayne@68 740 ctypedef int (*hts_name2id_f)(void*, const char*)
jpayne@68 741 ctypedef const char *(*hts_id2name_f)(void*, int)
jpayne@68 742 ctypedef hts_itr_t *hts_itr_query_func(
jpayne@68 743 const hts_idx_t *idx,
jpayne@68 744 int tid,
jpayne@68 745 int beg,
jpayne@68 746 int end,
jpayne@68 747 hts_readrec_func *readrec)
jpayne@68 748
jpayne@68 749 hts_itr_t *hts_itr_querys(
jpayne@68 750 const hts_idx_t *idx,
jpayne@68 751 const char *reg,
jpayne@68 752 hts_name2id_f getid,
jpayne@68 753 void *hdr,
jpayne@68 754 hts_itr_query_func *itr_query,
jpayne@68 755 hts_readrec_func *readrec)
jpayne@68 756
jpayne@68 757 int hts_itr_next(BGZF *fp, hts_itr_t *iter, void *r, void *data)
jpayne@68 758 const char **hts_idx_seqnames(const hts_idx_t *idx, int *n, hts_id2name_f getid, void *hdr) # free only the array, not the values
jpayne@68 759
jpayne@68 760 # hts_file_type() - Convenience function to determine file type
jpayne@68 761 # @fname: the file name
jpayne@68 762 #
jpayne@68 763 # Returns one of the FT_* defines.
jpayne@68 764 #
jpayne@68 765 # DEPRECATED: This function has been replaced by hts_detect_format().
jpayne@68 766 # It and these FT_* macros will be removed in a future HTSlib release.
jpayne@68 767 int FT_UNKN
jpayne@68 768 int FT_GZ
jpayne@68 769 int FT_VCF
jpayne@68 770 int FT_VCF_GZ
jpayne@68 771 int FT_BCF
jpayne@68 772 int FT_BCF_GZ
jpayne@68 773 int FT_STDIN
jpayne@68 774
jpayne@68 775 int hts_file_type(const char *fname)
jpayne@68 776
jpayne@68 777 # /***************************
jpayne@68 778 # * Revised MAQ error model *
jpayne@68 779 # ***************************/
jpayne@68 780
jpayne@68 781 ctypedef struct errmod_t
jpayne@68 782
jpayne@68 783 errmod_t *errmod_init(double depcorr)
jpayne@68 784 void errmod_destroy(errmod_t *em)
jpayne@68 785
jpayne@68 786 # /*
jpayne@68 787 # n: number of bases
jpayne@68 788 # m: maximum base
jpayne@68 789 # bases[i]: qual:6, strand:1, base:4
jpayne@68 790 # q[i*m+j]: phred-scaled likelihood of (i,j)
jpayne@68 791 # */
jpayne@68 792 int errmod_cal(const errmod_t *em, int n, int m, uint16_t *bases, float *Probabilistic)
jpayne@68 793
jpayne@68 794 # /*****************************************
jpayne@68 795 # * q banded glocal alignment *
jpayne@68 796 # *****************************************/
jpayne@68 797
jpayne@68 798 ctypedef struct probaln_par_t:
jpayne@68 799 float d, e
jpayne@68 800 int bw
jpayne@68 801
jpayne@68 802 int probaln_glocal(const uint8_t *ref,
jpayne@68 803 int l_ref,
jpayne@68 804 const uint8_t *query,
jpayne@68 805 int l_query, const uint8_t *iqual,
jpayne@68 806 const probaln_par_t *c,
jpayne@68 807 int *state, uint8_t *q)
jpayne@68 808
jpayne@68 809 # /**********************
jpayne@68 810 # * MD5 implementation *
jpayne@68 811 # **********************/
jpayne@68 812
jpayne@68 813 ctypedef struct hts_md5_context
jpayne@68 814
jpayne@68 815 # /*! @abstract Initialises an MD5 context.
jpayne@68 816 # * @discussion
jpayne@68 817 # * The expected use is to allocate an hts_md5_context using
jpayne@68 818 # * hts_md5_init(). This pointer is then passed into one or more calls
jpayne@68 819 # * of hts_md5_update() to compute successive internal portions of the
jpayne@68 820 # * MD5 sum, which can then be externalised as a full 16-byte MD5sum
jpayne@68 821 # * calculation by calling hts_md5_final(). This can then be turned
jpayne@68 822 # * into ASCII via hts_md5_hex().
jpayne@68 823 # *
jpayne@68 824 # * To dealloate any resources created by hts_md5_init() call the
jpayne@68 825 # * hts_md5_destroy() function.
jpayne@68 826 # *
jpayne@68 827 # * @return hts_md5_context pointer on success, NULL otherwise.
jpayne@68 828 # */
jpayne@68 829 hts_md5_context *hts_md5_init()
jpayne@68 830
jpayne@68 831 # /*! @abstract Updates the context with the MD5 of the data. */
jpayne@68 832 void hts_md5_update(hts_md5_context *ctx, const void *data, unsigned long size)
jpayne@68 833
jpayne@68 834 # /*! @abstract Computes the final 128-bit MD5 hash from the given context */
jpayne@68 835 void hts_md5_final(unsigned char *digest, hts_md5_context *ctx)
jpayne@68 836
jpayne@68 837 # /*! @abstract Resets an md5_context to the initial state, as returned
jpayne@68 838 # * by hts_md5_init().
jpayne@68 839 # */
jpayne@68 840 void hts_md5_reset(hts_md5_context *ctx)
jpayne@68 841
jpayne@68 842 # /*! @abstract Converts a 128-bit MD5 hash into a 33-byte nul-termninated
jpayne@68 843 # * hex string.
jpayne@68 844 # */
jpayne@68 845 void hts_md5_hex(char *hex, const unsigned char *digest)
jpayne@68 846
jpayne@68 847 # /*! @abstract Deallocates any memory allocated by hts_md5_init. */
jpayne@68 848 void hts_md5_destroy(hts_md5_context *ctx)
jpayne@68 849
jpayne@68 850 int hts_reg2bin(int64_t beg, int64_t end, int min_shift, int n_lvls)
jpayne@68 851 int hts_bin_bot(int bin, int n_lvls)
jpayne@68 852
jpayne@68 853 # * Endianness *
jpayne@68 854 int ed_is_big()
jpayne@68 855 uint16_t ed_swap_2(uint16_t v)
jpayne@68 856 void *ed_swap_2p(void *x)
jpayne@68 857 uint32_t ed_swap_4(uint32_t v)
jpayne@68 858 void *ed_swap_4p(void *x)
jpayne@68 859 uint64_t ed_swap_8(uint64_t v)
jpayne@68 860 void *ed_swap_8p(void *x)
jpayne@68 861
jpayne@68 862
jpayne@68 863 cdef extern from "htslib/sam.h" nogil:
jpayne@68 864 #**********************
jpayne@68 865 #*** SAM/BAM header ***
jpayne@68 866 #**********************
jpayne@68 867
jpayne@68 868 # @abstract Structure for the alignment header.
jpayne@68 869 # @field n_targets number of reference sequences
jpayne@68 870 # @field l_text length of the plain text in the header
jpayne@68 871 # @field target_len lengths of the reference sequences
jpayne@68 872 # @field target_name names of the reference sequences
jpayne@68 873 # @field text plain text
jpayne@68 874 # @field sdict header dictionary
jpayne@68 875
jpayne@68 876 ctypedef struct bam_hdr_t:
jpayne@68 877 int32_t n_targets, ignore_sam_err
jpayne@68 878 uint32_t l_text
jpayne@68 879 uint32_t *target_len
jpayne@68 880 uint8_t *cigar_tab
jpayne@68 881 char **target_name
jpayne@68 882 char *text
jpayne@68 883 void *sdict
jpayne@68 884
jpayne@68 885 #****************************
jpayne@68 886 #*** CIGAR related macros ***
jpayne@68 887 #****************************
jpayne@68 888
jpayne@68 889 int BAM_CMATCH
jpayne@68 890 int BAM_CINS
jpayne@68 891 int BAM_CDEL
jpayne@68 892 int BAM_CREF_SKIP
jpayne@68 893 int BAM_CSOFT_CLIP
jpayne@68 894 int BAM_CHARD_CLIP
jpayne@68 895 int BAM_CPAD
jpayne@68 896 int BAM_CEQUAL
jpayne@68 897 int BAM_CDIFF
jpayne@68 898 int BAM_CBACK
jpayne@68 899
jpayne@68 900 char *BAM_CIGAR_STR
jpayne@68 901 int BAM_CIGAR_SHIFT
jpayne@68 902 uint32_t BAM_CIGAR_MASK
jpayne@68 903 uint32_t BAM_CIGAR_TYPE
jpayne@68 904
jpayne@68 905 char bam_cigar_op(uint32_t c)
jpayne@68 906 uint32_t bam_cigar_oplen(uint32_t c)
jpayne@68 907 char bam_cigar_opchr(uint32_t)
jpayne@68 908 uint32_t bam_cigar_gen(char, uint32_t)
jpayne@68 909 int bam_cigar_type(char o)
jpayne@68 910
jpayne@68 911 # @abstract the read is paired in sequencing, no matter whether it is mapped in a pair
jpayne@68 912 int BAM_FPAIRED
jpayne@68 913 # @abstract the read is mapped in a proper pair
jpayne@68 914 int BAM_FPROPER_PAIR
jpayne@68 915 # @abstract the read itself is unmapped; conflictive with BAM_FPROPER_PAIR
jpayne@68 916 int BAM_FUNMAP
jpayne@68 917 # @abstract the mate is unmapped
jpayne@68 918 int BAM_FMUNMAP
jpayne@68 919 # @abstract the read is mapped to the reverse strand
jpayne@68 920 int BAM_FREVERSE
jpayne@68 921 # @abstract the mate is mapped to the reverse strand
jpayne@68 922 int BAM_FMREVERSE
jpayne@68 923 # @abstract this is read1
jpayne@68 924 int BAM_FREAD1
jpayne@68 925 # @abstract this is read2
jpayne@68 926 int BAM_FREAD2
jpayne@68 927 # @abstract not primary alignment
jpayne@68 928 int BAM_FSECONDARY
jpayne@68 929 # @abstract QC failure
jpayne@68 930 int BAM_FQCFAIL
jpayne@68 931 # @abstract optical or PCR duplicate
jpayne@68 932 int BAM_FDUP
jpayne@68 933 # @abstract supplementary alignment
jpayne@68 934 int BAM_FSUPPLEMENTARY
jpayne@68 935
jpayne@68 936 #*************************
jpayne@68 937 #*** Alignment records ***
jpayne@68 938 #*************************
jpayne@68 939
jpayne@68 940 # @abstract Structure for core alignment information.
jpayne@68 941 # @field tid chromosome ID, defined by bam_hdr_t
jpayne@68 942 # @field pos 0-based leftmost coordinate
jpayne@68 943 # @field bin bin calculated by bam_reg2bin()
jpayne@68 944 # @field qual mapping quality
jpayne@68 945 # @field l_qname length of the query name
jpayne@68 946 # @field flag bitwise flag
jpayne@68 947 # @field n_cigar number of CIGAR operations
jpayne@68 948 # @field l_qseq length of the query sequence (read)
jpayne@68 949 # @field mtid chromosome ID of next read in template, defined by bam_hdr_t
jpayne@68 950 # @field mpos 0-based leftmost coordinate of next read in template
jpayne@68 951
jpayne@68 952 ctypedef struct bam1_core_t:
jpayne@68 953 int32_t tid
jpayne@68 954 int32_t pos
jpayne@68 955 uint16_t bin
jpayne@68 956 uint8_t qual
jpayne@68 957 uint8_t l_qname
jpayne@68 958 uint16_t flag
jpayne@68 959 uint8_t unused1
jpayne@68 960 uint8_t l_extranul
jpayne@68 961 uint32_t n_cigar
jpayne@68 962 int32_t l_qseq
jpayne@68 963 int32_t mtid
jpayne@68 964 int32_t mpos
jpayne@68 965 int32_t isize
jpayne@68 966
jpayne@68 967 # @abstract Structure for one alignment.
jpayne@68 968 # @field core core information about the alignment
jpayne@68 969 # @field l_data current length of bam1_t::data
jpayne@68 970 # @field m_data maximum length of bam1_t::data
jpayne@68 971 # @field data all variable-length data, concatenated; structure: qname-cigar-seq-qual-aux
jpayne@68 972 #
jpayne@68 973 # @discussion Notes:
jpayne@68 974 #
jpayne@68 975 # 1. qname is zero tailing and core.l_qname includes the tailing '\0'.
jpayne@68 976 # 2. l_qseq is calculated from the total length of an alignment block
jpayne@68 977 # on reading or from CIGAR.
jpayne@68 978 # 3. cigar data is encoded 4 bytes per CIGAR operation.
jpayne@68 979 # 4. seq is nybble-encoded according to seq_nt16_table.
jpayne@68 980 ctypedef struct bam1_t:
jpayne@68 981 bam1_core_t core
jpayne@68 982 int l_data
jpayne@68 983 uint32_t m_data
jpayne@68 984 uint8_t *data
jpayne@68 985 uint64_t id
jpayne@68 986
jpayne@68 987 # @abstract Get whether the query is on the reverse strand
jpayne@68 988 # @param b pointer to an alignment
jpayne@68 989 # @return boolean true if query is on the reverse strand
jpayne@68 990 int bam_is_rev(bam1_t *b)
jpayne@68 991
jpayne@68 992 # @abstract Get whether the query's mate is on the reverse strand
jpayne@68 993 # @param b pointer to an alignment
jpayne@68 994 # @return boolean true if query's mate on the reverse strand
jpayne@68 995 int bam_is_mrev(bam1_t *b)
jpayne@68 996
jpayne@68 997 # @abstract Get the name of the query
jpayne@68 998 # @param b pointer to an alignment
jpayne@68 999 # @return pointer to the name string, null terminated
jpayne@68 1000 char *bam_get_qname(bam1_t *b)
jpayne@68 1001
jpayne@68 1002 # @abstract Get the CIGAR array
jpayne@68 1003 # @param b pointer to an alignment
jpayne@68 1004 # @return pointer to the CIGAR array
jpayne@68 1005 #
jpayne@68 1006 # @discussion In the CIGAR array, each element is a 32-bit integer. The
jpayne@68 1007 # lower 4 bits gives a CIGAR operation and the higher 28 bits keep the
jpayne@68 1008 # length of a CIGAR.
jpayne@68 1009 uint32_t *bam_get_cigar(bam1_t *b)
jpayne@68 1010
jpayne@68 1011 # @abstract Get query sequence
jpayne@68 1012 # @param b pointer to an alignment
jpayne@68 1013 # @return pointer to sequence
jpayne@68 1014 #
jpayne@68 1015 # @discussion Each base is encoded in 4 bits: 1 for A, 2 for C, 4 for G,
jpayne@68 1016 # 8 for T and 15 for N. Two bases are packed in one byte with the base
jpayne@68 1017 # at the higher 4 bits having smaller coordinate on the read. It is
jpayne@68 1018 # recommended to use bam_seqi() macro to get the base.
jpayne@68 1019 char *bam_get_seq(bam1_t *b)
jpayne@68 1020
jpayne@68 1021 # @abstract Get query quality
jpayne@68 1022 # @param b pointer to an alignment
jpayne@68 1023 # @return pointer to quality string
jpayne@68 1024 uint8_t *bam_get_qual(bam1_t *b)
jpayne@68 1025
jpayne@68 1026 # @abstract Get auxiliary data
jpayne@68 1027 # @param b pointer to an alignment
jpayne@68 1028 # @return pointer to the concatenated auxiliary data
jpayne@68 1029 uint8_t *bam_get_aux(bam1_t *b)
jpayne@68 1030
jpayne@68 1031 # @abstract Get length of auxiliary data
jpayne@68 1032 # @param b pointer to an alignment
jpayne@68 1033 # @return length of the concatenated auxiliary data
jpayne@68 1034 int bam_get_l_aux(bam1_t *b)
jpayne@68 1035
jpayne@68 1036 # @abstract Get a base on read
jpayne@68 1037 # @param s Query sequence returned by bam1_seq()
jpayne@68 1038 # @param i The i-th position, 0-based
jpayne@68 1039 # @return 4-bit integer representing the base.
jpayne@68 1040 char bam_seqi(char *s, int i)
jpayne@68 1041
jpayne@68 1042 #**************************
jpayne@68 1043 #*** Exported functions ***
jpayne@68 1044 #**************************
jpayne@68 1045
jpayne@68 1046 #***************
jpayne@68 1047 #*** BAM I/O ***
jpayne@68 1048 #***************
jpayne@68 1049
jpayne@68 1050 bam_hdr_t *bam_hdr_init()
jpayne@68 1051 bam_hdr_t *bam_hdr_read(BGZF *fp)
jpayne@68 1052 int bam_hdr_write(BGZF *fp, const bam_hdr_t *h)
jpayne@68 1053 void bam_hdr_destroy(bam_hdr_t *h)
jpayne@68 1054 int bam_name2id(bam_hdr_t *h, const char *ref)
jpayne@68 1055 bam_hdr_t* bam_hdr_dup(const bam_hdr_t *h0)
jpayne@68 1056
jpayne@68 1057 bam1_t *bam_init1()
jpayne@68 1058 void bam_destroy1(bam1_t *b)
jpayne@68 1059 int bam_read1(BGZF *fp, bam1_t *b)
jpayne@68 1060 int bam_write1(BGZF *fp, const bam1_t *b)
jpayne@68 1061 bam1_t *bam_copy1(bam1_t *bdst, const bam1_t *bsrc)
jpayne@68 1062 bam1_t *bam_dup1(const bam1_t *bsrc)
jpayne@68 1063
jpayne@68 1064 int bam_cigar2qlen(int n_cigar, const uint32_t *cigar)
jpayne@68 1065 int bam_cigar2rlen(int n_cigar, const uint32_t *cigar)
jpayne@68 1066
jpayne@68 1067 # @abstract Calculate the rightmost base position of an alignment on the
jpayne@68 1068 # reference genome.
jpayne@68 1069
jpayne@68 1070 # @param b pointer to an alignment
jpayne@68 1071 # @return the coordinate of the first base after the alignment, 0-based
jpayne@68 1072
jpayne@68 1073 # @discussion For a mapped read, this is just b->core.pos + bam_cigar2rlen.
jpayne@68 1074 # For an unmapped read (either according to its flags or if it has no cigar
jpayne@68 1075 # string), we return b->core.pos + 1 by convention.
jpayne@68 1076 int32_t bam_endpos(const bam1_t *b)
jpayne@68 1077
jpayne@68 1078 int bam_str2flag(const char *str) # returns negative value on error
jpayne@68 1079 char *bam_flag2str(int flag) # The string must be freed by the user
jpayne@68 1080
jpayne@68 1081 #*************************
jpayne@68 1082 #*** BAM/CRAM indexing ***
jpayne@68 1083 #*************************
jpayne@68 1084
jpayne@68 1085 # These BAM iterator functions work only on BAM files. To work with either
jpayne@68 1086 # BAM or CRAM files use the sam_index_load() & sam_itr_*() functions.
jpayne@68 1087 void bam_itr_destroy(hts_itr_t *iter)
jpayne@68 1088 hts_itr_t *bam_itr_queryi(const hts_idx_t *idx, int tid, int beg, int end)
jpayne@68 1089 hts_itr_t *bam_itr_querys(const hts_idx_t *idx, bam_hdr_t *hdr, const char *region)
jpayne@68 1090 int bam_itr_next(htsFile *htsfp, hts_itr_t *itr, void *r)
jpayne@68 1091
jpayne@68 1092 # Load/build .csi or .bai BAM index file. Does not work with CRAM.
jpayne@68 1093 # It is recommended to use the sam_index_* functions below instead.
jpayne@68 1094 hts_idx_t *bam_index_load(const char *fn)
jpayne@68 1095 int bam_index_build(const char *fn, int min_shift)
jpayne@68 1096
jpayne@68 1097 # Load a BAM (.csi or .bai) or CRAM (.crai) index file
jpayne@68 1098 # @param fp File handle of the data file whose index is being opened
jpayne@68 1099 # @param fn BAM/CRAM/etc filename to search alongside for the index file
jpayne@68 1100 # @return The index, or NULL if an error occurred.
jpayne@68 1101 hts_idx_t *sam_index_load(htsFile *fp, const char *fn)
jpayne@68 1102
jpayne@68 1103 # Load a specific BAM (.csi or .bai) or CRAM (.crai) index file
jpayne@68 1104 # @param fp File handle of the data file whose index is being opened
jpayne@68 1105 # @param fn BAM/CRAM/etc data file filename
jpayne@68 1106 # @param fnidx Index filename, or NULL to search alongside @a fn
jpayne@68 1107 # @return The index, or NULL if an error occurred.
jpayne@68 1108 hts_idx_t *sam_index_load2(htsFile *fp, const char *fn, const char *fnidx)
jpayne@68 1109
jpayne@68 1110 # Load or stream a BAM (.csi or .bai) or CRAM (.crai) index file
jpayne@68 1111 # @param fp File handle of the data file whose index is being opened
jpayne@68 1112 # @param fn BAM/CRAM/etc data file filename
jpayne@68 1113 # @param fnidx Index filename, or NULL to search alongside @a fn
jpayne@68 1114 # @param flags Flags to alter behaviour
jpayne@68 1115 # @return The index, or NULL if an error occurred.
jpayne@68 1116 hts_idx_t *sam_index_load3(htsFile *fp, const char *fn, const char *fnidx, int flags)
jpayne@68 1117
jpayne@68 1118 # Generate and save an index file
jpayne@68 1119 # @param fn Input BAM/etc filename, to which .csi/etc will be added
jpayne@68 1120 # @param min_shift Positive to generate CSI, or 0 to generate BAI
jpayne@68 1121 # @return 0 if successful, or negative if an error occurred (usually -1; or
jpayne@68 1122 # -2: opening fn failed; -3: format not indexable)
jpayne@68 1123 int sam_index_build(const char *fn, int min_shift)
jpayne@68 1124
jpayne@68 1125 # Generate and save an index to a specific file
jpayne@68 1126 # @param fn Input BAM/CRAM/etc filename
jpayne@68 1127 # @param fnidx Output filename, or NULL to add .bai/.csi/etc to @a fn
jpayne@68 1128 # @param min_shift Positive to generate CSI, or 0 to generate BAI
jpayne@68 1129 # @return 0 if successful, or negative if an error occurred.
jpayne@68 1130 int sam_index_build2(const char *fn, const char *fnidx, int min_shift)
jpayne@68 1131
jpayne@68 1132 void sam_itr_destroy(hts_itr_t *iter)
jpayne@68 1133 hts_itr_t *sam_itr_queryi(const hts_idx_t *idx, int tid, int beg, int end)
jpayne@68 1134 hts_itr_t *sam_itr_querys(const hts_idx_t *idx, bam_hdr_t *hdr, const char *region)
jpayne@68 1135 int sam_itr_next(htsFile *htsfp, hts_itr_t *itr, void *r)
jpayne@68 1136
jpayne@68 1137 #***************
jpayne@68 1138 #*** SAM I/O ***
jpayne@68 1139 #***************
jpayne@68 1140
jpayne@68 1141 htsFile *sam_open(const char *fn, const char *mode)
jpayne@68 1142 htsFile *sam_open_format(const char *fn, const char *mode, const htsFormat *fmt)
jpayne@68 1143 int sam_close(htsFile *fp)
jpayne@68 1144
jpayne@68 1145 int sam_open_mode(char *mode, const char *fn, const char *format)
jpayne@68 1146
jpayne@68 1147 # A version of sam_open_mode that can handle ,key=value options.
jpayne@68 1148 # The format string is allocated and returned, to be freed by the caller.
jpayne@68 1149 # Prefix should be "r" or "w",
jpayne@68 1150 char *sam_open_mode_opts(const char *fn, const char *mode, const char *format)
jpayne@68 1151
jpayne@68 1152 bam_hdr_t *sam_hdr_parse(int l_text, const char *text)
jpayne@68 1153 bam_hdr_t *sam_hdr_read(htsFile *fp)
jpayne@68 1154 int sam_hdr_write(htsFile *fp, const bam_hdr_t *h)
jpayne@68 1155
jpayne@68 1156 int sam_parse1(kstring_t *s, bam_hdr_t *h, bam1_t *b)
jpayne@68 1157 int sam_format1(const bam_hdr_t *h, const bam1_t *b, kstring_t *str)
jpayne@68 1158 int sam_read1(htsFile *fp, bam_hdr_t *h, bam1_t *b)
jpayne@68 1159 int sam_write1(htsFile *fp, const bam_hdr_t *h, const bam1_t *b)
jpayne@68 1160
jpayne@68 1161 #*************************************
jpayne@68 1162 #*** Manipulating auxiliary fields ***
jpayne@68 1163 #*************************************
jpayne@68 1164
jpayne@68 1165 uint8_t *bam_aux_get(const bam1_t *b, const char *tag)
jpayne@68 1166 int64_t bam_aux2i(const uint8_t *s)
jpayne@68 1167 double bam_aux2f(const uint8_t *s)
jpayne@68 1168 char bam_aux2A(const uint8_t *s)
jpayne@68 1169 char *bam_aux2Z(const uint8_t *s)
jpayne@68 1170
jpayne@68 1171 void bam_aux_append(bam1_t *b, const char *tag, char type, int len, uint8_t *data)
jpayne@68 1172 int bam_aux_del(bam1_t *b, uint8_t *s)
jpayne@68 1173
jpayne@68 1174 #**************************
jpayne@68 1175 #*** Pileup and Mpileup ***
jpayne@68 1176 #**************************
jpayne@68 1177
jpayne@68 1178 # @abstract Generic pileup 'client data'.
jpayne@68 1179 # @discussion The pileup iterator allows setting a constructor and
jpayne@68 1180 # destructor function, which will be called every time a sequence is
jpayne@68 1181 # fetched and discarded. This permits caching of per-sequence data in
jpayne@68 1182 # a tidy manner during the pileup process. This union is the cached
jpayne@68 1183 # data to be manipulated by the "client" (the caller of pileup).
jpayne@68 1184 #
jpayne@68 1185 union bam_pileup_cd:
jpayne@68 1186 void *p
jpayne@68 1187 int64_t i
jpayne@68 1188 double f
jpayne@68 1189
jpayne@68 1190 # @abstract Structure for one alignment covering the pileup position.
jpayne@68 1191 # @field b pointer to the alignment
jpayne@68 1192 # @field qpos position of the read base at the pileup site, 0-based
jpayne@68 1193 # @field indel indel length; 0 for no indel, positive for ins and negative for del
jpayne@68 1194 # @field level the level of the read in the "viewer" mode
jpayne@68 1195 # @field is_del 1 iff the base on the padded read is a deletion
jpayne@68 1196 # @field is_head ???
jpayne@68 1197 # @field is_tail ???
jpayne@68 1198 # @field is_refskip ???
jpayne@68 1199 # @field aux ???
jpayne@68 1200 #
jpayne@68 1201 # @discussion See also bam_plbuf_push() and bam_lplbuf_push(). The
jpayne@68 1202 # difference between the two functions is that the former does not
jpayne@68 1203 # set bam_pileup1_t::level, while the later does. Level helps the
jpayne@68 1204 # implementation of alignment viewers, but calculating this has some
jpayne@68 1205 # overhead.
jpayne@68 1206 #
jpayne@68 1207 # is_del, is_head, etc are a bit field, declaring as below should
jpayne@68 1208 # work as expected, see
jpayne@68 1209 # https://groups.google.com/forum/#!msg/cython-users/24tD1kwRY7A/pmoPuSmanM0J
jpayne@68 1210
jpayne@68 1211 ctypedef struct bam_pileup1_t:
jpayne@68 1212 bam1_t *b
jpayne@68 1213 int32_t qpos
jpayne@68 1214 int indel, level
jpayne@68 1215 uint32_t is_del
jpayne@68 1216 uint32_t is_head
jpayne@68 1217 uint32_t is_tail
jpayne@68 1218 uint32_t is_refskip
jpayne@68 1219 uint32_t aux
jpayne@68 1220 bam_pileup_cd cd
jpayne@68 1221
jpayne@68 1222 ctypedef int (*bam_plp_auto_f)(void *data, bam1_t *b)
jpayne@68 1223 ctypedef int (*bam_test_f)()
jpayne@68 1224
jpayne@68 1225 ctypedef struct __bam_plp_t
jpayne@68 1226 ctypedef __bam_plp_t *bam_plp_t
jpayne@68 1227
jpayne@68 1228 ctypedef struct __bam_mplp_t
jpayne@68 1229 ctypedef __bam_mplp_t *bam_mplp_t
jpayne@68 1230
jpayne@68 1231 # bam_plp_init() - sets an iterator over multiple
jpayne@68 1232 # @func: see mplp_func in bam_plcmd.c in samtools for an example. Expected return
jpayne@68 1233 # status: 0 on success, -1 on end, < -1 on non-recoverable errors
jpayne@68 1234 # @data: user data to pass to @func
jpayne@68 1235 bam_plp_t bam_plp_init(bam_plp_auto_f func, void *data)
jpayne@68 1236 void bam_plp_destroy(bam_plp_t iter)
jpayne@68 1237 int bam_plp_push(bam_plp_t iter, const bam1_t *b)
jpayne@68 1238 const bam_pileup1_t *bam_plp_next(bam_plp_t iter, int *_tid, int *_pos, int *_n_plp)
jpayne@68 1239 const bam_pileup1_t *bam_plp_auto(bam_plp_t iter, int *_tid, int *_pos, int *_n_plp)
jpayne@68 1240 void bam_plp_set_maxcnt(bam_plp_t iter, int maxcnt)
jpayne@68 1241 void bam_plp_reset(bam_plp_t iter)
jpayne@68 1242
jpayne@68 1243 bam_mplp_t bam_mplp_init(int n, bam_plp_auto_f func, void **data)
jpayne@68 1244
jpayne@68 1245 # bam_mplp_init_overlaps() - if called, mpileup will detect overlapping
jpayne@68 1246 # read pairs and for each base pair set the base quality of the
jpayne@68 1247 # lower-quality base to zero, thus effectively discarding it from
jpayne@68 1248 # calling. If the two bases are identical, the quality of the other base
jpayne@68 1249 # is increased to the sum of their qualities (capped at 200), otherwise
jpayne@68 1250 # it is multiplied by 0.8.
jpayne@68 1251 void bam_mplp_init_overlaps(bam_mplp_t iter)
jpayne@68 1252 void bam_mplp_destroy(bam_mplp_t iter)
jpayne@68 1253 void bam_mplp_set_maxcnt(bam_mplp_t iter, int maxcnt)
jpayne@68 1254 int bam_mplp_auto(bam_mplp_t iter, int *_tid, int *_pos, int *n_plp, const bam_pileup1_t **plp)
jpayne@68 1255 void bam_mplp_reset(bam_mplp_t iter)
jpayne@68 1256 void bam_mplp_constructor(bam_mplp_t iter,
jpayne@68 1257 int (*func)(void *data, const bam1_t *b, bam_pileup_cd *cd))
jpayne@68 1258 void bam_mplp_destructor(bam_mplp_t iter,
jpayne@68 1259 int (*func)(void *data, const bam1_t *b, bam_pileup_cd *cd))
jpayne@68 1260
jpayne@68 1261 # Added by AH
jpayne@68 1262 # ctypedef bam_pileup1_t * const_bam_pileup1_t_ptr "const bam_pileup1_t *"
jpayne@68 1263
jpayne@68 1264
jpayne@68 1265
jpayne@68 1266
jpayne@68 1267 # // ---------------------------
jpayne@68 1268 # // Base modification retrieval
jpayne@68 1269
jpayne@68 1270 # /*! @typedef
jpayne@68 1271 # @abstract Holds a single base modification.
jpayne@68 1272 # @field modified_base The short base code (m, h, etc) or -ChEBI (negative)
jpayne@68 1273 # @field canonical_base The canonical base referred to in the MM tag.
jpayne@68 1274 # One of A, C, G, T or N. Note this may not be the
jpayne@68 1275 # explicit base recorded in the SEQ column (esp. if N).
jpayne@68 1276 # @field strand 0 or 1, indicating + or - strand from MM tag.
jpayne@68 1277 # @field qual Quality code (256*probability), or -1 if unknown
jpayne@68 1278
jpayne@68 1279 # @discussion
jpayne@68 1280 # Note this doesn't hold any location data or information on which other
jpayne@68 1281 # modifications may be possible at this site.
jpayne@68 1282 ctypedef struct hts_base_mod:
jpayne@68 1283 int modified_base
jpayne@68 1284 int canonical_base
jpayne@68 1285 int strand
jpayne@68 1286 int qual
jpayne@68 1287
jpayne@68 1288 # /// Allocates an hts_base_mode_state.
jpayne@68 1289 # /**
jpayne@68 1290 # * @return An hts_base_mode_state pointer on success,
jpayne@68 1291 # * NULL on failure.
jpayne@68 1292 # *
jpayne@68 1293 # * This just allocates the memory. The initialisation of the contents is
jpayne@68 1294 # * done using bam_parse_basemod. Successive calls may be made to that
jpayne@68 1295 # * without the need to free and allocate a new state.
jpayne@68 1296 # *
jpayne@68 1297 # * The state be destroyed using the hts_base_mode_state_free function.
jpayne@68 1298 # */
jpayne@68 1299 ctypedef struct hts_base_mod_state
jpayne@68 1300 hts_base_mod_state *hts_base_mod_state_alloc()
jpayne@68 1301
jpayne@68 1302
jpayne@68 1303 # /// Destroys an hts_base_mode_state.
jpayne@68 1304 # /**
jpayne@68 1305 # * @param state The base modification state pointer.
jpayne@68 1306 # *
jpayne@68 1307 # * The should have previously been created by hts_base_mode_state_alloc.
jpayne@68 1308 # */
jpayne@68 1309 void hts_base_mod_state_free(hts_base_mod_state *state)
jpayne@68 1310
jpayne@68 1311 # /// Parses the Mm and Ml tags out of a bam record.
jpayne@68 1312 # /**
jpayne@68 1313 # * @param b BAM alignment record
jpayne@68 1314 # * @param state The base modification state pointer.
jpayne@68 1315 # * @return 0 on success,
jpayne@68 1316 # * -1 on failure.
jpayne@68 1317 # *
jpayne@68 1318 # * This fills out the contents of the modification state, resetting the
jpayne@68 1319 # * iterator location to the first sequence base.
jpayne@68 1320 # */
jpayne@68 1321 int bam_parse_basemod(const bam1_t *b, hts_base_mod_state *state)
jpayne@68 1322
jpayne@68 1323 # /// Finds the next location containing base modifications and returns them
jpayne@68 1324 # /**
jpayne@68 1325 # * @param b BAM alignment record
jpayne@68 1326 # * @param state The base modification state pointer.
jpayne@68 1327 # * @param mods A supplied array for returning base modifications
jpayne@68 1328 # * @param n_mods The size of the mods array
jpayne@68 1329 # * @return The number of modifications found on success,
jpayne@68 1330 # * 0 if no more modifications are present,
jpayne@68 1331 # * -1 on failure.
jpayne@68 1332 # *
jpayne@68 1333 # * Unlike bam_mods_at_next_pos this skips ahead to the next site
jpayne@68 1334 # * with modifications.
jpayne@68 1335 # *
jpayne@68 1336 # * If more than n_mods modifications are found, the total found is returned.
jpayne@68 1337 # * Note this means the caller needs to check whether this is higher than
jpayne@68 1338 # * n_mods.
jpayne@68 1339 # */
jpayne@68 1340
jpayne@68 1341 int bam_next_basemod(const bam1_t *b, hts_base_mod_state *state,hts_base_mod *mods, int n_mods, int *pos)
jpayne@68 1342
jpayne@68 1343 # ***********************************
jpayne@68 1344 # * BAQ calculation and realignment *
jpayne@68 1345 # ***********************************/
jpayne@68 1346 int sam_cap_mapq(bam1_t *b, const char *ref, int ref_len, int thres)
jpayne@68 1347 int sam_prob_realn(bam1_t *b, const char *ref, int ref_len, int flag)
jpayne@68 1348
jpayne@68 1349
jpayne@68 1350 cdef extern from "htslib/faidx.h" nogil:
jpayne@68 1351
jpayne@68 1352 ctypedef struct faidx_t:
jpayne@68 1353 pass
jpayne@68 1354
jpayne@68 1355 # /// Build index for a FASTA or bgzip-compressed FASTA file.
jpayne@68 1356 # /** @param fn FASTA file name
jpayne@68 1357 # @param fnfai Name of .fai file to build.
jpayne@68 1358 # @param fngzi Name of .gzi file to build (if fn is bgzip-compressed).
jpayne@68 1359 # @return 0 on success; or -1 on failure
jpayne@68 1360
jpayne@68 1361 # If fnfai is NULL, ".fai" will be appended to fn to make the FAI file name.
jpayne@68 1362 # If fngzi is NULL, ".gzi" will be appended to fn for the GZI file. The GZI
jpayne@68 1363 # file will only be built if fn is bgzip-compressed.
jpayne@68 1364 # */
jpayne@68 1365 int fai_build3(const char *fn,
jpayne@68 1366 const char *fnfai,
jpayne@68 1367 const char *fngzi)
jpayne@68 1368
jpayne@68 1369 # /// Build index for a FASTA or bgzip-compressed FASTA file.
jpayne@68 1370 # /** @param fn FASTA file name
jpayne@68 1371 # @return 0 on success; or -1 on failure
jpayne@68 1372 #
jpayne@68 1373 # File "fn.fai" will be generated. This function is equivalent to
jpayne@68 1374 # fai_build3(fn, NULL, NULL);
jpayne@68 1375 # */
jpayne@68 1376 int fai_build(char *fn)
jpayne@68 1377
jpayne@68 1378 # /// Destroy a faidx_t struct
jpayne@68 1379 void fai_destroy(faidx_t *fai)
jpayne@68 1380
jpayne@68 1381 # /// Load FASTA indexes.
jpayne@68 1382 # /** @param fn File name of the FASTA file (can be compressed with bgzip).
jpayne@68 1383 # @param fnfai File name of the FASTA index.
jpayne@68 1384 # @param fngzi File name of the bgzip index.
jpayne@68 1385 # @param flags Option flags to control index file caching and creation.
jpayne@68 1386 # @return Pointer to a faidx_t struct on success, NULL on failure.
jpayne@68 1387
jpayne@68 1388 # If fnfai is NULL, ".fai" will be appended to fn to make the FAI file name.
jpayne@68 1389 # If fngzi is NULL, ".gzi" will be appended to fn for the bgzip index name.
jpayne@68 1390 # The bgzip index is only needed if fn is compressed.
jpayne@68 1391
jpayne@68 1392 # If (flags & FAI_CREATE) is true, the index files will be built using
jpayne@68 1393 # fai_build3() if they are not already present.
jpayne@68 1394 # */
jpayne@68 1395 faidx_t *fai_load3(const char *fn,
jpayne@68 1396 const char *fnfai,
jpayne@68 1397 const char *fngzi,
jpayne@68 1398 int flags)
jpayne@68 1399
jpayne@68 1400 # /// Load index from "fn.fai".
jpayne@68 1401 # /** @param fn File name of the FASTA file
jpayne@68 1402 # @return Pointer to a faidx_t struct on success, NULL on failure.
jpayne@68 1403 # This function is equivalent to fai_load3(fn, NULL, NULL, FAI_CREATE|FAI_CACHE);
jpayne@68 1404 # */
jpayne@68 1405 faidx_t *fai_load(char *fn)
jpayne@68 1406
jpayne@68 1407 # /// Fetch the sequence in a region
jpayne@68 1408 # /** @param fai Pointer to the faidx_t struct
jpayne@68 1409 # @param reg Region in the format "chr2:20,000-30,000"
jpayne@68 1410 # @param len Length of the region; -2 if seq not present, -1 general error
jpayne@68 1411 # @return Pointer to the sequence; `NULL` on failure
jpayne@68 1412 # The returned sequence is allocated by `malloc()` family and should be destroyed
jpayne@68 1413 # by end users by calling `free()` on it.
jpayne@68 1414 # */
jpayne@68 1415 char *fai_fetch(faidx_t *fai,
jpayne@68 1416 char *reg,
jpayne@68 1417 int *len)
jpayne@68 1418
jpayne@68 1419 # /// Fetch the sequence in a region
jpayne@68 1420 # /** @param fai Pointer to the faidx_t struct
jpayne@68 1421 # @param c_name Region name
jpayne@68 1422 # @param p_beg_i Beginning position number (zero-based)
jpayne@68 1423 # @param p_end_i End position number (zero-based)
jpayne@68 1424 # @param len Length of the region; -2 if c_name not present, -1 general error
jpayne@68 1425 # @return Pointer to the sequence; null on failure
jpayne@68 1426 # The returned sequence is allocated by `malloc()` family and should be destroyed
jpayne@68 1427 # by end users by calling `free()` on it.
jpayne@68 1428 # */
jpayne@68 1429 char *faidx_fetch_seq(faidx_t *fai,
jpayne@68 1430 char *c_name,
jpayne@68 1431 int p_beg_i,
jpayne@68 1432 int p_end_i,
jpayne@68 1433 int *len)
jpayne@68 1434
jpayne@68 1435 # /// Query if sequence is present
jpayne@68 1436 # /** @param fai Pointer to the faidx_t struct
jpayne@68 1437 # @param seq Sequence name
jpayne@68 1438 # @return 1 if present or 0 if absent
jpayne@68 1439 # */
jpayne@68 1440 int faidx_has_seq(faidx_t *fai, const char *seq)
jpayne@68 1441
jpayne@68 1442 # /// Fetch the number of sequences
jpayne@68 1443 # /** @param fai Pointer to the faidx_t struct
jpayne@68 1444 # @return The number of sequences
jpayne@68 1445 # */
jpayne@68 1446 int faidx_nseq(const faidx_t *fai)
jpayne@68 1447
jpayne@68 1448 # /// Return name of i-th sequence
jpayne@68 1449 const char *faidx_iseq(const faidx_t *fai, int i)
jpayne@68 1450
jpayne@68 1451 # /// Return sequence length, -1 if not present
jpayne@68 1452 int faidx_seq_len(faidx_t *fai, const char *seq)
jpayne@68 1453
jpayne@68 1454 # tabix support
jpayne@68 1455 cdef extern from "htslib/tbx.h" nogil:
jpayne@68 1456
jpayne@68 1457 # tbx.h definitions
jpayne@68 1458 int8_t TBX_MAX_SHIFT
jpayne@68 1459 int32_t TBX_GENERIC
jpayne@68 1460 int32_t TBX_SAM
jpayne@68 1461 int32_t TBX_VCF
jpayne@68 1462 int32_t TBX_UCSC
jpayne@68 1463
jpayne@68 1464 ctypedef struct tbx_conf_t:
jpayne@68 1465 int32_t preset
jpayne@68 1466 int32_t sc, bc, ec # seq col., beg col. and end col.
jpayne@68 1467 int32_t meta_char, line_skip
jpayne@68 1468
jpayne@68 1469 ctypedef struct tbx_t:
jpayne@68 1470 tbx_conf_t conf
jpayne@68 1471 hts_idx_t *idx
jpayne@68 1472 void * dict
jpayne@68 1473
jpayne@68 1474 tbx_conf_t tbx_conf_gff
jpayne@68 1475 tbx_conf_t tbx_conf_bed
jpayne@68 1476 tbx_conf_t tbx_conf_psltbl
jpayne@68 1477 tbx_conf_t tbx_conf_sam
jpayne@68 1478 tbx_conf_t tbx_conf_vcf
jpayne@68 1479
jpayne@68 1480 void tbx_itr_destroy(hts_itr_t * iter)
jpayne@68 1481 hts_itr_t * tbx_itr_queryi(tbx_t * t, int tid, int bed, int end)
jpayne@68 1482 hts_itr_t * tbx_itr_querys(tbx_t * t, char * s)
jpayne@68 1483 int tbx_itr_next(htsFile * fp, tbx_t * t, hts_itr_t * iter, void * data)
jpayne@68 1484
jpayne@68 1485 int tbx_name2id(tbx_t *tbx, char *ss)
jpayne@68 1486
jpayne@68 1487 int tbx_index_build(char *fn, int min_shift, tbx_conf_t *conf)
jpayne@68 1488 int tbx_index_build2(const char *fn, const char *fnidx, int min_shift, const tbx_conf_t *conf)
jpayne@68 1489
jpayne@68 1490 tbx_t * tbx_index_load(char *fn)
jpayne@68 1491 tbx_t *tbx_index_load2(const char *fn, const char *fnidx)
jpayne@68 1492 tbx_t *tbx_index_load3(const char *fn, const char *fnidx, int flags)
jpayne@68 1493
jpayne@68 1494 # free the array but not the values
jpayne@68 1495 char **tbx_seqnames(tbx_t *tbx, int *n)
jpayne@68 1496
jpayne@68 1497 void tbx_destroy(tbx_t *tbx)
jpayne@68 1498
jpayne@68 1499
jpayne@68 1500 # VCF/BCF API
jpayne@68 1501 cdef extern from "htslib/vcf.h" nogil:
jpayne@68 1502
jpayne@68 1503 # Header struct
jpayne@68 1504
jpayne@68 1505 uint8_t BCF_HL_FLT # header line
jpayne@68 1506 uint8_t BCF_HL_INFO
jpayne@68 1507 uint8_t BCF_HL_FMT
jpayne@68 1508 uint8_t BCF_HL_CTG
jpayne@68 1509 uint8_t BCF_HL_STR # structured header line TAG=<A=..,B=..>
jpayne@68 1510 uint8_t BCF_HL_GEN # generic header line
jpayne@68 1511
jpayne@68 1512 uint8_t BCF_HT_FLAG # header type
jpayne@68 1513 uint8_t BCF_HT_INT
jpayne@68 1514 uint8_t BCF_HT_REAL
jpayne@68 1515 uint8_t BCF_HT_STR
jpayne@68 1516
jpayne@68 1517 uint8_t BCF_VL_FIXED # variable length
jpayne@68 1518 uint8_t BCF_VL_VAR
jpayne@68 1519 uint8_t BCF_VL_A
jpayne@68 1520 uint8_t BCF_VL_G
jpayne@68 1521 uint8_t BCF_VL_R
jpayne@68 1522
jpayne@68 1523 # === Dictionary ===
jpayne@68 1524 #
jpayne@68 1525 # The header keeps three dictionaries. The first keeps IDs in the
jpayne@68 1526 # "FILTER/INFO/FORMAT" lines, the second keeps the sequence names and lengths
jpayne@68 1527 # in the "contig" lines and the last keeps the sample names. bcf_hdr_t::dict[]
jpayne@68 1528 # is the actual hash table, which is opaque to the end users. In the hash
jpayne@68 1529 # table, the key is the ID or sample name as a C string and the value is a
jpayne@68 1530 # bcf_idinfo_t struct. bcf_hdr_t::id[] points to key-value pairs in the hash
jpayne@68 1531 # table in the order that they appear in the VCF header. bcf_hdr_t::n[] is the
jpayne@68 1532 # size of the hash table or, equivalently, the length of the id[] arrays.
jpayne@68 1533
jpayne@68 1534 uint8_t BCF_DT_ID # dictionary type
jpayne@68 1535 uint8_t BCF_DT_CTG
jpayne@68 1536 uint8_t BCF_DT_SAMPLE
jpayne@68 1537
jpayne@68 1538 # Complete textual representation of a header line
jpayne@68 1539 ctypedef struct bcf_hrec_t:
jpayne@68 1540 int type # One of the BCF_HL_* type
jpayne@68 1541 char *key # The part before '=', i.e. FILTER/INFO/FORMAT/contig/fileformat etc.
jpayne@68 1542 char *value # Set only for generic lines, NULL for FILTER/INFO, etc.
jpayne@68 1543 int nkeys # Number of structured fields
jpayne@68 1544 char **keys # The key=value pairs
jpayne@68 1545 char **vals
jpayne@68 1546
jpayne@68 1547 ctypedef struct bcf_idinfo_t:
jpayne@68 1548 uint32_t info[3] # stores Number:20, var:4, Type:4, ColType:4 in info[0..2]
jpayne@68 1549 bcf_hrec_t *hrec[3] # for BCF_HL_FLT,INFO,FMT and contig length in info[0] for BCF_HL_CTG
jpayne@68 1550 int id
jpayne@68 1551
jpayne@68 1552 ctypedef struct bcf_idpair_t:
jpayne@68 1553 const char *key
jpayne@68 1554 const bcf_idinfo_t *val
jpayne@68 1555
jpayne@68 1556 ctypedef struct bcf_hdr_t:
jpayne@68 1557 int32_t n[3] # n:the size of the dictionary block in use, (allocated size, m, is below to preserve ABI)
jpayne@68 1558 bcf_idpair_t *id[3]
jpayne@68 1559 void *dict[3] # ID dictionary, contig dict and sample dict
jpayne@68 1560 char **samples
jpayne@68 1561 bcf_hrec_t **hrec
jpayne@68 1562 int nhrec, dirty
jpayne@68 1563 int ntransl
jpayne@68 1564 int *transl[2] # for bcf_translate()
jpayne@68 1565 int nsamples_ori # for bcf_hdr_set_samples()
jpayne@68 1566 uint8_t *keep_samples
jpayne@68 1567 kstring_t mem
jpayne@68 1568 int32_t m[3] # m: allocated size of the dictionary block in use (see n above)
jpayne@68 1569
jpayne@68 1570 uint8_t bcf_type_shift[]
jpayne@68 1571
jpayne@68 1572 # * VCF record *
jpayne@68 1573
jpayne@68 1574 uint8_t BCF_BT_NULL
jpayne@68 1575 uint8_t BCF_BT_INT8
jpayne@68 1576 uint8_t BCF_BT_INT16
jpayne@68 1577 uint8_t BCF_BT_INT32
jpayne@68 1578 uint8_t BCF_BT_FLOAT
jpayne@68 1579 uint8_t BCF_BT_CHAR
jpayne@68 1580
jpayne@68 1581 uint8_t VCF_REF
jpayne@68 1582 uint8_t VCF_SNP
jpayne@68 1583 uint8_t VCF_MNP
jpayne@68 1584 uint8_t VCF_INDEL
jpayne@68 1585 uint8_t VCF_OTHER
jpayne@68 1586 uint8_t VCF_BND
jpayne@68 1587 uint8_t VCF_OVERLAP
jpayne@68 1588
jpayne@68 1589
jpayne@68 1590 ctypedef struct variant_t:
jpayne@68 1591 int type, n # variant type and the number of bases affected, negative for deletions
jpayne@68 1592
jpayne@68 1593 ctypedef struct bcf_fmt_t:
jpayne@68 1594 int id # id: numeric tag id, the corresponding string is bcf_hdr_t::id[BCF_DT_ID][$id].key
jpayne@68 1595 int n, size, type # n: number of values per-sample; size: number of bytes per-sample; type: one of BCF_BT_* types
jpayne@68 1596 uint8_t *p # same as vptr and vptr_* in bcf_info_t below
jpayne@68 1597 uint32_t p_len
jpayne@68 1598 uint32_t p_off
jpayne@68 1599 uint8_t p_free
jpayne@68 1600
jpayne@68 1601 union bcf_info_union_t:
jpayne@68 1602 int32_t i # integer value
jpayne@68 1603 float f # float value
jpayne@68 1604
jpayne@68 1605 ctypedef struct bcf_info_t:
jpayne@68 1606 int key # key: numeric tag id, the corresponding string is bcf_hdr_t::id[BCF_DT_ID][$key].key
jpayne@68 1607 int type, len # type: one of BCF_BT_* types; len: vector length, 1 for scalars
jpayne@68 1608
jpayne@68 1609 # v1 union only set if $len==1; for easier access
jpayne@68 1610 bcf_info_union_t v1
jpayne@68 1611 uint8_t *vptr # pointer to data array in bcf1_t->shared.s, excluding the size+type and tag id bytes
jpayne@68 1612 uint32_t vptr_len # length of the vptr block or, when set, of the vptr_mod block, excluding offset
jpayne@68 1613 uint32_t vptr_off # vptr offset, i.e., the size of the INFO key plus size+type bytes
jpayne@68 1614 uint8_t vptr_free # indicates that vptr-vptr_off must be freed; set only when modified and the new
jpayne@68 1615 # data block is bigger than the original
jpayne@68 1616
jpayne@68 1617 uint8_t BCF1_DIRTY_ID
jpayne@68 1618 uint8_t BCF1_DIRTY_ALS
jpayne@68 1619 uint8_t BCF1_DIRTY_FLT
jpayne@68 1620 uint8_t BCF1_DIRTY_INF
jpayne@68 1621
jpayne@68 1622 ctypedef struct bcf_dec_t:
jpayne@68 1623 int m_fmt, m_info, m_id, m_als, m_allele, m_flt # allocated size (high-water mark); do not change
jpayne@68 1624 int n_flt # Number of FILTER fields
jpayne@68 1625 int *flt # FILTER keys in the dictionary
jpayne@68 1626 char *id # ID
jpayne@68 1627 char *als # REF+ALT block (\0-seperated)
jpayne@68 1628 char **allele # allele[0] is the REF (allele[] pointers to the als block); all null terminated
jpayne@68 1629 bcf_info_t *info # INFO
jpayne@68 1630 bcf_fmt_t *fmt # FORMAT and individual sample
jpayne@68 1631 variant_t *var # $var and $var_type set only when set_variant_types called
jpayne@68 1632 int n_var, var_type
jpayne@68 1633 int shared_dirty # if set, shared.s must be recreated on BCF output
jpayne@68 1634 int indiv_dirty # if set, indiv.s must be recreated on BCF output
jpayne@68 1635
jpayne@68 1636 uint8_t BCF_ERR_CTG_UNDEF
jpayne@68 1637 uint8_t BCF_ERR_TAG_UNDEF
jpayne@68 1638 uint8_t BCF_ERR_NCOLS
jpayne@68 1639 uint8_t BCF_ERR_LIMITS
jpayne@68 1640 uint8_t BCF_ERR_CHAR
jpayne@68 1641 uint8_t BCF_ERR_CTG_INVALID
jpayne@68 1642 uint8_t BCF_ERR_TAG_INVALID
jpayne@68 1643
jpayne@68 1644 # The bcf1_t structure corresponds to one VCF/BCF line. Reading from VCF file
jpayne@68 1645 # is slower because the string is first to be parsed, packed into BCF line
jpayne@68 1646 # (done in vcf_parse), then unpacked into internal bcf1_t structure. If it
jpayne@68 1647 # is known in advance that some of the fields will not be required (notably
jpayne@68 1648 # the sample columns), parsing of these can be skipped by setting max_unpack
jpayne@68 1649 # appropriately.
jpayne@68 1650 # Similarly, it is fast to output a BCF line because the columns (kept in
jpayne@68 1651 # shared.s, indiv.s, etc.) are written directly by bcf_write, whereas a VCF
jpayne@68 1652 # line must be formatted in vcf_format.
jpayne@68 1653
jpayne@68 1654 ctypedef struct bcf1_t:
jpayne@68 1655 int32_t rid # CHROM
jpayne@68 1656 int32_t pos # POS
jpayne@68 1657 int32_t rlen # length of REF
jpayne@68 1658 float qual # QUAL
jpayne@68 1659 uint32_t n_info, n_allele
jpayne@68 1660 uint32_t n_fmt, n_sample
jpayne@68 1661 kstring_t shared, indiv
jpayne@68 1662 bcf_dec_t d # lazy evaluation: $d is not generated by bcf_read(), but by explicitly calling bcf_unpack()
jpayne@68 1663 int max_unpack # Set to BCF_UN_STR, BCF_UN_FLT, or BCF_UN_INFO to boost performance of vcf_parse when some of the fields won't be needed
jpayne@68 1664 int unpacked # remember what has been unpacked to allow calling bcf_unpack() repeatedly without redoing the work
jpayne@68 1665 int unpack_size[3] # the original block size of ID, REF+ALT and FILTER
jpayne@68 1666 int errcode # one of BCF_ERR_* codes
jpayne@68 1667
jpayne@68 1668 ####### API #######
jpayne@68 1669
jpayne@68 1670 # BCF and VCF I/O
jpayne@68 1671 #
jpayne@68 1672 # A note about naming conventions: htslib internally represents VCF
jpayne@68 1673 # records as bcf1_t data structures, therefore most functions are
jpayne@68 1674 # prefixed with bcf_. There are a few exceptions where the functions must
jpayne@68 1675 # be aware of both BCF and VCF worlds, such as bcf_parse vs vcf_parse. In
jpayne@68 1676 # these cases, functions prefixed with bcf_ are more general and work
jpayne@68 1677 # with both BCF and VCF.
jpayne@68 1678
jpayne@68 1679 # bcf_hdr_init() - create an empty BCF header.
jpayne@68 1680 # @param mode "r" or "w"
jpayne@68 1681 #
jpayne@68 1682 # When opened for writing, the mandatory fileFormat and
jpayne@68 1683 # FILTER=PASS lines are added automatically.
jpayne@68 1684 bcf_hdr_t *bcf_hdr_init(const char *mode)
jpayne@68 1685
jpayne@68 1686 # Destroy a BCF header struct
jpayne@68 1687 void bcf_hdr_destroy(bcf_hdr_t *h)
jpayne@68 1688
jpayne@68 1689 # Initialize a bcf1_t object; equivalent to calloc(1, sizeof(bcf1_t))
jpayne@68 1690 bcf1_t *bcf_init()
jpayne@68 1691
jpayne@68 1692 # Deallocate a bcf1_t object
jpayne@68 1693 void bcf_destroy(bcf1_t *v)
jpayne@68 1694
jpayne@68 1695 # Same as bcf_destroy() but frees only the memory allocated by bcf1_t,
jpayne@68 1696 # not the bcf1_t object itself.
jpayne@68 1697 void bcf_empty(bcf1_t *v)
jpayne@68 1698
jpayne@68 1699 # Make the bcf1_t object ready for next read. Intended mostly for
jpayne@68 1700 # internal use, the user should rarely need to call this function
jpayne@68 1701 # directly.
jpayne@68 1702 void bcf_clear(bcf1_t *v)
jpayne@68 1703
jpayne@68 1704 # Reads VCF or BCF header
jpayne@68 1705 bcf_hdr_t *bcf_hdr_read(htsFile *fp)
jpayne@68 1706
jpayne@68 1707 # bcf_hdr_set_samples() - for more efficient VCF parsing when only one/few samples are needed
jpayne@68 1708 # @samples: samples to include or exclude from file or as a comma-separated string.
jpayne@68 1709 # LIST|FILE .. select samples in list/file
jpayne@68 1710 # ^LIST|FILE .. exclude samples from list/file
jpayne@68 1711 # - .. include all samples
jpayne@68 1712 # NULL .. exclude all samples
jpayne@68 1713 # @is_file: @samples is a file (1) or a comma-separated list (0)
jpayne@68 1714 #
jpayne@68 1715 # The bottleneck of VCF reading is parsing of genotype fields. If the
jpayne@68 1716 # reader knows in advance that only subset of samples is needed (possibly
jpayne@68 1717 # no samples at all), the performance of bcf_read() can be significantly
jpayne@68 1718 # improved by calling bcf_hdr_set_samples after bcf_hdr_read().
jpayne@68 1719 # The function bcf_read() will subset the VCF/BCF records automatically
jpayne@68 1720 # with the notable exception when reading records via bcf_itr_next().
jpayne@68 1721 # In this case, bcf_subset_format() must be called explicitly, because
jpayne@68 1722 # bcf_readrec() does not see the header.
jpayne@68 1723 #
jpayne@68 1724 # Returns 0 on success, -1 on error or a positive integer if the list
jpayne@68 1725 # contains samples not present in the VCF header. In such a case, the
jpayne@68 1726 # return value is the index of the offending sample.
jpayne@68 1727 #
jpayne@68 1728 int bcf_hdr_set_samples(bcf_hdr_t *hdr, const char *samples, int is_file)
jpayne@68 1729 int bcf_subset_format(const bcf_hdr_t *hdr, bcf1_t *rec)
jpayne@68 1730
jpayne@68 1731 # Writes VCF or BCF header
jpayne@68 1732 int bcf_hdr_write(htsFile *fp, bcf_hdr_t *h)
jpayne@68 1733
jpayne@68 1734 # Parse VCF line contained in kstring and populate the bcf1_t struct
jpayne@68 1735 int vcf_parse(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v)
jpayne@68 1736
jpayne@68 1737 # The opposite of vcf_parse. It should rarely be called directly, see vcf_write
jpayne@68 1738 int vcf_format(const bcf_hdr_t *h, const bcf1_t *v, kstring_t *s)
jpayne@68 1739
jpayne@68 1740 # bcf_read() - read next VCF or BCF record
jpayne@68 1741 #
jpayne@68 1742 # Returns -1 on critical errors, 0 otherwise. On errors which are not
jpayne@68 1743 # critical for reading, such as missing header definitions, v->errcode is
jpayne@68 1744 # set to one of BCF_ERR* code and must be checked before calling
jpayne@68 1745 # vcf_write().
jpayne@68 1746 int bcf_read(htsFile *fp, const bcf_hdr_t *h, bcf1_t *v)
jpayne@68 1747
jpayne@68 1748 # bcf_unpack() - unpack/decode a BCF record (fills the bcf1_t::d field)
jpayne@68 1749 #
jpayne@68 1750 # Note that bcf_unpack() must be called even when reading VCF. It is safe
jpayne@68 1751 # to call the function repeatedly, it will not unpack the same field
jpayne@68 1752 # twice.
jpayne@68 1753 uint8_t BCF_UN_STR # up to ALT inclusive
jpayne@68 1754 uint8_t BCF_UN_FLT # up to FILTER
jpayne@68 1755 uint8_t BCF_UN_INFO # up to INFO
jpayne@68 1756 uint8_t BCF_UN_SHR # all shared information
jpayne@68 1757 uint8_t BCF_UN_FMT # unpack format and each sample
jpayne@68 1758 uint8_t BCF_UN_IND # a synonymo of BCF_UN_FMT
jpayne@68 1759 uint8_t BCF_UN_ALL # everything
jpayne@68 1760
jpayne@68 1761 int bcf_unpack(bcf1_t *b, int which)
jpayne@68 1762
jpayne@68 1763 # bcf_dup() - create a copy of BCF record.
jpayne@68 1764 #
jpayne@68 1765 # Note that bcf_unpack() must be called on the returned copy as if it was
jpayne@68 1766 # obtained from bcf_read(). Also note that bcf_dup() calls bcf_sync1(src)
jpayne@68 1767 # internally to reflect any changes made by bcf_update_* functions.
jpayne@68 1768 bcf1_t *bcf_dup(bcf1_t *src)
jpayne@68 1769 bcf1_t *bcf_copy(bcf1_t *dst, bcf1_t *src)
jpayne@68 1770
jpayne@68 1771 # bcf_write() - write one VCF or BCF record. The type is determined at the open() call.
jpayne@68 1772 int bcf_write(htsFile *fp, bcf_hdr_t *h, bcf1_t *v)
jpayne@68 1773
jpayne@68 1774 # The following functions work only with VCFs and should rarely be called
jpayne@68 1775 # directly. Usually one wants to use their bcf_* alternatives, which work
jpayne@68 1776 # transparently with both VCFs and BCFs.
jpayne@68 1777 bcf_hdr_t *vcf_hdr_read(htsFile *fp)
jpayne@68 1778 int vcf_hdr_write(htsFile *fp, const bcf_hdr_t *h)
jpayne@68 1779 int vcf_read(htsFile *fp, const bcf_hdr_t *h, bcf1_t *v)
jpayne@68 1780 int vcf_write(htsFile *fp, const bcf_hdr_t *h, bcf1_t *v)
jpayne@68 1781
jpayne@68 1782 #************************************************************************
jpayne@68 1783 # Header querying and manipulation routines
jpayne@68 1784 #************************************************************************
jpayne@68 1785
jpayne@68 1786 # Create a new header using the supplied template
jpayne@68 1787 bcf_hdr_t *bcf_hdr_dup(const bcf_hdr_t *hdr)
jpayne@68 1788
jpayne@68 1789 # Copy header lines from src to dst if not already present in dst. See also bcf_translate().
jpayne@68 1790 # Returns 0 on success or sets a bit on error:
jpayne@68 1791 # 1 .. conflicting definitions of tag length
jpayne@68 1792 # # todo
jpayne@68 1793 int bcf_hdr_combine(bcf_hdr_t *dst, const bcf_hdr_t *src)
jpayne@68 1794
jpayne@68 1795 # bcf_hdr_merge() - copy header lines from src to dst, see also bcf_translate()
jpayne@68 1796 # @param dst: the destination header to be merged into, NULL on the first pass
jpayne@68 1797 # @param src: the source header
jpayne@68 1798 #
jpayne@68 1799 # Notes:
jpayne@68 1800 # - use as:
jpayne@68 1801 # bcf_hdr_t *dst = NULL;
jpayne@68 1802 # for (i=0; i<nsrc; i++) dst = bcf_hdr_merge(dst,src[i]);
jpayne@68 1803 #
jpayne@68 1804 # - bcf_hdr_merge() replaces bcf_hdr_combine() which had a problem when
jpayne@68 1805 # combining multiple BCF headers. The current bcf_hdr_combine()
jpayne@68 1806 # does not have this problem, but became slow when used for many files.
jpayne@68 1807 bcf_hdr_t *bcf_hdr_merge(bcf_hdr_t *dst, const bcf_hdr_t *src)
jpayne@68 1808
jpayne@68 1809 # bcf_hdr_add_sample() - add a new sample.
jpayne@68 1810 # @param sample: sample name to be added
jpayne@68 1811 int bcf_hdr_add_sample(bcf_hdr_t *hdr, const char *sample)
jpayne@68 1812
jpayne@68 1813 # Read VCF header from a file and update the header
jpayne@68 1814 int bcf_hdr_set(bcf_hdr_t *hdr, const char *fname)
jpayne@68 1815
jpayne@68 1816 # Appends formatted header text to _str_.
jpayne@68 1817 # If _is_bcf_ is zero, `IDX` fields are discarded.
jpayne@68 1818 # @return 0 if successful, or negative if an error occurred
jpayne@68 1819 # @since 1.4
jpayne@68 1820 int bcf_hdr_format(const bcf_hdr_t *hdr, int is_bcf, kstring_t *str);
jpayne@68 1821
jpayne@68 1822 # Returns formatted header (newly allocated string) and its length,
jpayne@68 1823 # excluding the terminating \0. If is_bcf parameter is unset, IDX
jpayne@68 1824 # fields are discarded.
jpayne@68 1825 char *bcf_hdr_fmt_text(const bcf_hdr_t *hdr, int is_bcf, int *len)
jpayne@68 1826
jpayne@68 1827 # Append new VCF header line, returns 0 on success
jpayne@68 1828 int bcf_hdr_append(bcf_hdr_t *h, const char *line)
jpayne@68 1829 int bcf_hdr_printf(bcf_hdr_t *h, const char *format, ...)
jpayne@68 1830
jpayne@68 1831 # VCF version, e.g. VCFv4.2
jpayne@68 1832 const char *bcf_hdr_get_version(const bcf_hdr_t *hdr)
jpayne@68 1833 void bcf_hdr_set_version(bcf_hdr_t *hdr, const char *version)
jpayne@68 1834
jpayne@68 1835 # bcf_hdr_remove() - remove VCF header tag
jpayne@68 1836 # @param type: one of BCF_HL_*
jpayne@68 1837 # @param key: tag name or NULL to remove all tags of the given type
jpayne@68 1838 void bcf_hdr_remove(bcf_hdr_t *h, int type, const char *key)
jpayne@68 1839
jpayne@68 1840 # bcf_hdr_subset() - creates a new copy of the header removing unwanted samples
jpayne@68 1841 # @param n: number of samples to keep
jpayne@68 1842 # @param samples: names of the samples to keep
jpayne@68 1843 # @param imap: mapping from index in @samples to the sample index in the original file
jpayne@68 1844 #
jpayne@68 1845 # Sample names not present in h0 are ignored. The number of unmatched samples can be checked
jpayne@68 1846 # by comparing n and bcf_hdr_nsamples(out_hdr).
jpayne@68 1847 # This function can be used to reorder samples.
jpayne@68 1848 # See also bcf_subset() which subsets individual records.
jpayne@68 1849 #
jpayne@68 1850 bcf_hdr_t *bcf_hdr_subset(const bcf_hdr_t *h0, int n, char *const* samples, int *imap)
jpayne@68 1851
jpayne@68 1852 # Creates a list of sequence names. It is up to the caller to free the list (but not the sequence names)
jpayne@68 1853 const char **bcf_hdr_seqnames(const bcf_hdr_t *h, int *nseqs)
jpayne@68 1854
jpayne@68 1855 # Get number of samples
jpayne@68 1856 int32_t bcf_hdr_nsamples(const bcf_hdr_t *h)
jpayne@68 1857
jpayne@68 1858 # The following functions are for internal use and should rarely be called directly
jpayne@68 1859 int bcf_hdr_parse(bcf_hdr_t *hdr, char *htxt)
jpayne@68 1860 int bcf_hdr_sync(bcf_hdr_t *h)
jpayne@68 1861 bcf_hrec_t *bcf_hdr_parse_line(const bcf_hdr_t *h, const char *line, int *len)
jpayne@68 1862 void bcf_hrec_format(const bcf_hrec_t *hrec, kstring_t *str)
jpayne@68 1863 int bcf_hdr_add_hrec(bcf_hdr_t *hdr, bcf_hrec_t *hrec)
jpayne@68 1864
jpayne@68 1865 # bcf_hdr_get_hrec() - get header line info
jpayne@68 1866 # @param type: one of the BCF_HL_* types: FLT,INFO,FMT,CTG,STR,GEN
jpayne@68 1867 # @param key: the header key for generic lines (e.g. "fileformat"), any field
jpayne@68 1868 # for structured lines, typically "ID".
jpayne@68 1869 # @param value: the value which pairs with key. Can be be NULL for BCF_HL_GEN
jpayne@68 1870 # @param str_class: the class of BCF_HL_STR line (e.g. "ALT" or "SAMPLE"), otherwise NULL
jpayne@68 1871 #
jpayne@68 1872 bcf_hrec_t *bcf_hdr_get_hrec(const bcf_hdr_t *hdr, int type, const char *key, const char *value, const char *str_class)
jpayne@68 1873 bcf_hrec_t *bcf_hrec_dup(bcf_hrec_t *hrec)
jpayne@68 1874 void bcf_hrec_add_key(bcf_hrec_t *hrec, const char *str, int len)
jpayne@68 1875 void bcf_hrec_set_val(bcf_hrec_t *hrec, int i, const char *str, int len, int is_quoted)
jpayne@68 1876 int bcf_hrec_find_key(bcf_hrec_t *hrec, const char *key)
jpayne@68 1877 void hrec_add_idx(bcf_hrec_t *hrec, int idx)
jpayne@68 1878 void bcf_hrec_destroy(bcf_hrec_t *hrec)
jpayne@68 1879
jpayne@68 1880 #************************************************************************
jpayne@68 1881 # Individual record querying and manipulation routines
jpayne@68 1882 #************************************************************************
jpayne@68 1883
jpayne@68 1884 # See the description of bcf_hdr_subset()
jpayne@68 1885 int bcf_subset(const bcf_hdr_t *h, bcf1_t *v, int n, int *imap)
jpayne@68 1886
jpayne@68 1887 # bcf_translate() - translate tags ids to be consistent with different header. This function
jpayne@68 1888 # is useful when lines from multiple VCF need to be combined.
jpayne@68 1889 # @dst_hdr: the destination header, to be used in bcf_write(), see also bcf_hdr_combine()
jpayne@68 1890 # @src_hdr: the source header, used in bcf_read()
jpayne@68 1891 # @src_line: line obtained by bcf_read()
jpayne@68 1892 int bcf_translate(const bcf_hdr_t *dst_hdr, bcf_hdr_t *src_hdr, bcf1_t *src_line)
jpayne@68 1893
jpayne@68 1894 # bcf_get_variant_type[s]() - returns one of VCF_REF, VCF_SNP, etc
jpayne@68 1895 int bcf_get_variant_types(bcf1_t *rec)
jpayne@68 1896 int bcf_get_variant_type(bcf1_t *rec, int ith_allele)
jpayne@68 1897 int bcf_is_snp(bcf1_t *v)
jpayne@68 1898
jpayne@68 1899 # bcf_update_filter() - sets the FILTER column
jpayne@68 1900 # @flt_ids: The filter IDs to set, numeric IDs returned by bcf_hdr_id2int(hdr, BCF_DT_ID, "PASS")
jpayne@68 1901 # @n: Number of filters. If n==0, all filters are removed
jpayne@68 1902 int bcf_update_filter(const bcf_hdr_t *hdr, bcf1_t *line, int *flt_ids, int n)
jpayne@68 1903
jpayne@68 1904 # bcf_add_filter() - adds to the FILTER column
jpayne@68 1905 # @flt_id: The filter IDs to add, numeric IDs returned by bcf_hdr_id2int(hdr, BCF_DT_ID, "PASS")
jpayne@68 1906 #
jpayne@68 1907 # If flt_id is PASS, all existing filters are removed first. If other than PASS, existing PASS is removed.
jpayne@68 1908 int bcf_add_filter(const bcf_hdr_t *hdr, bcf1_t *line, int flt_id)
jpayne@68 1909
jpayne@68 1910 # bcf_remove_filter() - removes from the FILTER column
jpayne@68 1911 # @flt_id: filter ID to remove, numeric ID returned by bcf_hdr_id2int(hdr, BCF_DT_ID, "PASS")
jpayne@68 1912 # @pass: when set to 1 and no filters are present, set to PASS
jpayne@68 1913 int bcf_remove_filter(const bcf_hdr_t *hdr, bcf1_t *line, int flt_id, int set_pass)
jpayne@68 1914
jpayne@68 1915 # Returns 1 if present, 0 if absent, or -1 if filter does not exist. "PASS" and "." can be used interchangeably.
jpayne@68 1916 int bcf_has_filter(const bcf_hdr_t *hdr, bcf1_t *line, char *filter)
jpayne@68 1917
jpayne@68 1918 # bcf_update_alleles() and bcf_update_alleles_str() - update REF and ALT column
jpayne@68 1919 # @alleles: Array of alleles
jpayne@68 1920 # @nals: Number of alleles
jpayne@68 1921 # @alleles_string: Comma-separated alleles, starting with the REF allele
jpayne@68 1922 int bcf_update_alleles(const bcf_hdr_t *hdr, bcf1_t *line, const char **alleles, int nals)
jpayne@68 1923 int bcf_update_alleles_str(const bcf_hdr_t *hdr, bcf1_t *line, const char *alleles_string)
jpayne@68 1924
jpayne@68 1925 # bcf_update_id() - sets new ID string
jpayne@68 1926 # bcf_add_id() - adds to the ID string checking for duplicates
jpayne@68 1927 int bcf_update_id(const bcf_hdr_t *hdr, bcf1_t *line, const char *id)
jpayne@68 1928 int bcf_add_id(const bcf_hdr_t *hdr, bcf1_t *line, const char *id)
jpayne@68 1929
jpayne@68 1930 # bcf_update_info_*() - functions for updating INFO fields
jpayne@68 1931 # @hdr: the BCF header
jpayne@68 1932 # @line: VCF line to be edited
jpayne@68 1933 # @key: the INFO tag to be updated
jpayne@68 1934 # @values: pointer to the array of values. Pass NULL to remove the tag.
jpayne@68 1935 # @n: number of values in the array. When set to 0, the INFO tag is removed
jpayne@68 1936 #
jpayne@68 1937 # The @string in bcf_update_info_flag() is optional, @n indicates whether
jpayne@68 1938 # the flag is set or removed.
jpayne@68 1939 #
jpayne@68 1940 # Returns 0 on success or negative value on error.
jpayne@68 1941 #
jpayne@68 1942 int bcf_update_info_int32(const bcf_hdr_t *hdr, bcf1_t *line, const char *key, const int32_t *values, int n)
jpayne@68 1943 int bcf_update_info_float(const bcf_hdr_t *hdr, bcf1_t *line, const char *key, const float *values, int n)
jpayne@68 1944 int bcf_update_info_flag(const bcf_hdr_t *hdr, bcf1_t *line, const char *key, const char *values, int n)
jpayne@68 1945 int bcf_update_info_string(const bcf_hdr_t *hdr, bcf1_t *line, const char *key, const char *values, int n)
jpayne@68 1946 int bcf_update_info(const bcf_hdr_t *hdr, bcf1_t *line, const char *key, const void *values, int n, int type)
jpayne@68 1947
jpayne@68 1948 # bcf_update_format_*() - functions for updating FORMAT fields
jpayne@68 1949 # @values: pointer to the array of values, the same number of elements
jpayne@68 1950 # is expected for each sample. Missing values must be padded
jpayne@68 1951 # with bcf_*_missing or bcf_*_vector_end values.
jpayne@68 1952 # @n: number of values in the array. If n==0, existing tag is removed.
jpayne@68 1953 #
jpayne@68 1954 # The function bcf_update_format_string() is a higher-level (slower) variant of
jpayne@68 1955 # bcf_update_format_char(). The former accepts array of \0-terminated strings
jpayne@68 1956 # whereas the latter requires that the strings are collapsed into a single array
jpayne@68 1957 # of fixed-length strings. In case of strings with variable length, shorter strings
jpayne@68 1958 # can be \0-padded. Note that the collapsed strings passed to bcf_update_format_char()
jpayne@68 1959 # are not \0-terminated.
jpayne@68 1960 #
jpayne@68 1961 # Returns 0 on success or negative value on error.
jpayne@68 1962 #
jpayne@68 1963 int bcf_update_format_int32(const bcf_hdr_t *hdr, bcf1_t *line, const char *key, const int32_t *values, int n)
jpayne@68 1964 int bcf_update_format_float(const bcf_hdr_t *hdr, bcf1_t *line, const char *key, const float *values, int n)
jpayne@68 1965 int bcf_update_format_char(const bcf_hdr_t *hdr, bcf1_t *line, const char *key, const char *values, int n)
jpayne@68 1966 int bcf_update_genotypes(const bcf_hdr_t *hdr, bcf1_t *line, const int32_t *values, int n)
jpayne@68 1967 int bcf_update_format_string(const bcf_hdr_t *hdr, bcf1_t *line, const char *key, const char **values, int n)
jpayne@68 1968 int bcf_update_format(const bcf_hdr_t *hdr, bcf1_t *line, const char *key, const void *values, int n, int type)
jpayne@68 1969
jpayne@68 1970 # Macros for setting genotypes correctly, for use with bcf_update_genotypes only; idx corresponds
jpayne@68 1971 # to VCF's GT (1-based index to ALT or 0 for the reference allele) and val is the opposite, obtained
jpayne@68 1972 # from bcf_get_genotypes() below.
jpayne@68 1973 uint32_t bcf_gt_phased(uint32_t idx)
jpayne@68 1974 uint32_t bcf_gt_unphased(uint32_t idx)
jpayne@68 1975 uint32_t bcf_gt_missing
jpayne@68 1976 uint32_t bcf_gt_is_missing(uint32_t val)
jpayne@68 1977 uint32_t bcf_gt_is_phased(uint32_t idx)
jpayne@68 1978 uint32_t bcf_gt_allele(uint32_t val)
jpayne@68 1979
jpayne@68 1980 # Conversion between alleles indexes to Number=G genotype index (assuming diploid, all 0-based)
jpayne@68 1981 uint32_t bcf_alleles2gt(uint32_t a, uint32_t b)
jpayne@68 1982 void bcf_gt2alleles(int igt, int *a, int *b)
jpayne@68 1983
jpayne@68 1984 # bcf_get_fmt() - returns pointer to FORMAT's field data
jpayne@68 1985 # @header: for access to BCF_DT_ID dictionary
jpayne@68 1986 # @line: VCF line obtained from vcf_parse1
jpayne@68 1987 # @fmt: one of GT,PL,...
jpayne@68 1988 #
jpayne@68 1989 # Returns bcf_fmt_t* if the call succeeded, or returns NULL when the field
jpayne@68 1990 # is not available.
jpayne@68 1991 #
jpayne@68 1992 bcf_fmt_t *bcf_get_fmt(const bcf_hdr_t *hdr, bcf1_t *line, const char *key)
jpayne@68 1993 bcf_info_t *bcf_get_info(const bcf_hdr_t *hdr, bcf1_t *line, const char *key)
jpayne@68 1994
jpayne@68 1995 # bcf_get_*_id() - returns pointer to FORMAT/INFO field data given the header index instead of the string ID
jpayne@68 1996 # @line: VCF line obtained from vcf_parse1
jpayne@68 1997 # @id: The header index for the tag, obtained from bcf_hdr_id2int()
jpayne@68 1998 #
jpayne@68 1999 # Returns bcf_fmt_t* / bcf_info_t*. These functions do not check if the index is valid
jpayne@68 2000 # as their goal is to avoid the header lookup.
jpayne@68 2001 #
jpayne@68 2002 bcf_fmt_t *bcf_get_fmt_id(bcf1_t *line, const int id)
jpayne@68 2003 bcf_info_t *bcf_get_info_id(bcf1_t *line, const int id)
jpayne@68 2004
jpayne@68 2005 # bcf_get_info_*() - get INFO values, integers or floats
jpayne@68 2006 # @hdr: BCF header
jpayne@68 2007 # @line: BCF record
jpayne@68 2008 # @tag: INFO tag to retrieve
jpayne@68 2009 # @dst: *dst is pointer to a memory location, can point to NULL
jpayne@68 2010 # @ndst: pointer to the size of allocated memory
jpayne@68 2011 #
jpayne@68 2012 # Returns negative value on error or the number of written values on
jpayne@68 2013 # success. bcf_get_info_string() returns on success the number of
jpayne@68 2014 # characters written excluding the null-terminating byte. bcf_get_info_flag()
jpayne@68 2015 # returns 1 when flag is set or 0 if not.
jpayne@68 2016 #
jpayne@68 2017 # List of return codes:
jpayne@68 2018 # -1 .. no such INFO tag defined in the header
jpayne@68 2019 # -2 .. clash between types defined in the header and encountered in the VCF record
jpayne@68 2020 # -3 .. tag is not present in the VCF record
jpayne@68 2021 #
jpayne@68 2022 int bcf_get_info_int32(const bcf_hdr_t *hdr, bcf1_t *line, const char *tag, int32_t **dst, int *ndst)
jpayne@68 2023 int bcf_get_info_float(const bcf_hdr_t *hdr, bcf1_t *line, const char *tag, float **dst, int *ndst)
jpayne@68 2024 int bcf_get_info_string(const bcf_hdr_t *hdr, bcf1_t *line, const char *tag, char **dst, int *ndst)
jpayne@68 2025 int bcf_get_info_flag(const bcf_hdr_t *hdr, bcf1_t *line, const char *tag, int **dst, int *ndst)
jpayne@68 2026 int bcf_get_info_values(const bcf_hdr_t *hdr, bcf1_t *line, const char *tag, void **dst, int *ndst, int type)
jpayne@68 2027
jpayne@68 2028 # bcf_get_format_*() - same as bcf_get_info*() above
jpayne@68 2029 #
jpayne@68 2030 # The function bcf_get_format_string() is a higher-level (slower) variant of bcf_get_format_char().
jpayne@68 2031 # see the description of bcf_update_format_string() and bcf_update_format_char() above.
jpayne@68 2032 # Unlike other bcf_get_format__*() functions, bcf_get_format_string() allocates two arrays:
jpayne@68 2033 # a single block of \0-terminated strings collapsed into a single array and an array of pointers
jpayne@68 2034 # to these strings. Both arrays must be cleaned by the user.
jpayne@68 2035 #
jpayne@68 2036 # Returns negative value on error or the number of written values on success.
jpayne@68 2037 #
jpayne@68 2038 # Example:
jpayne@68 2039 # int ndst = 0; char **dst = NULL
jpayne@68 2040 # if ( bcf_get_format_string(hdr, line, "XX", &dst, &ndst) > 0 )
jpayne@68 2041 # for (i=0; i<bcf_hdr_nsamples(hdr); i++) printf("%s\n", dst[i])
jpayne@68 2042 # free(dst[0]); free(dst)
jpayne@68 2043 #
jpayne@68 2044 # Example:
jpayne@68 2045 # int ngt, *gt_arr = NULL, ngt_arr = 0
jpayne@68 2046 # ngt = bcf_get_genotypes(hdr, line, &gt_arr, &ngt_arr)
jpayne@68 2047 #
jpayne@68 2048 int bcf_get_format_int32(const bcf_hdr_t *hdr, bcf1_t *line, const char *tag, int32_t **dst, int *ndst)
jpayne@68 2049 int bcf_get_format_float(const bcf_hdr_t *hdr, bcf1_t *line, const char *tag, float **dst, int *ndst)
jpayne@68 2050 int bcf_get_format_char(const bcf_hdr_t *hdr, bcf1_t *line, const char *tag, char **dst, int *ndst)
jpayne@68 2051 int bcf_get_genotypes(const bcf_hdr_t *hdr, bcf1_t *line, int32_t **dst, int *ndst)
jpayne@68 2052 int bcf_get_format_string(const bcf_hdr_t *hdr, bcf1_t *line, const char *tag, char ***dst, int *ndst)
jpayne@68 2053 int bcf_get_format_values(const bcf_hdr_t *hdr, bcf1_t *line, const char *tag, void **dst, int *ndst, int type)
jpayne@68 2054
jpayne@68 2055 #************************************************************************
jpayne@68 2056 # Helper functions
jpayne@68 2057 #************************************************************************
jpayne@68 2058
jpayne@68 2059 #
jpayne@68 2060 # bcf_hdr_id2int() - Translates string into numeric ID
jpayne@68 2061 # bcf_hdr_int2id() - Translates numeric ID into string
jpayne@68 2062 # @type: one of BCF_DT_ID, BCF_DT_CTG, BCF_DT_SAMPLE
jpayne@68 2063 # @id: tag name, such as: PL, DP, GT, etc.
jpayne@68 2064 #
jpayne@68 2065 # Returns -1 if string is not in dictionary, otherwise numeric ID which identifies
jpayne@68 2066 # fields in BCF records.
jpayne@68 2067 #
jpayne@68 2068 int bcf_hdr_id2int(const bcf_hdr_t *hdr, int type, const char *id)
jpayne@68 2069 const char *bcf_hdr_int2id(const bcf_hdr_t *hdr, int type, int int_id)
jpayne@68 2070
jpayne@68 2071 # bcf_hdr_name2id() - Translates sequence names (chromosomes) into numeric ID
jpayne@68 2072 # bcf_hdr_id2name() - Translates numeric ID to sequence name
jpayne@68 2073 #
jpayne@68 2074 int bcf_hdr_name2id(const bcf_hdr_t *hdr, const char *id)
jpayne@68 2075 const char *bcf_hdr_id2name(const bcf_hdr_t *hdr, int rid)
jpayne@68 2076 const char *bcf_seqname(const bcf_hdr_t *hdr, bcf1_t *rec)
jpayne@68 2077
jpayne@68 2078 #
jpayne@68 2079 # bcf_hdr_id2*() - Macros for accessing bcf_idinfo_t
jpayne@68 2080 # @type: one of BCF_HL_FLT, BCF_HL_INFO, BCF_HL_FMT
jpayne@68 2081 # @int_id: return value of bcf_hdr_id2int, must be >=0
jpayne@68 2082 #
jpayne@68 2083 # The returned values are:
jpayne@68 2084 # bcf_hdr_id2length .. whether the number of values is fixed or variable, one of BCF_VL_*
jpayne@68 2085 # bcf_hdr_id2number .. the number of values, 0xfffff for variable length fields
jpayne@68 2086 # bcf_hdr_id2type .. the field type, one of BCF_HT_*
jpayne@68 2087 # bcf_hdr_id2coltype .. the column type, one of BCF_HL_*
jpayne@68 2088 #
jpayne@68 2089 # Notes: Prior to using the macros, the presence of the info should be
jpayne@68 2090 # tested with bcf_hdr_idinfo_exists().
jpayne@68 2091 #
jpayne@68 2092 int bcf_hdr_id2length(const bcf_hdr_t *hdr, int type, int int_id)
jpayne@68 2093 int bcf_hdr_id2number(const bcf_hdr_t *hdr, int type, int int_id)
jpayne@68 2094 int bcf_hdr_id2type(const bcf_hdr_t *hdr, int type, int int_id)
jpayne@68 2095 int bcf_hdr_id2coltype(const bcf_hdr_t *hdr, int type, int int_id)
jpayne@68 2096 int bcf_hdr_idinfo_exists(const bcf_hdr_t *hdr, int type, int int_id)
jpayne@68 2097 bcf_hrec_t *bcf_hdr_id2hrec(const bcf_hdr_t *hdr, int type, int col_type, int int_id)
jpayne@68 2098
jpayne@68 2099 void bcf_fmt_array(kstring_t *s, int n, int type, void *data)
jpayne@68 2100 uint8_t *bcf_fmt_sized_array(kstring_t *s, uint8_t *ptr)
jpayne@68 2101
jpayne@68 2102 void bcf_enc_vchar(kstring_t *s, int l, const char *a)
jpayne@68 2103 void bcf_enc_vint(kstring_t *s, int n, int32_t *a, int wsize)
jpayne@68 2104 void bcf_enc_vfloat(kstring_t *s, int n, float *a)
jpayne@68 2105
jpayne@68 2106 #************************************************************************
jpayne@68 2107 # BCF index
jpayne@68 2108 #
jpayne@68 2109 # Note that these functions work with BCFs only. See synced_bcf_reader.h
jpayne@68 2110 # which provides (amongst other things) an API to work transparently with
jpayne@68 2111 # both indexed BCFs and VCFs.
jpayne@68 2112 #************************************************************************
jpayne@68 2113
jpayne@68 2114 hts_idx_t *bcf_index_load2(const char *fn, const char *fnidx)
jpayne@68 2115 hts_idx_t *bcf_index_load3(const char *fn, const char *fnidx, int flags)
jpayne@68 2116 int bcf_index_build(const char *fn, int min_shift)
jpayne@68 2117 int bcf_index_build2(const char *fn, const char *fnidx, int min_shift)
jpayne@68 2118
jpayne@68 2119 #*******************
jpayne@68 2120 # Typed value I/O *
jpayne@68 2121 #******************
jpayne@68 2122
jpayne@68 2123 # Note that in contrast with BCFv2.1 specification, HTSlib implementation
jpayne@68 2124 # allows missing values in vectors. For integer types, the values 0x80,
jpayne@68 2125 # 0x8000, 0x80000000 are interpreted as missing values and 0x81, 0x8001,
jpayne@68 2126 # 0x80000001 as end-of-vector indicators. Similarly for floats, the value of
jpayne@68 2127 # 0x7F800001 is interpreted as a missing value and 0x7F800002 as an
jpayne@68 2128 # end-of-vector indicator.
jpayne@68 2129 # Note that the end-of-vector byte is not part of the vector.
jpayne@68 2130
jpayne@68 2131 # This trial BCF version (v2.2) is compatible with the VCF specification and
jpayne@68 2132 # enables to handle correctly vectors with different ploidy in presence of
jpayne@68 2133 # missing values.
jpayne@68 2134
jpayne@68 2135 int32_t bcf_int8_vector_end
jpayne@68 2136 int32_t bcf_int16_vector_end
jpayne@68 2137 int32_t bcf_int32_vector_end
jpayne@68 2138 int32_t bcf_str_vector_end
jpayne@68 2139 int32_t bcf_int8_missing
jpayne@68 2140 int32_t bcf_int16_missing
jpayne@68 2141 int32_t bcf_int32_missing
jpayne@68 2142 int32_t bcf_str_missing
jpayne@68 2143
jpayne@68 2144 uint32_t bcf_float_vector_end
jpayne@68 2145 uint32_t bcf_float_missing
jpayne@68 2146
jpayne@68 2147 void bcf_float_set(float *ptr, uint32_t value)
jpayne@68 2148 void bcf_float_set_vector_end(float *x)
jpayne@68 2149 void bcf_float_set_missing(float *x)
jpayne@68 2150
jpayne@68 2151 int bcf_float_is_missing(float f)
jpayne@68 2152 int bcf_float_is_vector_end(float f)
jpayne@68 2153 void bcf_format_gt(bcf_fmt_t *fmt, int isample, kstring_t *str)
jpayne@68 2154 void bcf_enc_size(kstring_t *s, int size, int type)
jpayne@68 2155 int bcf_enc_inttype(long x)
jpayne@68 2156 void bcf_enc_int1(kstring_t *s, int32_t x)
jpayne@68 2157 int32_t bcf_dec_int1(const uint8_t *p, int type, uint8_t **q)
jpayne@68 2158 int32_t bcf_dec_typed_int1(const uint8_t *p, uint8_t **q)
jpayne@68 2159 int32_t bcf_dec_size(const uint8_t *p, uint8_t **q, int *type)
jpayne@68 2160
jpayne@68 2161 # These trivial wrappers are defined only for consistency with other parts of htslib
jpayne@68 2162 bcf1_t *bcf_init1()
jpayne@68 2163 int bcf_read1(htsFile *fp, const bcf_hdr_t *h, bcf1_t *v)
jpayne@68 2164 int vcf_read1(htsFile *fp, const bcf_hdr_t *h, bcf1_t *v)
jpayne@68 2165 int bcf_write1(htsFile *fp, const bcf_hdr_t *h, bcf1_t *v)
jpayne@68 2166 int vcf_write1(htsFile *fp, const bcf_hdr_t *h, bcf1_t *v)
jpayne@68 2167 void bcf_destroy1(bcf1_t *v)
jpayne@68 2168 void bcf_empty1(bcf1_t *v)
jpayne@68 2169 int vcf_parse1(kstring_t *s, const bcf_hdr_t *h, bcf1_t *v)
jpayne@68 2170 void bcf_clear1(bcf1_t *v)
jpayne@68 2171 int vcf_format1(const bcf_hdr_t *h, const bcf1_t *v, kstring_t *s)
jpayne@68 2172
jpayne@68 2173 # Other nice wrappers
jpayne@68 2174 void bcf_itr_destroy(hts_itr_t *iter)
jpayne@68 2175 hts_itr_t *bcf_itr_queryi(const hts_idx_t *idx, int tid, int beg, int end)
jpayne@68 2176 hts_itr_t *bcf_itr_querys(const hts_idx_t *idx, const bcf_hdr_t *hdr, char *s)
jpayne@68 2177 int bcf_itr_next(htsFile *fp, hts_itr_t *iter, void *r)
jpayne@68 2178 hts_idx_t *bcf_index_load(const char *fn)
jpayne@68 2179 const char **bcf_index_seqnames(const hts_idx_t *idx, const bcf_hdr_t *hdr, int *nptr)
jpayne@68 2180
jpayne@68 2181
jpayne@68 2182 # VCF/BCF utility functions
jpayne@68 2183 cdef extern from "htslib/vcfutils.h" nogil:
jpayne@68 2184 struct kbitset_t
jpayne@68 2185
jpayne@68 2186 # bcf_trim_alleles() - remove ALT alleles unused in genotype fields
jpayne@68 2187 # @header: for access to BCF_DT_ID dictionary
jpayne@68 2188 # @line: VCF line obtain from vcf_parse1
jpayne@68 2189 #
jpayne@68 2190 # Returns the number of removed alleles on success or negative
jpayne@68 2191 # on error:
jpayne@68 2192 # -1 .. some allele index is out of bounds
jpayne@68 2193 int bcf_trim_alleles(const bcf_hdr_t *header, bcf1_t *line)
jpayne@68 2194
jpayne@68 2195 # bcf_remove_alleles() - remove ALT alleles according to bitmask @mask
jpayne@68 2196 # @header: for access to BCF_DT_ID dictionary
jpayne@68 2197 # @line: VCF line obtained from vcf_parse1
jpayne@68 2198 # @mask: alleles to remove
jpayne@68 2199 #
jpayne@68 2200 # If you have more than 31 alleles, then the integer bit mask will
jpayne@68 2201 # overflow, so use bcf_remove_allele_set instead
jpayne@68 2202 void bcf_remove_alleles(const bcf_hdr_t *header, bcf1_t *line, int mask)
jpayne@68 2203
jpayne@68 2204 # bcf_remove_allele_set() - remove ALT alleles according to bitset @rm_set
jpayne@68 2205 # @header: for access to BCF_DT_ID dictionary
jpayne@68 2206 # @line: VCF line obtained from vcf_parse1
jpayne@68 2207 # @rm_set: pointer to kbitset_t object with bits set for allele
jpayne@68 2208 # indexes to remove
jpayne@68 2209 #
jpayne@68 2210 # Number=A,R,G INFO and FORMAT fields will be updated accordingly.
jpayne@68 2211 void bcf_remove_allele_set(const bcf_hdr_t *header, bcf1_t *line, kbitset_t *rm_set)
jpayne@68 2212
jpayne@68 2213 # bcf_calc_ac() - calculate the number of REF and ALT alleles
jpayne@68 2214 # @header: for access to BCF_DT_ID dictionary
jpayne@68 2215 # @line: VCF line obtained from vcf_parse1
jpayne@68 2216 # @ac: array of length line->n_allele
jpayne@68 2217 # @which: determine if INFO/AN,AC and indv fields be used
jpayne@68 2218 #
jpayne@68 2219 # Returns 1 if the call succeeded, or 0 if the value could not
jpayne@68 2220 # be determined.
jpayne@68 2221 #
jpayne@68 2222 # The value of @which determines if existing INFO/AC,AN can be
jpayne@68 2223 # used (BCF_UN_INFO) and and if indv fields can be split (BCF_UN_FMT).
jpayne@68 2224 int bcf_calc_ac(const bcf_hdr_t *header, bcf1_t *line, int *ac, int which)
jpayne@68 2225
jpayne@68 2226 # bcf_gt_type() - determines type of the genotype
jpayne@68 2227 # @fmt_ptr: the GT format field as set for example by set_fmt_ptr
jpayne@68 2228 # @isample: sample index (starting from 0)
jpayne@68 2229 # @ial: index of the 1st non-reference allele (starting from 1)
jpayne@68 2230 # @jal: index of the 2nd non-reference allele (starting from 1)
jpayne@68 2231 #
jpayne@68 2232 # Returns the type of the genotype (one of GT_HOM_RR, GT_HET_RA,
jpayne@68 2233 # GT_HOM_AA, GT_HET_AA, GT_HAPL_R, GT_HAPL_A or GT_UNKN). If $ial
jpayne@68 2234 # is not NULL and the genotype has one or more non-reference
jpayne@68 2235 # alleles, $ial will be set. In case of GT_HET_AA, $ial is the
jpayne@68 2236 # position of the allele which appeared first in ALT. If $jal is
jpayne@68 2237 # not null and the genotype is GT_HET_AA, $jal will be set and is
jpayne@68 2238 # the position of the second allele in ALT.
jpayne@68 2239 uint8_t GT_HOM_RR # note: the actual value of GT_* matters, used in dosage r2 calculation
jpayne@68 2240 uint8_t GT_HOM_AA
jpayne@68 2241 uint8_t GT_HET_RA
jpayne@68 2242 uint8_t GT_HET_AA
jpayne@68 2243 uint8_t GT_HAPL_R
jpayne@68 2244 uint8_t GT_HAPL_A
jpayne@68 2245 uint8_t GT_UNKN
jpayne@68 2246 int bcf_gt_type(bcf_fmt_t *fmt_ptr, int isample, int *ial, int *jal)
jpayne@68 2247
jpayne@68 2248 int bcf_acgt2int(char c)
jpayne@68 2249 char bcf_int2acgt(int i)
jpayne@68 2250
jpayne@68 2251 # bcf_ij2G() - common task: allele indexes to Number=G index (diploid)
jpayne@68 2252 # @i,j: allele indexes, 0-based, i<=j
jpayne@68 2253 # Returns index to the Number=G diploid array
jpayne@68 2254 uint32_t bcf_ij2G(uint32_t i, uint32_t j)
jpayne@68 2255
jpayne@68 2256
jpayne@68 2257 cdef extern from "htslib/cram.h" nogil:
jpayne@68 2258
jpayne@68 2259 enum cram_block_method:
jpayne@68 2260 ERROR
jpayne@68 2261 RAW
jpayne@68 2262 GZIP
jpayne@68 2263 BZIP2
jpayne@68 2264 LZMA
jpayne@68 2265 RANS
jpayne@68 2266 RANS0
jpayne@68 2267 RANS1
jpayne@68 2268 GZIP_RLE
jpayne@68 2269
jpayne@68 2270 enum cram_content_type:
jpayne@68 2271 CT_ERROR
jpayne@68 2272 FILE_HEADER
jpayne@68 2273 COMPRESSION_HEADER
jpayne@68 2274 MAPPED_SLICE
jpayne@68 2275 UNMAPPED_SLICE
jpayne@68 2276 EXTERNAL
jpayne@68 2277 CORE
jpayne@68 2278
jpayne@68 2279 # Opaque data types, see cram_structs for the fully fledged versions.
jpayne@68 2280 ctypedef struct SAM_hdr
jpayne@68 2281 ctypedef struct cram_file_def
jpayne@68 2282 ctypedef struct cram_fd
jpayne@68 2283 ctypedef struct cram_container
jpayne@68 2284 ctypedef struct cram_block
jpayne@68 2285 ctypedef struct cram_slice
jpayne@68 2286 ctypedef struct cram_metrics
jpayne@68 2287 ctypedef struct cram_block_slice_hdr
jpayne@68 2288 ctypedef struct cram_block_compression_hdr
jpayne@68 2289 ctypedef struct refs_t
jpayne@68 2290
jpayne@68 2291 # Accessor functions
jpayne@68 2292
jpayne@68 2293 #
jpayne@68 2294 #-----------------------------------------------------------------------------
jpayne@68 2295 # cram_fd
jpayne@68 2296 #
jpayne@68 2297 SAM_hdr *cram_fd_get_header(cram_fd *fd)
jpayne@68 2298 void cram_fd_set_header(cram_fd *fd, SAM_hdr *hdr)
jpayne@68 2299
jpayne@68 2300 int cram_fd_get_version(cram_fd *fd)
jpayne@68 2301 void cram_fd_set_version(cram_fd *fd, int vers)
jpayne@68 2302
jpayne@68 2303 int cram_major_vers(cram_fd *fd)
jpayne@68 2304 int cram_minor_vers(cram_fd *fd)
jpayne@68 2305
jpayne@68 2306 hFILE *cram_fd_get_fp(cram_fd *fd)
jpayne@68 2307 void cram_fd_set_fp(cram_fd *fd, hFILE *fp)
jpayne@68 2308
jpayne@68 2309 #
jpayne@68 2310 #-----------------------------------------------------------------------------
jpayne@68 2311 # cram_container
jpayne@68 2312 #
jpayne@68 2313 int32_t cram_container_get_length(cram_container *c)
jpayne@68 2314 void cram_container_set_length(cram_container *c, int32_t length)
jpayne@68 2315 int32_t cram_container_get_num_blocks(cram_container *c)
jpayne@68 2316 void cram_container_set_num_blocks(cram_container *c, int32_t num_blocks)
jpayne@68 2317 int32_t *cram_container_get_landmarks(cram_container *c, int32_t *num_landmarks)
jpayne@68 2318 void cram_container_set_landmarks(cram_container *c, int32_t num_landmarks,
jpayne@68 2319 int32_t *landmarks)
jpayne@68 2320
jpayne@68 2321 # Returns true if the container is empty (EOF marker) */
jpayne@68 2322 int cram_container_is_empty(cram_fd *fd)
jpayne@68 2323
jpayne@68 2324
jpayne@68 2325 #
jpayne@68 2326 #-----------------------------------------------------------------------------
jpayne@68 2327 # cram_block
jpayne@68 2328 #
jpayne@68 2329 int32_t cram_block_get_content_id(cram_block *b)
jpayne@68 2330 int32_t cram_block_get_comp_size(cram_block *b)
jpayne@68 2331 int32_t cram_block_get_uncomp_size(cram_block *b)
jpayne@68 2332 int32_t cram_block_get_crc32(cram_block *b)
jpayne@68 2333 void * cram_block_get_data(cram_block *b)
jpayne@68 2334
jpayne@68 2335 cram_content_type cram_block_get_content_type(cram_block *b)
jpayne@68 2336
jpayne@68 2337 void cram_block_set_content_id(cram_block *b, int32_t id)
jpayne@68 2338 void cram_block_set_comp_size(cram_block *b, int32_t size)
jpayne@68 2339 void cram_block_set_uncomp_size(cram_block *b, int32_t size)
jpayne@68 2340 void cram_block_set_crc32(cram_block *b, int32_t crc)
jpayne@68 2341 void cram_block_set_data(cram_block *b, void *data)
jpayne@68 2342
jpayne@68 2343 int cram_block_append(cram_block *b, void *data, int size)
jpayne@68 2344 void cram_block_update_size(cram_block *b)
jpayne@68 2345
jpayne@68 2346 # Offset is known as "size" internally, but it can be confusing.
jpayne@68 2347 size_t cram_block_get_offset(cram_block *b)
jpayne@68 2348 void cram_block_set_offset(cram_block *b, size_t offset)
jpayne@68 2349
jpayne@68 2350 #
jpayne@68 2351 # Computes the size of a cram block, including the block
jpayne@68 2352 # header itself.
jpayne@68 2353 #
jpayne@68 2354 uint32_t cram_block_size(cram_block *b)
jpayne@68 2355
jpayne@68 2356 #
jpayne@68 2357 # Renumbers RG numbers in a cram compression header.
jpayne@68 2358 #
jpayne@68 2359 # CRAM stores RG as the Nth number in the header, rather than a
jpayne@68 2360 # string holding the ID: tag. This is smaller in space, but means
jpayne@68 2361 # "samtools cat" to join files together that contain single but
jpayne@68 2362 # different RG lines needs a way of renumbering them.
jpayne@68 2363 #
jpayne@68 2364 # The file descriptor is expected to be immediately after the
jpayne@68 2365 # cram_container structure (ie before the cram compression header).
jpayne@68 2366 # Due to the nature of the CRAM format, this needs to read and write
jpayne@68 2367 # the blocks itself. Note that there may be multiple slices within
jpayne@68 2368 # the container, meaning multiple compression headers to manipulate.
jpayne@68 2369 # Changing RG may change the size of the compression header and
jpayne@68 2370 # therefore the length field in the container. Hence we rewrite all
jpayne@68 2371 # blocks just in case and also emit the adjusted container.
jpayne@68 2372 #
jpayne@68 2373 # The current implementation can only cope with renumbering a single
jpayne@68 2374 # RG (and only then if it is using HUFFMAN or BETA codecs). In
jpayne@68 2375 # theory it *may* be possible to renumber multiple RGs if they use
jpayne@68 2376 # HUFFMAN to the CORE block or use an external block unshared by any
jpayne@68 2377 # other data series. So we have an API that can be upgraded to
jpayne@68 2378 # support this, but do not implement it for now. An example
jpayne@68 2379 # implementation of RG as an EXTERNAL block would be to find that
jpayne@68 2380 # block and rewrite it, returning the number of blocks consumed.
jpayne@68 2381 #
jpayne@68 2382 # Returns 0 on success;
jpayne@68 2383 # -1 if unable to edit;
jpayne@68 2384 # -2 on other errors (eg I/O).
jpayne@68 2385 #
jpayne@68 2386 int cram_transcode_rg(cram_fd *input, cram_fd *output,
jpayne@68 2387 cram_container *c,
jpayne@68 2388 int nrg, int *in_rg, int *out_rg)
jpayne@68 2389
jpayne@68 2390 #
jpayne@68 2391 # Copies the blocks representing the next num_slice slices from a
jpayne@68 2392 # container from 'in' to 'out'. It is expected that the file pointer
jpayne@68 2393 # is just after the read of the cram_container and cram compression
jpayne@68 2394 # header.
jpayne@68 2395 #
jpayne@68 2396 # Returns 0 on success
jpayne@68 2397 # -1 on failure
jpayne@68 2398 #
jpayne@68 2399 int cram_copy_slice(cram_fd *input, cram_fd *output, int32_t num_slice)
jpayne@68 2400
jpayne@68 2401 #
jpayne@68 2402 #-----------------------------------------------------------------------------
jpayne@68 2403 # SAM_hdr
jpayne@68 2404 #
jpayne@68 2405
jpayne@68 2406 # Tokenises a SAM header into a hash table.
jpayne@68 2407 #
jpayne@68 2408 # Also extracts a few bits on specific data types, such as @RG lines.
jpayne@68 2409 #
jpayne@68 2410 # @return
jpayne@68 2411 # Returns a SAM_hdr struct on success (free with sam_hdr_free())
jpayne@68 2412 # NULL on failure
jpayne@68 2413 #
jpayne@68 2414 SAM_hdr *sam_hdr_parse_(const char *hdr, int len)
jpayne@68 2415
jpayne@68 2416
jpayne@68 2417 #
jpayne@68 2418 #-----------------------------------------------------------------------------
jpayne@68 2419 # cram_io basics
jpayne@68 2420 #
jpayne@68 2421
jpayne@68 2422 # CRAM blocks - the dynamically growable data block. We have code to
jpayne@68 2423 # create, update, (un)compress and read/write.
jpayne@68 2424 #
jpayne@68 2425 # These are derived from the deflate_interlaced.c blocks, but with the
jpayne@68 2426 # CRAM extension of content types and IDs.
jpayne@68 2427 #
jpayne@68 2428
jpayne@68 2429 # Allocates a new cram_block structure with a specified content_type and
jpayne@68 2430 # id.
jpayne@68 2431 #
jpayne@68 2432 # @return
jpayne@68 2433 # Returns block pointer on success;
jpayne@68 2434 # NULL on failure
jpayne@68 2435 #
jpayne@68 2436 cram_block *cram_new_block(cram_content_type content_type,
jpayne@68 2437 int content_id)
jpayne@68 2438
jpayne@68 2439 # Reads a block from a cram file.
jpayne@68 2440 #
jpayne@68 2441 # @return
jpayne@68 2442 # Returns cram_block pointer on success;
jpayne@68 2443 # NULL on failure
jpayne@68 2444 #
jpayne@68 2445 cram_block *cram_read_block(cram_fd *fd)
jpayne@68 2446
jpayne@68 2447 # Writes a CRAM block.
jpayne@68 2448 #
jpayne@68 2449 # @return
jpayne@68 2450 # Returns 0 on success;
jpayne@68 2451 # -1 on failure
jpayne@68 2452 #
jpayne@68 2453 int cram_write_block(cram_fd *fd, cram_block *b)
jpayne@68 2454
jpayne@68 2455 # Frees a CRAM block, deallocating internal data too.
jpayne@68 2456 #
jpayne@68 2457 void cram_free_block(cram_block *b)
jpayne@68 2458
jpayne@68 2459 # Uncompresses a CRAM block, if compressed.
jpayne@68 2460 #
jpayne@68 2461 # @return
jpayne@68 2462 # Returns 0 on success;
jpayne@68 2463 # -1 on failure
jpayne@68 2464 #
jpayne@68 2465 int cram_uncompress_block(cram_block *b)
jpayne@68 2466
jpayne@68 2467 # Compresses a block.
jpayne@68 2468 #
jpayne@68 2469 # Compresses a block using one of two different zlib strategies. If we only
jpayne@68 2470 # want one choice set strat2 to be -1.
jpayne@68 2471 #
jpayne@68 2472 # The logic here is that sometimes Z_RLE does a better job than Z_FILTERED
jpayne@68 2473 # or Z_DEFAULT_STRATEGY on quality data. If so, we'd rather use it as it is
jpayne@68 2474 # significantly faster.
jpayne@68 2475 #
jpayne@68 2476 # @return
jpayne@68 2477 # Returns 0 on success;
jpayne@68 2478 # -1 on failure
jpayne@68 2479 #
jpayne@68 2480 int cram_compress_block(cram_fd *fd, cram_block *b, cram_metrics *metrics,
jpayne@68 2481 int method, int level)
jpayne@68 2482
jpayne@68 2483 # Containers
jpayne@68 2484 #
jpayne@68 2485
jpayne@68 2486 # Creates a new container, specifying the maximum number of slices
jpayne@68 2487 # and records permitted.
jpayne@68 2488 #
jpayne@68 2489 # @return
jpayne@68 2490 # Returns cram_container ptr on success;
jpayne@68 2491 # NULL on failure
jpayne@68 2492 #
jpayne@68 2493 cram_container *cram_new_container(int nrec, int nslice)
jpayne@68 2494 void cram_free_container(cram_container *c)
jpayne@68 2495
jpayne@68 2496 # Reads a container header.
jpayne@68 2497 #
jpayne@68 2498 # @return
jpayne@68 2499 # Returns cram_container on success;
jpayne@68 2500 # NULL on failure or no container left (fd->err == 0).
jpayne@68 2501 #
jpayne@68 2502 cram_container *cram_read_container(cram_fd *fd)
jpayne@68 2503
jpayne@68 2504 # Writes a container structure.
jpayne@68 2505 #
jpayne@68 2506 # @return
jpayne@68 2507 # Returns 0 on success;
jpayne@68 2508 # -1 on failure
jpayne@68 2509 #
jpayne@68 2510 int cram_write_container(cram_fd *fd, cram_container *h)
jpayne@68 2511
jpayne@68 2512 #
jpayne@68 2513 # Stores the container structure in dat and returns *size as the
jpayne@68 2514 # number of bytes written to dat[]. The input size of dat is also
jpayne@68 2515 # held in *size and should be initialised to cram_container_size(c).
jpayne@68 2516 #
jpayne@68 2517 # Returns 0 on success;
jpayne@68 2518 # -1 on failure
jpayne@68 2519 #
jpayne@68 2520 int cram_store_container(cram_fd *fd, cram_container *c, char *dat, int *size)
jpayne@68 2521
jpayne@68 2522 int cram_container_size(cram_container *c)
jpayne@68 2523
jpayne@68 2524 # The top-level cram opening, closing and option handling
jpayne@68 2525 #
jpayne@68 2526
jpayne@68 2527 # Opens a CRAM file for read (mode "rb") or write ("wb").
jpayne@68 2528 #
jpayne@68 2529 # The filename may be "-" to indicate stdin or stdout.
jpayne@68 2530 #
jpayne@68 2531 # @return
jpayne@68 2532 # Returns file handle on success;
jpayne@68 2533 # NULL on failure.
jpayne@68 2534 #
jpayne@68 2535 cram_fd *cram_open(const char *filename, const char *mode)
jpayne@68 2536
jpayne@68 2537 # Opens an existing stream for reading or writing.
jpayne@68 2538 #
jpayne@68 2539 # @return
jpayne@68 2540 # Returns file handle on success;
jpayne@68 2541 # NULL on failure.
jpayne@68 2542 #
jpayne@68 2543 cram_fd *cram_dopen(hFILE *fp, const char *filename, const char *mode)
jpayne@68 2544
jpayne@68 2545 # Closes a CRAM file.
jpayne@68 2546 #
jpayne@68 2547 # @return
jpayne@68 2548 # Returns 0 on success;
jpayne@68 2549 # -1 on failure
jpayne@68 2550 #
jpayne@68 2551 int cram_close(cram_fd *fd)
jpayne@68 2552
jpayne@68 2553 #
jpayne@68 2554 # Seek within a CRAM file.
jpayne@68 2555 #
jpayne@68 2556 # Returns 0 on success
jpayne@68 2557 # -1 on failure
jpayne@68 2558 #
jpayne@68 2559 int cram_seek(cram_fd *fd, off_t offset, int whence)
jpayne@68 2560
jpayne@68 2561 #
jpayne@68 2562 # Flushes a CRAM file.
jpayne@68 2563 # Useful for when writing to stdout without wishing to close the stream.
jpayne@68 2564 #
jpayne@68 2565 # Returns 0 on success
jpayne@68 2566 # -1 on failure
jpayne@68 2567 #
jpayne@68 2568 int cram_flush(cram_fd *fd)
jpayne@68 2569
jpayne@68 2570 # Checks for end of file on a cram_fd stream.
jpayne@68 2571 #
jpayne@68 2572 # @return
jpayne@68 2573 # Returns 0 if not at end of file
jpayne@68 2574 # 1 if we hit an expected EOF (end of range or EOF block)
jpayne@68 2575 # 2 for other EOF (end of stream without EOF block)
jpayne@68 2576 #
jpayne@68 2577 int cram_eof(cram_fd *fd)
jpayne@68 2578
jpayne@68 2579 # Sets options on the cram_fd.
jpayne@68 2580 #
jpayne@68 2581 # See CRAM_OPT_* definitions in hts.h.
jpayne@68 2582 # Use this immediately after opening.
jpayne@68 2583 #
jpayne@68 2584 # @return
jpayne@68 2585 # Returns 0 on success;
jpayne@68 2586 # -1 on failure
jpayne@68 2587 #
jpayne@68 2588 int cram_set_option(cram_fd *fd, hts_fmt_option opt, ...)
jpayne@68 2589
jpayne@68 2590 # Sets options on the cram_fd.
jpayne@68 2591 #
jpayne@68 2592 # See CRAM_OPT_* definitions in hts.h.
jpayne@68 2593 # Use this immediately after opening.
jpayne@68 2594 #
jpayne@68 2595 # @return
jpayne@68 2596 # Returns 0 on success;
jpayne@68 2597 # -1 on failure
jpayne@68 2598 #
jpayne@68 2599 int cram_set_voption(cram_fd *fd, hts_fmt_option opt, va_list args)
jpayne@68 2600
jpayne@68 2601 #
jpayne@68 2602 # Attaches a header to a cram_fd.
jpayne@68 2603 #
jpayne@68 2604 # This should be used when creating a new cram_fd for writing where
jpayne@68 2605 # we have an SAM_hdr already constructed (eg from a file we've read
jpayne@68 2606 # in).
jpayne@68 2607 #
jpayne@68 2608 # @return
jpayne@68 2609 # Returns 0 on success;
jpayne@68 2610 # -1 on failure
jpayne@68 2611 #
jpayne@68 2612 int cram_set_header(cram_fd *fd, SAM_hdr *hdr)
jpayne@68 2613
jpayne@68 2614 # Check if this file has a proper EOF block
jpayne@68 2615 #
jpayne@68 2616 # @return
jpayne@68 2617 # Returns 3 if the file is a version of CRAM that does not contain EOF blocks
jpayne@68 2618 # 2 if the file is a stream and thus unseekable
jpayne@68 2619 # 1 if the file contains an EOF block
jpayne@68 2620 # 0 if the file does not contain an EOF block
jpayne@68 2621 # -1 if an error occurred whilst reading the file or we could not seek back to where we were
jpayne@68 2622 #
jpayne@68 2623 #
jpayne@68 2624 int cram_check_EOF(cram_fd *fd)
jpayne@68 2625
jpayne@68 2626 # As int32_decoded/encode, but from/to blocks instead of cram_fd */
jpayne@68 2627 int int32_put_blk(cram_block *b, int32_t val)
jpayne@68 2628
jpayne@68 2629 # Deallocates all storage used by a SAM_hdr struct.
jpayne@68 2630 #
jpayne@68 2631 # This also decrements the header reference count. If after decrementing
jpayne@68 2632 # it is still non-zero then the header is assumed to be in use by another
jpayne@68 2633 # caller and the free is not done.
jpayne@68 2634 #
jpayne@68 2635 # This is a synonym for sam_hdr_dec_ref().
jpayne@68 2636 #
jpayne@68 2637 void sam_hdr_free(SAM_hdr *hdr)
jpayne@68 2638
jpayne@68 2639 # Returns the current length of the SAM_hdr in text form.
jpayne@68 2640 #
jpayne@68 2641 # Call sam_hdr_rebuild() first if editing has taken place.
jpayne@68 2642 #
jpayne@68 2643 int sam_hdr_length(SAM_hdr *hdr)
jpayne@68 2644
jpayne@68 2645 # Returns the string form of the SAM_hdr.
jpayne@68 2646 #
jpayne@68 2647 # Call sam_hdr_rebuild() first if editing has taken place.
jpayne@68 2648 #
jpayne@68 2649 char *sam_hdr_str(SAM_hdr *hdr)
jpayne@68 2650
jpayne@68 2651 # Appends a formatted line to an existing SAM header.
jpayne@68 2652 #
jpayne@68 2653 # Line is a full SAM header record, eg "@SQ\tSN:foo\tLN:100", with
jpayne@68 2654 # optional new-line. If it contains more than 1 line then multiple lines
jpayne@68 2655 # will be added in order.
jpayne@68 2656 #
jpayne@68 2657 # Len is the length of the text data, or 0 if unknown (in which case
jpayne@68 2658 # it should be null terminated).
jpayne@68 2659 #
jpayne@68 2660 # @return
jpayne@68 2661 # Returns 0 on success;
jpayne@68 2662 # -1 on failure
jpayne@68 2663 #
jpayne@68 2664
jpayne@68 2665 # Add an @PG line.
jpayne@68 2666 #
jpayne@68 2667 # If we wish complete control over this use sam_hdr_add() directly. This
jpayne@68 2668 # function uses that, but attempts to do a lot of tedious house work for
jpayne@68 2669 # you too.
jpayne@68 2670 #
jpayne@68 2671 # - It will generate a suitable ID if the supplied one clashes.
jpayne@68 2672 # - It will generate multiple @PG records if we have multiple PG chains.
jpayne@68 2673 #
jpayne@68 2674 # Call it as per sam_hdr_add() with a series of key,value pairs ending
jpayne@68 2675 # in NULL.
jpayne@68 2676 #
jpayne@68 2677 # @return
jpayne@68 2678 # Returns 0 on success;
jpayne@68 2679 # -1 on failure
jpayne@68 2680 #
jpayne@68 2681 int sam_hdr_add_PG(SAM_hdr *sh, const char *name, ...)
jpayne@68 2682
jpayne@68 2683 #
jpayne@68 2684 # A function to help with construction of CL tags in @PG records.
jpayne@68 2685 # Takes an argc, argv pair and returns a single space-separated string.
jpayne@68 2686 # This string should be deallocated by the calling function.
jpayne@68 2687 #
jpayne@68 2688 # @return
jpayne@68 2689 # Returns malloced char * on success;
jpayne@68 2690 # NULL on failure
jpayne@68 2691 #
jpayne@68 2692 char *stringify_argv(int argc, char *argv[])
jpayne@68 2693
jpayne@68 2694 #
jpayne@68 2695 # Returns the refs_t structure used by a cram file handle.
jpayne@68 2696 #
jpayne@68 2697 # This may be used in conjunction with option CRAM_OPT_SHARED_REF to
jpayne@68 2698 # share reference memory between multiple file handles.
jpayne@68 2699 #
jpayne@68 2700 # @return
jpayne@68 2701 # Returns NULL if none exists or the file handle is not a CRAM file.
jpayne@68 2702 #
jpayne@68 2703 refs_t *cram_get_refs(htsFile *fd)
jpayne@68 2704
jpayne@68 2705
jpayne@68 2706 cdef class HTSFile(object):
jpayne@68 2707 cdef htsFile *htsfile # pointer to htsFile structure
jpayne@68 2708 cdef int64_t start_offset # BGZF offset of first record
jpayne@68 2709
jpayne@68 2710 cdef readonly object filename # filename as supplied by user
jpayne@68 2711 cdef readonly object mode # file opening mode
jpayne@68 2712 cdef readonly object threads # number of threads to use
jpayne@68 2713 cdef readonly object index_filename # filename of index, if supplied by user
jpayne@68 2714
jpayne@68 2715 cdef readonly bint is_stream # Is htsfile a non-seekable stream
jpayne@68 2716 cdef readonly bint is_remote # Is htsfile a remote stream
jpayne@68 2717 cdef readonly bint duplicate_filehandle # Duplicate filehandle when opening via fh
jpayne@68 2718
jpayne@68 2719 cdef htsFile *_open_htsfile(self) except? NULL