pysam/csamtools.pyx

   1 # cython: embedsignature=True
   2 # cython: profile=True
   3 # adds doc-strings for sphinx
   4 import tempfile
   5 import os
   6 import sys
   7 import types
   8 import itertools
   9 import struct
  10 import ctypes
  11 import collections
  12 import re
  13 import platform
  14 import warnings
  15 from cpython cimport PyErr_SetString, PyBytes_Check, PyUnicode_Check, PyBytes_FromStringAndSize
  16 from cpython.version cimport PY_MAJOR_VERSION
  17
  18 ########################################################################
  19 ########################################################################
  20 ########################################################################
  21 ## Python 3 compatibility functions
  22 ########################################################################
  23 IS_PYTHON3 = PY_MAJOR_VERSION >= 3
  24 cdef from_string_and_size(char* s, size_t length):
  25     if PY_MAJOR_VERSION < 3:
  26         return s[:length]
  27     else:
  28         return s[:length].decode("ascii")
  29
  30 # filename encoding (copied from lxml.etree.pyx)
  31 cdef str _FILENAME_ENCODING
  32 _FILENAME_ENCODING = sys.getfilesystemencoding()
  33 if _FILENAME_ENCODING is None:
  34     _FILENAME_ENCODING = sys.getdefaultencoding()
  35 if _FILENAME_ENCODING is None:
  36     _FILENAME_ENCODING = 'ascii'
  37
  38 #cdef char* _C_FILENAME_ENCODING
  39 #_C_FILENAME_ENCODING = <char*>_FILENAME_ENCODING
  40
  41 cdef bytes _my_encodeFilename(object filename):
  42     u"""Make sure a filename is 8-bit encoded (or None).
  43     """
  44     if filename is None:
  45         return None
  46     elif PyBytes_Check(filename):
  47         return filename
  48     elif PyUnicode_Check(filename):
  49         return filename.encode(_FILENAME_ENCODING)
  50     else:
  51         raise TypeError, u"Argument must be string or unicode."
  52
  53
  54 cdef bytes _force_bytes(object s):
  55     u"""convert string or unicode object to bytes, assuming ascii encoding.
  56     """
  57     if PY_MAJOR_VERSION < 3:
  58         return s
  59     elif s is None:
  60         return None
  61     elif PyBytes_Check(s):
  62         return s
  63     elif PyUnicode_Check(s):
  64         return s.encode('ascii')
  65     else:
  66         raise TypeError, u"Argument must be string, bytes or unicode."
  67
  68 cdef inline bytes _force_cmdline_bytes(object s):
  69     return _force_bytes(s)
  70
  71 cdef _charptr_to_str(char* s):
  72     if PY_MAJOR_VERSION < 3:
  73         return s
  74     else:
  75         return s.decode("ascii")
  76
  77 cdef _force_str(object s):
  78     """Return s converted to str type of current Python (bytes in Py2, unicode in Py3)"""
  79     if s is None:
  80         return None
  81     if PY_MAJOR_VERSION < 3:
  82         return s
  83     elif PyBytes_Check(s):
  84         return s.decode('ascii')
  85     else:
  86         # assume unicode
  87         return s
  88 ########################################################################
  89 ########################################################################
  90 ########################################################################
  91 ## Constants and global variables
  92 ########################################################################
  93 # defines imported from samtools
  94 DEF SEEK_SET = 0
  95 DEF SEEK_CUR = 1
  96 DEF SEEK_END = 2
  97
  98 ## These are bits set in the flag.
  99 ## have to put these definitions here, in csamtools.pxd they got ignored
 100 ## @abstract the read is paired in sequencing, no matter whether it is mapped in a pair */
 101 DEF BAM_FPAIRED       =1
 102 ## @abstract the read is mapped in a proper pair */
 103 DEF BAM_FPROPER_PAIR  =2
 104 ## @abstract the read itself is unmapped; conflictive with BAM_FPROPER_PAIR */
 105 DEF BAM_FUNMAP        =4
 106 ## @abstract the mate is unmapped */
 107 DEF BAM_FMUNMAP       =8
 108 ## @abstract the read is mapped to the reverse strand */
 109 DEF BAM_FREVERSE      =16
 110 ## @abstract the mate is mapped to the reverse strand */
 111 DEF BAM_FMREVERSE     =32
 112 ## @abstract this is read1 */
 113 DEF BAM_FREAD1        =64
 114 ## @abstract this is read2 */
 115 DEF BAM_FREAD2       =128
 116 ## @abstract not primary alignment */
 117 DEF BAM_FSECONDARY   =256
 118 ## @abstract QC failure */
 119 DEF BAM_FQCFAIL      =512
 120 ## @abstract optical or PCR duplicate */
 121 DEF BAM_FDUP        =1024
 122
 123 #####################################################################
 124 # CIGAR operations
 125 DEF BAM_CIGAR_SHIFT=4
 126 DEF BAM_CIGAR_MASK=((1 << BAM_CIGAR_SHIFT) - 1)
 127
 128 DEF BAM_CMATCH     = 0
 129 DEF BAM_CINS       = 1
 130 DEF BAM_CDEL       = 2
 131 DEF BAM_CREF_SKIP  = 3
 132 DEF BAM_CSOFT_CLIP = 4
 133 DEF BAM_CHARD_CLIP = 5
 134 DEF BAM_CPAD       = 6
 135 DEF BAM_CEQUAL     = 7
 136 DEF BAM_CDIFF      = 8
 137
 138 cdef char* CODE2CIGAR= "MIDNSHP=X"
 139 if IS_PYTHON3:
 140     CIGAR2CODE = dict( [y,x] for x,y in enumerate( CODE2CIGAR) )
 141 else:
 142     CIGAR2CODE = dict( [ord(y),x] for x,y in enumerate( CODE2CIGAR) )
 143 CIGAR_REGEX = re.compile( "([MIDNSHP=X])(\d+)" )
 144
 145 #####################################################################
 146 ## set pysam stderr to /dev/null
 147 pysam_unset_stderr()
 148
 149 #####################################################################
 150 # hard-coded constants
 151 cdef char * bam_nt16_rev_table = "=ACMGRSVTWYHKDBN"
 152 cdef int max_pos = 2 << 29
 153
 154 #####################################################################
 155 #####################################################################
 156 #####################################################################
 157 ## private factory methods
 158 #####################################################################
 159 cdef class AlignedRead
 160 cdef makeAlignedRead(bam1_t * src):
 161     '''enter src into AlignedRead.'''
 162     cdef AlignedRead dest = AlignedRead.__new__(AlignedRead)
 163     dest._delegate = bam_dup1(src)
 164     return dest
 165
 166 cdef class PileupProxy
 167 cdef makePileupProxy( bam_pileup1_t ** plp, int tid, int pos, int n ):
 168      cdef PileupProxy dest = PileupProxy.__new__(PileupProxy)
 169      dest.plp = plp
 170      dest.tid = tid
 171      dest.pos = pos
 172      dest.n = n
 173      return dest
 174
 175 cdef class PileupRead
 176 cdef makePileupRead( bam_pileup1_t * src ):
 177     '''fill a  PileupRead object from a bam_pileup1_t * object.'''
 178     cdef PileupRead dest = PileupRead.__new__(PileupRead)
 179     dest._alignment = makeAlignedRead( src.b )
 180     dest._qpos = src.qpos
 181     dest._indel = src.indel
 182     dest._level = src.level
 183     dest._is_del = src.is_del
 184     dest._is_head = src.is_head
 185     dest._is_tail = src.is_tail
 186     return dest
 187
 188 cdef convertBinaryTagToList( uint8_t * s ):
 189     """return bytesize, number of values list of values in s."""
 190     cdef char auxtype
 191     cdef uint8_t byte_size
 192     cdef int32_t nvalues
 193
 194     # get byte size
 195     auxtype = s[0]
 196     byte_size = bam_aux_type2size( auxtype )
 197     s += 1
 198     # get number of values in array
 199     nvalues = (<int32_t*>s)[0]
 200     s += 4
 201     # get values
 202     values = []
 203     if auxtype == 'c':
 204         for x from 0 <= x < nvalues:
 205             values.append((<int8_t*>s)[0])
 206             s += 1
 207     elif auxtype == 'C':
 208         for x from 0 <= x < nvalues:
 209             values.append((<uint8_t*>s)[0])
 210             s += 1
 211     elif auxtype == 's':
 212         for x from 0 <= x < nvalues:
 213             values.append((<int16_t*>s)[0])
 214             s += 2
 215     elif auxtype == 'S':
 216         for x from 0 <= x < nvalues:
 217             values.append((<uint16_t*>s)[0])
 218             s += 2
 219     elif auxtype == 'i':
 220         for x from 0 <= x < nvalues:
 221             values.append((<int32_t*>s)[0])
 222             s += 4
 223     elif auxtype == 'I':
 224         for x from 0 <= x < nvalues:
 225             values.append((<uint32_t*>s)[0])
 226             s += 4
 227     elif auxtype == 'f':
 228         for x from 0 <= x < nvalues:
 229             values.append((<float*>s)[0])
 230             s += 4
 231
 232     return byte_size, nvalues, values
 233
 234 #####################################################################
 235 #####################################################################
 236 #####################################################################
 237 ## Generic callbacks for inserting python callbacks.
 238 #####################################################################
 239 cdef int fetch_callback( bam1_t *alignment, void *f):
 240     '''callback for bam_fetch.
 241
 242     calls function in *f* with a new :class:`AlignedRead` object as parameter.
 243     '''
 244     a = makeAlignedRead( alignment )
 245     (<object>f)(a)
 246
 247 class PileupColumn(object):
 248     '''A pileup column. A pileup column contains
 249     all the reads that map to a certain target base.
 250
 251     tid
 252         chromosome ID as is defined in the header
 253     pos
 254         the target base coordinate (0-based)
 255     n
 256         number of reads mapping to this column
 257     pileups
 258         list of reads (:class:`pysam.PileupRead`) aligned to this column
 259     '''
 260     def __str__(self):
 261         return "\t".join( map(str, (self.tid, self.pos, self.n))) +\
 262             "\n" + "\n".join( map(str, self.pileups) )
 263
 264 cdef int pileup_callback( uint32_t tid, uint32_t pos, int n, bam_pileup1_t *pl, void *f):
 265     '''callback for pileup.
 266
 267     calls function in *f* with a new :class:`Pileup` object as parameter.
 268
 269     tid
 270         chromosome ID as is defined in the header
 271     pos
 272         start coordinate of the alignment, 0-based
 273     n
 274         number of elements in pl array
 275     pl
 276         array of alignments
 277     data
 278         user provided data
 279     '''
 280
 281     p = PileupColumn()
 282     p.tid = tid
 283     p.pos = pos
 284     p.n = n
 285     pileups = []
 286
 287     cdef int x
 288     for x from 0 <= x < n:
 289         pileups.append( makePileupRead( &(pl[x]) ) )
 290     p.pileups = pileups
 291
 292     (<object>f)(p)
 293
 294 cdef int pileup_fetch_callback( bam1_t *b, void *data):
 295     '''callback for bam_fetch.
 296
 297     Fetches reads and submits them to pileup.
 298     '''
 299     cdef bam_plbuf_t * buf
 300     buf = <bam_plbuf_t*>data
 301     bam_plbuf_push(b, buf)
 302     return 0
 303
 304 class StderrStore():
 305     '''
 306     stderr is captured.
 307     '''
 308     def __init__(self):
 309         return
 310         self.stderr_h, self.stderr_f = tempfile.mkstemp()
 311         self.stderr_save = Outs( sys.stderr.fileno() )
 312         self.stderr_save.setfd( self.stderr_h )
 313
 314     def readAndRelease( self ):
 315         return []
 316         self.stderr_save.restore()
 317         lines = []
 318         if os.path.exists(self.stderr_f):
 319             lines = open( self.stderr_f, "r" ).readlines()
 320             os.remove( self.stderr_f )
 321         return lines
 322
 323     def release(self):
 324         return
 325         self.stderr_save.restore()
 326         if os.path.exists(self.stderr_f):
 327             os.remove( self.stderr_f )
 328
 329     def __del__(self):
 330         self.release()
 331
 332 class StderrStoreWindows():
 333     '''does nothing. stderr can't be redirected on windows'''
 334     def __init__(self): pass
 335     def readAndRelease(self): return []
 336     def release(self): pass
 337
 338 if platform.system()=='Windows':
 339     del StderrStore
 340     StderrStore = StderrStoreWindows
 341
 342
 343 ######################################################################
 344 ######################################################################
 345 ######################################################################
 346 # valid types for sam headers
 347 VALID_HEADER_TYPES = { "HD" : dict,
 348                        "SQ" : list,
 349                        "RG" : list,
 350                        "PG" : list,
 351                        "CO" : list }
 352
 353 # order of records within sam headers
 354 VALID_HEADERS = ("HD", "SQ", "RG", "PG", "CO" )
 355
 356 # type conversions within sam header records
 357 VALID_HEADER_FIELDS = { "HD" : { "VN" : str, "SO" : str, "GO" : str },
 358                         "SQ" : { "SN" : str, "LN" : int, "AS" : str, "M5" : str, "UR" : str, "SP" : str },
 359                         "RG" : { "ID" : str, "SM" : str, "LB" : str, "DS" : str, "PU" : str, "PI" : str,
 360                                  "CN" : str, "DT" : str, "PL" : str, "FO" : str, "KS" : str },
 361                         "PG" : { "PN" : str, "ID" : str, "VN" : str, "CL" : str }, }
 362
 363 # output order of fields within records
 364 VALID_HEADER_ORDER = { "HD" : ( "VN", "SO", "GO" ),
 365                        "SQ" : ( "SN", "LN", "AS", "M5" , "UR" , "SP" ),
 366                        "RG" : ( "ID", "SM", "LB", "DS" , "PU" , "PI" , "CN" , "DT", "PL", "FO", "KS" ),
 367                        "PG" : ( "PN", "ID", "VN", "CL" ), }
 368
 369
 370 ######################################################################
 371 ######################################################################
 372 ######################################################################
 373 ## Public methods
 374 ######################################################################
 375 cdef class Fastafile:
 376     '''*(filename)*
 377
 378     A *FASTA* file. The file is automatically opened.
 379
 380     The file expects an indexed fasta file.
 381
 382     TODO:
 383         add automatic indexing.
 384         add function to get sequence names.
 385     '''
 386
 387     def __cinit__(self, *args, **kwargs ):
 388         self.fastafile = NULL
 389         self._filename = NULL
 390         self._open( *args, **kwargs )
 391
 392     def _isOpen( self ):
 393         '''return true if samfile has been opened.'''
 394         return self.fastafile != NULL
 395
 396     def __len__(self):
 397         if self.fastafile == NULL:
 398             raise ValueError( "calling len() on closed file" )
 399
 400         return faidx_fetch_nseq(self.fastafile)
 401
 402     def _open( self,
 403                filename ):
 404         '''open an indexed fasta file.
 405
 406         This method expects an indexed fasta file.
 407         '''
 408
 409         # close a previously opened file
 410         if self.fastafile != NULL: self.close()
 411         if self._filename != NULL: free(self._filename)
 412         filename = _my_encodeFilename(filename)
 413         self._filename = strdup(filename)
 414         self.fastafile = fai_load( filename )
 415
 416         if self.fastafile == NULL:
 417             raise IOError("could not open file `%s`" % filename )
 418
 419     def close( self ):
 420         if self.fastafile != NULL:
 421             fai_destroy( self.fastafile )
 422             self.fastafile = NULL
 423
 424     def __dealloc__(self):
 425         self.close()
 426         if self._filename != NULL: free(self._filename)
 427
 428     property filename:
 429         '''number of :term:`filename` associated with this object.'''
 430         def __get__(self):
 431             if not self._isOpen(): raise ValueError( "I/O operation on closed file" )
 432             return self._filename
 433
 434     def fetch( self,
 435                reference = None,
 436                start = None,
 437                end = None,
 438                region = None):
 439
 440         '''*(reference = None, start = None, end = None, region = None)*
 441
 442         fetch :meth:`AlignedRead` objects in a :term:`region` using 0-based indexing.
 443
 444         The region is specified by :term:`reference`, *start* and *end*.
 445
 446         fetch returns an empty string if the region is out of range or addresses an unknown *reference*.
 447
 448         If *reference* is given and *start* is None, the sequence from the
 449         first base is returned. Similarly, if *end* is None, the sequence
 450         until the last base is returned.
 451
 452         Alternatively, a samtools :term:`region` string can be supplied.
 453         '''
 454
 455         if not self._isOpen():
 456             raise ValueError( "I/O operation on closed file" )
 457
 458         cdef int length
 459         cdef char * seq
 460
 461         if not region:
 462             if reference is None: raise ValueError( 'no sequence/region supplied.' )
 463             if start is None: start = 0
 464             if end is None: end = max_pos -1
 465
 466             if start > end: raise ValueError( 'invalid region: start (%i) > end (%i)' % (start, end) )
 467             if start == end: return b""
 468             # valid ranges are from 0 to 2^29-1
 469             if not 0 <= start < max_pos: raise ValueError( 'start out of range (%i)' % start )
 470             if not 0 <= end < max_pos: raise ValueError( 'end out of range (%i)' % end )
 471             # note: faidx_fetch_seq has a bug such that out-of-range access
 472             # always returns the last residue. Hence do not use faidx_fetch_seq,
 473             # but use fai_fetch instead
 474             # seq = faidx_fetch_seq(self.fastafile,
 475             #                       reference,
 476             #                       start,
 477             #                       end-1,
 478             #                       &length)
 479             region = "%s:%i-%i" % (reference, start+1, end)
 480             if PY_MAJOR_VERSION >= 3:
 481                 region = region.encode('ascii')
 482             seq = fai_fetch( self.fastafile,
 483                              region,
 484                              &length )
 485         else:
 486             # samtools adds a '\0' at the end
 487             seq = fai_fetch( self.fastafile, region, &length )
 488
 489         # copy to python
 490         if seq == NULL:
 491             return b""
 492         else:
 493             try:
 494                 py_seq = seq[:length]
 495             finally:
 496                 free(seq)
 497
 498         return py_seq
 499
 500     cdef char * _fetch( self, char * reference, int start, int end, int * length ):
 501         '''fetch sequence for reference, start and end'''
 502
 503         return faidx_fetch_seq(self.fastafile,
 504                                reference,
 505                                start,
 506                                end-1,
 507                                length )
 508
 509 #------------------------------------------------------------------------
 510 #------------------------------------------------------------------------
 511 #------------------------------------------------------------------------
 512 cdef int count_callback( bam1_t *alignment, void *f):
 513      '''callback for bam_fetch - count number of reads.
 514      '''
 515      cdef int* counter = (<int*>f)
 516      counter[0] += 1;
 517
 518 ctypedef struct MateData:
 519      char * name
 520      bam1_t * mate
 521      uint32_t flag
 522
 523 #------------------------------------------------------------------------
 524 #------------------------------------------------------------------------
 525 #------------------------------------------------------------------------
 526 cdef int mate_callback( bam1_t *alignment, void *f):
 527      '''callback for bam_fetch = filter mate
 528      '''
 529      cdef MateData * d = (<MateData*>f)
 530      # printf("mate = %p, name1 = %s, name2=%s\t%i\t%i\t%i\n",
 531      #        d.mate, d.name, bam1_qname(alignment),
 532      #        d.flag, alignment.core.flag, alignment.core.flag & d.flag)
 533
 534      if d.mate == NULL:
 535          # could be sped up by comparing the lengths of query strings first
 536          # using l_qname
 537          #
 538          # also, make sure that we get the other read by comparing
 539          # the flags
 540          if alignment.core.flag & d.flag != 0 and \
 541                  strcmp( bam1_qname( alignment ), d.name ) == 0:
 542              d.mate = bam_dup1( alignment )
 543
 544
 545 cdef class Samfile:
 546     '''*(filename, mode=None, template = None, referencenames = None, referencelengths = None, text = NULL, header = None,
 547          add_sq_text = False, check_header = True, check_sq = True )*
 548
 549     A :term:`SAM`/:term:`BAM` formatted file. The file is automatically opened.
 550
 551     *mode* should be ``r`` for reading or ``w`` for writing. The default is text mode (:term:`SAM`). For binary
 552     (:term:`BAM`) I/O you should append ``b`` for compressed or ``u`` for uncompressed :term:`BAM` output.
 553     Use ``h`` to output header information in text (:term:`TAM`)  mode.
 554
 555     If ``b`` is present, it must immediately follow ``r`` or ``w``.
 556     Valid modes are ``r``, ``w``, ``wh``, ``rb``, ``wb`` and ``wbu``. For instance, to open
 557     a :term:`BAM` formatted file for reading, type::
 558
 559         f = pysam.Samfile('ex1.bam','rb')
 560
 561     If mode is not specified, we will try to auto-detect in the order 'rb', 'r', thus both the following
 562     should work::
 563
 564         f1 = pysam.Samfile('ex1.bam' )
 565         f2 = pysam.Samfile('ex1.sam' )
 566
 567     If an index for a BAM file exists (.bai), it will be opened automatically. Without an index random
 568     access to reads via :meth:`fetch` and :meth:`pileup` is disabled.
 569
 570     For writing, the header of a :term:`SAM` file/:term:`BAM` file can be constituted from several
 571     sources (see also the samtools format specification):
 572
 573         1. If *template* is given, the header is copied from a another *Samfile*
 574            (*template* must be of type *Samfile*).
 575
 576         2. If *header* is given, the header is built from a multi-level dictionary. The first level
 577            are the four types ('HD', 'SQ', ...). The second level are a list of lines, with each line
 578            being a list of tag-value pairs. The header is constructed first from all the defined fields,
 579            followed by user tags in alphabetical order.
 580
 581         3. If *text* is given, new header text is copied from raw text.
 582
 583         4. The names (*referencenames*) and lengths (*referencelengths*) are supplied directly as lists.
 584            By default, 'SQ' and 'LN' tags will be added to the header text. This option can be
 585            changed by unsetting the flag *add_sq_text*.
 586
 587     By default, if file a file is opened in mode 'r', it is checked for a valid header
 588     (*check_header* = True) and a definition of chromosome names (*check_sq* = True).
 589
 590     '''
 591
 592     def __cinit__(self, *args, **kwargs ):
 593         self.samfile = NULL
 594         self._filename = NULL
 595         self.isbam = False
 596         self.isstream = False
 597         self._open( *args, **kwargs )
 598
 599         # allocate memory for iterator
 600         self.b = <bam1_t*>calloc(1, sizeof(bam1_t))
 601
 602     def _isOpen( self ):
 603         '''return true if samfile has been opened.'''
 604         return self.samfile != NULL
 605
 606     def _hasIndex( self ):
 607         '''return true if samfile has an existing (and opened) index.'''
 608         return self.index != NULL
 609
 610     def _open( self,
 611                filename,
 612                mode = None,
 613                Samfile template = None,
 614                referencenames = None,
 615                referencelengths = None,
 616                text = None,
 617                header = None,
 618                port = None,
 619                add_sq_text = True,
 620                check_header = True,
 621                check_sq = True,
 622               ):
 623         '''open a sam/bam file.
 624
 625         If _open is called on an existing bamfile, the current file will be
 626         closed and a new file will be opened.
 627         '''
 628
 629         # read mode autodetection
 630         if mode is None:
 631             try:
 632                 self._open(filename, 'rb', template=template,
 633                            referencenames=referencenames,
 634                            referencelengths=referencelengths,
 635                            text=text, header=header, port=port,
 636                            check_header=check_header,
 637                            check_sq=check_sq)
 638                 return
 639             except ValueError, msg:
 640                 pass
 641
 642             self._open(filename, 'r', template=template,
 643                        referencenames=referencenames,
 644                        referencelengths=referencelengths,
 645                        text=text, header=header, port=port,
 646                        check_header=check_header,
 647                        check_sq=check_sq)
 648             return
 649
 650         assert mode in ( "r","w","rb","wb", "wh", "wbu", "rU" ), "invalid file opening mode `%s`" % mode
 651
 652         #assert filename != NULL
 653
 654         # close a previously opened file
 655         if self.samfile != NULL: self.close()
 656         self.samfile = NULL
 657
 658         cdef bam_header_t * header_to_write
 659         header_to_write = NULL
 660
 661         if self._filename != NULL: free(self._filename )
 662         filename = _my_encodeFilename(filename)
 663         cdef bytes bmode = mode.encode('ascii')
 664         #cdef char* cfilename
 665         #cfilename = filename.encode(_FILENAME_ENCODING)
 666         self._filename = strdup(filename)
 667         self.isstream = strcmp( filename, "-" ) == 0
 668
 669         self.isbam = len(mode) > 1 and mode[1] == 'b'
 670
 671         self.isremote = strncmp(filename,"http:",5) == 0 or \
 672             strncmp(filename,"ftp:",4) == 0
 673
 674         cdef char * ctext
 675         ctext = NULL
 676
 677         if mode[0] == 'w':
 678             # open file for writing
 679
 680             # header structure (used for writing)
 681             if template:
 682                 # copy header from another file
 683                 header_to_write = template.samfile.header
 684
 685             elif header:
 686                 header_to_write = self._buildHeader( header )
 687
 688             else:
 689                 # build header from a target names and lengths
 690                 assert referencenames and referencelengths, "either supply options `template`, `header` or  both `referencenames` and `referencelengths` for writing"
 691                 assert len(referencenames) == len(referencelengths), "unequal names and lengths of reference sequences"
 692
 693                 # allocate and fill header
 694                 referencenames = [ _force_bytes(ref) for ref in referencenames ]
 695                 header_to_write = bam_header_init()
 696                 header_to_write.n_targets = len(referencenames)
 697                 n = 0
 698                 for x in referencenames: n += len(x) + 1
 699                 header_to_write.target_name = <char**>calloc(n, sizeof(char*))
 700                 header_to_write.target_len = <uint32_t*>calloc(n, sizeof(uint32_t))
 701                 for x from 0 <= x < header_to_write.n_targets:
 702                     header_to_write.target_len[x] = referencelengths[x]
 703                     name = referencenames[x]
 704                     header_to_write.target_name[x] = <char*>calloc(len(name)+1, sizeof(char))
 705                     strncpy( header_to_write.target_name[x], name, len(name) )
 706
 707                 # Optionally, if there is no text, add a SAM compatible header to output
 708                 # file.
 709                 if text is None and add_sq_text:
 710                     text = []
 711                     for x from 0 <= x < header_to_write.n_targets:
 712                         text.append( "@SQ\tSN:%s\tLN:%s\n" % (referencenames[x], referencelengths[x] ) )
 713                     text = ''.join(text)
 714
 715                 if text != None:
 716                     # copy without \0
 717                     text = _force_bytes(text)
 718                     ctext = text
 719                     header_to_write.l_text = strlen(ctext)
 720                     header_to_write.text = <char*>calloc( strlen(ctext), sizeof(char) )
 721                     memcpy( header_to_write.text, ctext, strlen(ctext) )
 722
 723                 header_to_write.hash = NULL
 724                 header_to_write.rg2lib = NULL
 725
 726             # open file. Header gets written to file at the same time for bam files
 727             # and sam files (in the latter case, the mode needs to be wh)
 728             store = StderrStore()
 729             self.samfile = samopen( filename, bmode, header_to_write )
 730             store.release()
 731
 732             # bam_header_destroy takes care of cleaning up of all the members
 733             if not template and header_to_write != NULL:
 734                 bam_header_destroy( header_to_write )
 735
 736         elif mode[0] == "r":
 737             # open file for reading
 738             if strncmp( filename, "-", 1) != 0 and \
 739                     not self.isremote and \
 740                     not os.path.exists( filename ):
 741                 raise IOError( "file `%s` not found" % filename)
 742
 743             # try to detect errors
 744             self.samfile = samopen( filename, bmode, NULL )
 745             if self.samfile == NULL:
 746                 raise ValueError( "could not open file (mode='%s') - is it SAM/BAM format?" % mode)
 747
 748             # bam files require a valid header
 749             if self.isbam:
 750                 if self.samfile.header == NULL:
 751                     raise ValueError( "file does not have valid header (mode='%s') - is it BAM format?" % mode )
 752             else:
 753                 # in sam files it is optional (samfile full of unmapped reads)
 754                 if check_header and self.samfile.header == NULL:
 755                     raise ValueError( "file does not have valid header (mode='%s') - is it SAM format?" % mode )
 756
 757             # disabled for autodetection to work
 758             # needs to be disabled so that reading from sam-files without headers works
 759             if check_sq and self.samfile.header.n_targets == 0:
 760                 raise ValueError( "file header is empty (mode='%s') - is it SAM/BAM format?" % mode)
 761
 762         if self.samfile == NULL:
 763             raise IOError("could not open file `%s`" % filename )
 764
 765         # check for index and open if present
 766         if mode[0] == "r" and self.isbam:
 767
 768             if not self.isremote:
 769                 if not os.path.exists(filename + b".bai"):
 770                     self.index = NULL
 771                 else:
 772                     # returns NULL if there is no index or index could not be opened
 773                     self.index = bam_index_load(filename)
 774                     if self.index == NULL:
 775                         raise IOError("error while opening index `%s` " % filename )
 776             else:
 777                 self.index = bam_index_load(filename)
 778                 if self.index == NULL:
 779                     raise IOError("error while opening index `%s` " % filename )
 780
 781             if not self.isstream:
 782                 self.start_offset = bam_tell( self.samfile.x.bam )
 783
 784     def gettid( self, reference ):
 785         '''
 786         convert :term:`reference` name into numerical :term:`tid`
 787
 788         returns -1 if reference is not known.
 789         '''
 790         if not self._isOpen(): raise ValueError( "I/O operation on closed file" )
 791         reference = _force_bytes(reference)
 792         return pysam_reference2tid( self.samfile.header, reference )
 793
 794     def getrname( self, tid ):
 795         '''
 796         convert numerical :term:`tid` into :term:`reference` name.'''
 797         if not self._isOpen(): raise ValueError( "I/O operation on closed file" )
 798         if not 0 <= tid < self.samfile.header.n_targets:
 799             raise ValueError( "tid %i out of range 0<=tid<%i" % (tid, self.samfile.header.n_targets ) )
 800         return _charptr_to_str(self.samfile.header.target_name[tid])
 801
 802     cdef char * _getrname( self, int tid ): # TODO unused
 803         '''
 804         convert numerical :term:`tid` into :term:`reference` name.'''
 805         if not self._isOpen(): raise ValueError( "I/O operation on closed file" )
 806         if not 0 <= tid < self.samfile.header.n_targets:
 807             raise ValueError( "tid %i out of range 0<=tid<%i" % (tid, self.samfile.header.n_targets ) )
 808         return self.samfile.header.target_name[tid]
 809
 810     def _parseRegion( self,
 811                       reference = None,
 812                       start = None,
 813                       end = None,
 814                       region = None ):
 815         '''
 816         parse region information.
 817
 818         raise ValueError for for invalid regions.
 819
 820         returns a tuple of flag, tid, start and end. Flag indicates
 821         whether some coordinates were supplied.
 822
 823         Note that regions are 1-based, while start,end are python coordinates.
 824         '''
 825         # This method's main objective is to translate from a reference to a tid.
 826         # For now, it calls bam_parse_region, which is clumsy. Might be worth
 827         # implementing it all in pysam (makes use of khash).
 828
 829         cdef int rtid
 830         cdef long long rstart
 831         cdef long long rend
 832
 833         rtid = -1
 834         rstart = 0
 835         rend = max_pos
 836         if start != None:
 837             try:
 838                 rstart = start
 839             except OverflowError:
 840                 raise ValueError( 'start out of range (%i)' % start )
 841
 842         if end != None:
 843             try:
 844                 rend = end
 845             except OverflowError:
 846                 raise ValueError( 'end out of range (%i)' % end )
 847
 848         if region:
 849             region = _force_str(region)
 850             parts = re.split( "[:-]", region )
 851             reference = parts[0]
 852             if len(parts) >= 2: rstart = int(parts[1]) - 1
 853             if len(parts) >= 3: rend = int(parts[2])
 854
 855         if not reference: return 0, 0, 0, 0
 856
 857         rtid = self.gettid( reference )
 858         if rtid < 0: raise ValueError( "invalid reference `%s`" % reference )
 859         if rstart > rend: raise ValueError( 'invalid coordinates: start (%i) > end (%i)' % (rstart, rend) )
 860         if not 0 <= rstart < max_pos: raise ValueError( 'start out of range (%i)' % rstart )
 861         if not 0 <= rend <= max_pos: raise ValueError( 'end out of range (%i)' % rend )
 862
 863         return 1, rtid, rstart, rend
 864
 865     def reset( self ):
 866         '''reset file position to beginning of read section.'''
 867         return self.seek( self.start_offset, 0 )
 868
 869     def seek( self, uint64_t offset, int where = 0):
 870         '''
 871         move file pointer to position *offset*, see :meth:`pysam.Samfile.tell`.
 872         '''
 873
 874         if not self._isOpen():
 875             raise ValueError( "I/O operation on closed file" )
 876         if not self.isbam:
 877             raise NotImplementedError("seek only available in bam files")
 878         if self.isstream:
 879             raise OSError("seek no available in streams")
 880
 881         return bam_seek( self.samfile.x.bam, offset, where )
 882
 883     def tell( self ):
 884         '''
 885         return current file position
 886         '''
 887         if not self._isOpen():
 888             raise ValueError( "I/O operation on closed file" )
 889         if not self.isbam:
 890             raise NotImplementedError("seek only available in bam files")
 891
 892         return bam_tell( self.samfile.x.bam )
 893
 894     def fetch( self,
 895                reference = None,
 896                start = None,
 897                end = None,
 898                region = None,
 899                callback = None,
 900                until_eof = False ):
 901         '''
 902         fetch aligned reads in a :term:`region` using 0-based indexing. The region is specified by
 903         :term:`reference`, *start* and *end*. Alternatively, a samtools :term:`region` string can
 904         be supplied.
 905
 906         Without *reference* or *region* all mapped reads will be fetched. The reads will be returned
 907         ordered by reference sequence, which will not necessarily be the order within the file.
 908
 909         If *until_eof* is given, all reads from the current file position will be returned
 910         in order as they are within the file. Using this option will also fetch unmapped reads.
 911
 912         If only *reference* is set, all reads aligned to *reference* will be fetched.
 913
 914         The method returns an iterator of type :class:`pysam.IteratorRow` unless
 915         a *callback is provided. If *callback* is given, the callback will be executed
 916         for each position within the :term:`region`. Note that callbacks currently work
 917         only, if *region* or *reference* is given.
 918
 919         Note that a :term:`SAM` file does not allow random access. If *region* or *reference* are given,
 920         an exception is raised.
 921         '''
 922         cdef int rtid, rstart, rend, has_coord
 923
 924         if not self._isOpen():
 925             raise ValueError( "I/O operation on closed file" )
 926
 927         has_coord, rtid, rstart, rend = self._parseRegion( reference, start, end, region )
 928
 929         if self.isstream: reopen = False
 930         else: reopen = True
 931
 932         if self.isbam:
 933             if not until_eof and not self._hasIndex() and not self.isremote:
 934                 raise ValueError( "fetch called on bamfile without index" )
 935
 936             if callback:
 937                 if not has_coord: raise ValueError( "callback functionality requires a region/reference" )
 938                 if not self._hasIndex(): raise ValueError( "no index available for fetch" )
 939                 return bam_fetch(self.samfile.x.bam,
 940                                  self.index,
 941                                  rtid,
 942                                  rstart,
 943                                  rend,
 944                                  <void*>callback,
 945                                  fetch_callback )
 946             else:
 947                 if has_coord:
 948                     return IteratorRowRegion( self, rtid, rstart, rend, reopen=reopen )
 949                 else:
 950                     if until_eof:
 951                         return IteratorRowAll( self, reopen=reopen )
 952                     else:
 953                         # AH: check - reason why no reopen for AllRefs?
 954                         return IteratorRowAllRefs(self ) # , reopen=reopen )
 955         else:
 956             if has_coord:
 957                 raise ValueError ("fetching by region is not available for sam files" )
 958
 959             if callback:
 960                 raise NotImplementedError( "callback not implemented yet" )
 961
 962             if self.samfile.header == NULL:
 963                 raise ValueError( "fetch called for samfile without header")
 964
 965             # check if targets are defined
 966             # give warning, sam_read1 segfaults
 967             if self.samfile.header.n_targets == 0:
 968                 warnings.warn( "fetch called for samfile without header")
 969
 970             return IteratorRowAll( self, reopen=reopen )
 971
 972     def mate( self,
 973               AlignedRead read ):
 974         '''return the mate of :class:`AlignedRead` *read*.
 975
 976         Throws a ValueError if read is unpaired or the mate
 977         is unmapped.
 978
 979         .. note::
 980             Calling this method will change the file position.
 981             This might interfere with any iterators that have
 982             not re-opened the file.
 983
 984         '''
 985         cdef uint32_t flag = read._delegate.core.flag
 986
 987         if flag & BAM_FPAIRED == 0:
 988             raise ValueError( "read %s: is unpaired" % (read.qname))
 989         if flag & BAM_FMUNMAP != 0:
 990             raise ValueError( "mate %s: is unmapped" % (read.qname))
 991
 992         cdef MateData mate_data
 993
 994         mate_data.name = <char *>bam1_qname(read._delegate)
 995         mate_data.mate = NULL
 996         # xor flags to get the other mate
 997         cdef int x = BAM_FREAD1 + BAM_FREAD2
 998         mate_data.flag = ( flag ^ x) & x
 999
1000         bam_fetch(self.samfile.x.bam,
1001                   self.index,
1002                   read._delegate.core.mtid,
1003                   read._delegate.core.mpos,
1004                   read._delegate.core.mpos + 1,
1005                   <void*>&mate_data,
1006                   mate_callback )
1007
1008         if mate_data.mate == NULL:
1009             raise ValueError( "mate not found" )
1010
1011         cdef AlignedRead dest = AlignedRead.__new__(AlignedRead)
1012         dest._delegate = mate_data.mate
1013         return dest
1014
1015     def count( self,
1016                reference = None,
1017                start = None,
1018                end = None,
1019                region = None,
1020                until_eof = False ):
1021         '''*(reference = None, start = None, end = None, region = None, callback = None, until_eof = False)*
1022
1023         count  reads :term:`region` using 0-based indexing. The region is specified by
1024         :term:`reference`, *start* and *end*. Alternatively, a samtools :term:`region` string can be supplied.
1025
1026         Note that a :term:`TAM` file does not allow random access. If *region* or *reference* are given,
1027         an exception is raised.
1028         '''
1029         cdef int rtid
1030         cdef int rstart
1031         cdef int rend
1032
1033         if not self._isOpen():
1034             raise ValueError( "I/O operation on closed file" )
1035
1036         region, rtid, rstart, rend = self._parseRegion( reference, start, end, region )
1037
1038         cdef int counter
1039         counter = 0;
1040
1041         if self.isbam:
1042             if not until_eof and not self._hasIndex() and not self.isremote:
1043                 raise ValueError( "fetch called on bamfile without index" )
1044
1045             if not region:
1046                 raise ValueError( "counting functionality requires a region/reference" )
1047             if not self._hasIndex(): raise ValueError( "no index available for fetch" )
1048             bam_fetch(self.samfile.x.bam,
1049                              self.index,
1050                              rtid,
1051                              rstart,
1052                              rend,
1053                              <void*>&counter,
1054                              count_callback )
1055             return counter
1056         else:
1057             raise ValueError ("count for a region is not available for sam files" )
1058
1059     def pileup( self,
1060                 reference = None,
1061                 start = None,
1062                 end = None,
1063                 region = None,
1064                 callback = None,
1065                 **kwargs ):
1066         '''
1067         perform a :term:`pileup` within a :term:`region`. The region is specified by
1068         :term:`reference`, *start* and *end* (using 0-based indexing).
1069         Alternatively, a samtools *region* string can be supplied.
1070
1071         Without *reference* or *region* all reads will be used for the pileup. The reads will be returned
1072         ordered by :term:`reference` sequence, which will not necessarily be the order within the file.
1073
1074         The method returns an iterator of type :class:`pysam.IteratorColumn` unless
1075         a *callback is provided. If a *callback* is given, the callback will be executed
1076         for each column within the :term:`region`.
1077
1078         Note that :term:`SAM` formatted files do not allow random access.
1079         In these files, if a *region* or *reference* are given an exception is raised.
1080
1081         Optional *kwargs* to the iterator:
1082
1083         stepper
1084            The stepper controlls how the iterator advances.
1085            Possible options for the stepper are
1086
1087            ``all``
1088               use all reads for pileup.
1089            ``samtools``
1090               same filter and read processing as in :term:`csamtools` pileup
1091
1092         fastafile
1093            A :class:`FastaFile` object
1094
1095          mask
1096            Skip all reads with bits set in mask if mask=True.
1097
1098          max_depth
1099            Maximum read depth permitted. The default limit is *8000*.
1100
1101          truncate
1102            By default, the samtools pileup engine outputs all reads overlapping a region (see note below).
1103            If truncate is True and a region is given, only output columns in the exact region
1104            specificied.
1105
1106         .. note::
1107
1108             *all* reads which overlap the region are returned. The first base returned will be the
1109             first base of the first read *not* necessarily the first base of the region used in the query.
1110
1111         '''
1112         cdef int rtid, rstart, rend, has_coord
1113         cdef bam_plbuf_t *buf
1114
1115         if not self._isOpen():
1116             raise ValueError( "I/O operation on closed file" )
1117
1118         has_coord, rtid, rstart, rend = self._parseRegion( reference, start, end, region )
1119
1120         if self.isbam:
1121             if not self._hasIndex(): raise ValueError( "no index available for pileup" )
1122
1123             if callback:
1124                 if not has_coord: raise ValueError( "callback functionality requires a region/reference" )
1125
1126                 buf = bam_plbuf_init( <bam_pileup_f>pileup_callback, <void*>callback )
1127                 bam_fetch(self.samfile.x.bam,
1128                           self.index, rtid, rstart, rend,
1129                           buf, pileup_fetch_callback )
1130
1131                 # finalize pileup
1132                 bam_plbuf_push( NULL, buf)
1133                 bam_plbuf_destroy(buf)
1134             else:
1135                 if has_coord:
1136                     return IteratorColumnRegion( self,
1137                                                  tid = rtid,
1138                                                  start = rstart,
1139                                                  end = rend,
1140                                                  **kwargs )
1141                 else:
1142                     return IteratorColumnAllRefs(self, **kwargs )
1143
1144         else:
1145             raise NotImplementedError( "pileup of samfiles not implemented yet" )
1146
1147     def close( self ):
1148         '''
1149         closes the :class:`pysam.Samfile`.'''
1150         if self.samfile != NULL:
1151             samclose( self.samfile )
1152             bam_index_destroy(self.index);
1153             self.samfile = NULL
1154
1155     def __dealloc__( self ):
1156         # remember: dealloc cannot call other methods
1157         # note: no doc string
1158         # note: __del__ is not called.
1159         self.close()
1160         bam_destroy1(self.b)
1161         if self._filename != NULL: free( self._filename )
1162
1163     cpdef int write( self, AlignedRead read ) except -1:
1164         '''
1165         write a single :class:`pysam.AlignedRead` to disk.
1166
1167         returns the number of bytes written.
1168         '''
1169         if not self._isOpen():
1170             return 0
1171
1172         return samwrite( self.samfile, read._delegate )
1173
1174     def __enter__(self):
1175         return self
1176
1177     def __exit__(self, exc_type, exc_value, traceback):
1178         self.close()
1179         return False
1180
1181     ###############################################################
1182     ###############################################################
1183     ###############################################################
1184     ## properties
1185     ###############################################################
1186     property filename:
1187         '''number of :term:`filename` associated with this object.'''
1188         def __get__(self):
1189             if not self._isOpen(): raise ValueError( "I/O operation on closed file" )
1190             return self._filename
1191
1192     property nreferences:
1193         '''number of :term:`reference` sequences in the file.'''
1194         def __get__(self):
1195             if not self._isOpen(): raise ValueError( "I/O operation on closed file" )
1196             return self.samfile.header.n_targets
1197
1198     property references:
1199         """tuple with the names of :term:`reference` sequences."""
1200         def __get__(self):
1201             if not self._isOpen(): raise ValueError( "I/O operation on closed file" )
1202             t = []
1203             for x from 0 <= x < self.samfile.header.n_targets:
1204                 t.append( _charptr_to_str(self.samfile.header.target_name[x]) )
1205             return tuple(t)
1206
1207     property lengths:
1208         """tuple of the lengths of the :term:`reference` sequences. The lengths are in the same order as
1209         :attr:`pysam.Samfile.references`
1210         """
1211         def __get__(self):
1212             if not self._isOpen(): raise ValueError( "I/O operation on closed file" )
1213             t = []
1214             for x from 0 <= x < self.samfile.header.n_targets:
1215                 t.append( self.samfile.header.target_len[x] )
1216             return tuple(t)
1217
1218     property mapped:
1219         """total number of mapped reads in file.
1220         """
1221         def __get__(self):
1222             if not self._isOpen(): raise ValueError( "I/O operation on closed file" )
1223             if not self.isbam: raise AttributeError( "Samfile.mapped only available in bam files" )
1224
1225             cdef int tid
1226             cdef uint32_t total = 0
1227             for tid from 0 <= tid < self.samfile.header.n_targets:
1228                 total += pysam_get_mapped( self.index, tid )
1229             return total
1230
1231     property unmapped:
1232         """total number of unmapped reads in file.
1233         """
1234         def __get__(self):
1235             if not self._isOpen(): raise ValueError( "I/O operation on closed file" )
1236             if not self.isbam: raise AttributeError( "Samfile.unmapped only available in bam files" )
1237             cdef int tid
1238             cdef uint32_t total = 0
1239             for tid from 0 <= tid < self.samfile.header.n_targets:
1240                 total += pysam_get_unmapped( self.index, tid )
1241             # get unmapped reads without coordinates
1242             total += pysam_get_unmapped( self.index, -1 )
1243             return total
1244
1245     property text:
1246         '''full contents of the :term:`sam file` header as a string.'''
1247         def __get__(self):
1248             if not self._isOpen(): raise ValueError( "I/O operation on closed file" )
1249             return from_string_and_size(self.samfile.header.text, self.samfile.header.l_text)
1250
1251     property header:
1252         '''header information within the :term:`sam file`. The records and fields are returned as
1253         a two-level dictionary.
1254         '''
1255         def __get__(self):
1256             if not self._isOpen(): raise ValueError( "I/O operation on closed file" )
1257
1258             result = {}
1259
1260             if self.samfile.header.text != NULL:
1261                 # convert to python string (note: call self.text to create 0-terminated string)
1262                 t = self.text
1263                 for line in t.split("\n"):
1264                     if not line.strip(): continue
1265                     assert line.startswith("@"), "header line without '@': '%s'" % line
1266                     fields = line[1:].split("\t")
1267                     record = fields[0]
1268                     assert record in VALID_HEADER_TYPES, "header line with invalid type '%s': '%s'" % (record, line)
1269
1270                     # treat comments
1271                     if record == "CO":
1272                         if record not in result: result[record] = []
1273                         result[record].append( "\t".join( fields[1:] ) )
1274                         continue
1275
1276                     # the following is clumsy as generators do not work?
1277                     x = {}
1278                     for field in fields[1:]:
1279                         key, value = field.split(":",1)
1280                         # uppercase keys must be valid
1281                         # lowercase are permitted for user fields
1282                         if key in VALID_HEADER_FIELDS[record]:
1283                             x[key] = VALID_HEADER_FIELDS[record][key](value)
1284                         elif not key.isupper():
1285                             x[key] = value
1286                         else:
1287                             raise ValueError( "unknown field code '%s' in record '%s'" % (key, record) )
1288
1289                     if VALID_HEADER_TYPES[record] == dict:
1290                         if record in result:
1291                             raise ValueError( "multiple '%s' lines are not permitted" % record )
1292                         result[record] = x
1293                     elif VALID_HEADER_TYPES[record] == list:
1294                         if record not in result: result[record] = []
1295                         result[record].append( x )
1296
1297                 # if there are no SQ lines in the header, add the reference names
1298                 # from the information in the bam file.
1299                 # Background: c-samtools keeps the textual part of the header separate from
1300                 # the list of reference names and lengths. Thus, if a header contains only
1301                 # SQ lines, the SQ information is not part of the textual header and thus
1302                 # are missing from the output. See issue 84.
1303                 if "SQ" not in result:
1304                     sq = []
1305                     for ref, length in zip( self.references, self.lengths ):
1306                         sq.append( {'LN': length, 'SN': ref } )
1307                     result["SQ"] = sq
1308
1309             return result
1310
1311     def _buildLine( self, fields, record ):
1312         '''build a header line from *fields* dictionary for *record*'''
1313
1314         # TODO: add checking for field and sort order
1315         line = ["@%s" % record ]
1316         # comment
1317         if record == "CO":
1318             line.append( fields )
1319         # user tags
1320         elif record.islower():
1321             for key in sorted(fields):
1322                 line.append( "%s:%s" % (key, str(fields[key])))
1323         # defined tags
1324         else:
1325             # write fields of the specification
1326             for key in VALID_HEADER_ORDER[record]:
1327                 if key in fields:
1328                     line.append( "%s:%s" % (key, str(fields[key])))
1329             # write user fields
1330             for key in fields:
1331                 if not key.isupper():
1332                     line.append( "%s:%s" % (key, str(fields[key])))
1333
1334         return "\t".join( line )
1335
1336     cdef bam_header_t * _buildHeader( self, new_header ):
1337         '''return a new header built from a dictionary in *new_header*.
1338
1339         This method inserts the text field, target_name and target_len.
1340         '''
1341
1342         lines = []
1343
1344         # check if hash exists
1345
1346         # create new header and copy old data
1347         cdef bam_header_t * dest
1348
1349         dest = bam_header_init()
1350
1351         # first: defined tags
1352         for record in VALID_HEADERS:
1353             if record in new_header:
1354                 ttype = VALID_HEADER_TYPES[record]
1355                 data = new_header[record]
1356                 if type( data ) != type( ttype() ):
1357                     raise ValueError( "invalid type for record %s: %s, expected %s" % (record, type(data), type(ttype()) ) )
1358                 if type( data ) is dict:
1359                     lines.append( self._buildLine( data, record ) )
1360                 else:
1361                     for fields in new_header[record]:
1362                         lines.append( self._buildLine( fields, record ) )
1363
1364         # then: user tags (lower case), sorted alphabetically
1365         for record, data in sorted(new_header.items()):
1366             if record in VALID_HEADERS: continue
1367             if type( data ) is dict:
1368                 lines.append( self._buildLine( data, record ) )
1369             else:
1370                 for fields in new_header[record]:
1371                     lines.append( self._buildLine( fields, record ) )
1372
1373         text = "\n".join(lines) + "\n"
1374         if dest.text != NULL: free( dest.text )
1375         dest.text = <char*>calloc( len(text), sizeof(char))
1376         dest.l_text = len(text)
1377         cdef bytes btext = text.encode('ascii')
1378         strncpy( dest.text, btext, dest.l_text )
1379
1380         cdef bytes bseqname
1381         # collect targets
1382         if "SQ" in new_header:
1383             seqs = []
1384             for fields in new_header["SQ"]:
1385                 try:
1386                     seqs.append( (fields["SN"], fields["LN"] ) )
1387                 except KeyError:
1388                     raise KeyError( "incomplete sequence information in '%s'" % str(fields))
1389
1390             dest.n_targets = len(seqs)
1391             dest.target_name = <char**>calloc( dest.n_targets, sizeof(char*) )
1392             dest.target_len = <uint32_t*>calloc( dest.n_targets, sizeof(uint32_t) )
1393
1394             for x from 0 <= x < dest.n_targets:
1395                 seqname, seqlen = seqs[x]
1396                 dest.target_name[x] = <char*>calloc( len( seqname ) + 1, sizeof(char) )
1397                 bseqname = seqname.encode('ascii')
1398                 strncpy( dest.target_name[x], bseqname, len(seqname) + 1 )
1399                 dest.target_len[x] = seqlen
1400
1401         return dest
1402
1403     ###############################################################
1404     ###############################################################
1405     ###############################################################
1406     ## file-object like iterator access
1407     ## note: concurrent access will cause errors (see IteratorRow
1408     ## and reopen)
1409     ## Possible solutions: deprecate or open new file handle
1410     ###############################################################
1411     def __iter__(self):
1412         if not self._isOpen(): raise ValueError( "I/O operation on closed file" )
1413         if not self.isbam and self.samfile.header.n_targets == 0:
1414                 raise NotImplementedError( "can not iterate over samfile without header")
1415         return self
1416
1417     cdef bam1_t * getCurrent( self ):
1418         return self.b
1419
1420     cdef int cnext(self):
1421         '''
1422         cversion of iterator. Used by :class:`pysam.Samfile.IteratorColumn`.
1423         '''
1424         cdef int ret
1425         return samread(self.samfile, self.b)
1426
1427     def __next__(self):
1428         """
1429         python version of next().
1430         """
1431         cdef int ret
1432         ret = samread(self.samfile, self.b)
1433         if (ret > 0):
1434             return makeAlignedRead( self.b )
1435         else:
1436             raise StopIteration
1437
1438 ##-------------------------------------------------------------------
1439 ##-------------------------------------------------------------------
1440 ##-------------------------------------------------------------------
1441 cdef class IteratorRow:
1442     '''abstract base class for iterators over mapped reads.
1443
1444     Various iterators implement different behaviours for wrapping around
1445     contig boundaries. Examples include:
1446
1447     :class:`pysam.IteratorRowRegion`
1448         iterate within a single contig and a defined region.
1449
1450     :class:`pysam.IteratorRowAll`
1451         iterate until EOF. This iterator will also include unmapped reads.
1452
1453     :class:`pysam.IteratorRowAllRefs`
1454         iterate over all reads in all reference sequences.
1455
1456     The method :meth:`Samfile.fetch` returns an IteratorRow.
1457     '''
1458     pass
1459
1460
1461 cdef class IteratorRowRegion(IteratorRow):
1462     """*(Samfile samfile, int tid, int beg, int end, int reopen = True )*
1463
1464     iterate over mapped reads in a region.
1465
1466     By default, the file is re-openend to avoid conflicts between
1467     multiple iterators working on the same file. Set *reopen* = False
1468     to not re-open *samfile*.
1469
1470     The samtools iterators assume that the file
1471     position between iterations do not change.
1472     As a consequence, no two iterators can work
1473     on the same file. To permit this, each iterator
1474     creates its own file handle by re-opening the
1475     file.
1476
1477     Note that the index will be shared between
1478     samfile and the iterator.
1479     """
1480
1481     def __cinit__(self, Samfile samfile, int tid, int beg, int end, int reopen = True ):
1482
1483         if not samfile._isOpen():
1484             raise ValueError( "I/O operation on closed file" )
1485
1486         if not samfile._hasIndex():
1487             raise ValueError( "no index available for iteration" )
1488
1489         # makes sure that samfile stays alive as long as the
1490         # iterator is alive
1491         self.samfile = samfile
1492
1493         if samfile.isbam: mode = b"rb"
1494         else: mode = b"r"
1495
1496         # reopen the file - note that this makes the iterator
1497         # slow and causes pileup to slow down significantly.
1498         if reopen:
1499             store = StderrStore()
1500             self.fp = samopen( samfile._filename, mode, NULL )
1501             store.release()
1502             assert self.fp != NULL
1503             self.owns_samfile = True
1504         else:
1505             self.fp = self.samfile.samfile
1506             self.owns_samfile = False
1507
1508         self.retval = 0
1509
1510         self.iter = bam_iter_query(self.samfile.index,
1511                                    tid,
1512                                    beg,
1513                                    end)
1514         self.b = bam_init1()
1515
1516     def __iter__(self):
1517         return self
1518
1519     cdef bam1_t * getCurrent( self ):
1520         return self.b
1521
1522     cdef int cnext(self):
1523         '''cversion of iterator. Used by IteratorColumn'''
1524         self.retval = bam_iter_read( self.fp.x.bam,
1525                                      self.iter,
1526                                      self.b)
1527
1528     def __next__(self):
1529         """python version of next().
1530         """
1531         self.cnext()
1532         if self.retval < 0: raise StopIteration
1533         return makeAlignedRead( self.b )
1534
1535     def __dealloc__(self):
1536         bam_destroy1(self.b)
1537         bam_iter_destroy( self.iter )
1538         if self.owns_samfile: samclose( self.fp )
1539
1540 cdef class IteratorRowAll(IteratorRow):
1541     """*(Samfile samfile, int reopen = True)*
1542
1543     iterate over all reads in *samfile*
1544
1545     By default, the file is re-openend to avoid conflicts between
1546     multiple iterators working on the same file. Set *reopen* = False
1547     to not re-open *samfile*.
1548     """
1549
1550     def __cinit__(self, Samfile samfile, int reopen = True ):
1551
1552         if not samfile._isOpen():
1553             raise ValueError( "I/O operation on closed file" )
1554
1555         if samfile.isbam: mode = b"rb"
1556         else: mode = b"r"
1557
1558         # reopen the file to avoid iterator conflict
1559         if reopen:
1560             store = StderrStore()
1561             self.fp = samopen( samfile._filename, mode, NULL )
1562             store.release()
1563             assert self.fp != NULL
1564             self.owns_samfile = True
1565         else:
1566             self.fp = samfile.samfile
1567             self.owns_samfile = False
1568
1569         # allocate memory for alignment
1570         self.b = <bam1_t*>calloc(1, sizeof(bam1_t))
1571
1572     def __iter__(self):
1573         return self
1574
1575     cdef bam1_t * getCurrent( self ):
1576         return self.b
1577
1578     cdef int cnext(self):
1579         '''cversion of iterator. Used by IteratorColumn'''
1580         return samread(self.fp, self.b)
1581
1582     def __next__(self):
1583         """python version of next().
1584
1585         pyrex uses this non-standard name instead of next()
1586         """
1587         cdef int ret
1588         ret = samread(self.fp, self.b)
1589         if (ret > 0):
1590             return makeAlignedRead( self.b )
1591         else:
1592             raise StopIteration
1593
1594     def __dealloc__(self):
1595         bam_destroy1(self.b)
1596         if self.owns_samfile: samclose( self.fp )
1597
1598 cdef class IteratorRowAllRefs(IteratorRow):
1599     """iterates over all mapped reads by chaining iterators over each reference
1600     """
1601
1602     def __cinit__(self, Samfile samfile):
1603         assert samfile._isOpen()
1604         if not samfile._hasIndex(): raise ValueError("no index available for fetch")
1605         self.samfile = samfile
1606         self.tid = -1
1607
1608     def nextiter(self):
1609         self.rowiter = IteratorRowRegion(self.samfile, self.tid, 0, 1<<29)
1610
1611     def __iter__(self):
1612         return self
1613
1614     def __next__(self):
1615         """python version of next().
1616
1617         pyrex uses this non-standard name instead of next()
1618         """
1619         # Create an initial iterator
1620         if self.tid==-1:
1621             if not self.samfile.nreferences:
1622                 raise StopIteration
1623             self.tid = 0
1624             self.nextiter()
1625
1626         while 1:
1627             self.rowiter.cnext()
1628
1629             # If current iterator is not exhausted, return aligned read
1630             if self.rowiter.retval>0:
1631                 return makeAlignedRead(self.rowiter.b)
1632
1633             self.tid += 1
1634
1635             # Otherwise, proceed to next reference or stop
1636             if self.tid<self.samfile.nreferences:
1637                 self.nextiter()
1638             else:
1639                 raise StopIteration
1640
1641 cdef class IteratorRowSelection(IteratorRow):
1642     """*(Samfile samfile)*
1643
1644     iterate over reads in *samfile* at a given list of file positions.
1645     """
1646
1647     def __cinit__(self, Samfile samfile, positions, int reopen = True ):
1648
1649         if not samfile._isOpen():
1650             raise ValueError( "I/O operation on closed file" )
1651
1652         if not samfile._isOpen():
1653             raise ValueError( "I/O operation on closed file" )
1654
1655         assert samfile.isbam, "can only use this iterator on bam files"
1656         mode = b"rb"
1657
1658         # reopen the file to avoid iterator conflict
1659         if reopen:
1660             store = StderrStore()
1661             self.fp = samopen( samfile._filename, mode, NULL )
1662             store.release()
1663             assert self.fp != NULL
1664             self.owns_samfile = True
1665         else:
1666             self.fp = samfile.samfile
1667             self.owns_samfile = False
1668
1669         # allocate memory for alignment
1670         self.b = <bam1_t*>calloc(1, sizeof(bam1_t))
1671
1672         self.positions = positions
1673         self.current_pos = 0
1674
1675     def __iter__(self):
1676         return self
1677
1678     cdef bam1_t * getCurrent( self ):
1679         return self.b
1680
1681     cdef int cnext(self):
1682         '''cversion of iterator'''
1683
1684         # end iteration if out of positions
1685         if self.current_pos >= len(self.positions): return -1
1686
1687         bam_seek( self.fp.x.bam, self.positions[self.current_pos], 0 )
1688         self.current_pos += 1
1689         return samread(self.fp, self.b)
1690
1691     def __next__(self):
1692         """python version of next().
1693
1694         pyrex uses this non-standard name instead of next()
1695         """
1696
1697         cdef int ret = self.cnext()
1698         if (ret > 0):
1699             return makeAlignedRead( self.b )
1700         else:
1701             raise StopIteration
1702
1703     def __dealloc__(self):
1704         bam_destroy1(self.b)
1705         if self.owns_samfile: samclose( self.fp )
1706
1707 ##-------------------------------------------------------------------
1708 ##-------------------------------------------------------------------
1709 ##-------------------------------------------------------------------
1710 cdef int __advance_all( void * data, bam1_t * b ):
1711     '''advance without any read filtering.
1712     '''
1713     cdef __iterdata * d
1714     d = <__iterdata*>data
1715     return bam_iter_read( d.samfile.x.bam, d.iter, b )
1716
1717 cdef int __advance_snpcalls( void * data, bam1_t * b ):
1718     '''advance using same filter and read processing as in
1719     the samtools pileup.
1720     '''
1721     cdef __iterdata * d
1722     d = <__iterdata*>data
1723
1724     cdef int ret = bam_iter_read( d.samfile.x.bam, d.iter, b )
1725     cdef int skip = 0
1726     cdef int q
1727     cdef int is_cns = 1
1728     cdef int is_nobaq = 0
1729     cdef int capQ_thres = 0
1730
1731     # reload sequence
1732     if d.fastafile != NULL and b.core.tid != d.tid:
1733         if d.seq != NULL: free(d.seq)
1734         d.tid = b.core.tid
1735         d.seq = faidx_fetch_seq(d.fastafile,
1736                                 d.samfile.header.target_name[d.tid],
1737                                 0, max_pos,
1738                                 &d.seq_len)
1739         if d.seq == NULL:
1740             raise ValueError( "reference sequence for '%s' (tid=%i) not found" % \
1741                                   (d.samfile.header.target_name[d.tid],
1742                                    d.tid))
1743
1744
1745     while ret >= 0:
1746
1747         skip = 0
1748
1749         # realign read - changes base qualities
1750         if d.seq != NULL and is_cns and not is_nobaq: bam_prob_realn( b, d.seq )
1751
1752         if d.seq != NULL and capQ_thres > 10:
1753             q = bam_cap_mapQ(b, d.seq, capQ_thres)
1754             if q < 0: skip = 1
1755             elif b.core.qual > q: b.core.qual = q
1756         if b.core.flag & BAM_FUNMAP: skip = 1
1757         elif b.core.flag & 1 and not b.core.flag & 2: skip = 1
1758
1759         if not skip: break
1760         # additional filters
1761
1762         ret = bam_iter_read( d.samfile.x.bam, d.iter, b )
1763
1764     return ret
1765
1766 cdef class IteratorColumn:
1767     '''abstract base class for iterators over columns.
1768
1769     IteratorColumn objects wrap the pileup functionality of samtools.
1770
1771     For reasons of efficiency, the iterator points to the current
1772     pileup buffer. The pileup buffer is updated at every iteration.
1773     This might cause some unexpected behavious. For example,
1774     consider the conversion to a list::
1775
1776        f = Samfile("file.bam", "rb")
1777        result = list( f.pileup() )
1778
1779     Here, ``result`` will contain ``n`` objects of type :class:`PileupProxy` for ``n`` columns,
1780     but each object in ``result`` will contain the same information.
1781
1782     The desired behaviour can be achieved by list comprehension::
1783
1784        result = [ x.pileups() for x in f.pileup() ]
1785
1786     ``result`` will be a list of ``n`` lists of objects of type :class:`PileupRead`.
1787
1788     If the iterator is associated with a :class:`Fastafile` using the :meth:`addReference`
1789     method, then the iterator will export the current sequence via the methods :meth:`getSequence`
1790     and :meth:`seq_len`.
1791
1792     Optional kwargs to the iterator
1793
1794     stepper
1795        The stepper controls how the iterator advances.
1796        Possible options for the stepper are
1797
1798        all
1799            use all reads for pileup.
1800        samtools
1801            same filter and read processing as in :term:`csamtools` pileup
1802
1803        The default is to use "all" if no stepper is given.
1804
1805     fastafile
1806        A :class:`FastaFile` object
1807     mask
1808        Skip all reads with bits set in mask.
1809     max_depth
1810        maximum read depth. The default is 8000.
1811     '''
1812
1813     def __cinit__( self, Samfile samfile, **kwargs ):
1814         self.samfile = samfile
1815         self.mask = kwargs.get("mask", BAM_DEF_MASK )
1816         self.fastafile = kwargs.get( "fastafile", None )
1817         self.stepper = kwargs.get( "stepper", None )
1818         self.max_depth = kwargs.get( "max_depth", 8000 )
1819         self.iterdata.seq = NULL
1820         self.tid = 0
1821         self.pos = 0
1822         self.n_plp = 0
1823         self.plp = NULL
1824         self.pileup_iter = <bam_plp_t>NULL
1825
1826
1827     def __iter__(self):
1828         return self
1829
1830     cdef int cnext(self):
1831         '''perform next iteration.
1832
1833         This method is analogous to the samtools bam_plp_auto method.
1834         It has been re-implemented to permit for filtering.
1835         '''
1836         self.plp = bam_plp_auto( self.pileup_iter,
1837                                  &self.tid,
1838                                  &self.pos,
1839                                  &self.n_plp )
1840
1841     cdef char * getSequence( self ):
1842         '''return current reference sequence underlying the iterator.
1843         '''
1844         return self.iterdata.seq
1845
1846     property seq_len:
1847         '''current sequence length.'''
1848         def __get__(self): return self.iterdata.seq_len
1849
1850     def addReference( self, Fastafile fastafile ):
1851        '''
1852        add reference sequences in *fastafile* to iterator.'''
1853        self.fastafile = fastafile
1854        if self.iterdata.seq != NULL: free(self.iterdata.seq)
1855        self.iterdata.tid = -1
1856        self.iterdata.fastafile = self.fastafile.fastafile
1857
1858     def hasReference( self ):
1859         '''
1860         return true if iterator is associated with a reference'''
1861         return self.fastafile
1862
1863     cdef setMask( self, mask ):
1864         '''set masking flag in iterator.
1865
1866         reads with bits set in *mask* will be skipped.
1867         '''
1868         self.mask = mask
1869         bam_plp_set_mask( self.pileup_iter, self.mask )
1870
1871     cdef setupIteratorData( self,
1872                             int tid,
1873                             int start,
1874                             int end,
1875                             int reopen = 0 ):
1876         '''setup the iterator structure'''
1877
1878         self.iter = IteratorRowRegion( self.samfile, tid, start, end, reopen )
1879         self.iterdata.samfile = self.samfile.samfile
1880         self.iterdata.iter = self.iter.iter
1881         self.iterdata.seq = NULL
1882         self.iterdata.tid = -1
1883
1884         if self.fastafile != None:
1885             self.iterdata.fastafile = self.fastafile.fastafile
1886         else:
1887             self.iterdata.fastafile = NULL
1888
1889         if self.stepper == None or self.stepper == "all":
1890             self.pileup_iter = bam_plp_init( &__advance_all, &self.iterdata )
1891         elif self.stepper == "samtools":
1892             self.pileup_iter = bam_plp_init( &__advance_snpcalls, &self.iterdata )
1893         else:
1894             raise ValueError( "unknown stepper option `%s` in IteratorColumn" % self.stepper)
1895
1896         if self.max_depth:
1897             bam_plp_set_maxcnt( self.pileup_iter, self.max_depth )
1898
1899         bam_plp_set_mask( self.pileup_iter, self.mask )
1900
1901     cdef reset( self, tid, start, end ):
1902         '''reset iterator position.
1903
1904         This permits using the iterator multiple times without
1905         having to incur the full set-up costs.
1906         '''
1907         self.iter = IteratorRowRegion( self.samfile, tid, start, end, reopen = 0 )
1908         self.iterdata.iter = self.iter.iter
1909
1910         # invalidate sequence if different tid
1911         if self.tid != tid:
1912             if self.iterdata.seq != NULL: free( self.iterdata.seq )
1913             self.iterdata.seq = NULL
1914             self.iterdata.tid = -1
1915
1916         # self.pileup_iter = bam_plp_init( &__advancepileup, &self.iterdata )
1917         bam_plp_reset(self.pileup_iter)
1918
1919     def __dealloc__(self):
1920         # reset in order to avoid memory leak messages for iterators that have
1921         # not been fully consumed
1922         if self.pileup_iter != <bam_plp_t>NULL:
1923             bam_plp_reset(self.pileup_iter)
1924             bam_plp_destroy(self.pileup_iter)
1925             self.pileup_iter = <bam_plp_t>NULL
1926
1927         if self.iterdata.seq != NULL:
1928             free(self.iterdata.seq)
1929             self.iterdata.seq = NULL
1930
1931 cdef class IteratorColumnRegion(IteratorColumn):
1932     '''iterates over a region only.
1933     '''
1934     def __cinit__(self, Samfile samfile,
1935                   int tid = 0,
1936                   int start = 0,
1937                   int end = max_pos,
1938                   int truncate = False,
1939                   **kwargs ):
1940
1941         # initialize iterator
1942         self.setupIteratorData( tid, start, end, 1 )
1943         self.start = start
1944         self.end = end
1945         self.truncate = truncate
1946
1947     def __next__(self):
1948         """python version of next().
1949         """
1950
1951         while 1:
1952             self.cnext()
1953             if self.n_plp < 0:
1954                 raise ValueError("error during iteration" )
1955
1956             if self.plp == NULL:
1957                 raise StopIteration
1958
1959             if self.truncate:
1960                 if self.start < self.pos: continue
1961                 if self.pos >= self.end: raise StopIteration
1962
1963             return makePileupProxy( &self.plp,
1964                                      self.tid,
1965                                      self.pos,
1966                                      self.n_plp )
1967
1968 cdef class IteratorColumnAllRefs(IteratorColumn):
1969     """iterates over all columns by chaining iterators over each reference
1970     """
1971
1972     def __cinit__(self,
1973                   Samfile samfile,
1974                   **kwargs ):
1975
1976         # no iteration over empty files
1977         if not samfile.nreferences: raise StopIteration
1978
1979         # initialize iterator
1980         self.setupIteratorData( self.tid, 0, max_pos, 1 )
1981
1982     def __next__(self):
1983         """python version of next().
1984         """
1985
1986         while 1:
1987             self.cnext()
1988
1989             if self.n_plp < 0:
1990                 raise ValueError("error during iteration" )
1991
1992             # return result, if within same reference
1993             if self.plp != NULL:
1994                 return makePileupProxy( &self.plp,
1995                                          self.tid,
1996                                          self.pos,
1997                                          self.n_plp )
1998
1999             # otherwise, proceed to next reference or stop
2000             self.tid += 1
2001             if self.tid < self.samfile.nreferences:
2002                 self.setupIteratorData( self.tid, 0, max_pos, 0 )
2003             else:
2004                 raise StopIteration
2005
2006 ##-------------------------------------------------------------------
2007 ##-------------------------------------------------------------------
2008 ##-------------------------------------------------------------------
2009 cdef inline int32_t query_start(bam1_t *src) except -1:
2010     cdef uint32_t * cigar_p, op
2011     cdef uint32_t k
2012     cdef uint32_t start_offset = 0
2013
2014     if src.core.n_cigar:
2015         cigar_p = bam1_cigar(src);
2016         for k from 0 <= k < src.core.n_cigar:
2017             op = cigar_p[k] & BAM_CIGAR_MASK
2018             if op==BAM_CHARD_CLIP:
2019                 if start_offset!=0 and start_offset!=src.core.l_qseq:
2020                     PyErr_SetString(ValueError, 'Invalid clipping in CIGAR string')
2021                     return -1
2022             elif op==BAM_CSOFT_CLIP:
2023                 start_offset += cigar_p[k] >> BAM_CIGAR_SHIFT
2024             else:
2025                 break
2026
2027     return start_offset
2028
2029 ##-------------------------------------------------------------------
2030 ##-------------------------------------------------------------------
2031 ##-------------------------------------------------------------------
2032 cdef inline int32_t query_end(bam1_t *src) except -1:
2033     cdef uint32_t * cigar_p, op
2034     cdef uint32_t k
2035     cdef uint32_t end_offset = src.core.l_qseq
2036
2037     if src.core.n_cigar>1:
2038         cigar_p = bam1_cigar(src);
2039         for k from src.core.n_cigar > k >= 1:
2040             op = cigar_p[k] & BAM_CIGAR_MASK
2041             if op==BAM_CHARD_CLIP:
2042                 if end_offset!=0 and end_offset!=src.core.l_qseq:
2043                     PyErr_SetString(ValueError, 'Invalid clipping in CIGAR string')
2044                     return -1
2045             elif op==BAM_CSOFT_CLIP:
2046                 end_offset -= cigar_p[k] >> BAM_CIGAR_SHIFT
2047             else:
2048                 break
2049
2050     if end_offset==0:
2051         end_offset = src.core.l_qseq
2052
2053     return end_offset
2054
2055
2056 cdef inline object get_seq_range(bam1_t *src, uint32_t start, uint32_t end):
2057     cdef uint8_t * p
2058     cdef uint32_t k
2059     cdef char * s
2060
2061     if not src.core.l_qseq:
2062         return None
2063
2064     seq = PyBytes_FromStringAndSize(NULL, end - start)
2065     s   = <char*>seq
2066     p   = bam1_seq(src)
2067
2068     for k from start <= k < end:
2069         # equivalent to bam_nt16_rev_table[bam1_seqi(s, i)] (see bam.c)
2070         # note: do not use string literal as it will be a python string
2071         s[k-start] = bam_nt16_rev_table[p[k/2] >> 4 * (1 - k%2) & 0xf]
2072
2073     return seq
2074
2075
2076 cdef inline object get_qual_range(bam1_t *src, uint32_t start, uint32_t end):
2077     cdef uint8_t * p
2078     cdef uint32_t k
2079     cdef char * q
2080
2081     p = bam1_qual(src)
2082     if p[0] == 0xff:
2083         return None
2084
2085     qual = PyBytes_FromStringAndSize(NULL, end - start)
2086     q    = <char*>qual
2087
2088     for k from start <= k < end:
2089         ## equivalent to t[i] + 33 (see bam.c)
2090         q[k-start] = p[k] + 33
2091
2092     return qual
2093
2094 cdef class AlignedRead:
2095     '''
2096     Class representing an aligned read. see SAM format specification for
2097     the meaning of fields (http://samtools.sourceforge.net/).
2098
2099     This class stores a handle to the samtools C-structure representing
2100     an aligned read. Member read access is forwarded to the C-structure
2101     and converted into python objects. This implementation should be fast,
2102     as only the data needed is converted.
2103
2104     For write access, the C-structure is updated in-place. This is
2105     not the most efficient way to build BAM entries, as the variable
2106     length data is concatenated and thus needs to resized if
2107     a field is updated. Furthermore, the BAM entry might be
2108     in an inconsistent state. The :meth:`~validate` method can
2109     be used to check if an entry is consistent.
2110
2111     One issue to look out for is that the sequence should always
2112     be set *before* the quality scores. Setting the sequence will
2113     also erase any quality scores that were set previously.
2114
2115     In Python 3, the fields containing sequence and quality
2116     (seq, query, qual and qqual) data are of type bytes. Other
2117     string data, such as the qname field and strings in the
2118     tags tuple, is represented as unicode strings. On assignment,
2119     both bytes and unicode objects are allowed, but unicode strings
2120     must contain only ASCII characters.
2121     '''
2122
2123     # Now only called when instances are created from Python
2124     def __init__(self):
2125         # see bam_init1
2126         self._delegate = <bam1_t*>calloc( 1, sizeof( bam1_t) )
2127         # allocate some memory
2128         # If size is 0, calloc does not return a pointer that can be passed to free()
2129         # so allocate 40 bytes for a new read
2130         self._delegate.m_data = 40
2131         self._delegate.data = <uint8_t *>calloc( self._delegate.m_data, 1 )
2132         self._delegate.data_len = 0
2133
2134     def __dealloc__(self):
2135         bam_destroy1(self._delegate)
2136
2137     def __str__(self):
2138         """return string representation of alignment.
2139
2140         The representation is an approximate :term:`sam` format.
2141
2142         An aligned read might not be associated with a :term:`Samfile`.
2143         As a result :term:`tid` is shown instead of the reference name.
2144
2145         Similarly, the tags field is returned in its parsed state.
2146         """
2147         # sam-parsing is done in sam.c/bam_format1_core which
2148         # requires a valid header.
2149         if sys.version_info[0] < 3:
2150             seq = self.seq
2151             qual = self.qual
2152         else:
2153             seq = self.seq.decode('ascii')
2154             qual = self.qual.decode('ascii')
2155         return "\t".join(map(str, (self.qname,
2156                                    self.flag,
2157                                    self.rname,
2158                                    self.pos,
2159                                    self.mapq,
2160                                    self.cigar,
2161                                    self.mrnm,
2162                                    self.mpos,
2163                                    self.rlen,
2164                                    seq,
2165                                    qual,
2166                                    self.tags )))
2167
2168     def compare(self, AlignedRead other):
2169         '''return -1,0,1, if contents in this are binary <,=,> to *other*'''
2170
2171         cdef int retval, x
2172         cdef bam1_t *t, *o
2173
2174         t = self._delegate
2175         o = other._delegate
2176
2177         # uncomment for debugging purposes
2178         # cdef unsigned char * oo, * tt
2179         # tt = <unsigned char*>(&t.core)
2180         # oo = <unsigned char*>(&o.core)
2181         # for x from 0 <= x < sizeof( bam1_core_t): print x, tt[x], oo[x]
2182         # tt = <unsigned char*>(t.data)
2183         # oo = <unsigned char*>(o.data)
2184         # for x from 0 <= x < max(t.data_len, o.data_len): print x, tt[x], oo[x], chr(tt[x]), chr(oo[x])
2185
2186         # Fast-path test for object identity
2187         if t==o:
2188             return 0
2189
2190         retval = memcmp(&t.core, &o.core, sizeof(bam1_core_t))
2191
2192         if retval: return retval
2193         retval = (t.data_len > o.data_len) - (t.data_len < o.data_len) # cmp(t.data_len, o.data_len)
2194         if retval: return retval
2195         return memcmp(t.data, o.data, t.data_len)
2196
2197     # Disabled so long as __cmp__ is a special method
2198     def __hash__(self):
2199         return _Py_HashPointer(<void *>self)
2200
2201     property qname:
2202         """the query name (None if not present)"""
2203         def __get__(self):
2204             cdef bam1_t * src
2205             src = self._delegate
2206             if src.core.l_qname == 0: return None
2207             return _charptr_to_str(<char *>bam1_qname( src ))
2208
2209         def __set__(self, qname ):
2210             if qname == None or len(qname) == 0: return
2211             qname = _force_bytes(qname)
2212             cdef bam1_t * src
2213             cdef int l
2214             cdef char * p
2215
2216             src = self._delegate
2217             p = bam1_qname( src )
2218
2219             # the qname is \0 terminated
2220             l = len(qname) + 1
2221             pysam_bam_update( src,
2222                               src.core.l_qname,
2223                               l,
2224                               <uint8_t*>p )
2225
2226             src.core.l_qname = l
2227
2228             # re-acquire pointer to location in memory
2229             # as it might have moved
2230             p = bam1_qname(src)
2231
2232             strncpy( p, qname, l )
2233
2234     property cigar:
2235         """the :term:`cigar` alignment (None if not present). The alignment
2236         is returned as a list of operations. The operations are:
2237
2238         +-----+--------------+-----+
2239         |M    |BAM_CMATCH    |0    |
2240         +-----+--------------+-----+
2241         |I    |BAM_CINS      |1    |
2242         +-----+--------------+-----+
2243         |D    |BAM_CDEL      |2    |
2244         +-----+--------------+-----+
2245         |N    |BAM_CREF_SKIP |3    |
2246         +-----+--------------+-----+
2247         |S    |BAM_CSOFT_CLIP|4    |
2248         +-----+--------------+-----+
2249         |H    |BAM_CHARD_CLIP|5    |
2250         +-----+--------------+-----+
2251         |P    |BAM_CPAD      |6    |
2252         +-----+--------------+-----+
2253         |=    |BAM_CEQUAL    |7    |
2254         +-----+--------------+-----+
2255         |X    |BAM_CDIFF     |8    |
2256         +-----+--------------+-----+
2257
2258         """
2259         def __get__(self):
2260             cdef uint32_t * cigar_p
2261             cdef bam1_t * src
2262             cdef op, l, cigar
2263             cdef int k
2264
2265             src = self._delegate
2266             if src.core.n_cigar == 0: return None
2267
2268             cigar = []
2269             cigar_p = bam1_cigar(src);
2270             for k from 0 <= k < src.core.n_cigar:
2271                 op = cigar_p[k] & BAM_CIGAR_MASK
2272                 l = cigar_p[k] >> BAM_CIGAR_SHIFT
2273                 cigar.append((op, l))
2274             return cigar
2275
2276         def __set__(self, values ):
2277             if values == None or len(values) == 0: return
2278             cdef uint32_t * p
2279             cdef bam1_t * src
2280             cdef op, l
2281             cdef int k
2282
2283             k = 0
2284
2285             src = self._delegate
2286
2287             # get location of cigar string
2288             p = bam1_cigar(src)
2289
2290             # create space for cigar data within src.data
2291             pysam_bam_update( src,
2292                               src.core.n_cigar * 4,
2293                               len(values) * 4,
2294                               <uint8_t*>p )
2295
2296             # length is number of cigar operations, not bytes
2297             src.core.n_cigar = len(values)
2298
2299             # re-acquire pointer to location in memory
2300             # as it might have moved
2301             p = bam1_cigar(src)
2302
2303             # insert cigar operations
2304             for op, l in values:
2305                 p[k] = l << BAM_CIGAR_SHIFT | op
2306                 k += 1
2307
2308             ## setting the cigar string also updates the "bin" attribute
2309             src.core.bin = bam_reg2bin( src.core.pos, bam_calend( &src.core, p))
2310
2311     property cigarstring:
2312         '''the :term:`cigar` alignment as a string.
2313
2314         Returns the empty string if not present.
2315         '''
2316         def __get__(self):
2317             c = self.cigar
2318             if c == None: return ""
2319             else: return "".join([ "%c%i" % (CODE2CIGAR[x],y) for x,y in c])
2320
2321         def __set__(self, cigar):
2322             if cigar == None or len(cigar) == 0: self.cigar = []
2323             parts = CIGAR_REGEX.findall( cigar )
2324             self.cigar = [ (CIGAR2CODE[ord(x)], int(y)) for x,y in parts ]
2325
2326     property seq:
2327         """read sequence bases, including :term:`soft clipped` bases (None if not present).
2328
2329         In Python 3, this property is of type bytes and assigning a unicode string to it consisting of ASCII characters only will work, but is inefficient."""
2330         def __get__(self):
2331             cdef bam1_t * src
2332             cdef char * s
2333             src = self._delegate
2334
2335             if src.core.l_qseq == 0: return None
2336
2337             return get_seq_range(src, 0, src.core.l_qseq)
2338
2339         def __set__(self,seq):
2340             # samtools manages sequence and quality length memory together
2341             # if no quality information is present, the first byte says 0xff.
2342
2343             if seq == None or len(seq) == 0: return
2344             seq = _force_bytes(seq)
2345             cdef bam1_t * src
2346             cdef uint8_t * p
2347             cdef char * s
2348             cdef int l, k, nbytes_new, nbytes_old
2349
2350             src = self._delegate
2351
2352             l = len(seq)
2353
2354             # as the sequence is stored in half-bytes, the total length (sequence
2355             # plus quality scores) is (l+1)/2 + l
2356             nbytes_new = (l+1)/2 + l
2357             nbytes_old = (src.core.l_qseq+1)/2 + src.core.l_qseq
2358             # acquire pointer to location in memory
2359             p = bam1_seq( src )
2360             src.core.l_qseq = l
2361
2362             pysam_bam_update( src,
2363                               nbytes_old,
2364                               nbytes_new,
2365                               p)
2366             # re-acquire pointer to location in memory
2367             # as it might have moved
2368             p = bam1_seq( src )
2369             for k from 0 <= k < nbytes_new: p[k] = 0
2370             # convert to C string
2371             s = seq
2372             for k from 0 <= k < l:
2373                 p[k/2] |= pysam_translate_sequence(s[k]) << 4 * (1 - k % 2)
2374
2375             # erase qualities
2376             p = bam1_qual( src )
2377             p[0] = 0xff
2378
2379
2380     property qual:
2381         """read sequence base qualities, including :term:`soft clipped` bases (None if not present).
2382
2383         In Python 3, this property is of type bytes and assigning a unicode string to it consisting of ASCII characters only will work, but is inefficient."""
2384         def __get__(self):
2385
2386             cdef bam1_t * src
2387             cdef char * q
2388
2389             src = self._delegate
2390
2391             if src.core.l_qseq == 0: return None
2392
2393             return get_qual_range(src, 0, src.core.l_qseq)
2394
2395         def __set__(self,qual):
2396             # note that space is already allocated via the sequences
2397             cdef bam1_t * src
2398             cdef uint8_t * p
2399             cdef char * q
2400             cdef int k
2401
2402             src = self._delegate
2403             p = bam1_qual( src )
2404             if qual == None or len(qual) == 0:
2405                 # if absent - set to 0xff
2406                 p[0] = 0xff
2407                 return
2408             qual = _force_bytes(qual)
2409             cdef int l
2410             # convert to C string
2411             q = qual
2412             l = len(qual)
2413             if src.core.l_qseq != l:
2414                 raise ValueError("quality and sequence mismatch: %i != %i" % (l, src.core.l_qseq))
2415             assert src.core.l_qseq == l
2416             for k from 0 <= k < l:
2417                 p[k] = <uint8_t>q[k] - 33
2418
2419     property query:
2420         """aligned portion of the read and excludes any flanking bases that were :term:`soft clipped` (None if not present).
2421
2422         In Python 3, this property is of type bytes. Assigning a unicode string to it consisting of ASCII characters only will work, but is inefficient.
2423
2424         SAM/BAM files may included extra flanking bases sequences that were
2425         not part of the alignment.  These bases may be the result of the
2426         Smith-Waterman or other algorithms, which may not require alignments
2427         that begin at the first residue or end at the last.  In addition,
2428         extra sequencing adapters, multiplex identifiers, and low-quality bases that
2429         were not considered for alignment may have been retained."""
2430
2431         def __get__(self):
2432             cdef bam1_t * src
2433             cdef uint32_t start, end
2434             cdef char * s
2435
2436             src = self._delegate
2437
2438             if src.core.l_qseq == 0: return None
2439
2440             start = query_start(src)
2441             end   = query_end(src)
2442
2443             return get_seq_range(src, start, end)
2444
2445     property qqual:
2446         """aligned query sequence quality values (None if not present). This property is read-only.
2447
2448         In Python 3, this property is of type bytes."""
2449         def __get__(self):
2450             cdef bam1_t * src
2451             cdef uint32_t start, end
2452
2453             src = self._delegate
2454
2455             if src.core.l_qseq == 0: return None
2456
2457             start = query_start(src)
2458             end   = query_end(src)
2459
2460             return get_qual_range(src, start, end)
2461
2462     property qstart:
2463         """start index of the aligned query portion of the sequence (0-based, inclusive)"""
2464         def __get__(self):
2465             return query_start(self._delegate)
2466
2467     property qend:
2468         """end index of the aligned query portion of the sequence (0-based, exclusive)"""
2469         def __get__(self):
2470             return query_end(self._delegate)
2471
2472     property qlen:
2473         """Length of the aligned query sequence"""
2474         def __get__(self):
2475             cdef bam1_t * src
2476             src = self._delegate
2477             return query_end(src)-query_start(src)
2478
2479     property tags:
2480         """the tags in the AUX field.
2481
2482         This property permits convenience access to
2483         the tags. Changes it the returned list will
2484         not update the tags automatically. Instead,
2485         the following is required for adding a
2486         new tag::
2487
2488             read.tags = read.tags + [("RG",0)]
2489
2490
2491         This method will happily write the same tag
2492         multiple times.
2493         """
2494         def __get__(self):
2495             cdef char * ctag
2496             cdef bam1_t * src
2497             cdef uint8_t * s
2498             cdef char auxtag[3]
2499             cdef char auxtype
2500             cdef uint8_t byte_size
2501             cdef int32_t nvalues
2502
2503             src = self._delegate
2504             if src.l_aux == 0: return []
2505             s = bam1_aux( src )
2506             result = []
2507             auxtag[2] = 0
2508             while s < (src.data + src.data_len):
2509                 # get tag
2510                 auxtag[0] = s[0]
2511                 auxtag[1] = s[1]
2512                 s += 2
2513                 auxtype = s[0]
2514                 if auxtype in ('c', 'C'):
2515                     value = <int>bam_aux2i(s)
2516                     s += 1
2517                 elif auxtype in ('s', 'S'):
2518                     value = <int>bam_aux2i(s)
2519                     s += 2
2520                 elif auxtype in ('i', 'I'):
2521                     value = <int32_t>bam_aux2i(s)
2522                     s += 4
2523                 elif auxtype == 'f':
2524                     value = <float>bam_aux2f(s)
2525                     s += 4
2526                 elif auxtype == 'd':
2527                     value = <double>bam_aux2d(s)
2528                     s += 8
2529                 elif auxtype == 'A':
2530                     value = "%c" % <char>bam_aux2A(s)
2531                     s += 1
2532                 elif auxtype in ('Z', 'H'):
2533                     value = _charptr_to_str(<char*>bam_aux2Z(s))
2534                     # +1 for NULL terminated string
2535                     s += len(value) + 1
2536                 elif auxtype == 'B':
2537                     s += 1
2538                     byte_size, nvalues, value = convertBinaryTagToList( s )
2539                     # 5 for 1 char and 1 int
2540                     s += 5 + ( nvalues * byte_size) - 1
2541
2542                 s += 1
2543
2544                 result.append( (_charptr_to_str(auxtag), value) )
2545
2546             return result
2547
2548         def __set__(self, tags):
2549             cdef bam1_t * src
2550             cdef uint8_t * s
2551             cdef uint8_t * new_data
2552             cdef char * temp
2553
2554             src = self._delegate
2555
2556             fmts, args = ["<"], []
2557
2558             if tags != None:
2559
2560                 # map samtools code to python.struct code and byte size
2561                 for pytag, value in tags:
2562                     if not type(pytag) is bytes:
2563                         pytag = pytag.encode('ascii')
2564                     t = type(value)
2565
2566                     if t is tuple or t is list:
2567                         # binary tags - treat separately
2568                         pytype = 'B'
2569                         # get data type - first value determines type
2570                         if type(value[0]) is float:
2571                             datafmt, datatype = "f", "f"
2572                         else:
2573                             mi, ma = min(value), max(value)
2574                             absmax = max( abs(mi), abs(ma) )
2575                             # signed ints
2576                             if mi < 0:
2577                                 if mi >= -127: datafmt, datatype = "b", 'c'
2578                                 elif mi >= -32767: datafmt, datatype = "h", 's'
2579                                 elif absmax < -2147483648: raise ValueError( "integer %i out of range of BAM/SAM specification" % value )
2580                                 else: datafmt, datatype = "i", 'i'
2581
2582                             # unsigned ints
2583                             else:
2584                                 if absmax <= 255: datafmt, datatype = "B", 'C'
2585                                 elif absmax <= 65535: datafmt, datatype = "H", 'S'
2586                                 elif absmax > 4294967295: raise ValueError( "integer %i out of range of BAM/SAM specification" % value )
2587                                 else: datafmt, datatype = "I", 'I'
2588
2589                         datafmt = "2sccI%i%s" % (len(value), datafmt)
2590                         args.extend( [pytag[:2],
2591                                       pytype.encode('ascii'),
2592                                       datatype.encode('ascii'),
2593                                       len(value)] + list(value) )
2594                         fmts.append( datafmt )
2595                         continue
2596
2597                     if t is float:
2598                         fmt, pytype = "2scf", 'f'
2599                     elif t is int:
2600                         # negative values
2601                         if value < 0:
2602                             if value >= -127: fmt, pytype = "2scb", 'c'
2603                             elif value >= -32767: fmt, pytype = "2sch", 's'
2604                             elif value < -2147483648: raise ValueError( "integer %i out of range of BAM/SAM specification" % value )
2605                             else: fmt, pytype = "2sci", 'i'
2606                         # positive values
2607                         else:
2608                             if value <= 255: fmt, pytype = "2scB", 'C'
2609                             elif value <= 65535: fmt, pytype = "2scH", 'S'
2610                             elif value > 4294967295: raise ValueError( "integer %i out of range of BAM/SAM specification" % value )
2611                             else: fmt, pytype = "2scI", 'I'
2612                     else:
2613                         # Note: hex strings (H) are not supported yet
2614                         if t is not bytes:
2615                             value = value.encode('ascii')
2616                         if len(value) == 1:
2617                             fmt, pytype = "2scc", 'A'
2618                         else:
2619                             fmt, pytype = "2sc%is" % (len(value)+1), 'Z'
2620
2621                     args.extend( [pytag[:2],
2622                                   pytype.encode('ascii'),
2623                                   value ] )
2624
2625                     fmts.append( fmt )
2626
2627                 fmt = "".join(fmts)
2628                 total_size = struct.calcsize(fmt)
2629                 buffer = ctypes.create_string_buffer(total_size)
2630                 struct.pack_into( fmt,
2631                                   buffer,
2632                                   0,
2633                                   *args )
2634
2635             # delete the old data and allocate new space.
2636             # If total_size == 0, the aux field will be
2637             # empty
2638             pysam_bam_update( src,
2639                               src.l_aux,
2640                               total_size,
2641                               bam1_aux( src ) )
2642
2643             src.l_aux = total_size
2644
2645             # copy data only if there is any
2646             if total_size != 0:
2647
2648                 # get location of new data
2649                 s = bam1_aux( src )
2650
2651                 # check if there is direct path from buffer.raw to tmp
2652                 temp = buffer.raw
2653                 memcpy( s, temp, total_size )
2654
2655     property flag:
2656         """properties flag"""
2657         def __get__(self): return self._delegate.core.flag
2658         def __set__(self, flag): self._delegate.core.flag = flag
2659
2660     property rname:
2661         """
2662         :term:`target` ID
2663
2664         DEPRECATED from pysam-0.4 - use tid in the future.
2665         The rname field caused a lot of confusion as it returns
2666         the :term:`target` ID instead of the reference sequence
2667         name.
2668
2669         .. note::
2670
2671             This field contains the index of the reference sequence
2672             in the sequence dictionary. To obtain the name
2673             of the reference sequence, use :meth:`pysam.Samfile.getrname()`
2674
2675         """
2676         def __get__(self): return self._delegate.core.tid
2677         def __set__(self, tid): self._delegate.core.tid = tid
2678
2679     property tid:
2680         """
2681         :term:`target` ID
2682
2683         .. note::
2684
2685             This field contains the index of the reference sequence
2686             in the sequence dictionary. To obtain the name
2687             of the reference sequence, use :meth:`pysam.Samfile.getrname()`
2688
2689         """
2690         def __get__(self): return self._delegate.core.tid
2691         def __set__(self, tid): self._delegate.core.tid = tid
2692
2693     property pos:
2694         """0-based leftmost coordinate"""
2695         def __get__(self): return self._delegate.core.pos
2696         def __set__(self, pos):
2697             ## setting the cigar string also updates the "bin" attribute
2698             cdef bam1_t * src
2699             src = self._delegate
2700             if src.core.n_cigar:
2701                 src.core.bin = bam_reg2bin( src.core.pos, bam_calend( &src.core, bam1_cigar(src)) )
2702             else:
2703                 src.core.bin = bam_reg2bin( src.core.pos, src.core.pos + 1)
2704             self._delegate.core.pos = pos
2705     property bin:
2706         """properties bin"""
2707         def __get__(self): return self._delegate.core.bin
2708         def __set__(self, bin): self._delegate.core.bin = bin
2709     property rlen:
2710         '''length of the read (read only). Returns 0 if not given.'''
2711         def __get__(self): return self._delegate.core.l_qseq
2712     property aend:
2713         '''aligned reference position of the read on the reference genome.
2714
2715         aend points to one past the last aligned residue.
2716         Returns None if not available.'''
2717         def __get__(self):
2718             cdef bam1_t * src
2719             src = self._delegate
2720             if (self.flag & BAM_FUNMAP) or src.core.n_cigar == 0:
2721                 return None
2722             return bam_calend(&src.core, bam1_cigar(src))
2723
2724     property alen:
2725         '''aligned length of the read on the reference genome.  Returns None if
2726         not available.'''
2727         def __get__(self):
2728             cdef bam1_t * src
2729             src = self._delegate
2730             if (self.flag & BAM_FUNMAP) or src.core.n_cigar == 0:
2731                 return None
2732             return bam_calend(&src.core,
2733                                bam1_cigar(src)) - \
2734                                self._delegate.core.pos
2735
2736     property mapq:
2737         """mapping quality"""
2738         def __get__(self): return self._delegate.core.qual
2739         def __set__(self, qual): self._delegate.core.qual = qual
2740
2741     property mrnm:
2742         """the :term:`reference` id of the mate
2743         deprecated, use RNEXT instead.
2744         """
2745         def __get__(self): return self._delegate.core.mtid
2746         def __set__(self, mtid): self._delegate.core.mtid = mtid
2747     property rnext:
2748         """the :term:`reference` id of the mate """
2749         def __get__(self): return self._delegate.core.mtid
2750         def __set__(self, mtid): self._delegate.core.mtid = mtid
2751     property mpos:
2752         """the position of the mate
2753         deprecated, use PNEXT instead."""
2754         def __get__(self): return self._delegate.core.mpos
2755         def __set__(self, mpos): self._delegate.core.mpos = mpos
2756     property pnext:
2757         """the position of the mate"""
2758         def __get__(self): return self._delegate.core.mpos
2759         def __set__(self, mpos): self._delegate.core.mpos = mpos
2760     property isize:
2761         """the insert size
2762         deprecated: use tlen instead"""
2763         def __get__(self): return self._delegate.core.isize
2764         def __set__(self, isize): self._delegate.core.isize = isize
2765     property tlen:
2766         """the insert size"""
2767         def __get__(self): return self._delegate.core.isize
2768         def __set__(self, isize): self._delegate.core.isize = isize
2769     property is_paired:
2770         """true if read is paired in sequencing"""
2771         def __get__(self): return (self._delegate.core.flag & BAM_FPAIRED) != 0
2772         def __set__(self,val):
2773             if val: self._delegate.core.flag |= BAM_FPAIRED
2774             else: self._delegate.core.flag &= ~BAM_FPAIRED
2775     property is_proper_pair:
2776         """true if read is mapped in a proper pair"""
2777         def __get__(self): return (self.flag & BAM_FPROPER_PAIR) != 0
2778         def __set__(self,val):
2779             if val: self._delegate.core.flag |= BAM_FPROPER_PAIR
2780             else: self._delegate.core.flag &= ~BAM_FPROPER_PAIR
2781     property is_unmapped:
2782         """true if read itself is unmapped"""
2783         def __get__(self): return (self.flag & BAM_FUNMAP) != 0
2784         def __set__(self,val):
2785             if val: self._delegate.core.flag |= BAM_FUNMAP
2786             else: self._delegate.core.flag &= ~BAM_FUNMAP
2787     property mate_is_unmapped:
2788         """true if the mate is unmapped"""
2789         def __get__(self): return (self.flag & BAM_FMUNMAP) != 0
2790         def __set__(self,val):
2791             if val: self._delegate.core.flag |= BAM_FMUNMAP
2792             else: self._delegate.core.flag &= ~BAM_FMUNMAP
2793     property is_reverse:
2794         """true if read is mapped to reverse strand"""
2795         def __get__(self): return (self.flag & BAM_FREVERSE) != 0
2796         def __set__(self,val):
2797             if val: self._delegate.core.flag |= BAM_FREVERSE
2798             else: self._delegate.core.flag &= ~BAM_FREVERSE
2799     property mate_is_reverse:
2800         """true is read is mapped to reverse strand"""
2801         def __get__(self): return (self.flag & BAM_FMREVERSE) != 0
2802         def __set__(self,val):
2803             if val: self._delegate.core.flag |= BAM_FMREVERSE
2804             else: self._delegate.core.flag &= ~BAM_FMREVERSE
2805     property is_read1:
2806         """true if this is read1"""
2807         def __get__(self): return (self.flag & BAM_FREAD1) != 0
2808         def __set__(self,val):
2809             if val: self._delegate.core.flag |= BAM_FREAD1
2810             else: self._delegate.core.flag &= ~BAM_FREAD1
2811     property is_read2:
2812         """true if this is read2"""
2813         def __get__(self): return (self.flag & BAM_FREAD2) != 0
2814         def __set__(self,val):
2815             if val: self._delegate.core.flag |= BAM_FREAD2
2816             else: self._delegate.core.flag &= ~BAM_FREAD2
2817     property is_secondary:
2818         """true if not primary alignment"""
2819         def __get__(self): return (self.flag & BAM_FSECONDARY) != 0
2820         def __set__(self,val):
2821             if val: self._delegate.core.flag |= BAM_FSECONDARY
2822             else: self._delegate.core.flag &= ~BAM_FSECONDARY
2823     property is_qcfail:
2824         """true if QC failure"""
2825         def __get__(self): return (self.flag & BAM_FQCFAIL) != 0
2826         def __set__(self,val):
2827             if val: self._delegate.core.flag |= BAM_FQCFAIL
2828             else: self._delegate.core.flag &= ~BAM_FQCFAIL
2829     property is_duplicate:
2830         """true if optical or PCR duplicate"""
2831         def __get__(self): return (self.flag & BAM_FDUP) != 0
2832         def __set__(self,val):
2833             if val: self._delegate.core.flag |= BAM_FDUP
2834             else: self._delegate.core.flag &= ~BAM_FDUP
2835     property positions:
2836         """a list of reference positions that this read aligns to."""
2837         def __get__(self):
2838             cdef uint32_t k, i, pos
2839             cdef int op
2840             cdef uint32_t * cigar_p
2841             cdef bam1_t * src
2842
2843             src = self._delegate
2844             if src.core.n_cigar == 0: return []
2845
2846             result = []
2847             pos = src.core.pos
2848             cigar_p = bam1_cigar(src)
2849
2850             for k from 0 <= k < src.core.n_cigar:
2851                 op = cigar_p[k] & BAM_CIGAR_MASK
2852                 l = cigar_p[k] >> BAM_CIGAR_SHIFT
2853                 if op == BAM_CMATCH:
2854                     for i from pos <= i < pos + l:
2855                         result.append( i )
2856
2857                 if op == BAM_CMATCH or op == BAM_CDEL or op == BAM_CREF_SKIP:
2858                     pos += l
2859
2860             return result
2861
2862     property aligned_pairs:
2863        """a list of aligned read and reference positions.
2864
2865        Unaligned position are marked by None.
2866        """
2867        def __get__(self):
2868            cdef uint32_t k, i, pos, qpos
2869            cdef int op
2870            cdef uint32_t * cigar_p
2871            cdef bam1_t * src
2872
2873            src = self._delegate
2874            if src.core.n_cigar == 0: return []
2875
2876            result = []
2877            pos = src.core.pos
2878            qpos = 0
2879            cigar_p = bam1_cigar(src)
2880
2881            for k from 0 <= k < src.core.n_cigar:
2882                op = cigar_p[k] & BAM_CIGAR_MASK
2883                l = cigar_p[k] >> BAM_CIGAR_SHIFT
2884
2885                if op == BAM_CMATCH:
2886                    for i from pos <= i < pos + l:
2887                        result.append( (qpos, i) )
2888                        qpos += 1
2889                    pos += l
2890
2891                elif op == BAM_CINS:
2892                    for i from pos <= i < pos + l:
2893                        result.append( (qpos, None) )
2894                        qpos += 1
2895
2896                elif op == BAM_CDEL or op == BAM_CREF_SKIP:
2897                    for i from pos <= i < pos + l:
2898                        result.append( (None, i) )
2899                    pos += l
2900
2901            return result
2902
2903
2904     def overlap( self, uint32_t start, uint32_t end ):
2905         """return number of aligned bases of read overlapping the interval *start* and *end*
2906         on the reference sequence.
2907         """
2908         cdef uint32_t k, i, pos, overlap
2909         cdef int op, o
2910         cdef uint32_t * cigar_p
2911         cdef bam1_t * src
2912
2913         overlap = 0
2914
2915         src = self._delegate
2916         if src.core.n_cigar == 0: return 0
2917         pos = src.core.pos
2918         o = 0
2919
2920         cigar_p = bam1_cigar(src)
2921         for k from 0 <= k < src.core.n_cigar:
2922             op = cigar_p[k] & BAM_CIGAR_MASK
2923             l = cigar_p[k] >> BAM_CIGAR_SHIFT
2924
2925             if op == BAM_CMATCH:
2926                 o = min( pos + l, end) - max( pos, start )
2927                 if o > 0: overlap += o
2928
2929             if op == BAM_CMATCH or op == BAM_CDEL or op == BAM_CREF_SKIP:
2930                 pos += l
2931
2932         return overlap
2933
2934     def opt(self, tag):
2935         """retrieves optional data given a two-letter *tag*"""
2936         #see bam_aux.c: bam_aux_get() and bam_aux2i() etc
2937         cdef uint8_t * v
2938         cdef int nvalues
2939         btag = _force_bytes(tag)
2940         v = bam_aux_get(self._delegate, btag)
2941         if v == NULL: raise KeyError( "tag '%s' not present" % tag )
2942         auxtype = chr(v[0])
2943         if auxtype == 'c' or auxtype == 'C' or auxtype == 's' or auxtype == 'S':
2944             return <int>bam_aux2i(v)
2945         elif auxtype == 'i' or auxtype == 'I':
2946             return <int32_t>bam_aux2i(v)
2947         elif auxtype == 'f' or auxtype == 'F':
2948             return <float>bam_aux2f(v)
2949         elif auxtype == 'd' or auxtype == 'D':
2950             return <double>bam_aux2d(v)
2951         elif auxtype == 'A':
2952             # there might a more efficient way
2953             # to convert a char into a string
2954             return '%c' % <char>bam_aux2A(v)
2955         elif auxtype == 'Z':
2956             return _charptr_to_str(<char*>bam_aux2Z(v))
2957         elif auxtype == 'B':
2958             bytesize, nvalues, values = convertBinaryTagToList( v + 1 )
2959             return values
2960         else:
2961             raise ValueError("unknown auxilliary type '%s'" % auxtype)
2962
2963
2964     def fancy_str (self):
2965         """returns list of fieldnames/values in pretty format for debugging
2966         """
2967         ret_string = []
2968         field_names = {
2969            "tid":           "Contig index",
2970            "pos":           "Mapped position on contig",
2971            "mtid":          "Contig index for mate pair",
2972            "mpos":          "Position of mate pair",
2973            "isize":         "Insert size",
2974            "flag":          "Binary flag",
2975            "n_cigar":       "Count of cigar entries",
2976            "cigar":         "Cigar entries",
2977            "qual":          "Mapping quality",
2978            "bin":           "Bam index bin number",
2979            "l_qname":       "Length of query name",
2980            "qname":         "Query name",
2981            "l_qseq":        "Length of query sequence",
2982            "qseq":          "Query sequence",
2983            "bqual":         "Quality scores",
2984            "l_aux":         "Length of auxilary data",
2985            "m_data":        "Maximum data length",
2986            "data_len":      "Current data length",
2987            }
2988         fields_names_in_order = ["tid", "pos", "mtid", "mpos", "isize", "flag",
2989                                  "n_cigar", "cigar", "qual", "bin", "l_qname", "qname",
2990                                  "l_qseq", "qseq", "bqual", "l_aux", "m_data", "data_len"]
2991
2992         for f in fields_names_in_order:
2993             if not f in self.__dict__:
2994                 continue
2995             ret_string.append("%-30s %-10s= %s" % (field_names[f], "(" + f + ")", self.__getattribute__(f)))
2996
2997         for f in self.__dict__:
2998             if not f in field_names:
2999                 ret_string.append("%-30s %-10s= %s" % (f, "", self.__getattribute__(f)))
3000         return ret_string
3001
3002 cdef class PileupProxy:
3003     '''A pileup column. A pileup column contains
3004     all the reads that map to a certain target base.
3005
3006     tid
3007         chromosome ID as is defined in the header
3008     pos
3009         the target base coordinate (0-based)
3010     n
3011         number of reads mapping to this column
3012     pileups
3013         list of reads (:class:`pysam.PileupRead`) aligned to this column
3014
3015     This class is a proxy for results returned by the samtools pileup engine.
3016     If the underlying engine iterator advances, the results of this column
3017     will change.
3018     '''
3019     def __init__(self):
3020         raise TypeError("This class cannot be instantiated from Python")
3021
3022     def __str__(self):
3023         return "\t".join( map(str, (self.tid, self.pos, self.n))) +\
3024             "\n" +\
3025             "\n".join( map(str, self.pileups) )
3026
3027     property tid:
3028         '''the chromosome ID as is defined in the header'''
3029         def __get__(self): return self.tid
3030
3031     property n:
3032         '''number of reads mapping to this column.'''
3033         def __get__(self): return self.n_pu
3034         def __set__(self, n): self.n_pu = n
3035
3036     property pos:
3037         def __get__(self): return self.pos
3038
3039     property pileups:
3040         '''list of reads (:class:`pysam.PileupRead`) aligned to this column'''
3041         def __get__(self):
3042             cdef int x
3043             pileups = []
3044
3045             if self.plp[0] == NULL:
3046                 raise ValueError("PileupProxy accessed after iterator finished")
3047
3048             # warning: there could be problems if self.n and self.buf are
3049             # out of sync.
3050             for x from 0 <= x < self.n_pu:
3051                 pileups.append( makePileupRead( &(self.plp[0][x])) )
3052             return pileups
3053
3054 cdef class PileupRead:
3055     '''A read aligned to a column.
3056     '''
3057
3058     def __init__(self):
3059         raise TypeError("This class cannot be instantiated from Python")
3060
3061     def __str__(self):
3062         return "\t".join( map(str, (self.alignment, self.qpos, self.indel, self.level, self.is_del, self.is_head, self.is_tail ) ) )
3063
3064     property alignment:
3065         """a :class:`pysam.AlignedRead` object of the aligned read"""
3066         def __get__(self):
3067             return self._alignment
3068     property qpos:
3069         """position of the read base at the pileup site, 0-based"""
3070         def __get__(self):
3071             return self._qpos
3072     property indel:
3073         """indel length; 0 for no indel, positive for ins and negative for del"""
3074         def __get__(self):
3075             return self._indel
3076     property is_del:
3077         """1 iff the base on the padded read is a deletion"""
3078         def __get__(self):
3079             return self._is_del
3080     property is_head:
3081         def __get__(self):
3082             return self._is_head
3083     property is_tail:
3084         def __get__(self):
3085             return self._is_tail
3086     property level:
3087         def __get__(self):
3088             return self._level
3089
3090 class Outs:
3091     '''http://mail.python.org/pipermail/python-list/2000-June/038406.html'''
3092     def __init__(self, id = 1):
3093         self.streams = []
3094         self.id = id
3095
3096     def setdevice(self, filename):
3097         '''open an existing file, like "/dev/null"'''
3098         fd = os.open(filename, os.O_WRONLY)
3099         self.setfd(fd)
3100
3101     def setfile(self, filename):
3102         '''open a new file.'''
3103         fd = os.open(filename, os.O_WRONLY|os.O_CREAT, 0660);
3104         self.setfd(fd)
3105
3106     def setfd(self, fd):
3107         ofd = os.dup(self.id)      #  Save old stream on new unit.
3108         self.streams.append(ofd)
3109         sys.stdout.flush()          #  Buffered data goes to old stream.
3110         sys.stderr.flush()          #  Buffered data goes to old stream.
3111         os.dup2(fd, self.id)        #  Open unit 1 on new stream.
3112         os.close(fd)                #  Close other unit (look out, caller.)
3113
3114     def restore(self):
3115         '''restore previous output stream'''
3116         if self.streams:
3117             # the following was not sufficient, hence flush both stderr and stdout
3118             # os.fsync( self.id )
3119             sys.stdout.flush()
3120             sys.stderr.flush()
3121             os.dup2(self.streams[-1], self.id)
3122             os.close(self.streams[-1])
3123             del self.streams[-1]
3124
3125 def _samtools_dispatch( method,
3126                         args = (),
3127                         catch_stdout = True ):
3128     '''call ``method`` in samtools providing arguments in args.
3129
3130     .. note::
3131        This method redirects stdout to capture it
3132        from samtools. If for some reason stdout disappears
3133        the reason might be in this method.
3134
3135     .. note::
3136        The current implementation might only work on linux.
3137
3138     .. note::
3139        This method captures stdout and stderr using temporary files,
3140        which are then read into memory in their entirety. This method
3141        is slow and might cause large memory overhead.
3142
3143     See http://bytes.com/topic/c/answers/487231-how-capture-stdout-temporarily
3144     on the topic of redirecting stderr/stdout.
3145     '''
3146
3147     # note that debugging this module can be a problem
3148     # as stdout/stderr will not appear on the terminal
3149
3150     # some special cases
3151     if method == "index":
3152         if not os.path.exists( args[0] ):
3153             raise IOError( "No such file or directory: '%s'" % args[0] )
3154
3155     # redirect stderr and stdout to file
3156     stderr_h, stderr_f = tempfile.mkstemp()
3157     pysam_set_stderr( stderr_h )
3158
3159     if catch_stdout:
3160         stdout_h, stdout_f = tempfile.mkstemp()
3161         try:
3162             stdout_save = Outs( sys.stdout.fileno() )
3163             stdout_save.setfd( stdout_h )
3164         except AttributeError:
3165             # stdout has already been redirected
3166             catch_stdout = False
3167
3168         # patch for `samtools view`
3169         # samtools `view` closes stdout, from which I can not
3170         # recover. Thus redirect output to file with -o option.
3171         if method == "view":
3172             if "-o" in args: raise ValueError("option -o is forbidden in samtools view")
3173             args = ( "-o", stdout_f ) + args
3174
3175     # do the function call to samtools
3176     cdef char ** cargs
3177     cdef int i, n, retval
3178
3179     n = len(args)
3180     method = _force_cmdline_bytes(method)
3181     args = [ _force_cmdline_bytes(a) for a in args ]
3182
3183     # allocate two more for first (dummy) argument (contains command)
3184     cargs = <char**>calloc( n+2, sizeof( char *) )
3185     cargs[0] = "samtools"
3186     cargs[1] = method
3187     for i from 0 <= i < n: cargs[i+2] = args[i]
3188
3189     retval = pysam_dispatch(n+2, cargs)
3190     free( cargs )
3191
3192     # restore stdout/stderr. This will also flush, so
3193     # needs to be before reading back the file contents
3194     if catch_stdout:
3195         stdout_save.restore()
3196         try:
3197             with open( stdout_f, "r") as inf:
3198                 out_stdout = inf.readlines()
3199         except UnicodeDecodeError:
3200             with open( stdout_f, "rb") as inf:
3201                 # read binary output
3202                 out_stdout = inf.read()
3203         os.remove( stdout_f )
3204     else:
3205         out_stdout = []
3206
3207     # get error messages
3208     pysam_unset_stderr()
3209     try:
3210         with open( stderr_f, "r") as inf:
3211             out_stderr = inf.readlines()
3212     except UnicodeDecodeError:
3213         with open( stderr_f, "rb") as inf:
3214             # read binary output
3215             out_stderr = inf.read()
3216         os.remove( stderr_f )
3217     else:
3218         out_stderr = []
3219
3220     return retval, out_stderr, out_stdout
3221
3222 cdef class SNPCall:
3223     '''the results of a SNP call.'''
3224     cdef int _tid
3225     cdef int _pos
3226     cdef char _reference_base
3227     cdef char _genotype
3228     cdef int _consensus_quality
3229     cdef int _snp_quality
3230     cdef int _rms_mapping_quality
3231     cdef int _coverage
3232
3233     property tid:
3234         '''the chromosome ID as is defined in the header'''
3235         def __get__(self):
3236             return self._tid
3237
3238     property pos:
3239        '''nucleotide position of SNP.'''
3240        def __get__(self): return self._pos
3241
3242     property reference_base:
3243        '''reference base at pos. ``N`` if no reference sequence supplied.'''
3244        def __get__(self): return from_string_and_size( &self._reference_base, 1 )
3245
3246     property genotype:
3247        '''the genotype called.'''
3248        def __get__(self): return from_string_and_size( &self._genotype, 1 )
3249
3250     property consensus_quality:
3251        '''the genotype quality (Phred-scaled).'''
3252        def __get__(self): return self._consensus_quality
3253
3254     property snp_quality:
3255        '''the snp quality (Phred scaled) - probability of consensus being identical to reference sequence.'''
3256        def __get__(self): return self._snp_quality
3257
3258     property mapping_quality:
3259        '''the root mean square (rms) of the mapping quality of all reads involved in the call.'''
3260        def __get__(self): return self._rms_mapping_quality
3261
3262     property coverage:
3263        '''coverage or read depth - the number of reads involved in the call.'''
3264        def __get__(self): return self._coverage
3265
3266     def __str__(self):
3267
3268         return "\t".join( map(str, (
3269                     self.tid,
3270                     self.pos,
3271                     self.reference_base,
3272                     self.genotype,
3273                     self.consensus_quality,
3274                     self.snp_quality,
3275                     self.mapping_quality,
3276                     self.coverage ) ) )
3277
3278
3279 # cdef class SNPCallerBase:
3280 #     '''Base class for SNP callers.
3281
3282 #     *min_baseQ*
3283 #        minimum base quality (possibly capped by BAQ)
3284 #     *capQ_threshold*
3285 #        coefficient for adjusting mapQ of poor mappings
3286 #     *theta*
3287 #        theta in maq consensus calling model
3288 #     *n_haplotypes*
3289 #        number of haplotypes in the sample
3290 #     *het_rate*
3291 #        prior of a difference between two haplotypes
3292 #     '''
3293
3294 #     cdef bam_maqcns_t * c
3295 #     cdef IteratorColumn iter
3296
3297 #     def __cinit__(self,
3298 #                   IteratorColumn iterator_column,
3299 #                   **kwargs ):
3300
3301 #         self.iter = iterator_column
3302 #         self.c =  bam_maqcns_init()
3303
3304 #         # set the default parameterization according to
3305 #         # samtools
3306
3307 #         # new default mode for samtools >0.1.10
3308 #         self.c.errmod = kwargs.get( "errmod", BAM_ERRMOD_MAQ2 )
3309
3310 #         self.c.min_baseQ = kwargs.get( "min_baseQ", 13 )
3311 #         # self.c.capQ_thres = kwargs.get( "capQ_threshold", 60 )
3312 #         self.c.n_hap = kwargs.get( "n_haplotypes", 2 )
3313 #         self.c.het_rate = kwargs.get( "het_rate", 0.001 )
3314 #         self.c.theta = kwargs.get( "theta", 0.83 )
3315
3316 #         if self.c.errmod != BAM_ERRMOD_MAQ2:
3317 #             self.c.theta += 0.02
3318
3319 #         # call prepare AFTER setting parameters
3320 #         bam_maqcns_prepare( self.c )
3321
3322 #     def __dealloc__(self):
3323 #         bam_maqcns_destroy( self.c )
3324
3325     # cdef __dump( self, glf1_t * g, uint32_t cns, int rb ):
3326     #     '''debugging output.'''
3327
3328     #     pysam_dump_glf( g, self.c );
3329     #     print ""
3330     #     for x in range(self.iter.n_plp):
3331     #         print "--> read %i %s %i" % (x,
3332     #                                      bam1_qname(self.iter.plp[x].b),
3333     #                                      self.iter.plp[x].qpos,
3334     #                                      )
3335
3336     #     print "pos=%i, cns=%i, q_r = %f, depth=%i, n=%i, rb=%i, cns-cq=%i %i %i %i" \
3337     #         % (self.iter.pos,
3338     #            cns,
3339     #            self.c.q_r,
3340     #            self.iter.n_plp,
3341     #            self.iter.n_plp,
3342     #            rb,
3343     #            cns >> 8 & 0xff,
3344     #            cns >> 16 & 0xff,
3345     #            cns & 0xff,
3346     #            cns >> 28,
3347     #            )
3348
3349     #     printf("-------------------------------------\n");
3350     #     sys.stdout.flush()
3351
3352 # cdef class IteratorSNPCalls( SNPCallerBase ):
3353 #     """*(IteratorColumn iterator)*
3354
3355 #     call SNPs within a region.
3356
3357 #     *iterator* is a pileup iterator. SNPs will be called
3358 #     on all positions returned by this iterator.
3359
3360 #     This caller is fast if SNPs are called over large continuous
3361 #     regions. It is slow, if instantiated frequently and in random
3362 #     order as the sequence will have to be reloaded.
3363
3364 #     """
3365
3366 #     def __cinit__(self,
3367 #                   IteratorColumn iterator_column,
3368 #                   **kwargs ):
3369
3370 #         assert self.iter.hasReference(), "IteratorSNPCalls requires an pileup iterator with reference sequence"
3371
3372 #     def __iter__(self):
3373 #         return self
3374
3375 #     def __next__(self):
3376 #         """python version of next().
3377 #         """
3378
3379 #         # the following code was adapted from bam_plcmd.c:pileup_func()
3380 #         self.iter.cnext()
3381
3382 #         if self.iter.n_plp < 0:
3383 #             raise ValueError("error during iteration" )
3384
3385 #         if self.iter.plp == NULL:
3386 #            raise StopIteration
3387
3388 #         cdef char * seq = self.iter.getSequence()
3389 #         cdef int seq_len = self.iter.seq_len
3390
3391 #         assert seq != NULL
3392
3393 #         # reference base
3394 #         if self.iter.pos >= seq_len:
3395 #             raise ValueError( "position %i out of bounds on reference sequence (len=%i)" % (self.iter.pos, seq_len) )
3396
3397 #         cdef int rb = seq[self.iter.pos]
3398 #         cdef uint32_t cns
3399 #        cdef glf1_t * g
3400
3401 #        g = bam_maqcns_glfgen( self.iter.n_plp,
3402 #                               self.iter.plp,
3403 #                               bam_nt16_table[rb],
3404 #                               self.c )
3405
3406 #        if pysam_glf_depth( g ) == 0:
3407 #            cns = 0xfu << 28 | 0xf << 24
3408 #        else:
3409 #            cns = glf2cns(g, <int>(self.c.q_r + .499))
3410
3411 #        free(g)
3412
3413 #         cdef SNPCall call
3414
3415 #         call = SNPCall()
3416 #         call._tid = self.iter.tid
3417 #         call._pos = self.iter.pos
3418 #         call._reference_base = rb
3419 #         call._genotype = bam_nt16_rev_table[cns>>28]
3420 #         call._consensus_quality = cns >> 8 & 0xff
3421 #         call._snp_quality = cns & 0xff
3422 #         call._rms_mapping_quality = cns >> 16&0xff
3423 #         call._coverage = self.iter.n_plp
3424
3425 #         return call
3426
3427 # cdef class SNPCaller( SNPCallerBase ):
3428 #     '''*(IteratorColumn iterator_column )*
3429
3430 #     The samtools SNP caller.
3431
3432 #     This object will call SNPs in *samfile* against the reference
3433 #     sequence in *fasta*.
3434
3435 #     This caller is fast for calling few SNPs in selected regions.
3436
3437 #     It is slow, if called over large genomic regions.
3438 #     '''
3439
3440
3441 #     def __cinit__(self,
3442 #                   IteratorColumn iterator_column,
3443 #                   **kwargs ):
3444
3445 #         pass
3446
3447 #     def call(self, reference, int pos ):
3448 #         """call a snp on chromosome *reference*
3449 #         and position *pos*.
3450
3451 #         returns a :class:`SNPCall` object.
3452 #         """
3453
3454 #         cdef int tid = self.iter.samfile.gettid( reference )
3455
3456 #         self.iter.reset( tid, pos, pos + 1 )
3457
3458 #         while 1:
3459 #             self.iter.cnext()
3460
3461 #             if self.iter.n_plp < 0:
3462 #                 raise ValueError("error during iteration" )
3463
3464 #             if self.iter.plp == NULL:
3465 #                 raise ValueError( "no reads in region - no call" )
3466
3467 #             if self.iter.pos == pos: break
3468
3469 #         cdef char * seq = self.iter.getSequence()
3470 #         cdef int seq_len = self.iter.seq_len
3471
3472 #         assert seq != NULL
3473
3474 #         # reference base
3475 #         if self.iter.pos >= seq_len:
3476 #             raise ValueError( "position %i out of bounds on reference sequence (len=%i)" % (self.iter.pos, seq_len) )
3477
3478 #         cdef int rb = seq[self.iter.pos]
3479 #         cdef uint32_t cns
3480 # #        cdef glf1_t * g
3481 # #
3482 # #        g = bam_maqcns_glfgen( self.iter.n_plp,
3483 # #                               self.iter.plp,
3484 # #                               bam_nt16_table[rb],
3485 # #                               self.c )
3486 # ##
3487 # #
3488 # #        if pysam_glf_depth( g ) == 0:
3489 # #            cns = 0xfu << 28 | 0xf << 24
3490 # #        else:
3491 # #            cns = glf2cns(g, <int>(self.c.q_r + .499))
3492 # #
3493 # #        free(g)
3494
3495 #         cdef SNPCall call
3496
3497 #         call = SNPCall()
3498 #         call._tid = self.iter.tid
3499 #         call._pos = self.iter.pos
3500 #         call._reference_base = rb
3501 #         call._genotype = bam_nt16_rev_table[cns>>28]
3502 #         call._consensus_quality = cns >> 8 & 0xff
3503 #         call._snp_quality = cns & 0xff
3504 #         call._rms_mapping_quality = cns >> 16&0xff
3505 #         call._coverage = self.iter.n_plp
3506
3507 #         return call
3508
3509 # cdef class IndelCall:
3510 #     '''the results of an indel call.'''
3511 #     cdef int _tid
3512 #     cdef int _pos
3513 #     cdef int _coverage
3514 #     cdef int _rms_mapping_quality
3515 #     cdef bam_maqindel_ret_t * _r
3516
3517 #     def __cinit__(self):
3518 #         #assert r != NULL
3519 #         #self._r = r
3520 #         pass
3521
3522 #     property tid:
3523 #         '''the chromosome ID as is defined in the header'''
3524 #         def __get__(self):
3525 #             return self._tid
3526
3527 #     property pos:
3528 #        '''nucleotide position of SNP.'''
3529 #        def __get__(self): return self._pos
3530
3531 #     property genotype:
3532 #        '''the genotype called.'''
3533 #        def __get__(self):
3534 #            if self._r.gt == 0:
3535 #                s = PyString_FromStringAndSize( self._r.s[0], self._r.indel1 + 1)
3536 #                return "%s/%s" % (s,s)
3537 #            elif self._r.gt == 1:
3538 #                s = PyString_FromStringAndSize( self._r.s[1], self._r.indel2 + 1)
3539 #                return "%s/%s" % (s,s)
3540 #            else:
3541 #                return "%s/%s" % (self.first_allele, self.second_allele )
3542
3543 #     property consensus_quality:
3544 #        '''the genotype quality (Phred-scaled).'''
3545 #        def __get__(self): return self._r.q_cns
3546
3547 #     property snp_quality:
3548 #        '''the snp quality (Phred scaled) - probability of consensus being identical to reference sequence.'''
3549 #        def __get__(self): return self._r.q_ref
3550
3551 #     property mapping_quality:
3552 #        '''the root mean square (rms) of the mapping quality of all reads involved in the call.'''
3553 #        def __get__(self): return self._rms_mapping_quality
3554
3555 #     property coverage:
3556 #        '''coverage or read depth - the number of reads involved in the call.'''
3557 #        def __get__(self): return self._coverage
3558
3559 #     property first_allele:
3560 #        '''sequence of first allele.'''
3561 #        def __get__(self): return PyString_FromStringAndSize( self._r.s[0], self._r.indel1 + 1)
3562
3563 #     property second_allele:
3564 #        '''sequence of second allele.'''
3565 #        def __get__(self): return PyString_FromStringAndSize( self._r.s[1], self._r.indel2 + 1)
3566
3567 #     property reads_first:
3568 #        '''reads supporting first allele.'''
3569 #        def __get__(self): return self._r.cnt1
3570
3571 #     property reads_second:
3572 #        '''reads supporting first allele.'''
3573 #        def __get__(self): return self._r.cnt2
3574
3575 #     property reads_diff:
3576 #        '''reads supporting first allele.'''
3577 #        def __get__(self): return self._r.cnt_anti
3578
3579 #     def __str__(self):
3580
3581 #         return "\t".join( map(str, (
3582 #                     self.tid,
3583 #                     self.pos,
3584 #                     self.genotype,
3585 #                     self.consensus_quality,
3586 #                     self.snp_quality,
3587 #                     self.mapping_quality,
3588 #                     self.coverage,
3589 #                     self.first_allele,
3590 #                     self.second_allele,
3591 #                     self.reads_first,
3592 #                     self.reads_second,
3593 #                     self.reads_diff ) ) )
3594
3595 #     def __dealloc__(self ):
3596 #         bam_maqindel_ret_destroy(self._r)
3597
3598 # cdef class IndelCallerBase:
3599 #     '''Base class for SNP callers.
3600
3601 #     *min_baseQ*
3602 #        minimum base quality (possibly capped by BAQ)
3603 #     *capQ_threshold*
3604 #        coefficient for adjusting mapQ of poor mappings
3605 #     *theta*
3606 #        theta in maq consensus calling model
3607 #     *n_haplotypes*
3608 #        number of haplotypes in the sample
3609 #     *het_rate*
3610 #        prior of a difference between two haplotypes
3611 #     '''
3612
3613 #     cdef bam_maqindel_opt_t * options
3614 #     cdef IteratorColumn iter
3615 #     cdef int cap_mapQ
3616 #     cdef int max_depth
3617
3618 #     def __cinit__(self,
3619 #                   IteratorColumn iterator_column,
3620 #                   **kwargs ):
3621
3622
3623 #         self.iter = iterator_column
3624
3625 #         assert iterator_column.hasReference(), "IndelCallerBase requires an pileup iterator with reference sequence"
3626
3627 #         self.options = bam_maqindel_opt_init()
3628
3629 #         # set the default parameterization according to
3630 #         # samtools
3631
3632 #         self.options.r_indel = kwargs.get( "r_indel", 0.00015 )
3633 #         self.options.q_indel = kwargs.get( "q_indel", 40 )
3634 #         self.cap_mapQ = kwargs.get( "cap_mapQ", 60 )
3635 #         self.max_depth = kwargs.get( "max_depth", 1024 )
3636
3637 #     def __dealloc__(self):
3638 #         free( self.options )
3639
3640 #     def _call( self ):
3641
3642 #         cdef char * seq = self.iter.getSequence()
3643 #         cdef int seq_len = self.iter.seq_len
3644
3645 #         assert seq != NULL
3646
3647 #         # reference base
3648 #         if self.iter.pos >= seq_len:
3649 #             raise ValueError( "position %i out of bounds on reference sequence (len=%i)" % (self.iter.pos, seq_len) )
3650
3651 #         cdef bam_maqindel_ret_t * r
3652
3653 #         cdef int m = min( self.max_depth, self.iter.n_plp )
3654
3655 #         # printf("pysam: m=%i, q_indel=%i, r_indel=%f, r_snp=%i, mm_penalty=%i, indel_err=%i, ambi_thres=%i\n",
3656 #         #        m, self.options.q_indel, self.options.r_indel, self.options.r_snp, self.options.mm_penalty,
3657 #         #        self.options.indel_err, self.options.ambi_thres );
3658
3659 #         r = bam_maqindel(m,
3660 #                          self.iter.pos,
3661 #                          self.options,
3662 #                          self.iter.plp,
3663 #                          seq,
3664 #                          0,
3665 #                          NULL)
3666
3667 #         if r == NULL: return None
3668
3669 #         cdef IndelCall call
3670 #         call = IndelCall()
3671 #         call._r = r
3672 #         call._tid = self.iter.tid
3673 #         call._pos = self.iter.pos
3674 #         call._coverage = self.iter.n_plp
3675
3676 #         cdef uint64_t rms_aux = 0
3677 #         cdef int i = 0
3678 #         cdef bam_pileup1_t * p
3679 #         cdef int tmp
3680
3681 #         for i from 0 <= i < self.iter.n_plp:
3682 #             p = self.iter.plp + i
3683 #             if p.b.core.qual < self.cap_mapQ:
3684 #                 tmp = p.b.core.qual
3685 #             else:
3686 #                 tmp = self.cap_mapQ
3687 #             rms_aux += tmp * tmp
3688
3689 #         call._rms_mapping_quality = <uint64_t>(sqrt(<double>rms_aux / self.iter.n_plp) + .499)
3690
3691 #         return call
3692
3693 # cdef class IndelCaller( IndelCallerBase ):
3694 #     '''*(IteratorColumn iterator_column )*
3695
3696 #     The samtools SNP caller.
3697
3698 #     This object will call SNPs in *samfile* against the reference
3699 #     sequence in *fasta*.
3700
3701 #     This caller is fast for calling few SNPs in selected regions.
3702
3703 #     It is slow, if called over large genomic regions.
3704 #     '''
3705
3706 #     def __cinit__(self,
3707 #                   IteratorColumn iterator_column,
3708 #                   **kwargs ):
3709
3710 #         pass
3711
3712 #     def call(self, reference, int pos ):
3713 #         """call a snp on chromosome *reference*
3714 #         and position *pos*.
3715
3716 #         returns a :class:`SNPCall` object or None, if no indel call could be made.
3717 #         """
3718
3719 #         cdef int tid = self.iter.samfile.gettid( reference )
3720
3721 #         self.iter.reset( tid, pos, pos + 1 )
3722
3723 #         while 1:
3724 #             self.iter.cnext()
3725
3726 #             if self.iter.n_plp < 0:
3727 #                 raise ValueError("error during iteration" )
3728
3729 #             if self.iter.plp == NULL:
3730 #                 raise ValueError( "no reads in region - no call" )
3731
3732 #             if self.iter.pos == pos: break
3733
3734 #         return self._call()
3735
3736 # cdef class IteratorIndelCalls( IndelCallerBase ):
3737 #     """*(IteratorColumn iterator)*
3738
3739 #     call indels within a region.
3740
3741 #     *iterator* is a pileup iterator. SNPs will be called
3742 #     on all positions returned by this iterator.
3743
3744 #     This caller is fast if SNPs are called over large continuous
3745 #     regions. It is slow, if instantiated frequently and in random
3746 #     order as the sequence will have to be reloaded.
3747
3748 #     """
3749
3750 #     def __cinit__(self,
3751 #                   IteratorColumn iterator_column,
3752 #                   **kwargs ):
3753 #         pass
3754
3755
3756 #     def __iter__(self):
3757 #         return self
3758
3759 #     def __next__(self):
3760 #         """python version of next().
3761 #         """
3762
3763 #         # the following code was adapted from bam_plcmd.c:pileup_func()
3764 #         self.iter.cnext()
3765
3766 #         if self.iter.n_plp < 0:
3767 #             raise ValueError("error during iteration" )
3768
3769 #         if self.iter.plp == NULL:
3770 #            raise StopIteration
3771
3772 #         return self._call()
3773
3774
3775
3776 cdef class IndexedReads:
3777     """index a bamfile by read.
3778
3779     The index is kept in memory.
3780
3781     By default, the file is re-openend to avoid conflicts if
3782     multiple operators work on the same file. Set *reopen* = False
3783     to not re-open *samfile*.
3784     """
3785
3786     def __init__(self, Samfile samfile, int reopen = True ):
3787         self.samfile = samfile
3788
3789         if samfile.isbam: mode = b"rb"
3790         else: mode = b"r"
3791
3792         # reopen the file - note that this makes the iterator
3793         # slow and causes pileup to slow down significantly.
3794         if reopen:
3795             store = StderrStore()
3796             self.fp = samopen( samfile._filename, mode, NULL )
3797             store.release()
3798             assert self.fp != NULL
3799             self.owns_samfile = True
3800         else:
3801             self.fp = samfile.samfile
3802             self.owns_samfile = False
3803
3804         assert samfile.isbam, "can only IndexReads on bam files"
3805
3806     def build( self ):
3807         '''build index.'''
3808
3809         self.index = collections.defaultdict( list )
3810
3811         # this method will start indexing from the current file position
3812         # if you decide
3813         cdef int ret = 1
3814         cdef bam1_t * b = <bam1_t*> calloc(1, sizeof( bam1_t) )
3815
3816         cdef uint64_t pos
3817
3818         while ret > 0:
3819             pos = bam_tell( self.fp.x.bam )
3820             ret = samread( self.fp, b)
3821             if ret > 0:
3822                 qname = _charptr_to_str(bam1_qname( b ))
3823                 self.index[qname].append( pos )
3824
3825         bam_destroy1( b )
3826
3827     def find( self, qname ):
3828         if qname in self.index:
3829             return IteratorRowSelection( self.samfile, self.index[qname], reopen = False )
3830         else:
3831             raise KeyError( "read %s not found" % qname )
3832
3833     def __dealloc__(self):
3834         if self.owns_samfile: samclose( self.fp )
3835
3836 __all__ = ["Samfile",
3837            "Fastafile",
3838            "IteratorRow",
3839            "IteratorColumn",
3840            "AlignedRead",
3841            "PileupColumn",
3842            "PileupProxy",
3843            "PileupRead",
3844            # "IteratorSNPCalls",
3845            # "SNPCaller",
3846            # "IndelCaller",
3847            # "IteratorIndelCalls",
3848            "IndexedReads" ]
3849
3850
3851