# cython: embedsignature=True # cython: profile=True ######################################################## ######################################################## # Cython wrapper for SAM/BAM/CRAM files based on htslib ######################################################## # The principal classes defined in this module are: # # class AlignmentFile read/write access to SAM/BAM/CRAM formatted files # # class AlignmentHeader manage SAM/BAM/CRAM header data # # class IndexedReads index a SAM/BAM/CRAM file by query name while keeping # the original sort order intact # # Additionally this module defines numerous additional classes that # are part of the internal API. These are: # # Various iterator classes to iterate over alignments in sequential # (IteratorRow) or in a stacked fashion (IteratorColumn): # # class IteratorRow # class IteratorRowRegion # class IteratorRowHead # class IteratorRowAll # class IteratorRowAllRefs # class IteratorRowSelection # class IteratorColumn # class IteratorColumnRegion # class IteratorColumnAll # class IteratorColumnAllRefs # ######################################################## # # The MIT License # # Copyright (c) 2015 Andreas Heger # # Permission is hereby granted, free of charge, to any person obtaining a # copy of this software and associated documentation files (the "Software"), # to deal in the Software without restriction, including without limitation # the rights to use, copy, modify, merge, publish, distribute, sublicense, # and/or sell copies of the Software, and to permit persons to whom the # Software is furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included in # all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER # DEALINGS IN THE SOFTWARE. # ######################################################## import os import collections try: from collections.abc import Sequence, Mapping # noqa except ImportError: from collections import Sequence, Mapping # noqa import re import warnings import array from libc.errno cimport errno, EPIPE from libc.string cimport strcmp, strpbrk, strerror from libc.stdint cimport INT32_MAX from cpython cimport array as c_array from cpython.version cimport PY_MAJOR_VERSION from pysam.libcutils cimport force_bytes, force_str, charptr_to_str from pysam.libcutils cimport encode_filename, from_string_and_size from pysam.libcalignedsegment cimport makeAlignedSegment, makePileupColumn from pysam.libchtslib cimport HTSFile, hisremote if PY_MAJOR_VERSION >= 3: from io import StringIO else: from StringIO import StringIO cimport cython __all__ = [ "AlignmentFile", "AlignmentHeader", "IteratorRow", "IteratorColumn", "IndexedReads"] IndexStats = collections.namedtuple("IndexStats", ("contig", "mapped", "unmapped", "total")) ######################################################## ## global variables # maximum genomic coordinace # for some reason, using 'int' causes overflow cdef int MAX_POS = (1 << 31) - 1 # valid types for SAM headers VALID_HEADER_TYPES = {"HD" : Mapping, "SQ" : Sequence, "RG" : Sequence, "PG" : Sequence, "CO" : Sequence} # order of records within SAM headers VALID_HEADERS = ("HD", "SQ", "RG", "PG", "CO") # default type conversions within SAM header records KNOWN_HEADER_FIELDS = {"HD" : {"VN" : str, "SO" : str, "GO" : str}, "SQ" : {"SN" : str, "LN" : int, "AS" : str, "M5" : str, "SP" : str, "UR" : str, "AH" : str,}, "RG" : {"ID" : str, "CN" : str, "DS" : str, "DT" : str, "FO" : str, "KS" : str, "LB" : str, "PG" : str, "PI" : str, "PL" : str, "PM" : str, "PU" : str, "SM" : str,}, "PG" : {"ID" : str, "PN" : str, "CL" : str, "PP" : str, "DS" : str, "VN" : str,},} # output order of fields within records. Ensure that CL is at # the end as parsing a CL will ignore any subsequent records. VALID_HEADER_ORDER = {"HD" : ("VN", "SO", "GO"), "SQ" : ("SN", "LN", "AS", "M5", "UR", "SP", "AH"), "RG" : ("ID", "CN", "SM", "LB", "PU", "PI", "DT", "DS", "PL", "FO", "KS", "PG", "PM"), "PG" : ("PN", "ID", "VN", "PP", "DS", "CL"),} def build_header_line(fields, record): '''build a header line from `fields` dictionary for `record`''' # TODO: add checking for field and sort order line = ["@%s" % record] # comment if record == "CO": line.append(fields) # user tags elif record.islower(): for key in sorted(fields): line.append("%s:%s" % (key, str(fields[key]))) # defined tags else: # write fields of the specification for key in VALID_HEADER_ORDER[record]: if key in fields: line.append("%s:%s" % (key, str(fields[key]))) # write user fields for key in fields: if not key.isupper(): line.append("%s:%s" % (key, str(fields[key]))) return "\t".join(line) cdef AlignmentHeader makeAlignmentHeader(bam_hdr_t *hdr): if not hdr: raise ValueError('cannot create AlignmentHeader, received NULL pointer') # check: is AlignmetHeader.__cinit__ called? cdef AlignmentHeader header = AlignmentHeader.__new__(AlignmentHeader) header.ptr = hdr return header def read_failure_reason(code): if code == -2: return 'truncated file' else: return "error {} while reading file".format(code) # the following should be class-method for VariantHeader, but cdef @classmethods # are not implemented in cython. cdef int fill_AlignmentHeader_from_list(bam_hdr_t *dest, reference_names, reference_lengths, add_sq_text=True, text=None) except -1: """build header from list of reference names and lengths. """ cdef class AlignmentHeader(object): """header information for a :class:`AlignmentFile` object Parameters ---------- header_dict : dict build header from a multi-level dictionary. The first level are the four types ('HD', 'SQ', ...). The second level are a list of lines, with each line being a list of tag-value pairs. The header is constructed first from all the defined fields, followed by user tags in alphabetical order. Alternatively, an :class:`~pysam.AlignmentHeader` object can be passed directly. text : string use the string provided as the header reference_names : list see reference_lengths reference_lengths : list build header from list of chromosome names and lengths. By default, 'SQ' and 'LN' tags will be added to the header text. This option can be changed by unsetting the flag `add_sq_text`. add_sq_text : bool do not add 'SQ' and 'LN' tags to header. This option permits construction :term:`SAM` formatted files without a header. """ # See makeVariantHeader for C constructor def __cinit__(self): self.ptr = NULL # Python constructor def __init__(self): self.ptr = bam_hdr_init() if self.ptr is NULL: raise MemoryError("could not create header") @classmethod def _from_text_and_lengths(cls, text, reference_names, reference_lengths): cdef AlignmentHeader self = AlignmentHeader() cdef char *ctext cdef int l_text cdef int n, x if text is not None: btext = force_bytes(text) ctext = btext l_text = len(btext) self.ptr.text = calloc(l_text + 1, sizeof(char)) if self.ptr.text == NULL: raise MemoryError("could not allocate {} bytes".format(l_text + 1), sizeof(char)) self.ptr.l_text = l_text memcpy(self.ptr.text, ctext, l_text + 1) if reference_names and reference_lengths: reference_names = [force_bytes(ref) for ref in reference_names] self.ptr.n_targets = len(reference_names) n = sum([len(reference_names) + 1]) self.ptr.target_name = calloc(n, sizeof(char*)) if self.ptr.target_name == NULL: raise MemoryError("could not allocate {} bytes".format(n, sizeof(char *))) self.ptr.target_len = calloc(n, sizeof(uint32_t)) if self.ptr.target_len == NULL: raise MemoryError("could not allocate {} bytes".format(n, sizeof(uint32_t))) for x from 0 <= x < self.ptr.n_targets: self.ptr.target_len[x] = reference_lengths[x] name = reference_names[x] self.ptr.target_name[x] = calloc(len(name) + 1, sizeof(char)) if self.ptr.target_name[x] == NULL: raise MemoryError("could not allocate {} bytes".format(len(name) + 1, sizeof(char))) strncpy(self.ptr.target_name[x], name, len(name)) return self @classmethod def from_text(cls, text): reference_names, reference_lengths = [], [] for line in text.splitlines(): if line.startswith("@SQ"): fields = dict([x.split(":", 1) for x in line.split("\t")[1:]]) try: reference_names.append(fields["SN"]) reference_lengths.append(int(fields["LN"])) except KeyError: raise KeyError("incomplete sequence information in '%s'" % str(fields)) except ValueError: raise ValueError("wrong sequence information in '%s'" % str(fields)) return cls._from_text_and_lengths(text, reference_names, reference_lengths) @classmethod def from_dict(cls, header_dict): cdef list lines = [] # first: defined tags for record in VALID_HEADERS: if record in header_dict: data = header_dict[record] if not isinstance(data, VALID_HEADER_TYPES[record]): raise ValueError( "invalid type for record %s: %s, expected %s".format( record, type(data), VALID_HEADER_TYPES[record])) if isinstance(data, Mapping): lines.append(build_header_line(data, record)) else: for fields in header_dict[record]: lines.append(build_header_line(fields, record)) # then: user tags (lower case), sorted alphabetically for record, data in sorted(header_dict.items()): if record in VALID_HEADERS: continue if isinstance(data, Mapping): lines.append(build_header_line(data, record)) else: for fields in header_dict[record]: lines.append(build_header_line(fields, record)) text = "\n".join(lines) + "\n" reference_names, reference_lengths = [], [] if "SQ" in header_dict: for fields in header_dict["SQ"]: try: reference_names.append(fields["SN"]) reference_lengths.append(fields["LN"]) except KeyError: raise KeyError("incomplete sequence information in '%s'" % str(fields)) return cls._from_text_and_lengths(text, reference_names, reference_lengths) @classmethod def from_references(cls, reference_names, reference_lengths, text=None, add_sq_text=True): if len(reference_names) != len(reference_lengths): raise ValueError("number of reference names and lengths do not match") # optionally, if there is no text, add a SAM compatible header to output file. if text is None and add_sq_text: text = "".join(["@SQ\tSN:{}\tLN:{}\n".format(x, y) for x, y in zip( reference_names, reference_lengths)]) return cls._from_text_and_lengths(text, reference_names, reference_lengths) def __dealloc__(self): bam_hdr_destroy(self.ptr) self.ptr = NULL def __bool__(self): return self.ptr != NULL def copy(self): return makeAlignmentHeader(bam_hdr_dup(self.ptr)) property nreferences: """int with the number of :term:`reference` sequences in the file. This is a read-only attribute.""" def __get__(self): return self.ptr.n_targets property references: """tuple with the names of :term:`reference` sequences. This is a read-only attribute""" def __get__(self): t = [] cdef int x for x in range(self.ptr.n_targets): t.append(charptr_to_str(self.ptr.target_name[x])) return tuple(t) property lengths: """tuple of the lengths of the :term:`reference` sequences. This is a read-only attribute. The lengths are in the same order as :attr:`pysam.AlignmentFile.references` """ def __get__(self): t = [] cdef int x for x in range(self.ptr.n_targets): t.append(self.ptr.target_len[x]) return tuple(t) def _build_sequence_section(self): """return sequence section of header. The sequence section is built from the list of reference names and lengths stored in the BAM-file and not from any @SQ entries that are part of the header's text section. """ cdef int x text = [] for x in range(self.ptr.n_targets): text.append("@SQ\tSN:{}\tLN:{}\n".format( force_str(self.ptr.target_name[x]), self.ptr.target_len[x])) return "".join(text) def to_dict(self): """return two-level dictionary with header information from the file. The first level contains the record (``HD``, ``SQ``, etc) and the second level contains the fields (``VN``, ``LN``, etc). The parser is validating and will raise an AssertionError if if encounters any record or field tags that are not part of the SAM specification. Use the :attr:`pysam.AlignmentFile.text` attribute to get the unparsed header. The parsing follows the SAM format specification with the exception of the ``CL`` field. This option will consume the rest of a header line irrespective of any additional fields. This behaviour has been added to accommodate command line options that contain characters that are not valid field separators. If no @SQ entries are within the text section of the header, this will be automatically added from the reference names and lengths stored in the binary part of the header. """ result = collections.OrderedDict() # convert to python string t = self.__str__() for line in t.split("\n"): line = line.strip(' \0') if not line: continue assert line.startswith("@"), \ "header line without '@': '%s'" % line fields = line[1:].split("\t") record = fields[0] assert record in VALID_HEADER_TYPES, \ "header line with invalid type '%s': '%s'" % (record, line) # treat comments if record == "CO": if record not in result: result[record] = [] result[record].append("\t".join( fields[1:])) continue # the following is clumsy as generators do not work? x = {} for idx, field in enumerate(fields[1:]): if ":" not in field: raise ValueError("malformatted header: no ':' in field" ) key, value = field.split(":", 1) if key in ("CL",): # special treatment for command line # statements (CL). These might contain # characters that are non-conformant with # the valid field separators in the SAM # header. Thus, in contravention to the # SAM API, consume the rest of the line. key, value = "\t".join(fields[idx+1:]).split(":", 1) x[key] = KNOWN_HEADER_FIELDS[record][key](value) break # interpret type of known header record tags, default to str x[key] = KNOWN_HEADER_FIELDS[record].get(key, str)(value) if VALID_HEADER_TYPES[record] == Mapping: if record in result: raise ValueError( "multiple '%s' lines are not permitted" % record) result[record] = x elif VALID_HEADER_TYPES[record] == Sequence: if record not in result: result[record] = [] result[record].append(x) # if there are no SQ lines in the header, add the # reference names from the information in the bam # file. # # Background: c-samtools keeps the textual part of the # header separate from the list of reference names and # lengths. Thus, if a header contains only SQ lines, # the SQ information is not part of the textual header # and thus are missing from the output. See issue 84. if "SQ" not in result: sq = [] for ref, length in zip(self.references, self.lengths): sq.append({'LN': length, 'SN': ref }) result["SQ"] = sq return result def as_dict(self): """deprecated: use :meth:`to_dict()`""" return self.to_dict() def get_reference_name(self, tid): if tid == -1: return None if not 0 <= tid < self.ptr.n_targets: raise ValueError("reference_id %i out of range 0<=tid<%i" % (tid, self.ptr.n_targets)) return charptr_to_str(self.ptr.target_name[tid]) def get_reference_length(self, reference): cdef int tid = self.get_tid(reference) if tid < 0: raise KeyError("unknown reference {}".format(reference)) else: return self.ptr.target_len[tid] def is_valid_tid(self, int tid): """ return True if the numerical :term:`tid` is valid; False otherwise. Note that the unmapped tid code (-1) counts as an invalid. """ return 0 <= tid < self.ptr.n_targets def get_tid(self, reference): """ return the numerical :term:`tid` corresponding to :term:`reference` returns -1 if reference is not known. """ reference = force_bytes(reference) tid = bam_name2id(self.ptr, reference) if tid < -1: raise ValueError('could not parse header') return tid def __str__(self): '''string with the full contents of the :term:`sam file` header as a string. If no @SQ entries are within the text section of the header, this will be automatically added from the reference names and lengths stored in the binary part of the header. See :attr:`pysam.AlignmentFile.header.to_dict()` to get a parsed representation of the header. ''' text = from_string_and_size(self.ptr.text, self.ptr.l_text) if "@SQ" not in text: text += "\n" + self._build_sequence_section() return text # dictionary access methods, for backwards compatibility. def __setitem__(self, key, value): raise TypeError("AlignmentHeader does not support item assignment (use header.to_dict()") def __getitem__(self, key): return self.to_dict().__getitem__(key) def items(self): return self.to_dict().items() # PY2 compatibility def iteritems(self): return self.to_dict().items() def keys(self): return self.to_dict().keys() def values(self): return self.to_dict().values() def get(self, *args): return self.to_dict().get(*args) def __len__(self): return self.to_dict().__len__() def __contains__(self, key): return self.to_dict().__contains__(key) cdef class AlignmentFile(HTSFile): """AlignmentFile(filepath_or_object, mode=None, template=None, reference_names=None, reference_lengths=None, text=NULL, header=None, add_sq_text=False, check_header=True, check_sq=True, reference_filename=None, filename=None, index_filename=None, filepath_index=None, require_index=False, duplicate_filehandle=True, ignore_truncation=False, threads=1) A :term:`SAM`/:term:`BAM`/:term:`CRAM` formatted file. If `filepath_or_object` is a string, the file is automatically opened. If `filepath_or_object` is a python File object, the already opened file will be used. If the file is opened for reading and an index exists (if file is BAM, a .bai file or if CRAM a .crai file), it will be opened automatically. `index_filename` may be specified explicitly. If the index is not named in the standard manner, not located in the same directory as the BAM/CRAM file, or is remote. Without an index, random access via :meth:`~pysam.AlignmentFile.fetch` and :meth:`~pysam.AlignmentFile.pileup` is disabled. For writing, the header of a :term:`SAM` file/:term:`BAM` file can be constituted from several sources (see also the samtools format specification): 1. If `template` is given, the header is copied from a another `AlignmentFile` (`template` must be a :class:`~pysam.AlignmentFile`). 2. If `header` is given, the header is built from a multi-level dictionary. 3. If `text` is given, new header text is copied from raw text. 4. The names (`reference_names`) and lengths (`reference_lengths`) are supplied directly as lists. When reading or writing a CRAM file, the filename of a FASTA-formatted reference can be specified with `reference_filename`. By default, if a file is opened in mode 'r', it is checked for a valid header (`check_header` = True) and a definition of chromosome names (`check_sq` = True). Parameters ---------- mode : string `mode` should be ``r`` for reading or ``w`` for writing. The default is text mode (:term:`SAM`). For binary (:term:`BAM`) I/O you should append ``b`` for compressed or ``u`` for uncompressed :term:`BAM` output. Use ``h`` to output header information in text (:term:`TAM`) mode. Use ``c`` for :term:`CRAM` formatted files. If ``b`` is present, it must immediately follow ``r`` or ``w``. Valid modes are ``r``, ``w``, ``wh``, ``rb``, ``wb``, ``wbu``, ``wb0``, ``rc`` and ``wc``. For instance, to open a :term:`BAM` formatted file for reading, type:: f = pysam.AlignmentFile('ex1.bam','rb') If mode is not specified, the method will try to auto-detect in the order 'rb', 'r', thus both the following should work:: f1 = pysam.AlignmentFile('ex1.bam') f2 = pysam.AlignmentFile('ex1.sam') template : AlignmentFile when writing, copy header from file `template`. header : dict or AlignmentHeader when writing, build header from a multi-level dictionary. The first level are the four types ('HD', 'SQ', ...). The second level are a list of lines, with each line being a list of tag-value pairs. The header is constructed first from all the defined fields, followed by user tags in alphabetical order. Alternatively, an :class:`~pysam.AlignmentHeader` object can be passed directly. text : string when writing, use the string provided as the header reference_names : list see reference_lengths reference_lengths : list when writing or opening a SAM file without header build header from list of chromosome names and lengths. By default, 'SQ' and 'LN' tags will be added to the header text. This option can be changed by unsetting the flag `add_sq_text`. add_sq_text : bool do not add 'SQ' and 'LN' tags to header. This option permits construction :term:`SAM` formatted files without a header. add_sam_header : bool when outputting SAM the default is to output a header. This is equivalent to opening the file in 'wh' mode. If this option is set to False, no header will be output. To read such a file, set `check_header=False`. check_header : bool obsolete: when reading a SAM file, check if header is present (default=True) check_sq : bool when reading, check if SQ entries are present in header (default=True) reference_filename : string Path to a FASTA-formatted reference file. Valid only for CRAM files. When reading a CRAM file, this overrides both ``$REF_PATH`` and the URL specified in the header (``UR`` tag), which are normally used to find the reference. index_filename : string Explicit path to the index file. Only needed if the index is not named in the standard manner, not located in the same directory as the BAM/CRAM file, or is remote. An IOError is raised if the index cannot be found or is invalid. filepath_index : string Alias for `index_filename`. require_index : bool When reading, require that an index file is present and is valid or raise an IOError. (default=False) filename : string Alternative to filepath_or_object. Filename of the file to be opened. duplicate_filehandle: bool By default, file handles passed either directly or through File-like objects will be duplicated before passing them to htslib. The duplication prevents issues where the same stream will be closed by htslib and through destruction of the high-level python object. Set to False to turn off duplication. ignore_truncation: bool Issue a warning, instead of raising an error if the current file appears to be truncated due to a missing EOF marker. Only applies to bgzipped formats. (Default=False) format_options: list A list of key=value strings, as accepted by --input-fmt-option and --output-fmt-option in samtools. threads: integer Number of threads to use for compressing/decompressing BAM/CRAM files. Setting threads to > 1 cannot be combined with `ignore_truncation`. (Default=1) """ def __cinit__(self, *args, **kwargs): self.htsfile = NULL self.filename = None self.mode = None self.threads = 1 self.is_stream = False self.is_remote = False self.index = NULL if "filename" in kwargs: args = [kwargs["filename"]] del kwargs["filename"] self._open(*args, **kwargs) # allocate memory for iterator self.b = calloc(1, sizeof(bam1_t)) if self.b == NULL: raise MemoryError("could not allocate memory of size {}".format(sizeof(bam1_t))) def has_index(self): """return true if htsfile has an existing (and opened) index. """ return self.index != NULL def check_index(self): """return True if index is present. Raises ------ AttributeError if htsfile is :term:`SAM` formatted and thus has no index. ValueError if htsfile is closed or index could not be opened. """ if not self.is_open: raise ValueError("I/O operation on closed file") if not self.is_bam and not self.is_cram: raise AttributeError( "AlignmentFile.mapped only available in bam files") if self.index == NULL: raise ValueError( "mapping information not recorded in index " "or index not available") return True def _open(self, filepath_or_object, mode=None, AlignmentFile template=None, reference_names=None, reference_lengths=None, reference_filename=None, text=None, header=None, port=None, add_sq_text=True, add_sam_header=True, check_header=True, check_sq=True, index_filename=None, filepath_index=None, require_index=False, referencenames=None, referencelengths=None, duplicate_filehandle=True, ignore_truncation=False, format_options=None, threads=1): '''open a sam, bam or cram formatted file. If _open is called on an existing file, the current file will be closed and a new file will be opened. ''' cdef char *cfilename = NULL cdef char *creference_filename = NULL cdef char *cindexname = NULL cdef char *cmode = NULL cdef bam_hdr_t * hdr = NULL if threads > 1 and ignore_truncation: # This won't raise errors if reaching a truncated alignment, # because bgzf_mt_reader in htslib does not deal with # bgzf_mt_read_block returning non-zero values, contrary # to bgzf_read (https://github.com/samtools/htslib/blob/1.7/bgzf.c#L888) # Better to avoid this (for now) than to produce seemingly correct results. raise ValueError('Cannot add extra threads when "ignore_truncation" is True') self.threads = threads # for backwards compatibility: if referencenames is not None: reference_names = referencenames if referencelengths is not None: reference_lengths = referencelengths # close a previously opened file if self.is_open: self.close() # autodetection for read if mode is None: mode = "r" if add_sam_header and mode == "w": mode = "wh" assert mode in ("r", "w", "rb", "wb", "wh", "wbu", "rU", "wb0", "rc", "wc"), \ "invalid file opening mode `%s`" % mode self.duplicate_filehandle = duplicate_filehandle # StringIO not supported if isinstance(filepath_or_object, StringIO): raise NotImplementedError( "access from StringIO objects not supported") # reading from a file descriptor elif isinstance(filepath_or_object, int): self.filename = filepath_or_object filename = None self.is_remote = False self.is_stream = True # reading from a File object or other object with fileno elif hasattr(filepath_or_object, "fileno"): if filepath_or_object.closed: raise ValueError('I/O operation on closed file') self.filename = filepath_or_object # .name can be TextIOWrapper try: filename = encode_filename(str(filepath_or_object.name)) cfilename = filename except AttributeError: filename = None self.is_remote = False self.is_stream = True # what remains is a filename else: self.filename = filename = encode_filename(filepath_or_object) cfilename = filename self.is_remote = hisremote(cfilename) self.is_stream = self.filename == b'-' # for htslib, wbu seems to not work if mode == "wbu": mode = "wb0" self.mode = force_bytes(mode) self.reference_filename = reference_filename = encode_filename( reference_filename) if mode[0] == 'w': # open file for writing if not (template or header or text or (reference_names and reference_lengths)): raise ValueError( "either supply options `template`, `header`, `text` or both `reference_names` " "and `reference_lengths` for writing") if template: # header is copied, though at the moment not strictly # necessary as AlignmentHeader is immutable. self.header = template.header.copy() elif isinstance(header, AlignmentHeader): self.header = header.copy() elif isinstance(header, Mapping): self.header = AlignmentHeader.from_dict(header) elif reference_names and reference_lengths: self.header = AlignmentHeader.from_references( reference_names, reference_lengths, add_sq_text=add_sq_text, text=text) elif text: self.header = AlignmentHeader.from_text(text) else: raise ValueError("not enough information to construct header. Please provide template, " "header, text or reference_names/reference_lengths") self.htsfile = self._open_htsfile() if self.htsfile == NULL: if errno: raise IOError(errno, "could not open alignment file `{}`: {}".format( force_str(filename), force_str(strerror(errno)))) else: raise ValueError("could not open alignment file `{}`".format(force_str(filename))) if format_options and len(format_options): self.add_hts_options(format_options) # set filename with reference sequences. If no filename # is given, the CRAM reference arrays will be built from # the @SQ header in the header if "c" in mode and reference_filename: if (hts_set_fai_filename(self.htsfile, self.reference_filename) != 0): raise ValueError("failure when setting reference filename") # write header to htsfile if "b" in mode or "c" in mode or "h" in mode: hdr = self.header.ptr with nogil: sam_hdr_write(self.htsfile, hdr) elif mode[0] == "r": # open file for reading self.htsfile = self._open_htsfile() if self.htsfile == NULL: if errno: raise IOError(errno, "could not open alignment file `{}`: {}".format(force_str(filename), force_str(strerror(errno)))) else: raise ValueError("could not open alignment file `{}`".format(force_str(filename))) if self.htsfile.format.category != sequence_data: raise ValueError("file does not contain alignment data") if format_options and len(format_options): self.add_hts_options(format_options) self.check_truncation(ignore_truncation) # bam/cram files require a valid header if self.is_bam or self.is_cram: with nogil: hdr = sam_hdr_read(self.htsfile) if hdr == NULL: raise ValueError( "file does not have a valid header (mode='%s') " "- is it BAM/CRAM format?" % mode) self.header = makeAlignmentHeader(hdr) else: # in sam files a header is optional. If not given, # user may provide reference names and lengths to built # an on-the-fly header. if reference_names and reference_lengths: # build header from a target names and lengths self.header = AlignmentHeader.from_references( reference_names=reference_names, reference_lengths=reference_lengths, add_sq_text=add_sq_text, text=text) else: with nogil: hdr = sam_hdr_read(self.htsfile) if hdr == NULL: raise ValueError( "SAM? file does not have a valid header (mode='%s'), " "please provide reference_names and reference_lengths") self.header = makeAlignmentHeader(hdr) # set filename with reference sequences if self.is_cram and reference_filename: creference_filename = self.reference_filename hts_set_opt(self.htsfile, CRAM_OPT_REFERENCE, creference_filename) if check_sq and self.header.nreferences == 0: raise ValueError( ("file has no sequences defined (mode='%s') - " "is it SAM/BAM format? Consider opening with " "check_sq=False") % mode) if self.is_bam or self.is_cram: self.index_filename = index_filename or filepath_index if self.index_filename: cindexname = bfile_name = encode_filename(self.index_filename) if cfilename or cindexname: with nogil: self.index = sam_index_load2(self.htsfile, cfilename, cindexname) if not self.index and (cindexname or require_index): if errno: raise IOError(errno, force_str(strerror(errno))) else: raise IOError('unable to open index file `%s`' % self.index_filename) elif require_index: raise IOError('unable to open index file') # save start of data section if not self.is_stream: self.start_offset = self.tell() def fetch(self, contig=None, start=None, stop=None, region=None, tid=None, until_eof=False, multiple_iterators=False, reference=None, end=None): """fetch reads aligned in a :term:`region`. See :meth:`~pysam.HTSFile.parse_region` for more information on how genomic regions can be specified. :term:`reference` and `end` are also accepted for backward compatibility as synonyms for :term:`contig` and `stop`, respectively. Without a `contig` or `region` all mapped reads in the file will be fetched. The reads will be returned ordered by reference sequence, which will not necessarily be the order within the file. This mode of iteration still requires an index. If there is no index, use `until_eof=True`. If only `contig` is set, all reads aligned to `contig` will be fetched. A :term:`SAM` file does not allow random access. If `region` or `contig` are given, an exception is raised. Parameters ---------- until_eof : bool If `until_eof` is True, all reads from the current file position will be returned in order as they are within the file. Using this option will also fetch unmapped reads. multiple_iterators : bool If `multiple_iterators` is True, multiple iterators on the same file can be used at the same time. The iterator returned will receive its own copy of a filehandle to the file effectively re-opening the file. Re-opening a file creates some overhead, so beware. Returns ------- An iterator over a collection of reads. Raises ------ ValueError if the genomic coordinates are out of range or invalid or the file does not permit random access to genomic coordinates. """ cdef int rtid, rstart, rstop, has_coord if not self.is_open: raise ValueError( "I/O operation on closed file" ) has_coord, rtid, rstart, rstop = self.parse_region( contig, start, stop, region, tid, end=end, reference=reference) # Turn of re-opening if htsfile is a stream if self.is_stream: multiple_iterators = False if self.is_bam or self.is_cram: if not until_eof and not self.is_remote: if not self.has_index(): raise ValueError( "fetch called on bamfile without index") if has_coord: return IteratorRowRegion( self, rtid, rstart, rstop, multiple_iterators=multiple_iterators) else: if until_eof: return IteratorRowAll( self, multiple_iterators=multiple_iterators) else: # AH: check - reason why no multiple_iterators for # AllRefs? return IteratorRowAllRefs( self, multiple_iterators=multiple_iterators) else: if has_coord: raise ValueError( "fetching by region is not available for SAM files") if multiple_iterators == True: raise ValueError( "multiple iterators not implemented for SAM files") return IteratorRowAll(self, multiple_iterators=multiple_iterators) def head(self, n, multiple_iterators=True): '''return an iterator over the first n alignments. This iterator is is useful for inspecting the bam-file. Parameters ---------- multiple_iterators : bool is set to True by default in order to avoid changing the current file position. Returns ------- an iterator over a collection of reads ''' return IteratorRowHead(self, n, multiple_iterators=multiple_iterators) def mate(self, AlignedSegment read): '''return the mate of :class:`~pysam.AlignedSegment` `read`. .. note:: Calling this method will change the file position. This might interfere with any iterators that have not re-opened the file. .. note:: This method is too slow for high-throughput processing. If a read needs to be processed with its mate, work from a read name sorted file or, better, cache reads. Returns ------- :class:`~pysam.AlignedSegment` : the mate Raises ------ ValueError if the read is unpaired or the mate is unmapped ''' cdef uint32_t flag = read._delegate.core.flag if flag & BAM_FPAIRED == 0: raise ValueError("read %s: is unpaired" % (read.query_name)) if flag & BAM_FMUNMAP != 0: raise ValueError("mate %s: is unmapped" % (read.query_name)) # xor flags to get the other mate cdef int x = BAM_FREAD1 + BAM_FREAD2 flag = (flag ^ x) & x # Make sure to use a separate file to jump around # to mate as otherwise the original file position # will be lost # The following code is not using the C API and # could thus be made much quicker, for example # by using tell and seek. for mate in self.fetch( read._delegate.core.mpos, read._delegate.core.mpos + 1, tid=read._delegate.core.mtid, multiple_iterators=True): if mate.flag & flag != 0 and \ mate.query_name == read.query_name: break else: raise ValueError("mate not found") return mate def pileup(self, contig=None, start=None, stop=None, region=None, reference=None, end=None, **kwargs): """perform a :term:`pileup` within a :term:`region`. The region is specified by :term:`contig`, `start` and `stop` (using 0-based indexing). :term:`reference` and `end` are also accepted for backward compatibility as synonyms for :term:`contig` and `stop`, respectively. Alternatively, a samtools 'region' string can be supplied. Without 'contig' or 'region' all reads will be used for the pileup. The reads will be returned ordered by :term:`contig` sequence, which will not necessarily be the order within the file. Note that :term:`SAM` formatted files do not allow random access. In these files, if a 'region' or 'contig' are given an exception is raised. .. note:: 'all' reads which overlap the region are returned. The first base returned will be the first base of the first read 'not' necessarily the first base of the region used in the query. Parameters ---------- truncate : bool By default, the samtools pileup engine outputs all reads overlapping a region. If truncate is True and a region is given, only columns in the exact region specified are returned. max_depth : int Maximum read depth permitted. The default limit is '8000'. stepper : string The stepper controls how the iterator advances. Possible options for the stepper are ``all`` skip reads in which any of the following flags are set: BAM_FUNMAP, BAM_FSECONDARY, BAM_FQCFAIL, BAM_FDUP ``nofilter`` uses every single read turning off any filtering. ``samtools`` same filter and read processing as in :term:`csamtools` pileup. For full compatibility, this requires a 'fastafile' to be given. The following options all pertain to filtering of the ``samtools`` stepper. fastafile : :class:`~pysam.FastaFile` object. This is required for some of the steppers. ignore_overlaps: bool If set to True, detect if read pairs overlap and only take the higher quality base. This is the default. flag_filter : int ignore reads where any of the bits in the flag are set. The default is BAM_FUNMAP | BAM_FSECONDARY | BAM_FQCFAIL | BAM_FDUP. flag_require : int only use reads where certain flags are set. The default is 0. ignore_orphans: bool ignore orphans (paired reads that are not in a proper pair). The default is to ignore orphans. min_base_quality: int Minimum base quality. Bases below the minimum quality will not be output. The default is 13. adjust_capq_threshold: int adjust mapping quality. The default is 0 for no adjustment. The recommended value for adjustment is 50. min_mapping_quality : int only use reads above a minimum mapping quality. The default is 0. compute_baq: bool re-alignment computing per-Base Alignment Qualities (BAQ). The default is to do re-alignment. Realignment requires a reference sequence. If none is present, no realignment will be performed. redo_baq: bool recompute per-Base Alignment Quality on the fly ignoring existing base qualities. The default is False (use existing base qualities). Returns ------- an iterator over genomic positions. """ cdef int rtid, has_coord cdef int32_t rstart, rstop if not self.is_open: raise ValueError("I/O operation on closed file") has_coord, rtid, rstart, rstop = self.parse_region( contig, start, stop, region, reference=reference, end=end) if has_coord: if not self.has_index(): raise ValueError("no index available for pileup") return IteratorColumnRegion(self, tid=rtid, start=rstart, stop=rstop, **kwargs) else: if self.has_index(): return IteratorColumnAllRefs(self, **kwargs) else: return IteratorColumnAll(self, **kwargs) def count(self, contig=None, start=None, stop=None, region=None, until_eof=False, read_callback="nofilter", reference=None, end=None): '''count the number of reads in :term:`region` The region is specified by :term:`contig`, `start` and `stop`. :term:`reference` and `end` are also accepted for backward compatibility as synonyms for :term:`contig` and `stop`, respectively. Alternatively, a :term:`samtools` :term:`region` string can be supplied. A :term:`SAM` file does not allow random access and if `region` or `contig` are given, an exception is raised. Parameters ---------- contig : string reference_name of the genomic region (chromosome) start : int start of the genomic region (0-based inclusive) stop : int end of the genomic region (0-based exclusive) region : string a region string in samtools format. until_eof : bool count until the end of the file, possibly including unmapped reads as well. read_callback: string or function select a call-back to ignore reads when counting. It can be either a string with the following values: ``all`` skip reads in which any of the following flags are set: BAM_FUNMAP, BAM_FSECONDARY, BAM_FQCFAIL, BAM_FDUP ``nofilter`` uses every single read Alternatively, `read_callback` can be a function ``check_read(read)`` that should return True only for those reads that shall be included in the counting. reference : string backward compatible synonym for `contig` end : int backward compatible synonym for `stop` Raises ------ ValueError if the genomic coordinates are out of range or invalid. ''' cdef AlignedSegment read cdef long counter = 0 if not self.is_open: raise ValueError("I/O operation on closed file") cdef int filter_method = 0 if read_callback == "all": filter_method = 1 elif read_callback == "nofilter": filter_method = 2 for read in self.fetch(contig=contig, start=start, stop=stop, reference=reference, end=end, region=region, until_eof=until_eof): # apply filter if filter_method == 1: # filter = "all" if (read.flag & (0x4 | 0x100 | 0x200 | 0x400)): continue elif filter_method == 2: # filter = "nofilter" pass else: if not read_callback(read): continue counter += 1 return counter @cython.boundscheck(False) # we do manual bounds checking def count_coverage(self, contig, start=None, stop=None, region=None, quality_threshold=15, read_callback='all', reference=None, end=None): """count the coverage of genomic positions by reads in :term:`region`. The region is specified by :term:`contig`, `start` and `stop`. :term:`reference` and `end` are also accepted for backward compatibility as synonyms for :term:`contig` and `stop`, respectively. Alternatively, a :term:`samtools` :term:`region` string can be supplied. The coverage is computed per-base [ACGT]. Parameters ---------- contig : string reference_name of the genomic region (chromosome) start : int start of the genomic region (0-based inclusive). If not given, count from the start of the chromosome. stop : int end of the genomic region (0-based exclusive). If not given, count to the end of the chromosome. region : string a region string. quality_threshold : int quality_threshold is the minimum quality score (in phred) a base has to reach to be counted. read_callback: string or function select a call-back to ignore reads when counting. It can be either a string with the following values: ``all`` skip reads in which any of the following flags are set: BAM_FUNMAP, BAM_FSECONDARY, BAM_FQCFAIL, BAM_FDUP ``nofilter`` uses every single read Alternatively, `read_callback` can be a function ``check_read(read)`` that should return True only for those reads that shall be included in the counting. reference : string backward compatible synonym for `contig` end : int backward compatible synonym for `stop` Raises ------ ValueError if the genomic coordinates are out of range or invalid. Returns ------- four array.arrays of the same length in order A C G T : tuple """ cdef uint32_t contig_length = self.get_reference_length(contig) cdef int _start = start if start is not None else 0 cdef int _stop = stop if stop is not None else contig_length _stop = _stop if _stop < contig_length else contig_length if _stop == _start: raise ValueError("interval of size 0") if _stop < _start: raise ValueError("interval of size less than 0") cdef int length = _stop - _start cdef c_array.array int_array_template = array.array('L', []) cdef c_array.array count_a cdef c_array.array count_c cdef c_array.array count_g cdef c_array.array count_t count_a = c_array.clone(int_array_template, length, zero=True) count_c = c_array.clone(int_array_template, length, zero=True) count_g = c_array.clone(int_array_template, length, zero=True) count_t = c_array.clone(int_array_template, length, zero=True) cdef AlignedSegment read cdef cython.str seq cdef c_array.array quality cdef int qpos cdef int refpos cdef int c = 0 cdef int filter_method = 0 if read_callback == "all": filter_method = 1 elif read_callback == "nofilter": filter_method = 2 cdef int _threshold = quality_threshold or 0 for read in self.fetch(contig=contig, reference=reference, start=start, stop=stop, end=end, region=region): # apply filter if filter_method == 1: # filter = "all" if (read.flag & (0x4 | 0x100 | 0x200 | 0x400)): continue elif filter_method == 2: # filter = "nofilter" pass else: if not read_callback(read): continue # count seq = read.seq if seq is None: continue quality = read.query_qualities for qpos, refpos in read.get_aligned_pairs(True): if qpos is not None and refpos is not None and \ _start <= refpos < _stop: # only check base quality if _threshold > 0 if (_threshold and quality and quality[qpos] >= _threshold) or not _threshold: if seq[qpos] == 'A': count_a.data.as_ulongs[refpos - _start] += 1 if seq[qpos] == 'C': count_c.data.as_ulongs[refpos - _start] += 1 if seq[qpos] == 'G': count_g.data.as_ulongs[refpos - _start] += 1 if seq[qpos] == 'T': count_t.data.as_ulongs[refpos - _start] += 1 return count_a, count_c, count_g, count_t def find_introns_slow(self, read_iterator): """Return a dictionary {(start, stop): count} Listing the intronic sites in the reads (identified by 'N' in the cigar strings), and their support ( = number of reads ). read_iterator can be the result of a .fetch(...) call. Or it can be a generator filtering such reads. Example samfile.find_introns((read for read in samfile.fetch(...) if read.is_reverse) """ res = collections.Counter() for r in read_iterator: if 'N' in r.cigarstring: last_read_pos = False for read_loc, genome_loc in r.get_aligned_pairs(): if read_loc is None and last_read_pos: start = genome_loc elif read_loc and last_read_pos is None: stop = genome_loc # we are right exclusive ,so this is correct res[(start, stop)] += 1 del start del stop last_read_pos = read_loc return res def find_introns(self, read_iterator): """Return a dictionary {(start, stop): count} Listing the intronic sites in the reads (identified by 'N' in the cigar strings), and their support ( = number of reads ). read_iterator can be the result of a .fetch(...) call. Or it can be a generator filtering such reads. Example samfile.find_introns((read for read in samfile.fetch(...) if read.is_reverse) """ cdef: uint32_t base_position, junc_start, nt int op AlignedSegment r int BAM_CREF_SKIP = 3 #BAM_CREF_SKIP res = collections.Counter() match_or_deletion = {0, 2, 7, 8} # only M/=/X (0/7/8) and D (2) are related to genome position for r in read_iterator: base_position = r.pos for op, nt in r.cigartuples: if op in match_or_deletion: base_position += nt elif op == BAM_CREF_SKIP: junc_start = base_position base_position += nt res[(junc_start, base_position)] += 1 return res def close(self): '''closes the :class:`pysam.AlignmentFile`.''' if self.htsfile == NULL: return cdef int ret = hts_close(self.htsfile) self.htsfile = NULL if self.index != NULL: hts_idx_destroy(self.index) self.index = NULL self.header = None if ret < 0: global errno if errno == EPIPE: errno = 0 else: raise IOError(errno, force_str(strerror(errno))) def __dealloc__(self): cdef int ret = 0 if self.htsfile != NULL: ret = hts_close(self.htsfile) self.htsfile = NULL if self.index != NULL: hts_idx_destroy(self.index) self.index = NULL self.header = None if self.b: bam_destroy1(self.b) self.b = NULL if ret < 0: global errno if errno == EPIPE: errno = 0 else: raise IOError(errno, force_str(strerror(errno))) cpdef int write(self, AlignedSegment read) except -1: ''' write a single :class:`pysam.AlignedSegment` to disk. Raises: ValueError if the writing failed Returns: int : the number of bytes written. If the file is closed, this will be 0. ''' if not self.is_open: return 0 if self.header.ptr.n_targets <= read._delegate.core.tid: raise ValueError( "AlignedSegment refers to reference number {} that " "is larger than the number of references ({}) in the header".format( read._delegate.core.tid, self.header.ptr.n_targets)) cdef int ret with nogil: ret = sam_write1(self.htsfile, self.header.ptr, read._delegate) # kbj: Still need to raise an exception with except -1. Otherwise # when ret == -1 we get a "SystemError: error return without # exception set". if ret < 0: raise IOError( "sam_write1 failed with error code {}".format(ret)) return ret # context manager interface def __enter__(self): return self def __exit__(self, exc_type, exc_value, traceback): self.close() return False ############################################################### ############################################################### ############################################################### ## properties ############################################################### property mapped: """int with total number of mapped alignments according to the statistics recorded in the index. This is a read-only attribute. (This will be 0 for a CRAM file indexed by a .crai index, as that index format does not record these statistics.) """ def __get__(self): self.check_index() cdef int tid cdef uint64_t total = 0 cdef uint64_t mapped, unmapped for tid from 0 <= tid < self.header.nreferences: with nogil: hts_idx_get_stat(self.index, tid, &mapped, &unmapped) total += mapped return total property unmapped: """int with total number of unmapped reads according to the statistics recorded in the index. This number of reads includes the number of reads without coordinates. This is a read-only attribute. (This will be 0 for a CRAM file indexed by a .crai index, as that index format does not record these statistics.) """ def __get__(self): self.check_index() cdef int tid cdef uint64_t total = hts_idx_get_n_no_coor(self.index) cdef uint64_t mapped, unmapped for tid from 0 <= tid < self.header.nreferences: with nogil: hts_idx_get_stat(self.index, tid, &mapped, &unmapped) total += unmapped return total property nocoordinate: """int with total number of reads without coordinates according to the statistics recorded in the index, i.e., the statistic printed for "*" by the ``samtools idxstats`` command. This is a read-only attribute. (This will be 0 for a CRAM file indexed by a .crai index, as that index format does not record these statistics.) """ def __get__(self): self.check_index() cdef uint64_t n with nogil: n = hts_idx_get_n_no_coor(self.index) return n def get_index_statistics(self): """return statistics about mapped/unmapped reads per chromosome as they are stored in the index, similarly to the statistics printed by the ``samtools idxstats`` command. CRAI indexes do not record these statistics, so for a CRAM file with a .crai index the returned statistics will all be 0. Returns: list : a list of records for each chromosome. Each record has the attributes 'contig', 'mapped', 'unmapped' and 'total'. """ self.check_index() cdef int tid cdef uint64_t mapped, unmapped results = [] # TODO: use header for tid from 0 <= tid < self.nreferences: with nogil: hts_idx_get_stat(self.index, tid, &mapped, &unmapped) results.append( IndexStats._make(( self.get_reference_name(tid), mapped, unmapped, mapped + unmapped))) return results ############################################################### ## file-object like iterator access ## note: concurrent access will cause errors (see IteratorRow ## and multiple_iterators) ## Possible solutions: deprecate or open new file handle def __iter__(self): if not self.is_open: raise ValueError("I/O operation on closed file") if not self.is_bam and self.header.nreferences == 0: raise NotImplementedError( "can not iterate over samfile without header") return self cdef bam1_t * getCurrent(self): return self.b cdef int cnext(self): ''' cversion of iterator. Used by :class:`pysam.AlignmentFile.IteratorColumn`. ''' cdef int ret cdef bam_hdr_t * hdr = self.header.ptr with nogil: ret = sam_read1(self.htsfile, hdr, self.b) return ret def __next__(self): cdef int ret = self.cnext() if ret >= 0: return makeAlignedSegment(self.b, self.header) elif ret == -1: raise StopIteration else: raise IOError(read_failure_reason(ret)) ########################################### # methods/properties referencing the header def is_valid_tid(self, int tid): """ return True if the numerical :term:`tid` is valid; False otherwise. Note that the unmapped tid code (-1) counts as an invalid. """ if self.header is None: raise ValueError("header not available in closed files") return self.header.is_valid_tid(tid) def get_tid(self, reference): """ return the numerical :term:`tid` corresponding to :term:`reference` returns -1 if reference is not known. """ if self.header is None: raise ValueError("header not available in closed files") return self.header.get_tid(reference) def get_reference_name(self, tid): """ return :term:`reference` name corresponding to numerical :term:`tid` """ if self.header is None: raise ValueError("header not available in closed files") return self.header.get_reference_name(tid) def get_reference_length(self, reference): """ return :term:`reference` length corresponding to numerical :term:`tid` """ if self.header is None: raise ValueError("header not available in closed files") return self.header.get_reference_length(reference) property nreferences: """int with the number of :term:`reference` sequences in the file. This is a read-only attribute.""" def __get__(self): if self.header: return self.header.nreferences else: raise ValueError("header not available in closed files") property references: """tuple with the names of :term:`reference` sequences. This is a read-only attribute""" def __get__(self): if self.header: return self.header.references else: raise ValueError("header not available in closed files") property lengths: """tuple of the lengths of the :term:`reference` sequences. This is a read-only attribute. The lengths are in the same order as :attr:`pysam.AlignmentFile.references` """ def __get__(self): if self.header: return self.header.lengths else: raise ValueError("header not available in closed files") # Compatibility functions for pysam < 0.14 property text: """deprecated, use .header directly""" def __get__(self): if self.header: return self.header.__str__() else: raise ValueError("header not available in closed files") # Compatibility functions for pysam < 0.8.3 def gettid(self, reference): """deprecated, use get_tid() instead""" return self.get_tid(reference) def getrname(self, tid): """deprecated, use get_reference_name() instead""" return self.get_reference_name(tid) cdef class IteratorRow: '''abstract base class for iterators over mapped reads. Various iterators implement different behaviours for wrapping around contig boundaries. Examples include: :class:`pysam.IteratorRowRegion` iterate within a single contig and a defined region. :class:`pysam.IteratorRowAll` iterate until EOF. This iterator will also include unmapped reads. :class:`pysam.IteratorRowAllRefs` iterate over all reads in all reference sequences. The method :meth:`AlignmentFile.fetch` returns an IteratorRow. .. note:: It is usually not necessary to create an object of this class explicitly. It is returned as a result of call to a :meth:`AlignmentFile.fetch`. ''' def __init__(self, AlignmentFile samfile, int multiple_iterators=False): cdef char *cfilename cdef char *creference_filename cdef char *cindexname = NULL if not samfile.is_open: raise ValueError("I/O operation on closed file") # makes sure that samfile stays alive as long as the # iterator is alive self.samfile = samfile # reopen the file - note that this makes the iterator # slow and causes pileup to slow down significantly. if multiple_iterators: cfilename = samfile.filename with nogil: self.htsfile = hts_open(cfilename, 'r') assert self.htsfile != NULL if samfile.has_index(): if samfile.index_filename: cindexname = bindex_filename = encode_filename(samfile.index_filename) with nogil: self.index = sam_index_load2(self.htsfile, cfilename, cindexname) else: self.index = NULL # need to advance in newly opened file to position after header # better: use seek/tell? with nogil: hdr = sam_hdr_read(self.htsfile) if hdr is NULL: raise IOError("unable to read header information") self.header = makeAlignmentHeader(hdr) self.owns_samfile = True # options specific to CRAM files if samfile.is_cram and samfile.reference_filename: creference_filename = samfile.reference_filename hts_set_opt(self.htsfile, CRAM_OPT_REFERENCE, creference_filename) else: self.htsfile = samfile.htsfile self.index = samfile.index self.owns_samfile = False self.header = samfile.header self.retval = 0 self.b = bam_init1() def __dealloc__(self): bam_destroy1(self.b) if self.owns_samfile: hts_close(self.htsfile) hts_idx_destroy(self.index) cdef class IteratorRowRegion(IteratorRow): """*(AlignmentFile samfile, int tid, int beg, int stop, int multiple_iterators=False)* iterate over mapped reads in a region. .. note:: It is usually not necessary to create an object of this class explicitly. It is returned as a result of call to a :meth:`AlignmentFile.fetch`. """ def __init__(self, AlignmentFile samfile, int tid, int beg, int stop, int multiple_iterators=False): if not samfile.has_index(): raise ValueError("no index available for iteration") IteratorRow.__init__(self, samfile, multiple_iterators=multiple_iterators) with nogil: self.iter = sam_itr_queryi( self.index, tid, beg, stop) def __iter__(self): return self cdef bam1_t * getCurrent(self): return self.b cdef int cnext(self): '''cversion of iterator. Used by IteratorColumn''' with nogil: self.retval = hts_itr_next(hts_get_bgzfp(self.htsfile), self.iter, self.b, self.htsfile) def __next__(self): self.cnext() if self.retval >= 0: return makeAlignedSegment(self.b, self.header) elif self.retval == -1: raise StopIteration elif self.retval == -2: # Note: it is currently not the case that hts_iter_next # returns -2 for a truncated file. # See https://github.com/pysam-developers/pysam/pull/50#issuecomment-64928625 raise IOError('truncated file') else: raise IOError("error while reading file {}: {}".format(self.samfile.filename, self.retval)) def __dealloc__(self): hts_itr_destroy(self.iter) cdef class IteratorRowHead(IteratorRow): """*(AlignmentFile samfile, n, int multiple_iterators=False)* iterate over first n reads in `samfile` .. note:: It is usually not necessary to create an object of this class explicitly. It is returned as a result of call to a :meth:`AlignmentFile.head`. """ def __init__(self, AlignmentFile samfile, int n, int multiple_iterators=False): IteratorRow.__init__(self, samfile, multiple_iterators=multiple_iterators) self.max_rows = n self.current_row = 0 def __iter__(self): return self cdef bam1_t * getCurrent(self): return self.b cdef int cnext(self): '''cversion of iterator. Used by IteratorColumn''' cdef int ret cdef bam_hdr_t * hdr = self.header.ptr with nogil: ret = sam_read1(self.htsfile, hdr, self.b) return ret def __next__(self): if self.current_row >= self.max_rows: raise StopIteration cdef int ret = self.cnext() if ret >= 0: self.current_row += 1 return makeAlignedSegment(self.b, self.header) elif ret == -1: raise StopIteration else: raise IOError(read_failure_reason(ret)) cdef class IteratorRowAll(IteratorRow): """*(AlignmentFile samfile, int multiple_iterators=False)* iterate over all reads in `samfile` .. note:: It is usually not necessary to create an object of this class explicitly. It is returned as a result of call to a :meth:`AlignmentFile.fetch`. """ def __init__(self, AlignmentFile samfile, int multiple_iterators=False): IteratorRow.__init__(self, samfile, multiple_iterators=multiple_iterators) def __iter__(self): return self cdef bam1_t * getCurrent(self): return self.b cdef int cnext(self): '''cversion of iterator. Used by IteratorColumn''' cdef int ret cdef bam_hdr_t * hdr = self.header.ptr with nogil: ret = sam_read1(self.htsfile, hdr, self.b) return ret def __next__(self): cdef int ret = self.cnext() if ret >= 0: return makeAlignedSegment(self.b, self.header) elif ret == -1: raise StopIteration else: raise IOError(read_failure_reason(ret)) cdef class IteratorRowAllRefs(IteratorRow): """iterates over all mapped reads by chaining iterators over each reference .. note:: It is usually not necessary to create an object of this class explicitly. It is returned as a result of call to a :meth:`AlignmentFile.fetch`. """ def __init__(self, AlignmentFile samfile, multiple_iterators=False): IteratorRow.__init__(self, samfile, multiple_iterators=multiple_iterators) if not samfile.has_index(): raise ValueError("no index available for fetch") self.tid = -1 def nextiter(self): # get a new iterator for a chromosome. The file # will not be re-opened. self.rowiter = IteratorRowRegion(self.samfile, self.tid, 0, MAX_POS) # set htsfile and header of the rowiter # to the values in this iterator to reflect multiple_iterators self.rowiter.htsfile = self.htsfile self.rowiter.header = self.header # make sure the iterator understand that IteratorRowAllRefs # has ownership self.rowiter.owns_samfile = False def __iter__(self): return self def __next__(self): # Create an initial iterator if self.tid == -1: if not self.samfile.nreferences: raise StopIteration self.tid = 0 self.nextiter() while 1: self.rowiter.cnext() # If current iterator is not exhausted, return aligned read if self.rowiter.retval > 0: return makeAlignedSegment(self.rowiter.b, self.header) self.tid += 1 # Otherwise, proceed to next reference or stop if self.tid < self.samfile.nreferences: self.nextiter() else: raise StopIteration cdef class IteratorRowSelection(IteratorRow): """*(AlignmentFile samfile)* iterate over reads in `samfile` at a given list of file positions. .. note:: It is usually not necessary to create an object of this class explicitly. It is returned as a result of call to a :meth:`AlignmentFile.fetch`. """ def __init__(self, AlignmentFile samfile, positions, int multiple_iterators=True): IteratorRow.__init__(self, samfile, multiple_iterators=multiple_iterators) self.positions = positions self.current_pos = 0 def __iter__(self): return self cdef bam1_t * getCurrent(self): return self.b cdef int cnext(self): '''cversion of iterator''' # end iteration if out of positions if self.current_pos >= len(self.positions): return -1 cdef uint64_t pos = self.positions[self.current_pos] with nogil: bgzf_seek(hts_get_bgzfp(self.htsfile), pos, 0) self.current_pos += 1 cdef int ret cdef bam_hdr_t * hdr = self.header.ptr with nogil: ret = sam_read1(self.htsfile, hdr, self.b) return ret def __next__(self): cdef int ret = self.cnext() if ret >= 0: return makeAlignedSegment(self.b, self.header) elif ret == -1: raise StopIteration else: raise IOError(read_failure_reason(ret)) cdef int __advance_nofilter(void *data, bam1_t *b): '''advance without any read filtering. ''' cdef __iterdata * d = <__iterdata*>data cdef int ret with nogil: ret = sam_itr_next(d.htsfile, d.iter, b) return ret cdef int __advance_raw_nofilter(void *data, bam1_t *b): '''advance (without iterator) without any read filtering. ''' cdef __iterdata * d = <__iterdata*>data cdef int ret with nogil: ret = sam_read1(d.htsfile, d.header, b) return ret cdef int __advance_all(void *data, bam1_t *b): '''only use reads for pileup passing basic filters such as BAM_FUNMAP, BAM_FSECONDARY, BAM_FQCFAIL, BAM_FDUP ''' cdef __iterdata * d = <__iterdata*>data cdef mask = BAM_FUNMAP | BAM_FSECONDARY | BAM_FQCFAIL | BAM_FDUP cdef int ret while 1: with nogil: ret = sam_itr_next(d.htsfile, d.iter, b) if ret < 0: break if b.core.flag & d.flag_filter: continue break return ret cdef int __advance_raw_all(void *data, bam1_t *b): '''only use reads for pileup passing basic filters such as BAM_FUNMAP, BAM_FSECONDARY, BAM_FQCFAIL, BAM_FDUP ''' cdef __iterdata * d = <__iterdata*>data cdef int ret while 1: with nogil: ret = sam_read1(d.htsfile, d.header, b) if ret < 0: break if b.core.flag & d.flag_filter: continue break return ret cdef int __advance_samtools(void * data, bam1_t * b): '''advance using same filter and read processing as in the samtools pileup. ''' cdef __iterdata * d = <__iterdata*>data cdef int ret cdef int q while 1: with nogil: ret = sam_itr_next(d.htsfile, d.iter, b) if d.iter else sam_read1(d.htsfile, d.header, b) if ret < 0: break if b.core.flag & d.flag_filter: continue if d.flag_require and not (b.core.flag & d.flag_require): continue # reload sequence if d.fastafile != NULL and b.core.tid != d.tid: if d.seq != NULL: free(d.seq) d.tid = b.core.tid with nogil: d.seq = faidx_fetch_seq( d.fastafile, d.header.target_name[d.tid], 0, MAX_POS, &d.seq_len) if d.seq == NULL: raise ValueError( "reference sequence for '{}' (tid={}) not found".format( d.header.target_name[d.tid], d.tid)) # realign read - changes base qualities if d.seq != NULL and d.compute_baq: # 4th option to realign is flag: # apply_baq = flag&1, extend_baq = flag&2, redo_baq = flag&4 if d.redo_baq: sam_prob_realn(b, d.seq, d.seq_len, 7) else: sam_prob_realn(b, d.seq, d.seq_len, 3) if d.seq != NULL and d.adjust_capq_threshold > 10: q = sam_cap_mapq(b, d.seq, d.seq_len, d.adjust_capq_threshold) if q < 0: continue elif b.core.qual > q: b.core.qual = q if b.core.qual < d.min_mapping_quality: continue if d.ignore_orphans and b.core.flag & BAM_FPAIRED and not (b.core.flag & BAM_FPROPER_PAIR): continue break return ret cdef class IteratorColumn: '''abstract base class for iterators over columns. IteratorColumn objects wrap the pileup functionality of samtools. For reasons of efficiency, the iterator points to the current pileup buffer. The pileup buffer is updated at every iteration. This might cause some unexpected behaviour. For example, consider the conversion to a list:: f = AlignmentFile("file.bam", "rb") result = list(f.pileup()) Here, ``result`` will contain ``n`` objects of type :class:`~pysam.PileupColumn` for ``n`` columns, but each object in ``result`` will contain the same information. The desired behaviour can be achieved by list comprehension:: result = [x.pileups() for x in f.pileup()] ``result`` will be a list of ``n`` lists of objects of type :class:`~pysam.PileupRead`. If the iterator is associated with a :class:`~pysam.Fastafile` using the :meth:`add_reference` method, then the iterator will export the current sequence via the methods :meth:`get_sequence` and :meth:`seq_len`. See :class:`~AlignmentFile.pileup` for kwargs to the iterator. ''' def __cinit__( self, AlignmentFile samfile, **kwargs): self.samfile = samfile self.fastafile = kwargs.get("fastafile", None) self.stepper = kwargs.get("stepper", "samtools") self.max_depth = kwargs.get("max_depth", 8000) self.ignore_overlaps = kwargs.get("ignore_overlaps", True) self.min_base_quality = kwargs.get("min_base_quality", 13) self.iterdata.seq = NULL self.iterdata.min_mapping_quality = kwargs.get("min_mapping_quality", 0) self.iterdata.flag_require = kwargs.get("flag_require", 0) self.iterdata.flag_filter = kwargs.get("flag_filter", BAM_FUNMAP | BAM_FSECONDARY | BAM_FQCFAIL | BAM_FDUP) self.iterdata.adjust_capq_threshold = kwargs.get("adjust_capq_threshold", 0) self.iterdata.compute_baq = kwargs.get("compute_baq", True) self.iterdata.redo_baq = kwargs.get("redo_baq", False) self.iterdata.ignore_orphans = kwargs.get("ignore_orphans", True) self.tid = 0 self.pos = 0 self.n_plp = 0 self.plp = NULL self.pileup_iter = NULL def __iter__(self): return self cdef int cnext(self): '''perform next iteration. ''' # do not release gil here because of call-backs cdef int ret = bam_mplp_auto(self.pileup_iter, &self.tid, &self.pos, &self.n_plp, &self.plp) return ret cdef char * get_sequence(self): '''return current reference sequence underlying the iterator. ''' return self.iterdata.seq property seq_len: '''current sequence length.''' def __get__(self): return self.iterdata.seq_len def add_reference(self, FastaFile fastafile): ''' add reference sequences in `fastafile` to iterator.''' self.fastafile = fastafile if self.iterdata.seq != NULL: free(self.iterdata.seq) self.iterdata.tid = -1 self.iterdata.fastafile = self.fastafile.fastafile def has_reference(self): ''' return true if iterator is associated with a reference''' return self.fastafile cdef _setup_iterator(self, int tid, int start, int stop, int multiple_iterators=0): '''setup the iterator structure''' self.iter = IteratorRowRegion(self.samfile, tid, start, stop, multiple_iterators) self.iterdata.htsfile = self.samfile.htsfile self.iterdata.iter = self.iter.iter self.iterdata.seq = NULL self.iterdata.tid = -1 self.iterdata.header = self.samfile.header.ptr if self.fastafile is not None: self.iterdata.fastafile = self.fastafile.fastafile else: self.iterdata.fastafile = NULL # Free any previously allocated memory before reassigning # pileup_iter self._free_pileup_iter() cdef void * data[1] data[0] = &self.iterdata if self.stepper is None or self.stepper == "all": with nogil: self.pileup_iter = bam_mplp_init(1, &__advance_all, data) elif self.stepper == "nofilter": with nogil: self.pileup_iter = bam_mplp_init(1, &__advance_nofilter, data) elif self.stepper == "samtools": with nogil: self.pileup_iter = bam_mplp_init(1, &__advance_samtools, data) else: raise ValueError( "unknown stepper option `%s` in IteratorColumn" % self.stepper) if self.max_depth: with nogil: bam_mplp_set_maxcnt(self.pileup_iter, self.max_depth) if self.ignore_overlaps: with nogil: bam_mplp_init_overlaps(self.pileup_iter) cdef _setup_raw_rest_iterator(self): '''set up an "iterator" that just uses sam_read1(), similar to HTS_IDX_REST''' self.iter = None self.iterdata.iter = NULL self.iterdata.htsfile = self.samfile.htsfile self.iterdata.seq = NULL self.iterdata.tid = -1 self.iterdata.header = self.samfile.header.ptr if self.fastafile is not None: self.iterdata.fastafile = self.fastafile.fastafile else: self.iterdata.fastafile = NULL # Free any previously allocated memory before reassigning # pileup_iter self._free_pileup_iter() cdef void * data[1] data[0] = &self.iterdata if self.stepper is None or self.stepper == "all": with nogil: self.pileup_iter = bam_mplp_init(1, &__advance_raw_all, data) elif self.stepper == "nofilter": with nogil: self.pileup_iter = bam_mplp_init(1, &__advance_raw_nofilter, data) elif self.stepper == "samtools": with nogil: self.pileup_iter = bam_mplp_init(1, &__advance_samtools, data) else: raise ValueError( "unknown stepper option `%s` in IteratorColumn" % self.stepper) if self.max_depth: with nogil: bam_mplp_set_maxcnt(self.pileup_iter, self.max_depth) if self.ignore_overlaps: with nogil: bam_mplp_init_overlaps(self.pileup_iter) cdef reset(self, tid, start, stop): '''reset iterator position. This permits using the iterator multiple times without having to incur the full set-up costs. ''' if self.iter is None: raise TypeError("Raw iterator set up without region cannot be reset") self.iter = IteratorRowRegion(self.samfile, tid, start, stop, multiple_iterators=0) self.iterdata.iter = self.iter.iter # invalidate sequence if different tid if self.tid != tid: if self.iterdata.seq != NULL: free(self.iterdata.seq) self.iterdata.seq = NULL self.iterdata.tid = -1 # self.pileup_iter = bam_mplp_init(1 # &__advancepileup, # &self.iterdata) with nogil: bam_mplp_reset(self.pileup_iter) cdef _free_pileup_iter(self): '''free the memory alloc'd by bam_plp_init. This is needed before setup_iterator allocates another pileup_iter, or else memory will be lost. ''' if self.pileup_iter != NULL: with nogil: bam_mplp_reset(self.pileup_iter) bam_mplp_destroy(self.pileup_iter) self.pileup_iter = NULL def __dealloc__(self): # reset in order to avoid memory leak messages for iterators # that have not been fully consumed self._free_pileup_iter() self.plp = NULL if self.iterdata.seq != NULL: free(self.iterdata.seq) self.iterdata.seq = NULL # backwards compatibility def hasReference(self): return self.has_reference() cdef char * getSequence(self): return self.get_sequence() def addReference(self, FastaFile fastafile): return self.add_reference(fastafile) cdef class IteratorColumnRegion(IteratorColumn): '''iterates over a region only. ''' def __cinit__(self, AlignmentFile samfile, int tid = 0, int start = 0, int stop = MAX_POS, int truncate = False, int multiple_iterators = True, **kwargs ): # initialize iterator. Multiple iterators not available # for CRAM. if multiple_iterators and samfile.is_cram: warnings.warn("multiple_iterators not implemented for CRAM") multiple_iterators = False self._setup_iterator(tid, start, stop, multiple_iterators) self.start = start self.stop = stop self.truncate = truncate def __next__(self): cdef int n while 1: n = self.cnext() if n < 0: raise ValueError("error during iteration" ) if n == 0: raise StopIteration if self.truncate: if self.start > self.pos: continue if self.pos >= self.stop: raise StopIteration return makePileupColumn(&self.plp, self.tid, self.pos, self.n_plp, self.min_base_quality, self.iterdata.seq, self.samfile.header) cdef class IteratorColumnAllRefs(IteratorColumn): """iterates over all columns by chaining iterators over each reference """ def __cinit__(self, AlignmentFile samfile, **kwargs): # no iteration over empty files if not samfile.nreferences: raise StopIteration # initialize iterator self._setup_iterator(self.tid, 0, MAX_POS, 1) def __next__(self): cdef int n while 1: n = self.cnext() if n < 0: raise ValueError("error during iteration") # proceed to next reference or stop if n == 0: self.tid += 1 if self.tid < self.samfile.nreferences: self._setup_iterator(self.tid, 0, MAX_POS, 0) else: raise StopIteration continue # return result, if within same reference return makePileupColumn(&self.plp, self.tid, self.pos, self.n_plp, self.min_base_quality, self.iterdata.seq, self.samfile.header) cdef class IteratorColumnAll(IteratorColumn): """iterates over all columns, without using an index """ def __cinit__(self, AlignmentFile samfile, **kwargs): self._setup_raw_rest_iterator() def __next__(self): cdef int n n = self.cnext() if n < 0: raise ValueError("error during iteration") if n == 0: raise StopIteration return makePileupColumn(&self.plp, self.tid, self.pos, self.n_plp, self.min_base_quality, self.iterdata.seq, self.samfile.header) cdef class SNPCall: '''the results of a SNP call.''' cdef int _tid cdef int _pos cdef char _reference_base cdef char _genotype cdef int _consensus_quality cdef int _snp_quality cdef int _rms_mapping_quality cdef int _coverage property tid: '''the chromosome ID as is defined in the header''' def __get__(self): return self._tid property pos: '''nucleotide position of SNP.''' def __get__(self): return self._pos property reference_base: '''reference base at pos. ``N`` if no reference sequence supplied.''' def __get__(self): return from_string_and_size( &self._reference_base, 1 ) property genotype: '''the genotype called.''' def __get__(self): return from_string_and_size( &self._genotype, 1 ) property consensus_quality: '''the genotype quality (Phred-scaled).''' def __get__(self): return self._consensus_quality property snp_quality: '''the snp quality (Phred scaled) - probability of consensus being identical to reference sequence.''' def __get__(self): return self._snp_quality property mapping_quality: '''the root mean square (rms) of the mapping quality of all reads involved in the call.''' def __get__(self): return self._rms_mapping_quality property coverage: '''coverage or read depth - the number of reads involved in the call.''' def __get__(self): return self._coverage def __str__(self): return "\t".join( map(str, ( self.tid, self.pos, self.reference_base, self.genotype, self.consensus_quality, self.snp_quality, self.mapping_quality, self.coverage ) ) ) cdef class IndexedReads: """Index a Sam/BAM-file by query name while keeping the original sort order intact. The index is kept in memory and can be substantial. By default, the file is re-openend to avoid conflicts if multiple operators work on the same file. Set `multiple_iterators` = False to not re-open `samfile`. Parameters ---------- samfile : AlignmentFile File to be indexed. multiple_iterators : bool Flag indicating whether the file should be reopened. Reopening prevents existing iterators being affected by the indexing. """ def __init__(self, AlignmentFile samfile, int multiple_iterators=True): cdef char *cfilename # makes sure that samfile stays alive as long as this # object is alive. self.samfile = samfile cdef bam_hdr_t * hdr = NULL assert samfile.is_bam, "can only apply IndexReads on bam files" # multiple_iterators the file - note that this makes the iterator # slow and causes pileup to slow down significantly. if multiple_iterators: cfilename = samfile.filename with nogil: self.htsfile = hts_open(cfilename, 'r') if self.htsfile == NULL: raise OSError("unable to reopen htsfile") # need to advance in newly opened file to position after header # better: use seek/tell? with nogil: hdr = sam_hdr_read(self.htsfile) if hdr == NULL: raise OSError("unable to read header information") self.header = makeAlignmentHeader(hdr) self.owns_samfile = True else: self.htsfile = self.samfile.htsfile self.header = samfile.header self.owns_samfile = False def build(self): '''build the index.''' self.index = collections.defaultdict(list) # this method will start indexing from the current file position cdef int ret = 1 cdef bam1_t * b = calloc(1, sizeof( bam1_t)) if b == NULL: raise MemoryError("could not allocate {} bytes".format(sizeof(bam1_t))) cdef uint64_t pos cdef bam_hdr_t * hdr = self.header.ptr while ret > 0: with nogil: pos = bgzf_tell(hts_get_bgzfp(self.htsfile)) ret = sam_read1(self.htsfile, hdr, b) if ret > 0: qname = charptr_to_str(pysam_bam_get_qname(b)) self.index[qname].append(pos) bam_destroy1(b) def find(self, query_name): '''find `query_name` in index. Returns ------- IteratorRowSelection Returns an iterator over all reads with query_name. Raises ------ KeyError if the `query_name` is not in the index. ''' if query_name in self.index: return IteratorRowSelection( self.samfile, self.index[query_name], multiple_iterators = False) else: raise KeyError("read %s not found" % query_name) def __dealloc__(self): if self.owns_samfile: hts_close(self.htsfile)