pysam/csamtools.pyx

   1 # cython: embedsignature=True
   2 # adds doc-strings for sphinx
   3
   4 import tempfile, os, sys, types, itertools, struct, ctypes
   5
   6 # defines imported from samtools
   7 DEF SEEK_SET = 0
   8 DEF SEEK_CUR = 1
   9 DEF SEEK_END = 2
  10
  11 ## These are bits set in the flag.
  12 ## have to put these definitions here, in csamtools.pxd they got ignored
  13 ## @abstract the read is paired in sequencing, no matter whether it is mapped in a pair */
  14 DEF BAM_FPAIRED       =1
  15 ## @abstract the read is mapped in a proper pair */
  16 DEF BAM_FPROPER_PAIR  =2
  17 ## @abstract the read itself is unmapped; conflictive with BAM_FPROPER_PAIR */
  18 DEF BAM_FUNMAP        =4
  19 ## @abstract the mate is unmapped */
  20 DEF BAM_FMUNMAP       =8
  21 ## @abstract the read is mapped to the reverse strand */
  22 DEF BAM_FREVERSE      =16
  23 ## @abstract the mate is mapped to the reverse strand */
  24 DEF BAM_FMREVERSE     =32
  25 ## @abstract this is read1 */
  26 DEF BAM_FREAD1        =64
  27 ## @abstract this is read2 */
  28 DEF BAM_FREAD2       =128
  29 ## @abstract not primary alignment */
  30 DEF BAM_FSECONDARY   =256
  31 ## @abstract QC failure */
  32 DEF BAM_FQCFAIL      =512
  33 ## @abstract optical or PCR duplicate */
  34 DEF BAM_FDUP        =1024
  35
  36 DEF BAM_CIGAR_SHIFT=4
  37 DEF BAM_CIGAR_MASK=((1 << BAM_CIGAR_SHIFT) - 1)
  38
  39 #####################################################################
  40 #####################################################################
  41 #####################################################################
  42 ## private factory methods
  43 #####################################################################
  44 cdef class AlignedRead
  45 cdef makeAlignedRead( bam1_t * src):
  46     '''enter src into AlignedRead.'''
  47     cdef AlignedRead dest
  48     dest = AlignedRead()
  49     # destroy dummy delegate created in constructor
  50     # to prevent memory leak.
  51     pysam_bam_destroy1(dest._delegate)
  52     dest._delegate = bam_dup1(src)
  53     return dest
  54
  55 cdef class PileupProxy
  56 cdef makePileupProxy( bam_plbuf_t * buf, int n ):
  57      cdef PileupProxy dest
  58      dest = PileupProxy()
  59      dest.buf = buf
  60      dest.n = n
  61      return dest
  62
  63 cdef class PileupRead
  64 cdef makePileupRead( bam_pileup1_t * src ):
  65     '''fill a  PileupRead object from a bam_pileup1_t * object.'''
  66     cdef PileupRead dest
  67     dest = PileupRead()
  68     dest._alignment = makeAlignedRead( src.b )
  69     dest._qpos = src.qpos
  70     dest._indel = src.indel
  71     dest._level = src.level
  72     dest._is_del = src.is_del
  73     dest._is_head = src.is_head
  74     dest._is_tail = src.is_tail
  75     return dest
  76
  77 #####################################################################
  78 #####################################################################
  79 #####################################################################
  80 ## Generic callbacks for inserting python callbacks.
  81 #####################################################################
  82 cdef int fetch_callback( bam1_t *alignment, void *f):
  83     '''callback for bam_fetch.
  84
  85     calls function in *f* with a new :class:`AlignedRead` object as parameter.
  86     '''
  87     a = makeAlignedRead( alignment )
  88     (<object>f)(a)
  89
  90 class PileupColumn(object):
  91     '''A pileup column. A pileup column contains
  92     all the reads that map to a certain target base.
  93
  94     tid
  95         chromosome ID as is defined in the header
  96     pos
  97         the target base coordinate (0-based)
  98     n
  99         number of reads mapping to this column
 100     pileups
 101         list of reads (:class:`pysam.PileupRead`) aligned to this column
 102     '''
 103     def __str__(self):
 104         return "\t".join( map(str, (self.tid, self.pos, self.n))) +\
 105             "\n" + "\n".join( map(str, self.pileups) )
 106
 107 cdef int pileup_callback( uint32_t tid, uint32_t pos, int n, bam_pileup1_t *pl, void *f):
 108     '''callback for pileup.
 109
 110     calls function in *f* with a new :class:`Pileup` object as parameter.
 111
 112     tid
 113         chromosome ID as is defined in the header
 114     pos
 115         start coordinate of the alignment, 0-based
 116     n
 117         number of elements in pl array
 118     pl
 119         array of alignments
 120     data
 121         user provided data
 122     '''
 123
 124     p = PileupColumn()
 125     p.tid = tid
 126     p.pos = pos
 127     p.n = n
 128     pileups = []
 129
 130     for x from 0 <= x < n:
 131         pileups.append( makePileupRead( &(pl[x]) ) )
 132     p.pileups = pileups
 133
 134     (<object>f)(p)
 135
 136 cdef int pileup_fetch_callback( bam1_t *b, void *data):
 137     '''callback for bam_fetch.
 138
 139     Fetches reads and submits them to pileup.
 140     '''
 141     cdef bam_plbuf_t * buf
 142     buf = <bam_plbuf_t*>data
 143     bam_plbuf_push(b, buf)
 144     return 0
 145
 146 class StderrStore():
 147     '''
 148     stderr is captured.
 149     '''
 150     def __init__(self):
 151         self.stderr_h, self.stderr_f = tempfile.mkstemp()
 152         self.stderr_save = Outs( sys.stderr.fileno() )
 153         self.stderr_save.setfd( self.stderr_h )
 154
 155     def release(self):
 156         self.stderr_save.restore()
 157         if os.path.exists(self.stderr_f):
 158             os.remove( self.stderr_f )
 159
 160     def __del__(self):
 161         self.release()
 162
 163 ######################################################################
 164 ######################################################################
 165 ######################################################################
 166 # valid types for sam headers
 167 VALID_HEADER_TYPES = { "HD" : dict,
 168                        "SQ" : list,
 169                        "RG" : list,
 170                        "PG" : list,
 171                        "CO" : list }
 172
 173 # order of records within sam headers
 174 VALID_HEADERS = ("HD", "SQ", "RG", "PG", "CO" )
 175
 176 # type conversions within sam header records
 177 VALID_HEADER_FIELDS = { "HD" : { "VN" : str, "SO" : str, "GO" : str },
 178                         "SQ" : { "SN" : str, "LN" : int, "AS" : str, "M5" : str, "UR" : str, "SP" : str },
 179                         "RG" : { "ID" : str, "SM" : str, "LB" : str, "DS" : str, "PU" : str, "PI" : str, "CN" : str, "DT" : str, "PL" : str, },
 180                         "PG" : { "ID" : str, "VN" : str, "CL" : str }, }
 181
 182 # output order of fields within records
 183 VALID_HEADER_ORDER = { "HD" : ( "VN", "SO", "GO" ),
 184                        "SQ" : ( "SN", "LN", "AS", "M5" , "UR" , "SP" ),
 185                        "RG" : ( "ID", "SM", "LB", "DS" , "PU" , "PI" , "CN" , "DT", "PL" ),
 186                        "PG" : ( "ID", "VN", "CL" ), }
 187
 188 ######################################################################
 189 ######################################################################
 190 ######################################################################
 191 ## Public methods
 192 ######################################################################
 193 cdef class Samfile:
 194     '''*(filename, mode='r', template = None, referencenames = None, referencelengths = None, text = NULL, header = None)*
 195
 196     A *SAM* file. The file is automatically opened.
 197
 198     *mode* should be ``r`` for reading or ``w`` for writing. The default is text mode so for binary
 199     (:term:`BAM`) I/O you should append ``b`` for compressed or ``u`` for uncompressed :term:`BAM` output.
 200     Use ``h`` to output header information  in text (:term:`TAM`)  mode.
 201
 202     If ``b`` is present, it must immediately follow ``r`` or ``w``.
 203     Currently valid modes are ``r``, ``w``, ``wh``, ``rb``, ``wb`` and ``wbu``.
 204
 205     so to open a :term:`BAM` file for reading::
 206
 207         f=Samfile('ex1.bam','rb')
 208
 209
 210     For writing, the header of a :term:`TAM` file/:term:`BAM` file can be constituted from several
 211     sources:
 212
 213         1. If *template* is given, the header is copied from a another *Samfile* (*template* must be of type *Samfile*).
 214
 215         2. If *header* is given, the header is build from a multi-level dictionary. The first level are the four types ('HD', 'SQ', ...). The second level is then a list of lines, with each line being a list of tag-value pairs.
 216
 217         3. If *text* is given, new header text is copied from raw text.
 218
 219         4. The names (*referencenames*) and lengths (*referencelengths*) are supplied directly as lists.
 220
 221     If an index for a BAM file exists (.bai), it will be opened automatically. Without an index random
 222     access to reads via :meth:`fetch` and :meth:`pileup` is disabled.
 223     '''
 224
 225     cdef char * filename
 226     # pointer to samfile
 227     cdef samfile_t * samfile
 228     # pointer to index
 229     cdef bam_index_t *index
 230     # true if file is a bam file
 231     cdef int isbam
 232
 233     # current read within iteration
 234     cdef bam1_t * b
 235
 236     def __cinit__(self, *args, **kwargs ):
 237         self.samfile = NULL
 238         self.isbam = False
 239         self._open( *args, **kwargs )
 240
 241         # allocate memory for iterator
 242         self.b = <bam1_t*>calloc(1, sizeof(bam1_t))
 243
 244     def _isOpen( self ):
 245         '''return true if samfile has been opened.'''
 246         return self.samfile != NULL
 247
 248     def _hasIndex( self ):
 249         '''return true if samfile has an existing (and opened) index.'''
 250         return self.index != NULL
 251
 252     def _open( self,
 253                char * filename,
 254                mode ='r',
 255                Samfile template = None,
 256                referencenames = None,
 257                referencelengths = None,
 258                char * text = NULL,
 259                header = None,
 260               ):
 261         '''open a sam/bam file.
 262
 263         If _open is called on an existing bamfile, the current file will be
 264         closed and a new file will be opened.
 265         '''
 266
 267         assert mode in ( "r","w","rb","wb", "wh", "wbu" ), "invalid file opening mode `%s`" % mode
 268
 269         # close a previously opened file
 270         if self.samfile != NULL: self.close()
 271         self.samfile = NULL
 272
 273         cdef bam_header_t * header_to_write
 274         header_to_write = NULL
 275
 276         self.filename = filename
 277
 278         self.isbam = len(mode) > 1 and mode[1] == 'b'
 279
 280         if mode[0] == 'w':
 281             # open file for writing
 282
 283             # header structure (used for writing)
 284             if template:
 285                 # copy header from another file
 286                 header_to_write = template.samfile.header
 287
 288             elif header:
 289                 header_to_write = self._buildHeader( header )
 290
 291             else:
 292                 # build header from a target names and lengths
 293                 assert referencenames and referencelengths, "either supply options `template`, `header` or  both `refernencenames` and `referencelengths` for writing"
 294                 assert len(referencenames) == len(referencelengths), "unequal names and lengths of reference sequences"
 295
 296                 # allocate and fill header
 297                 header_to_write = bam_header_init()
 298                 header_to_write.n_targets = len(referencenames)
 299                 n = 0
 300                 for x in referencenames: n += len(x) + 1
 301                 header_to_write.target_name = <char**>calloc(n, sizeof(char*))
 302                 header_to_write.target_len = <uint32_t*>calloc(n, sizeof(uint32_t))
 303                 for x from 0 <= x < header_to_write.n_targets:
 304                     header_to_write.target_len[x] = referencelengths[x]
 305                     name = referencenames[x]
 306                     header_to_write.target_name[x] = <char*>calloc(len(name)+1, sizeof(char))
 307                     strncpy( header_to_write.target_name[x], name, len(name) )
 308
 309                 if text != NULL:
 310                     # copy without \0
 311                     header_to_write.l_text = strlen(text)
 312                     header_to_write.text = <char*>calloc( strlen(text), sizeof(char) )
 313                     memcpy( header_to_write.text, text, strlen(text) )
 314
 315                 header_to_write.hash = NULL
 316                 header_to_write.rg2lib = NULL
 317
 318             # open file. Header gets written to file at the same time for bam files
 319             # and sam files (in the latter case, the mode needs to be wh)
 320             store = StderrStore()
 321             self.samfile = samopen( filename, mode, header_to_write )
 322             store.release()
 323
 324             # bam_header_destroy takes care of cleaning up of all the members
 325             if not template and header_to_write != NULL:
 326                 bam_header_destroy( header_to_write )
 327
 328         elif mode[0] == "r":
 329             # open file for reading
 330             if strncmp( filename, "-", 1) != 0 and not os.path.exists( filename ):
 331                 raise IOError( "file `%s` not found" % filename)
 332
 333             store = StderrStore()
 334             self.samfile = samopen( filename, mode, NULL )
 335             store.release()
 336
 337         if self.samfile == NULL:
 338             raise IOError("could not open file `%s`" % filename )
 339
 340         if mode[0] == "r" and self.isbam:
 341             if not os.path.exists(filename + ".bai"):
 342                 self.index = NULL
 343             else:
 344                 # returns NULL if there is no index or index could not be opened
 345                 self.index = bam_index_load(filename)
 346                 if self.index == NULL:
 347                     raise IOError("error while opening index `%s` " % filename )
 348
 349     def getrname( self, tid ):
 350         '''(tid )
 351         convert numerical :term:`tid` into :ref:`reference` name.'''
 352         if not 0 <= tid < self.samfile.header.n_targets:
 353             raise ValueError( "tid out of range 0<=tid<%i" % self.samfile.header.n_targets )
 354         return self.samfile.header.target_name[tid]
 355
 356     def _parseRegion( self,
 357                       reference = None,
 358                       start = None,
 359                       end = None,
 360                       region = None ):
 361         '''parse region information.
 362
 363         raise Value for for invalid regions.
 364
 365         returns a tuple of region, tid, start and end. Region
 366         is a valid samtools :term:`region` or None if the region
 367         extends over the whole file.
 368
 369         Note that regions are 1-based, while start,end are python coordinates.
 370         '''
 371
 372         cdef int rtid
 373         cdef int rstart
 374         cdef int rend
 375         cdef int max_pos
 376         max_pos = 2 << 29
 377
 378         rtid = rstart = rend = 0
 379
 380         # translate to a region
 381         if reference:
 382             if start != None and end != None:
 383                 region = "%s:%i-%i" % (reference, start+1, end)
 384             else:
 385                 region = reference
 386
 387         if region:
 388             store = StderrStore()
 389             bam_parse_region( self.samfile.header, region, &rtid, &rstart, &rend)
 390             store.release()
 391             if rtid < 0: raise ValueError( "invalid region `%s`" % region )
 392             if rstart > rend: raise ValueError( 'invalid region: start (%i) > end (%i)' % (rstart, rend) )
 393             if not 0 <= rstart < max_pos: raise ValueError( 'start out of range (%i)' % rstart )
 394             if not 0 <= rend < max_pos: raise ValueError( 'end out of range (%i)' % rend )
 395
 396         return region, rtid, rstart, rend
 397
 398     def fetch( self,
 399                reference = None,
 400                start = None,
 401                end = None,
 402                region = None,
 403                callback = None,
 404                until_eof = False ):
 405         '''*(reference = None, start = None, end = None, region = None, callback = None, until_eof = False)*
 406
 407         fetch :meth:`AlignedRead` objects in a :term:`region` using 0-based indexing. The region is specified by
 408         :term:`reference`, *start* and *end*. Alternatively, a samtools :term:`region` string can be supplied.
 409
 410         Without *reference* or *region* all reads will be fetched. The reads will be returned
 411         ordered by reference sequence, which will not necessarily be the order within the file.
 412         If *until_eof* is given, all reads from the current file position will be returned
 413         *as they are sorted within the file*.
 414
 415         If only *reference* is set, all reads matching on *reference* will be fetched.
 416
 417         The method returns an iterator of type :class:`pysam.IteratorRow` unless
 418         a *callback is provided. If *callback* is given, the callback will be executed
 419         for each position within the :term:`region`. Note that callbacks currently work
 420         only, if *region* or *reference* is given.
 421
 422         Note that a :term:`TAM` file does not allow random access. If *region* or *reference* are given,
 423         an exception is raised.
 424         '''
 425         cdef int rtid
 426         cdef int rstart
 427         cdef int rend
 428
 429         if not self._isOpen():
 430             raise ValueError( "I/O operation on closed file" )
 431
 432         region, rtid, rstart, rend = self._parseRegion( reference, start, end, region )
 433
 434         if self.isbam:
 435             if callback:
 436                 if not region:
 437                     raise ValueError( "callback functionality requires a region/reference" )
 438                 if not self._hasIndex(): raise ValueError( "no index available for fetch" )
 439                 return bam_fetch(self.samfile.x.bam,
 440                                  self.index, rtid, rstart, rend, <void*>callback, fetch_callback )
 441             else:
 442                 if region:
 443                     return IteratorRow( self, rtid, rstart, rend )
 444                 else:
 445                     if until_eof:
 446                         return IteratorRowAll( self )
 447                     else:
 448                         # return all targets by chaining the individual targets together.
 449                         if not self._hasIndex(): raise ValueError( "no index available for fetch" )
 450                         i = []
 451                         rstart = 0
 452                         rend = 1<<29
 453                         for rtid from 0 <= rtid < self.nreferences:
 454                             i.append( IteratorRow( self, rtid, rstart, rend))
 455                         return itertools.chain( *i )
 456         else:
 457             if region != None:
 458                 raise ValueError ("fetch for a region is not available for sam files" )
 459             if callback:
 460                 raise NotImplementedError( "callback not implemented yet" )
 461             else:
 462                 return IteratorRowAll( self )
 463
 464     def pileup( self, reference = None, start = None, end = None, region = None, callback = None ):
 465         '''run a pileup within a :term:`region` using 0-based indexing. The region is specified by
 466         :term:`reference`, *start* and *end*. Alternatively, a samtools *region* string can be supplied.
 467
 468         Without *reference* or *region* all reads will be fetched. The reads will be returned
 469         ordered by :term:`reference` sequence, which will not necessarily be the order within the file.
 470
 471         The method returns an iterator of type :class:`pysam.IteratorColumn` unless
 472         a *callback is provided. If *callback* is given, the callback will be executed
 473         for each position within the :term:`region`.
 474
 475         Note that samfiles do not allow random access. If *region* or *reference* are given,
 476         an exception is raised.
 477
 478         .. Note::
 479
 480             *all* reads which overlap the region are returned. The first base returned will be the
 481             first base of the first read *not* necessarily the first base of the region used in the query.
 482         '''
 483         cdef int rtid
 484         cdef int rstart
 485         cdef int rend
 486         cdef bam_plbuf_t *buf
 487
 488         if not self._isOpen():
 489             raise ValueError( "I/O operation on closed file" )
 490
 491         region, rtid, rstart, rend = self._parseRegion( reference, start, end, region )
 492
 493         if self.isbam:
 494             if not self._hasIndex(): raise ValueError( "no index available for pileup" )
 495
 496             if callback:
 497                 if not region:
 498                     raise ValueError( "callback functionality requires a region/reference" )
 499
 500                 buf = bam_plbuf_init( <bam_pileup_f>pileup_callback, <void*>callback )
 501                 bam_fetch(self.samfile.x.bam,
 502                           self.index, rtid, rstart, rend,
 503                           buf, pileup_fetch_callback )
 504
 505                 # finalize pileup
 506                 bam_plbuf_push( NULL, buf)
 507                 bam_plbuf_destroy(buf)
 508             else:
 509                 if region:
 510                     return IteratorColumn( self, rtid, rstart, rend )
 511                 else:
 512                     # return all targets by chaining the individual targets together.
 513                     i = []
 514                     rstart = 0
 515                     rend = 1<<29
 516                     for rtid from 0 <= rtid < self.nreferences:
 517                         i.append( IteratorColumn( self, rtid, rstart, rend))
 518                     return itertools.chain( *i )
 519
 520         else:
 521             raise NotImplementedError( "pileup of samfiles not implemented yet" )
 522
 523     def close( self ):
 524         '''closes file.'''
 525         if self.samfile != NULL:
 526             samclose( self.samfile )
 527             bam_index_destroy(self.index);
 528             self.samfile = NULL
 529
 530     def __dealloc__( self ):
 531         '''clean up.'''
 532         # remember: dealloc cannot call other methods
 533         # Note that __del__ is not called.
 534         self.close()
 535         pysam_bam_destroy1(self.b)
 536
 537     def write( self, AlignedRead read ):
 538         '''(AlignedRead read )
 539         write a single :class:`pysam.AlignedRead`..
 540
 541         return the number of bytes written.
 542         '''
 543         return samwrite( self.samfile, read._delegate )
 544
 545     property nreferences:
 546         '''number of :term:`reference` sequences in the file.'''
 547         def __get__(self):
 548             return self.samfile.header.n_targets
 549
 550     property references:
 551         """tuple with the names of :term:`reference` sequences."""
 552         def __get__(self):
 553             t = []
 554             for x from 0 <= x < self.samfile.header.n_targets:
 555                 t.append( self.samfile.header.target_name[x] )
 556             return tuple(t)
 557
 558     property lengths:
 559         """tuple of the lengths of the :term:`reference` sequences. The lengths are in the same order as :attr:`pysam.Samfile.reference`
 560         """
 561         def __get__(self):
 562             t = []
 563             for x from 0 <= x < self.samfile.header.n_targets:
 564                 t.append( self.samfile.header.target_len[x] )
 565             return tuple(t)
 566
 567     property text:
 568         '''full contents of the :term:`sam file` header as a string.'''
 569         def __get__(self):
 570             # create a temporary 0-terminated copy
 571             cdef char * t
 572             t = <char*>calloc( self.samfile.header.l_text + 1, sizeof(char) )
 573             memcpy( t, self.samfile.header.text, self.samfile.header.l_text )
 574             result = t
 575             free(t)
 576             return result
 577
 578     property header:
 579         '''header information within the :term:`sam file`. The records and fields are returned as
 580         a two-level dictionary.
 581         '''
 582         def __get__(self):
 583             result = {}
 584
 585             if self.samfile.header.text != NULL:
 586                 # convert to python string (note: call self.text to create 0-terminated string)
 587                 t = self.text
 588                 for line in t.split("\n"):
 589                     if not line.strip(): continue
 590                     assert line.startswith("@"), "header line without '@': '%s'" % line
 591                     fields = line[1:].split("\t")
 592                     record = fields[0]
 593                     assert record in VALID_HEADER_TYPES, "header line with invalid type '%s': '%s'" % (record, line)
 594
 595                     # treat comments
 596                     if record == "CO":
 597                         if record not in result: result[record] = []
 598                         result[record].append( "\t".join( fields[1:] ) )
 599                         continue
 600
 601                     # the following is clumsy as generators do not work?
 602                     x = {}
 603                     for field in fields[1:]:
 604                         key, value = field.split(":",1)
 605                         if key not in VALID_HEADER_FIELDS[record]:
 606                             raise ValueError( "unknown field code '%s' in record '%s'" % (key, record) )
 607                         x[key] = VALID_HEADER_FIELDS[record][key](value)
 608
 609                     if VALID_HEADER_TYPES[record] == dict:
 610                         if record in result:
 611                             raise ValueError( "multiple '%s' lines are not permitted" % record )
 612                         result[record] = x
 613                     elif VALID_HEADER_TYPES[record] == list:
 614                         if record not in result: result[record] = []
 615                         result[record].append( x )
 616
 617             return result
 618
 619     def _buildLine( self, fields, record ):
 620         '''build a header line from *fields* dictionary for *record*'''
 621
 622         # TODO: add checking for field and sort order
 623         line = ["@%s" % record ]
 624         if record == "CO":
 625             line.append( fields )
 626         else:
 627             for key in VALID_HEADER_ORDER[record]:
 628                 if key in fields:
 629                     line.append( "%s:%s" % (key, str(fields[key])))
 630         return "\t".join( line )
 631
 632     cdef bam_header_t * _buildHeader( self, new_header ):
 633         '''return a new header built from a dictionary in *new_header*.
 634
 635         This method inserts the text field, target_name and target_len.
 636         '''
 637
 638         lines = []
 639
 640         # check if hash exists
 641
 642         # create new header and copy old data
 643         cdef bam_header_t * dest
 644
 645         dest = bam_header_init()
 646
 647         for record in VALID_HEADERS:
 648             if record in new_header:
 649                 ttype = VALID_HEADER_TYPES[record]
 650                 data = new_header[record]
 651                 if type( data ) != type( ttype() ):
 652                     raise ValueError( "invalid type for record %s: %s, expected %s" % (record, type(data), type(ttype()) ) )
 653                 if type( data ) == types.DictType:
 654                     lines.append( self._buildLine( data, record ) )
 655                 else:
 656                     for fields in new_header[record]:
 657                         lines.append( self._buildLine( fields, record ) )
 658
 659         text = "\n".join(lines) + "\n"
 660         if dest.text != NULL: free( dest.text )
 661         dest.text = <char*>calloc( len(text), sizeof(char))
 662         dest.l_text = len(text)
 663         strncpy( dest.text, text, dest.l_text )
 664
 665         # collect targets
 666         if "SQ" in new_header:
 667             seqs = []
 668             for fields in new_header["SQ"]:
 669                 try:
 670                     seqs.append( (fields["SN"], fields["LN"] ) )
 671                 except KeyError:
 672                     raise KeyError( "incomplete sequence information in '%s'" % str(fields))
 673
 674             dest.n_targets = len(seqs)
 675             dest.target_name = <char**>calloc( dest.n_targets, sizeof(char*) )
 676             dest.target_len = <uint32_t*>calloc( dest.n_targets, sizeof(uint32_t) )
 677
 678             for x from 0 <= x < dest.n_targets:
 679                 seqname, seqlen = seqs[x]
 680                 dest.target_name[x] = <char*>calloc( len( seqname ) + 1, sizeof(char) )
 681                 strncpy( dest.target_name[x], seqname, len(seqname) + 1 )
 682                 dest.target_len[x] = seqlen
 683
 684         return dest
 685
 686     def __iter__(self):
 687         return self
 688
 689     cdef bam1_t * getCurrent( self ):
 690         return self.b
 691
 692     cdef int cnext(self):
 693         '''cversion of iterator. Used by IteratorColumn'''
 694         cdef int ret
 695         return samread(self.samfile, self.b)
 696
 697     def __next__(self):
 698         """python version of next().
 699
 700         pyrex uses this non-standard name instead of next()
 701         """
 702         cdef int ret
 703         ret = samread(self.samfile, self.b)
 704         if (ret > 0):
 705             return makeAlignedRead( self.b )
 706         else:
 707             raise StopIteration
 708
 709 cdef class Fastafile:
 710     '''*(filename)*
 711
 712     A *FASTA* file. The file is automatically opened.
 713
 714     The file expects an indexed fasta file.
 715
 716     TODO:
 717         add automatic indexing.
 718         add function to get sequence names.
 719     '''
 720
 721     cdef char * filename
 722     # pointer to fastafile
 723     cdef faidx_t * fastafile
 724
 725     def __cinit__(self, *args, **kwargs ):
 726         self.fastafile = NULL
 727         self._open( *args, **kwargs )
 728
 729     def _isOpen( self ):
 730         '''return true if samfile has been opened.'''
 731         return self.fastafile != NULL
 732
 733     def _open( self,
 734                char * filename ):
 735         '''open an indexed fasta file.
 736
 737         This method expects an indexed fasta file.
 738         '''
 739
 740         # close a previously opened file
 741         if self.fastafile != NULL: self.close()
 742         self.filename = filename
 743         self.fastafile = fai_load( filename )
 744
 745         if self.fastafile == NULL:
 746             raise IOError("could not open file `%s`" % filename )
 747
 748     def close( self ):
 749         if self.fastafile != NULL:
 750             fai_destroy( self.fastafile )
 751             self.fastafile = NULL
 752
 753     def fetch( self,
 754                reference = None,
 755                start = None,
 756                end = None,
 757                region = None):
 758
 759         '''*(reference = None, start = None, end = None, region = None)*
 760
 761         fetch :meth:`AlignedRead` objects in a :term:`region` using 0-based indexing. The region is specified by
 762         :term:`reference`, *start* and *end*. Alternatively, a samtools :term:`region` string can be supplied.
 763         '''
 764
 765         if not self._isOpen():
 766             raise ValueError( "I/O operation on closed file" )
 767
 768         cdef int len, max_pos
 769         cdef char * seq
 770         max_pos = 2 << 29
 771
 772         if not region:
 773             if reference == None: raise ValueError( 'no sequence/region supplied.' )
 774             if start == None and end == None:
 775                 region = "%s" % str(reference)
 776             elif start == None or end == None:
 777                 raise ValueError( 'only start or only end of region supplied' )
 778             else:
 779                 if start > end: raise ValueError( 'invalid region: start (%i) > end (%i)' % (start, end) )
 780                 # valid ranges are from 0 to 2^29-1
 781                 if not 0 <= start < max_pos: raise ValueError( 'start out of range (%i)' % start )
 782                 if not 0 <= end < max_pos: raise ValueError( 'end out of range (%i)' % end )
 783                 region = "%s:%i-%i" % (reference, start+1, end )
 784
 785         # samtools adds a '\0' at the end
 786         seq = fai_fetch( self.fastafile, region, &len )
 787         # copy to python
 788         result = seq
 789         # clean up
 790         free(seq)
 791
 792         return result
 793
 794 ## turning callbacks elegantly into iterators is an unsolved problem, see the following threads:
 795 ## http://groups.google.com/group/comp.lang.python/browse_frm/thread/0ce55373f128aa4e/1d27a78ca6408134?hl=en&pli=1
 796 ## http://www.velocityreviews.com/forums/t359277-turning-a-callback-function-into-a-generator.html
 797 ## Thus I chose to rewrite the functions requiring callbacks. The downside is that if the samtools C-API or code
 798 ## changes, the changes have to be manually entered.
 799
 800 cdef class IteratorRow:
 801     """iterates over mapped reads in a region.
 802     """
 803
 804     cdef bam_fetch_iterator_t*  bam_iter # iterator state object
 805     cdef bam1_t *               b
 806     cdef                        error_msg
 807     cdef int                    error_state
 808     cdef Samfile                samfile
 809     def __cinit__(self, Samfile samfile, int tid, int beg, int end ):
 810         self.bam_iter = NULL
 811
 812         assert samfile._isOpen()
 813         assert samfile._hasIndex()
 814
 815         # makes sure that samfile stays alive as long as the
 816         # iterator is alive.
 817         self.samfile = samfile
 818
 819         # parse the region
 820         self.error_state = 0
 821         self.error_msg = None
 822
 823         cdef bamFile  fp
 824         fp = samfile.samfile.x.bam
 825         self.bam_iter = bam_init_fetch_iterator(fp, samfile.index, tid, beg, end)
 826
 827     def __iter__(self):
 828         return self
 829
 830     cdef bam1_t * getCurrent( self ):
 831         return self.b
 832
 833     cdef int cnext(self):
 834         '''cversion of iterator. Used by IteratorColumn'''
 835         self.b = bam_fetch_iterate(self.bam_iter)
 836         if self.b == NULL: return 0
 837         return 1
 838
 839     def __next__(self):
 840         """python version of next().
 841
 842         pyrex uses this non-standard name instead of next()
 843         """
 844         if self.error_state:
 845             raise ValueError( self.error_msg)
 846
 847         self.b = bam_fetch_iterate(self.bam_iter)
 848         if self.b != NULL:
 849             return makeAlignedRead( self.b )
 850         else:
 851             raise StopIteration
 852
 853     def __dealloc__(self):
 854         '''remember: dealloc cannot call other methods!'''
 855         if self.bam_iter:
 856             bam_cleanup_fetch_iterator(self.bam_iter)
 857
 858 cdef class IteratorRowAll:
 859     """iterates over all mapped reads
 860     """
 861
 862     cdef bam1_t * b
 863     cdef samfile_t * fp
 864
 865     def __cinit__(self, Samfile samfile):
 866
 867         assert samfile._isOpen()
 868
 869         self.fp = samfile.samfile
 870
 871         # allocate memory for alignment
 872         self.b = <bam1_t*>calloc(1, sizeof(bam1_t))
 873
 874     def __iter__(self):
 875         return self
 876
 877     cdef bam1_t * getCurrent( self ):
 878         return self.b
 879
 880     cdef int cnext(self):
 881         '''cversion of iterator. Used by IteratorColumn'''
 882         cdef int ret
 883         return samread(self.fp, self.b)
 884
 885     def __next__(self):
 886         """python version of next().
 887
 888         pyrex uses this non-standard name instead of next()
 889         """
 890         cdef int ret
 891         ret = samread(self.fp, self.b)
 892         if (ret > 0):
 893             return makeAlignedRead( self.b )
 894         else:
 895             raise StopIteration
 896
 897     def __dealloc__(self):
 898         '''remember: dealloc cannot call other methods!'''
 899         pysam_bam_destroy1(self.b)
 900
 901 cdef class IteratorColumn:
 902     '''iterates over columns.
 903
 904     This iterator wraps the pileup functionality of samtools.
 905
 906     For reasons of efficiency, the iterator returns the current
 907     pileup buffer. As this buffer is updated at every iteration,
 908     the contents of this iterator will change accordingly. Hence the conversion to
 909     a list will not produce the expected result::
 910
 911        f = Samfile("file.bam", "rb")
 912        result = list( f.pileup() )
 913
 914     Here, result will contain ``n`` objects of type :class:`PileupProxy` for ``n`` columns,
 915     but each object will contain the same information.
 916
 917     If the results of several columns are required at the same time, the results
 918     need to be stored explicitely::
 919
 920        result = [ x.pileups() for x in f.pileup() ]
 921
 922     Here, result will be a list of ``n`` lists of objects of type :class:`PileupRead`.
 923
 924     '''
 925     cdef bam_plbuf_t *buf
 926
 927     # check if first iteration
 928     cdef int notfirst
 929     # result of the last plbuf_push
 930     cdef int n_pu
 931     cdef int eof
 932     cdef IteratorRow iter
 933
 934     def __cinit__(self, Samfile samfile, int tid, int start, int end ):
 935
 936         self.iter = IteratorRow( samfile, tid, start, end )
 937         self.buf = bam_plbuf_init(NULL, NULL )
 938         self.n_pu = 0
 939         self.eof = 0
 940
 941     def __iter__(self):
 942         return self
 943
 944     cdef int cnext(self):
 945         '''perform next iteration.
 946
 947         return 1 if there is a buffer to emit. Return 0 for end of iteration.
 948         '''
 949
 950         cdef int retval1, retval2
 951
 952         # pysam bam_plbuf_push returns:
 953         # 1: if buf is full and can be emitted
 954         # 0: if b has been added
 955         # -1: if there was an error
 956
 957         # check if previous plbuf was incomplete. If so, continue within
 958         # the loop and yield if necessary
 959         if self.n_pu > 0:
 960             self.n_pu = pysam_bam_plbuf_push( self.iter.getCurrent(), self.buf, 1)
 961             if self.n_pu > 0: return 1
 962
 963         if self.eof: return 0
 964
 965         # get next alignments and submit until plbuf indicates that
 966         # an new column has finished
 967         while self.n_pu == 0:
 968             retval1 = self.iter.cnext()
 969             # wrap up if no more input
 970             if retval1 == 0:
 971                 self.n_pu = pysam_bam_plbuf_push( NULL, self.buf, 0)
 972                 self.eof = 1
 973                 return self.n_pu
 974
 975             # submit to plbuf
 976             self.n_pu = pysam_bam_plbuf_push( self.iter.getCurrent(), self.buf, 0)
 977             if self.n_pu < 0: raise ValueError( "error while iterating" )
 978
 979         # plbuf has yielded
 980         return 1
 981
 982     def __next__(self):
 983         """python version of next().
 984
 985         pyrex uses this non-standard name instead of next()
 986         """
 987         cdef int ret
 988         ret = self.cnext()
 989         cdef bam_pileup1_t * pl
 990
 991         if ret > 0 :
 992             return makePileupProxy( self.buf, self.n_pu )
 993         else:
 994             raise StopIteration
 995
 996     def __dealloc__(self):
 997         bam_plbuf_destroy(self.buf);
 998
 999 cdef class AlignedRead:
1000     '''
1001     Class representing an aligned read. see SAM format specification for meaning of fields (http://samtools.sourceforge.net/).
1002
1003     This class stores a handle to the samtools C-structure representing
1004     an aligned read. Member read access is forwarded to the C-structure
1005     and converted into python objects. This implementation should be fast,
1006     as only the data needed is converted.
1007
1008     For write access, the C-structure is updated in-place. This is
1009     not the most efficient way to build BAM entries, as the variable
1010     length data is concatenated and thus needs to resized if
1011     a field is updated. Furthermore, the BAM entry might be
1012     in an inconsistent state. The :meth:`~validate` method can
1013     be used to check if an entry is consistent.
1014
1015     One issue to look out for is that the sequence should always
1016     be set *before* the quality scores. Setting the sequence will
1017     also erase any quality scores that were set previously.
1018     '''
1019     cdef:
1020          bam1_t * _delegate
1021
1022     def __cinit__( self ):
1023         # see bam_init1
1024         self._delegate = <bam1_t*>calloc( 1, sizeof( bam1_t) )
1025         # allocate some memory
1026         # If size is 0, calloc does not return a pointer that can be passed to free()
1027         # so allocate 40 bytes for a new read
1028         self._delegate.m_data = 40
1029         self._delegate.data = <uint8_t *>calloc( self._delegate.m_data, 1 )
1030         self._delegate.data_len = 0
1031
1032     def __dealloc__(self):
1033         '''clear up memory.'''
1034         pysam_bam_destroy1(self._delegate)
1035
1036     def __str__(self):
1037         """todo"""
1038         return "\t".join(map(str, (self.qname,
1039                                    self.rname,
1040                                    self.pos,
1041                                    self.cigar,
1042                                    self.qual,
1043                                    self.flag,
1044                                    self.seq,
1045                                    self.mapq,
1046                                    self.tags)))
1047
1048
1049     def __cmp__(self, AlignedRead other):
1050         '''return true, if contents in this are binary equal to ``other``.'''
1051         cdef int retval, x
1052         cdef bam1_t *t, *o
1053         t = self._delegate
1054         o = other._delegate
1055
1056         # uncomment for debugging purposes
1057         # cdef unsigned char * oo, * tt
1058         # tt = <unsigned char*>(&t.core)
1059         # oo = <unsigned char*>(&o.core)
1060         # for x from 0 <= x < sizeof( bam1_core_t): print x, tt[x], oo[x]
1061         # tt = <unsigned char*>(t.data)
1062         # oo = <unsigned char*>(o.data)
1063         # for x from 0 <= x < max(t.data_len, o.data_len): print x, tt[x], oo[x], chr(tt[x]), chr(oo[x])
1064
1065         retval = memcmp( &t.core,
1066                           &o.core,
1067                           sizeof( bam1_core_t ))
1068
1069         if retval: return retval
1070         retval = cmp( t.data_len, o.data_len)
1071         if retval: return retval
1072         return memcmp( t.data,
1073                        o.data,
1074                        sizeof( t.data_len ))
1075
1076     property qname:
1077         """the query name (None if not present)"""
1078         def __get__(self):
1079             cdef bam1_t * src
1080             src = self._delegate
1081             if src.core.l_qname == 0: return None
1082             return <char *>pysam_bam1_qname( src )
1083
1084         def __set__(self, qname ):
1085             if qname == None or len(qname) == 0: return
1086             cdef bam1_t * src
1087             cdef int l
1088             cdef char * p
1089
1090             src = self._delegate
1091             p = pysam_bam1_qname( src )
1092
1093             # the qname is \0 terminated
1094             l = len(qname) + 1
1095             pysam_bam_update( src,
1096                               src.core.l_qname,
1097                               l,
1098                               <uint8_t*>p )
1099
1100             src.core.l_qname = l
1101
1102             # re-acquire pointer to location in memory
1103             # as it might have moved
1104             p = pysam_bam1_qname(src)
1105
1106             strncpy( p, qname, l )
1107
1108     property cigar:
1109         """the :term:`cigar` alignment (None if not present).
1110         """
1111         def __get__(self):
1112             cdef uint32_t * cigar_p
1113             cdef bam1_t * src
1114             cdef op, l, cigar
1115             src = self._delegate
1116             if src.core.n_cigar == 0: return None
1117
1118             cigar = []
1119             cigar_p = pysam_bam1_cigar(src);
1120             for k from 0 <= k < src.core.n_cigar:
1121                 op = cigar_p[k] & BAM_CIGAR_MASK
1122                 l = cigar_p[k] >> BAM_CIGAR_SHIFT
1123                 cigar.append((op, l))
1124             return cigar
1125
1126         def __set__(self, values ):
1127             if values == None or len(values) == 0: return
1128             cdef uint32_t * p
1129             cdef bam1_t * src
1130             cdef op, l
1131             cdef int k
1132
1133             k = 0
1134
1135             src = self._delegate
1136
1137             # get location of cigar string
1138             p = pysam_bam1_cigar(src)
1139
1140             # create space for cigar data within src.data
1141             pysam_bam_update( src,
1142                               src.core.n_cigar * 4,
1143                               len(values) * 4,
1144                               p )
1145
1146             # length is number of cigar operations, not bytes
1147             src.core.n_cigar = len(values)
1148
1149             # re-acquire pointer to location in memory
1150             # as it might have moved
1151             p = pysam_bam1_cigar(src)
1152
1153             # insert cigar operations
1154             for op, l in values:
1155                 p[k] = l << BAM_CIGAR_SHIFT | op
1156                 k += 1
1157
1158             ## setting the cigar string also updates the "bin" attribute
1159             src.core.bin = bam_reg2bin( src.core.pos, bam_calend( &src.core, p))
1160
1161     property seq:
1162         """the query sequence (None if not present)"""
1163         def __get__(self):
1164             cdef bam1_t * src
1165             cdef uint8_t * p
1166             cdef char * s
1167             src = self._delegate
1168             bam_nt16_rev_table = "=ACMGRSVTWYHKDBN"
1169             ## parse qseq (bam1_seq)
1170             if src.core.l_qseq == 0: return None
1171
1172             s = < char *> calloc(src.core.l_qseq + 1 , sizeof(char))
1173             p = pysam_bam1_seq( src )
1174             for k from 0 <= k < src.core.l_qseq:
1175             ## equivalent to bam_nt16_rev_table[bam1_seqi(s, i)] (see bam.c)
1176                 s[k] = "=ACMGRSVTWYHKDBN"[((p)[(k) / 2] >> 4 * (1 - (k) % 2) & 0xf)]
1177             retval=s
1178             free(s)
1179             return retval
1180
1181         def __set__(self,seq):
1182             # samtools manages sequence and quality length memory together
1183             # if no quality information is present, the first byte says 0xff.
1184
1185             if seq == None or len(seq) == 0: return
1186             cdef bam1_t * src
1187             cdef uint8_t * p
1188             cdef char * s
1189             src = self._delegate
1190             cdef int l, k, nbytes_new, nbytes_old
1191
1192             l = len(seq)
1193
1194             # as the sequence is stored in half-bytes, the total length (sequence
1195             # plus quality scores) is (l+1)/2 + l
1196             nbytes_new = (l+1)/2 + l
1197             nbytes_old = (src.core.l_qseq+1)/2 + src.core.l_qseq
1198             # acquire pointer to location in memory
1199             p = pysam_bam1_seq( src )
1200             src.core.l_qseq = l
1201
1202             pysam_bam_update( src,
1203                               nbytes_old,
1204                               nbytes_new,
1205                               p)
1206             # re-acquire pointer to location in memory
1207             # as it might have moved
1208             p = pysam_bam1_seq( src )
1209             for k from 0 <= k < nbytes_new: p[k] = 0
1210             # convert to C string
1211             s = seq
1212             for k from 0 <= k < l:
1213                 p[k/2] |= pysam_translate_sequence(s[k]) << 4 * (1 - k % 2)
1214
1215             # erase qualities
1216             p = pysam_bam1_qual( src )
1217             p[0] = 0xff
1218
1219     property qual:
1220         """the base quality (None if not present)"""
1221         def __get__(self):
1222             cdef bam1_t * src
1223             cdef uint8_t * p
1224             cdef char * q
1225             src = self._delegate
1226             if src.core.l_qseq == 0: return None
1227
1228             p = pysam_bam1_qual( src )
1229             if p[0] == 0xff: return None
1230
1231             q = < char *>calloc(src.core.l_qseq + 1 , sizeof(char))
1232             for k from 0 <= k < src.core.l_qseq:
1233             ## equivalent to t[i] + 33 (see bam.c)
1234                 q[k] = p[k] + 33
1235             # convert to python string
1236             retval=q
1237             # clean up
1238             free(q)
1239             return retval
1240
1241         def __set__(self,qual):
1242             # note that space is already allocated via the sequences
1243             cdef bam1_t * src
1244             cdef uint8_t * p
1245             cdef char * q
1246             src = self._delegate
1247             p = pysam_bam1_qual( src )
1248             if qual == None or len(qual) == 0:
1249                 # if absent - set to 0xff
1250                 p[0] = 0xff
1251                 return
1252             cdef int l
1253             # convert to C string
1254             q = qual
1255             l = len(qual)
1256             if src.core.l_qseq != l:
1257                 raise ValueError("quality and sequence mismatch: %i != %i" % (l, src.core.l_qseq))
1258             assert src.core.l_qseq == l
1259             for k from 0 <= k < l:
1260                 p[k] = <uint8_t>q[k] - 33
1261
1262     property tags:
1263         """the tags in the AUX field."""
1264         def __get__(self):
1265             cdef char * ctag
1266             cdef bam1_t * src
1267             cdef uint8_t * s
1268             cdef char tpe
1269
1270             src = self._delegate
1271             if src.l_aux == 0: return None
1272
1273             s = pysam_bam1_aux( src )
1274             result = []
1275             ctag = <char*>calloc( 3, sizeof(char) )
1276             cdef int x
1277             while s < (src.data + src.data_len):
1278                 # get tag
1279                 ctag[0] = s[0]
1280                 ctag[1] = s[1]
1281                 pytag = ctag
1282
1283                 s += 2
1284
1285                 # convert type - is there a better way?
1286                 ctag[0] = s[0]
1287                 ctag[1] = 0
1288                 pytype = ctag
1289                 # get type and value
1290                 # how do I do char literal comparison in cython?
1291                 # the code below works (i.e, is C comparison)
1292                 tpe = toupper(s[0])
1293                 if tpe == 'S'[0]:
1294                     value = <int>bam_aux2i(s)
1295                     s += 2
1296                 elif tpe == 'I'[0]:
1297                     value = <int>bam_aux2i(s)
1298                     s += 4
1299                 elif tpe == 'F'[0]:
1300                     value = <float>bam_aux2f(s)
1301                     s += 4
1302                 elif tpe == 'D'[0]:
1303                     value = <double>bam_aux2d(s)
1304                     s += 8
1305                 elif tpe == 'C'[0]:
1306                     value = <int>bam_aux2i(s)
1307                     s += 1
1308                 elif tpe == 'A'[0]:
1309                     # there might a more efficient way
1310                     # to convert a char into a string
1311                     value = "%c" % <char>bam_aux2A(s)
1312                     s += 1
1313                 elif tpe == 'Z'[0]:
1314                     value = <char*>bam_aux2Z(s)
1315                     # +1 for NULL terminated string
1316                     s += len(value) + 1
1317
1318                 # skip over type
1319                 s += 1
1320
1321                 # ignore pytype
1322                 result.append( (pytag, value) )
1323
1324             free( ctag )
1325             return result
1326
1327         def __set__(self, tags):
1328             cdef char * ctag
1329             cdef bam1_t * src
1330             cdef uint8_t * s
1331             cdef uint8_t * new_data
1332             cdef int guessed_size, control_size
1333             src = self._delegate
1334             cdef int max_size, size
1335             max_size = 4000
1336
1337             # map samtools code to python.struct code and byte size
1338             buffer = ctypes.create_string_buffer(max_size)
1339
1340             offset = 0
1341             for pytag, value in tags:
1342                 t = type(value)
1343                 if t == types.FloatType:
1344                     fmt = "<cccf"
1345                 elif t == types.IntType:
1346                     if value < 0:
1347                         if value >= -127: fmt, pytype = "<cccb", 'c'
1348                         elif value >= -32767: fmt, pytype = "<ccch", 's'
1349                         elif value < -2147483648: raise ValueError( "integer %i out of range of BAM/SAM specification" % value )
1350                         else: fmt, ctype = "<ccci", 'i'[0]
1351                     else:
1352                         if value <= 255: fmt, pytype = "<cccB", 'C'
1353                         elif value <= 65535: fmt, pytype = "<cccH", 'S'
1354                         elif value > 4294967295: raise ValueError( "integer %i out of range of BAM/SAM specification" % value )
1355                         else: fmt, pytype = "<cccI", 'I'
1356                 else:
1357                     # Note: hex strings (H) are not supported yet
1358                     if len(value) == 1:
1359                         fmt, pytype = "<cccc", 'A'
1360                     else:
1361                         fmt, pytype = "<ccc%is" % (len(value)+1), 'Z'
1362
1363                 size = struct.calcsize(fmt)
1364                 if offset + size > max_size:
1365                     raise NotImplementedError("tags field too large")
1366
1367                 struct.pack_into( fmt,
1368                                   buffer,
1369                                   offset,
1370                                   pytag[0],
1371                                   pytag[1],
1372                                   pytype,
1373                                   value )
1374                 offset += size
1375
1376             # delete the old data and allocate new
1377             pysam_bam_update( src,
1378                               src.l_aux,
1379                               offset,
1380                               pysam_bam1_aux( src ) )
1381
1382             src.l_aux = offset
1383
1384             if offset == 0: return
1385
1386             # get location of new data
1387             s = pysam_bam1_aux( src )
1388
1389             # check if there is direct path from buffer.raw to tmp
1390             cdef char * temp
1391             temp = buffer.raw
1392             memcpy( s, temp, offset )
1393
1394     property flag:
1395         """properties flag"""
1396         def __get__(self): return self._delegate.core.flag
1397         def __set__(self, flag): self._delegate.core.flag = flag
1398     property rname:
1399         """
1400         :term:`target` ID
1401
1402         .. note::
1403
1404             This field contains the index of the reference sequence
1405             in the sequence dictionary. To obtain the name
1406             of the reference sequence, use :meth:`pysam.Samfile.getrname()`
1407
1408         """
1409         def __get__(self): return self._delegate.core.tid
1410         def __set__(self, tid): self._delegate.core.tid = tid
1411     property pos:
1412         """0-based leftmost coordinate"""
1413         def __get__(self): return self._delegate.core.pos
1414         def __set__(self, pos):
1415             ## setting the cigar string also updates the "bin" attribute
1416             cdef bam1_t * src
1417             src = self._delegate
1418             if src.core.n_cigar:
1419                 src.core.bin = bam_reg2bin( src.core.pos, bam_calend( &src.core, pysam_bam1_cigar(src)) )
1420             else:
1421                 src.core.bin = bam_reg2bin( src.core.pos, src.core.pos + 1)
1422             self._delegate.core.pos = pos
1423     property bin:
1424         """properties bin"""
1425         def __get__(self): return self._delegate.core.bin
1426         def __set__(self, bin): self._delegate.core.bin = bin
1427     property rlen:
1428         '''length of the read (read only). Returns 0 if not given.'''
1429         def __get__(self): return self._delegate.core.l_qseq
1430     property mapq:
1431         """mapping quality"""
1432         def __get__(self): return self._delegate.core.qual
1433         def __set__(self, qual): self._delegate.core.qual = qual
1434     property mrnm:
1435         """the :term:`reference` id of the mate """
1436         def __get__(self): return self._delegate.core.mtid
1437         def __set__(self, mtid): self._delegate.core.mtid = mtid
1438     property mpos:
1439         """the position of the mate"""
1440         def __get__(self): return self._delegate.core.mpos
1441         def __set__(self, mpos): self._delegate.core.mpos = mpos
1442     property isize:
1443         """the insert size"""
1444         def __get__(self): return self._delegate.core.isize
1445         def __set__(self, isize): self._delegate.core.isize = isize
1446     property is_paired:
1447         """true if read is paired in sequencing"""
1448         def __get__(self): return (self._delegate.core.flag & BAM_FPAIRED) != 0
1449         def __set__(self,val):
1450             if val: self._delegate.core.flag |= BAM_FPAIRED
1451             else: self._delegate.core.flag &= ~BAM_FPAIRED
1452     property is_proper_pair:
1453         """true if read is mapped in a proper pair"""
1454         def __get__(self): return (self.flag & BAM_FPROPER_PAIR) != 0
1455         def __set__(self,val):
1456             if val: self._delegate.core.flag |= BAM_FPROPER_PAIR
1457             else: self._delegate.core.flag &= ~BAM_FPROPER_PAIR
1458     property is_unmapped:
1459         """true if read itself is unmapped"""
1460         def __get__(self): return (self.flag & BAM_FUNMAP) != 0
1461         def __set__(self,val):
1462             if val: self._delegate.core.flag |= BAM_FUNMAP
1463             else: self._delegate.core.flag &= ~BAM_FUNMAP
1464     property mate_is_unmapped:
1465         """true if the mate is unmapped"""
1466         def __get__(self): return (self.flag & BAM_FMUNMAP) != 0
1467         def __set__(self,val):
1468             if val: self._delegate.core.flag |= BAM_FMUNMAP
1469             else: self._delegate.core.flag &= ~BAM_FMUNMAP
1470     property is_reverse:
1471         """true if read is mapped to reverse strand"""
1472         def __get__(self):return (self.flag & BAM_FREVERSE) != 0
1473         def __set__(self,val):
1474             if val: self._delegate.core.flag |= BAM_FREVERSE
1475             else: self._delegate.core.flag &= ~BAM_FREVERSE
1476     property mate_is_reverse:
1477         """true is read is mapped to reverse strand"""
1478         def __get__(self): return (self.flag & BAM_FMREVERSE) != 0
1479         def __set__(self,val):
1480             if val: self._delegate.core.flag |= BAM_FMREVERSE
1481             else: self._delegate.core.flag &= ~BAM_FMREVERSE
1482     property is_read1:
1483         """true if this is read1"""
1484         def __get__(self): return (self.flag & BAM_FREAD1) != 0
1485         def __set__(self,val):
1486             if val: self._delegate.core.flag |= BAM_FREAD1
1487             else: self._delegate.core.flag &= ~BAM_FREAD1
1488     property is_read2:
1489         """true if this is read2"""
1490         def __get__(self): return (self.flag & BAM_FREAD2) != 0
1491         def __set__(self,val):
1492             if val: self._delegate.core.flag |= BAM_FREAD2
1493             else: self._delegate.core.flag &= ~BAM_FREAD2
1494     property is_secondary:
1495         """true if not primary alignment"""
1496         def __get__(self): return (self.flag & BAM_FSECONDARY) != 0
1497         def __set__(self,val):
1498             if val: self._delegate.core.flag |= BAM_FSECONDARY
1499             else: self._delegate.core.flag &= ~BAM_FSECONDARY
1500     property is_qcfail:
1501         """true if QC failure"""
1502         def __get__(self): return (self.flag & BAM_FQCFAIL) != 0
1503         def __set__(self,val):
1504             if val: self._delegate.core.flag |= BAM_FQCFAIL
1505             else: self._delegate.core.flag &= ~BAM_FQCFAIL
1506     property is_duplicate:
1507         """ true if optical or PCR duplicate"""
1508         def __get__(self): return (self.flag & BAM_FDUP) != 0
1509         def __set__(self,val):
1510             if val: self._delegate.core.flag |= BAM_FDUP
1511             else: self._delegate.core.flag &= ~BAM_FDUP
1512
1513     def opt(self, tag):
1514         """retrieves optional data given a two-letter *tag*"""
1515         #see bam_aux.c: bam_aux_get() and bam_aux2i() etc
1516         cdef uint8_t * v
1517         v = bam_aux_get(self._delegate, tag)
1518         if v == NULL: raise KeyError( "tag '%s' not present" % tag )
1519         type = chr(v[0])
1520         if type == 'c' or type == 'C' or type == 's' or type == 'S' or type == 'i':
1521             return <int>bam_aux2i(v)
1522         elif type == 'f':
1523             return <float>bam_aux2f(v)
1524         elif type == 'd':
1525             return <double>bam_aux2d(v)
1526         elif type == 'A':
1527             # there might a more efficient way
1528             # to convert a char into a string
1529             return '%c' % <char>bam_aux2A(v)
1530         elif type == 'Z':
1531             return <char*>bam_aux2Z(v)
1532
1533     def fancy_str (self):
1534         """returns list of fieldnames/values in pretty format for debugging
1535         """
1536         ret_string = []
1537         field_names = {
1538            "tid":           "Contig index",
1539            "pos":           "Mapped position on contig",
1540            "mtid":          "Contig index for mate pair",
1541            "mpos":          "Position of mate pair",
1542            "isize":         "Insert size",
1543            "flag":          "Binary flag",
1544            "n_cigar":       "Count of cigar entries",
1545            "cigar":         "Cigar entries",
1546            "qual":          "Mapping quality",
1547            "bin":           "Bam index bin number",
1548            "l_qname":       "Length of query name",
1549            "qname":         "Query name",
1550            "l_qseq":        "Length of query sequence",
1551            "qseq":          "Query sequence",
1552            "bqual":         "Quality scores",
1553            "l_aux":         "Length of auxilary data",
1554            "m_data":        "Maximum data length",
1555            "data_len":      "Current data length",
1556            }
1557         fields_names_in_order = ["tid", "pos", "mtid", "mpos", "isize", "flag",
1558                                  "n_cigar", "cigar", "qual", "bin", "l_qname", "qname",
1559                                  "l_qseq", "qseq", "bqual", "l_aux", "m_data", "data_len"]
1560
1561         for f in fields_names_in_order:
1562             if not f in self.__dict__:
1563                 continue
1564             ret_string.append("%-30s %-10s= %s" % (field_names[f], "(" + f + ")", self.__getattribute__(f)))
1565
1566         for f in self.__dict__:
1567             if not f in field_names:
1568                 ret_string.append("%-30s %-10s= %s" % (f, "", self.__getattribute__(f)))
1569         return ret_string
1570
1571 cdef class PileupProxy:
1572     '''A pileup column. A pileup column contains
1573     all the reads that map to a certain target base.
1574
1575     tid
1576         chromosome ID as is defined in the header
1577     pos
1578         the target base coordinate (0-based)
1579     n
1580         number of reads mapping to this column
1581     pileups
1582         list of reads (:class:`pysam.PileupRead`) aligned to this column
1583
1584     This class is a proxy for results returned by the samtools pileup engine.
1585     If the underlying engine iterator advances, the results of this column
1586     will change.
1587     '''
1588     cdef bam_plbuf_t * buf
1589     cdef int n_pu
1590
1591     def __cinit__(self ):
1592         pass
1593
1594     def __str__(self):
1595         return "\t".join( map(str, (self.tid, self.pos, self.n))) +\
1596             "\n" +\
1597             "\n".join( map(str, self.pileups) )
1598
1599     property tid:
1600         '''the chromosome ID as is defined in the header'''
1601         def __get__(self): return pysam_get_tid( self.buf )
1602
1603     property n:
1604         '''number of reads mapping to this column.'''
1605         def __get__(self): return self.n_pu
1606         def __set__(self, n): self.n_pu = n
1607
1608     property pos:
1609         def __get__(self): return pysam_get_pos( self.buf )
1610
1611     property pileups:
1612         '''list of reads (:class:`pysam.PileupRead`) aligned to this column'''
1613         def __get__(self):
1614             cdef bam_pileup1_t * pl
1615             pl = pysam_get_pileup( self.buf )
1616             pileups = []
1617             # warning: there could be problems if self.n and self.buf are
1618             # out of sync.
1619             for x from 0 <= x < self.n_pu:
1620                 pileups.append( makePileupRead( &pl[x]) )
1621             return pileups
1622
1623 cdef class PileupRead:
1624     '''A read aligned to a column.
1625     '''
1626
1627     cdef:
1628          AlignedRead _alignment
1629          int32_t  _qpos
1630          int _indel
1631          int _level
1632          uint32_t _is_del
1633          uint32_t _is_head
1634          uint32_t _is_tail
1635
1636     def __cinit__( self ):
1637         pass
1638
1639     def __str__(self):
1640         return "\t".join( map(str, (self.alignment, self.qpos, self.indel, self.level, self.is_del, self.is_head, self.is_tail ) ) )
1641
1642     property alignment:
1643         """a :class:`pysam.AlignedRead` object of the aligned read"""
1644         def __get__(self):
1645             return self._alignment
1646     property qpos:
1647         """position of the read base at the pileup site, 0-based"""
1648         def __get__(self):
1649             return self._qpos
1650     property indel:
1651         """indel length; 0 for no indel, positive for ins and negative for del"""
1652         def __get__(self):
1653             return self._indel
1654     property is_del:
1655         """1 iff the base on the padded read is a deletion"""
1656         def __get__(self):
1657             return self._is_del
1658     property is_head:
1659         def __get__(self):
1660             return self._is_head
1661     property is_tail:
1662         def __get__(self):
1663             return self._is_tail
1664     property level:
1665         def __get__(self):
1666             return self._level
1667
1668 class Outs:
1669     '''http://mail.python.org/pipermail/python-list/2000-June/038406.html'''
1670     def __init__(self, id = 1):
1671         self.streams = []
1672         self.id = id
1673
1674     def setdevice(self, filename):
1675         '''open an existing file, like "/dev/null"'''
1676         fd = os.open(filename, os.O_WRONLY)
1677         self.setfd(fd)
1678
1679     def setfile(self, filename):
1680         '''open a new file.'''
1681         fd = os.open(filename, os.O_WRONLY|os.O_CREAT, 0660);
1682         self.setfd(fd)
1683
1684     def setfd(self, fd):
1685         ofd = os.dup(self.id)      #  Save old stream on new unit.
1686         self.streams.append(ofd)
1687         sys.stdout.flush()          #  Buffered data goes to old stream.
1688         os.dup2(fd, self.id)        #  Open unit 1 on new stream.
1689         os.close(fd)                #  Close other unit (look out, caller.)
1690
1691     def restore(self):
1692         '''restore previous output stream'''
1693         if self.streams:
1694             # the following was not sufficient, hence flush both stderr and stdout
1695             # os.fsync( self.id )
1696             sys.stdout.flush()
1697             sys.stderr.flush()
1698             os.dup2(self.streams[-1], self.id)
1699             os.close(self.streams[-1])
1700             del self.streams[-1]
1701
1702 def _samtools_dispatch( method, args = () ):
1703     '''call ``method`` in samtools providing arguments in args.
1704
1705     .. note::
1706        This method redirects stdout and stderr to capture it
1707        from samtools. If for some reason stdout/stderr disappears
1708        the reason might be in this method.
1709
1710     .. note::
1711        The current implementation might only work on linux.
1712
1713     .. note::
1714        This method captures stdout and stderr using temporary files,
1715        which are then read into memory in their entirety. This method
1716        is slow and might cause large memory overhead.
1717
1718     See http://bytes.com/topic/c/answers/487231-how-capture-stdout-temporarily
1719     on the topic of redirecting stderr/stdout.
1720     '''
1721
1722     # note that debugging this module can be a problem
1723     # as stdout/stderr will not appear
1724
1725     # redirect stderr and stdout to file
1726
1727     # open files and redirect into it
1728     stderr_h, stderr_f = tempfile.mkstemp()
1729     stdout_h, stdout_f = tempfile.mkstemp()
1730
1731     # patch for `samtools view`
1732     # samtools `view` closes stdout, from which I can not
1733     # recover. Thus redirect output to file with -o option.
1734     if method == "view":
1735         if "-o" in args: raise ValueError("option -o is forbidden in samtools view")
1736         args = ( "-o", stdout_f ) + args
1737
1738     stdout_save = Outs( sys.stdout.fileno() )
1739     stdout_save.setfd( stdout_h )
1740     stderr_save = Outs( sys.stderr.fileno() )
1741     stderr_save.setfd( stderr_h )
1742
1743     # do the function call to samtools
1744     cdef char ** cargs
1745     cdef int i, n, retval
1746
1747     n = len(args)
1748     # allocate two more for first (dummy) argument (contains command)
1749     cargs = <char**>calloc( n+2, sizeof( char *) )
1750     cargs[0] = "samtools"
1751     cargs[1] = method
1752     for i from 0 <= i < n: cargs[i+2] = args[i]
1753     retval = pysam_dispatch(n+2, cargs)
1754     free( cargs )
1755
1756     # restore stdout/stderr. This will also flush, so
1757     # needs to be before reading back the file contents
1758     stdout_save.restore()
1759     stderr_save.restore()
1760
1761     # capture stderr/stdout.
1762     out_stderr = open( stderr_f, "r").readlines()
1763     out_stdout = open( stdout_f, "r").readlines()
1764
1765     # clean up files
1766     os.remove( stderr_f )
1767     os.remove( stdout_f )
1768
1769     return retval, out_stderr, out_stdout
1770
1771 __all__ = ["Samfile",
1772            "Fastafile",
1773            "IteratorRow",
1774            "IteratorRowAll",
1775            "IteratorColumn",
1776            "AlignedRead",
1777            "PileupColumn",
1778            "PileupProxy",
1779            "PileupRead" ]
1780
1781
1782