pysam/csamtools.pyx

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