pysam/ctabix.pyx

   1 # cython: embedsignature=True
   2 # adds doc-strings for sphinx
   3
   4 import tempfile, os, sys, types, itertools, struct, ctypes, gzip
   5 from cpython cimport PyString_FromStringAndSize, PyString_AS_STRING
   6 cimport TabProxies
   7
   8 cdef class Tabixfile:
   9     '''*(filename, mode='r')*
  10
  11     opens a :term:`tabix file` for reading. A missing
  12     index (*filename* + ".tbi") will raise an exception.
  13     '''
  14
  15     def __cinit__(self, *args, **kwargs ):
  16         self.tabixfile = NULL
  17         self._open( *args, **kwargs )
  18
  19     def _isOpen( self ):
  20         '''return true if samfile has been opened.'''
  21         return self.tabixfile != NULL
  22
  23     def _open( self,
  24                char * filename,
  25                mode ='r',
  26               ):
  27         '''open a :term:`tabix file` for reading.
  28         '''
  29
  30         assert mode in ( "r",), "invalid file opening mode `%s`" % mode
  31
  32         # close a previously opened file
  33         if self.tabixfile != NULL: self.close()
  34         self.tabixfile = NULL
  35
  36         if self._filename != NULL: free(self._filename )
  37         self._filename = strdup( filename )
  38
  39         filename_index = filename + ".tbi"
  40
  41         if mode[0] == 'w':
  42             # open file for writing
  43             pass
  44
  45         elif mode[0] == "r":
  46             # open file for reading
  47             if not os.path.exists( self._filename ):
  48                 raise IOError( "file `%s` not found" % self._filename)
  49
  50             if not os.path.exists( filename_index ):
  51                 raise IOError( "index `%s` not found" % filename_index)
  52
  53             # open file and load index
  54             self.tabixfile = ti_open( self._filename, filename_index )
  55
  56         if self.tabixfile == NULL:
  57             raise IOError("could not open file `%s`" % filename )
  58
  59     def _parseRegion( self,
  60                       reference = None,
  61                       start = None,
  62                       end = None,
  63                       region = None ):
  64         '''parse region information.
  65
  66         raise ValueError for for invalid regions.
  67
  68         returns a tuple of region, tid, start and end. Region
  69         is a valid samtools :term:`region` or None if the region
  70         extends over the whole file.
  71
  72         Note that regions are 1-based, while start,end are python coordinates.
  73         '''
  74         ti_lazy_index_load( self.tabixfile )
  75
  76         cdef int rtid
  77         cdef int rstart
  78         cdef int rend
  79         cdef int max_pos
  80         max_pos = 2 << 29
  81
  82         rtid = rstart = rend = 0
  83
  84         # translate to a region
  85         if reference:
  86             if start != None and end != None:
  87                 region = "%s:%i-%i" % (reference, start+1, end)
  88             elif start == None and end != None:
  89                 region = "%s:%i-%i" % (reference, 1, end)
  90             elif end == None and start != None:
  91                 region = "%s:%i-%i" % (reference, start+1, max_pos-1)
  92             else:
  93                 region = reference
  94
  95         if region:
  96             ti_parse_region( self.tabixfile.idx, region, &rtid, &rstart, &rend)
  97             if rtid < 0: raise ValueError( "invalid region `%s`" % region )
  98             if rstart > rend: raise ValueError( 'invalid region: start (%i) > end (%i)' % (rstart, rend) )
  99             if not 0 <= rstart < max_pos: raise ValueError( 'start out of range (%i)' % rstart )
 100             if not 0 <= rend < max_pos: raise ValueError( 'end out of range (%i)' % rend )
 101
 102         return region, rtid, rstart, rend
 103
 104     def fetch( self,
 105                reference = None,
 106                start = None,
 107                end = None,
 108                region = None,
 109                parser = None ):
 110         '''
 111
 112         fetch one or more rows in a :term:`region` using 0-based indexing. The region is specified by
 113         :term:`reference`, *start* and *end*. Alternatively, a samtools :term:`region` string can be supplied.
 114
 115         Without *reference* or *region* all entries will be fetched.
 116
 117         If only *reference* is set, all reads matching on *reference* will be fetched.
 118
 119         If *parser* is None, the results are returned as an unparsed string.
 120         Otherwise, *parser* is assumed to be a functor that will return parsed
 121         data (see for example :meth:`asTuple` and :meth:`asGTF`).
 122         '''
 123         ti_lazy_index_load( self.tabixfile )
 124
 125         if not self._isOpen():
 126             raise ValueError( "I/O operation on closed file" )
 127
 128         region, rtid, rstart, rend = self._parseRegion( reference, start, end, region )
 129
 130         if parser == None:
 131             if region:
 132                 return TabixIterator( self, rtid, rstart, rend )
 133             else:
 134                 return TabixIterator( self, -1, 0, 0 )
 135         else:
 136             if region:
 137                 return TabixIteratorParsed( self, rtid, rstart, rend, parser )
 138             else:
 139                 return TabixIteratorParsed( self, -1, 0, 0, parser )
 140
 141     ###############################################################
 142     ###############################################################
 143     ###############################################################
 144     ## properties
 145     ###############################################################
 146     property filename:
 147         '''filename associated with this object.'''
 148         def __get__(self):
 149             if not self._isOpen(): raise ValueError( "I/O operation on closed file" )
 150             return self._filename
 151
 152     property header:
 153         '''the file header.
 154
 155         .. note::
 156             The header is returned as an iterator over lines without the
 157             newline character.
 158         '''
 159
 160         def __get__( self ):
 161             return TabixHeaderIterator( self )
 162
 163     property contigs:
 164         '''chromosome names'''
 165         def __get__(self):
 166             cdef char ** sequences
 167             cdef int nsequences
 168
 169             ti_lazy_index_load( self.tabixfile )
 170             sequences = ti_seqname( self.tabixfile.idx, &nsequences )
 171             cdef int x
 172             result = []
 173             for x from 0 <= x < nsequences:
 174                 result.append( sequences[x] )
 175             return result
 176
 177     def close( self ):
 178         '''
 179         closes the :class:`pysam.Tabixfile`.'''
 180         if self.tabixfile != NULL:
 181             ti_close( self.tabixfile )
 182             self.tabixfile = NULL
 183
 184     def __dealloc__( self ):
 185         # remember: dealloc cannot call other python methods
 186         # note: no doc string
 187         # note: __del__ is not called.
 188         if self.tabixfile != NULL:
 189             ti_close( self.tabixfile )
 190             self.tabixfile = NULL
 191         if self._filename != NULL: free( self._filename )
 192
 193 cdef class TabixIterator:
 194     """iterates over rows in *tabixfile* in region
 195     given by *tid*, *start* and *end*.
 196     """
 197
 198     cdef ti_iter_t iterator
 199     cdef tabix_t * tabixfile
 200
 201     def __cinit__(self, Tabixfile tabixfile,
 202                   int tid, int start, int end ):
 203
 204         assert tabixfile._isOpen()
 205
 206         # makes sure that samfile stays alive as long as the
 207         # iterator is alive.
 208         self.tabixfile = tabixfile.tabixfile
 209
 210         if tid < 0:
 211             # seek to start of file to ensure iteration is over
 212             # all entries.
 213             bgzf_seek( self.tabixfile.fp, 0, 0)
 214             self.iterator = ti_iter_first()
 215         else:
 216             self.iterator = ti_queryi(self.tabixfile, tid, start, end)
 217
 218         if <void*>self.iterator == NULL:
 219             raise ValueError("malformatted query or wrong sequence name.\n")
 220
 221     def __iter__(self):
 222         return self
 223
 224     def __next__(self):
 225         """python version of next().
 226
 227         pyrex uses this non-standard name instead of next()
 228         """
 229
 230         cdef char * s
 231         cdef int len
 232         # metachar filtering does not work within tabix
 233         # though it should. Getting the metachar is a pain
 234         # as ti_index_t is incomplete type.
 235
 236         # simply use '#' for now.
 237         while 1:
 238             s = ti_read(self.tabixfile, self.iterator, &len)
 239             if s == NULL: raise StopIteration
 240             if s[0] != '#': break
 241
 242         return s
 243
 244     def __dealloc__(self):
 245         if <void*>self.iterator != NULL:
 246             ti_iter_destroy(self.iterator)
 247
 248 cdef class TabixHeaderIterator:
 249     """return header lines.
 250     """
 251
 252     cdef ti_iter_t iterator
 253     cdef tabix_t * tabixfile
 254
 255     def __cinit__(self, Tabixfile tabixfile ):
 256
 257         assert tabixfile._isOpen()
 258
 259         # makes sure that samfile stays alive as long as the
 260         # iterator is alive.
 261         self.tabixfile = tabixfile.tabixfile
 262
 263         self.iterator = ti_query(self.tabixfile, NULL, 0, 0)
 264
 265         if <void*>self.iterator == NULL:
 266             raise ValueError("can't open header.\n")
 267
 268     def __iter__(self):
 269         return self
 270
 271     def __next__(self):
 272         """python version of next().
 273
 274         pyrex uses this non-standard name instead of next()
 275         """
 276
 277         cdef char * s
 278         cdef int len
 279
 280         # Getting the metachar is a pain as ti_index_t is incomplete type.
 281         # simply use '#' for now.
 282         s = ti_read(self.tabixfile, self.iterator, &len)
 283         if s == NULL: raise StopIteration
 284         # stop at first non-header line
 285         if s[0] != '#': raise StopIteration
 286
 287         return s
 288
 289     def __dealloc__(self):
 290         if <void*>self.iterator != NULL:
 291             ti_iter_destroy(self.iterator)
 292
 293 #########################################################
 294 #########################################################
 295 #########################################################
 296 cdef class Parser:
 297     pass
 298
 299 cdef class asTuple(Parser):
 300     '''converts a :term:`tabix row` into a python tuple.
 301
 302     Access is by numeric index.
 303     '''
 304     def __call__(self, char * buffer, int len):
 305         cdef TabProxies.TupleProxy r
 306         r = TabProxies.TupleProxy()
 307         # need to copy - there were some
 308         # persistence issues with "present"
 309         r.copy( buffer, len )
 310         return r
 311
 312 cdef class asGTF(Parser):
 313     '''converts a :term:`tabix row` into a GTF record with the following
 314     fields:
 315
 316     contig
 317        contig
 318     feature
 319        feature
 320     source
 321        source
 322     start
 323        genomic start coordinate (0-based)
 324     end
 325        genomic end coordinate plus one (0-based)
 326     score
 327        feature score
 328     strand
 329        strand
 330     frame
 331        frame
 332     attributes
 333        attribute string.
 334
 335     GTF formatted entries also defined the attributes:
 336
 337     gene_id
 338        the gene identifier
 339     transcript_ind
 340        the transcript identifier
 341
 342     '''
 343     def __call__(self, char * buffer, int len):
 344         cdef TabProxies.GTFProxy r
 345         r = TabProxies.GTFProxy()
 346         r.copy( buffer, len )
 347         return r
 348
 349 cdef class asBed( Parser ):
 350     '''converts a :term:`tabix row` into a bed record
 351     with the following fields:
 352
 353     contig
 354        contig
 355     start
 356        genomic start coordinate (zero-based)
 357     end
 358        genomic end coordinate plus one (zero-based)
 359     name
 360        name of feature.
 361     score
 362        score of feature
 363     strand
 364        strand of feature
 365     thickStart
 366        thickStart
 367     thickEnd
 368        thickEnd
 369     itemRGB
 370        itemRGB
 371     blockCount
 372        number of bocks
 373     blockSizes
 374        ',' separated string of block sizes
 375     blockStarts
 376        ',' separated string of block genomic start positions
 377
 378     Only the first three fields are required. Additional
 379     fields are optional, but if one is defined, all the preceeding
 380     need to be defined as well.
 381
 382     '''
 383     def __call__(self, char * buffer, int len):
 384         cdef TabProxies.BedProxy r
 385         r = TabProxies.BedProxy()
 386         r.copy( buffer, len )
 387         return r
 388
 389 cdef class asVCF( Parser ):
 390     '''converts a :term:`tabix row` into a VCF record with
 391     the following fields:
 392
 393     contig
 394        contig
 395     pos
 396        chromosomal position, zero-based
 397     id
 398        id
 399     ref
 400        reference
 401     alt
 402        alt
 403     qual
 404        qual
 405     filter
 406        filter
 407     info
 408        info
 409     format
 410        format specifier.
 411
 412     Access to genotypes is via index::
 413
 414         contig = vcf.contig
 415         first_sample_genotype = vcf[0]
 416         second_sample_genotype = vcf[1]
 417
 418     '''
 419     def __call__(self, char * buffer, int len ):
 420         cdef TabProxies.VCFProxy r
 421         r = TabProxies.VCFProxy()
 422         r.copy( buffer, len )
 423         return r
 424
 425 #########################################################
 426 #########################################################
 427 #########################################################
 428 cdef class TabixIteratorParsed:
 429     """iterates over mapped reads in a region.
 430
 431     Returns parsed data.
 432     """
 433
 434     cdef ti_iter_t iterator
 435     cdef tabix_t * tabixfile
 436     cdef Parser parser
 437
 438     def __cinit__(self,
 439                   Tabixfile tabixfile,
 440                   int tid,
 441                   int start,
 442                   int end,
 443                   Parser parser ):
 444
 445         assert tabixfile._isOpen()
 446         self.parser = parser
 447
 448         # makes sure that samfile stays alive as long as the
 449         # iterator is alive.
 450         self.tabixfile = tabixfile.tabixfile
 451
 452         if tid < 0:
 453             # seek to start of file to ensure iteration is over
 454             # all entries.
 455             bgzf_seek( self.tabixfile.fp, 0, 0)
 456             self.iterator = ti_iter_first()
 457         else:
 458             self.iterator = ti_queryi(self.tabixfile, tid, start, end)
 459
 460         if <void*>self.iterator == NULL:
 461             raise ValueError("malformatted query or wrong sequence name.\n")
 462
 463     def __iter__(self):
 464         return self
 465
 466     def __next__(self):
 467         """python version of next().
 468
 469         pyrex uses this non-standard name instead of next()
 470         """
 471
 472         cdef char * s
 473         cdef int len
 474         while 1:
 475             s = ti_read(self.tabixfile, self.iterator, &len)
 476             if s == NULL: raise StopIteration
 477             # todo: read metachar from configuration
 478             if s[0] != '#': break
 479
 480         return self.parser(s, len)
 481
 482     def __dealloc__(self):
 483         if <void*>self.iterator != NULL:
 484             ti_iter_destroy(self.iterator)
 485
 486 def tabix_compress( filename_in,
 487               filename_out,
 488               force = False ):
 489
 490     '''
 491     compress *filename_in* writing the output to *filename_out*.
 492
 493     Raise an IOError if *filename_out* already exists, unless *force* is set.
 494     '''
 495
 496     if not force and os.path.exists(filename_out ):
 497         raise IOError( "Filename '%s' already exists, use *force* to overwrite" % filename_out)
 498
 499     cdef int WINDOW_SIZE
 500     cdef int c, r
 501     cdef void * buffer
 502     cdef BGZF * fp
 503     cdef int fd_src
 504
 505     cdef int O_RDONLY
 506     O_RDONLY = os.O_RDONLY
 507
 508     WINDOW_SIZE = 64 * 1024
 509
 510     fp = bgzf_open( filename_out, "w")
 511     if fp == NULL:
 512         raise IOError( "could not open '%s' for writing" )
 513
 514     fd_src = open(filename_in, O_RDONLY)
 515     if fd_src == 0:
 516         raise IOError( "could not open '%s' for reading" )
 517
 518     buffer = malloc(WINDOW_SIZE)
 519
 520     while c > 0:
 521         c = read(fd_src, buffer, WINDOW_SIZE)
 522         r = bgzf_write(fp, buffer, c)
 523         if r < 0:
 524             free( buffer )
 525             raise OSError("writing failed")
 526
 527     free( buffer )
 528     r = bgzf_close(fp)
 529     if r < 0: raise OSError("writing failed")
 530
 531 def tabix_index( filename,
 532                  force = False,
 533                  seq_col = None,
 534                  start_col = None,
 535                  end_col = None,
 536                  preset = None,
 537                  meta_char = "#",
 538                  zerobased = False,
 539                 ):
 540     '''
 541     index tab-separated *filename* using tabix.
 542
 543     An existing index will not be overwritten unless
 544     *force* is set.
 545
 546     The index will be built from coordinates
 547     in columns *seq_col*, *start_col* and *end_col*.
 548
 549     The contents of *filename* have to be sorted by
 550     contig and position - the method does not check
 551     if the file is sorted.
 552
 553     Column indices are 0-based. Coordinates in the file
 554     are assumed to be 1-based.
 555
 556     If *preset* is provided, the column coordinates
 557     are taken from a preset. Valid values for preset
 558     are "gff", "bed", "sam", "vcf", psltbl", "pileup".
 559
 560     Lines beginning with *meta_char* and the first
 561     *line_skip* lines will be skipped.
 562
 563     If *filename* does not end in ".gz", it will be automatically
 564     compressed. The original file will be removed and only the
 565     compressed file will be retained.
 566
 567     If *filename* ends in *gz*, the file is assumed to be already
 568     compressed with bgzf.
 569
 570     returns the filename of the compressed data
 571     '''
 572
 573     if not os.path.exists(filename): raise IOError("No such file '%s'" % filename)
 574
 575     if not filename.endswith(".gz"):
 576
 577         tabix_compress( filename, filename + ".gz", force = force )
 578         os.unlink( filename )
 579         filename += ".gz"
 580
 581     if not force and os.path.exists(filename + ".tbi" ):
 582         raise IOError( "Filename '%s.tbi' already exists, use *force* to overwrite" )
 583
 584     # columns (1-based)
 585     # preset-code, contig, start, end, metachar for commends, lines to ignore at beginning
 586     # 0 is a missing column
 587     preset2conf = {
 588         'gff' : ( 0, 1, 4, 5, ord('#'), 0 ),
 589         'bed' : ( 0x10000, 1, 2, 3, ord('#'), 0 ),
 590         'psltbl' : ( 0x10000, 15, 17, 18, ord('#'), 0 ),
 591         'sam' : ( 1, 3, 4, 0, ord('#'), 0 ),
 592         'vcf' : ( 2, 1, 2, 0, ord('#'), 0 ),
 593         'pileup': (3, 1, 2, 0, ord('#'), 0 ),
 594         }
 595
 596     if preset:
 597         try:
 598             conf_data = preset2conf[preset]
 599         except KeyError:
 600             raise KeyError( "unknown preset '%s', valid presets are '%s'" % (preset, ",".join(preset2conf.keys() )))
 601     else:
 602         if end_col == None: end_col = -1
 603         preset = 0
 604
 605         # note that tabix internally works with 0-based coordinates and open/closed intervals.
 606         # When using a preset, conversion is automatically taken care of.
 607         # Otherwise, the coordinates are assumed to be 1-based closed intervals and
 608         # -1 is subtracted from the start coordinate. To avoid doing this, set
 609         # the TI_FLAG_UCSC=0x10000 flag:
 610         if zerobased: preset = preset | 0x10000
 611
 612         conf_data = (preset, seq_col+1, start_col+1, end_col+1, ord(meta_char), 0)
 613
 614     cdef ti_conf_t conf
 615     conf.preset, conf.sc, conf.bc, conf.ec, conf.meta_char, conf.line_skip = conf_data
 616
 617     ti_index_build( filename, &conf)
 618
 619     return filename
 620
 621 __all__ = ["tabix_index",
 622            "tabix_compress",
 623            "Tabixfile",
 624            "asTuple",
 625            "asGTF",
 626            "asVCF",
 627            "asBed",
 628            ]