pysam/TabProxies.pyx

   1 import types, sys
   2
   3 from cpython.version cimport PY_MAJOR_VERSION
   4
   5 from cpython cimport PyErr_SetString, PyBytes_Check, PyUnicode_Check, PyBytes_FromStringAndSize
   6
   7 cdef from_string_and_size(char* s, size_t length):
   8     if PY_MAJOR_VERSION < 3:
   9         return s[:length]
  10     else:
  11         return s[:length].decode("ascii")
  12
  13 # filename encoding (copied from lxml.etree.pyx)
  14 cdef str _FILENAME_ENCODING
  15 _FILENAME_ENCODING = sys.getfilesystemencoding()
  16 if _FILENAME_ENCODING is None:
  17     _FILENAME_ENCODING = sys.getdefaultencoding()
  18 if _FILENAME_ENCODING is None:
  19     _FILENAME_ENCODING = 'ascii'
  20
  21 cdef bytes _my_encodeFilename(object filename):
  22     u"""Make sure a filename is 8-bit encoded (or None).
  23     """
  24     if filename is None:
  25         return None
  26     elif PyBytes_Check(filename):
  27         return filename
  28     elif PyUnicode_Check(filename):
  29         return filename.encode(_FILENAME_ENCODING)
  30     else:
  31         raise TypeError, u"Argument must be string or unicode."
  32
  33 cdef bytes _force_bytes(object s):
  34     u"""convert string or unicode object to bytes, assuming ascii encoding.
  35     """
  36     if PY_MAJOR_VERSION < 3:
  37         return s
  38     elif s is None:
  39         return None
  40     elif PyBytes_Check(s):
  41         return s
  42     elif PyUnicode_Check(s):
  43         return s.encode('ascii')
  44     else:
  45         raise TypeError, u"Argument must be string, bytes or unicode."
  46
  47 cdef inline bytes _force_cmdline_bytes(object s):
  48     return _force_bytes(s)
  49
  50 cdef _charptr_to_str(char* s):
  51     if PY_MAJOR_VERSION < 3:
  52         return s
  53     else:
  54         return s.decode("ascii")
  55
  56 cdef _force_str(object s):
  57     """Return s converted to str type of current Python (bytes in Py2, unicode in Py3)"""
  58     if s is None:
  59         return None
  60     if PY_MAJOR_VERSION < 3:
  61         return s
  62     elif PyBytes_Check(s):
  63         return s.decode('ascii')
  64     else:
  65         # assume unicode
  66         return s
  67
  68 cdef char * nextItem( char * buffer ):
  69     cdef char * pos
  70     pos = strchr( buffer, '\t' )
  71     if pos == NULL: raise ValueError( "malformatted entry at %s" % buffer )
  72     pos[0] = '\0'
  73     pos += 1
  74     return pos
  75
  76 cdef char *StrOrEmpty( char * buffer ):
  77      if buffer == NULL: return ""
  78      else: return buffer
  79
  80 cdef int isNew( char * p, char * buffer, size_t nbytes ):
  81      if p == NULL: return 0
  82      return not (buffer <= p < buffer + nbytes )
  83
  84 cdef class TupleProxy:
  85     '''Proxy class for access to parsed row as a tuple.
  86
  87     This class represents a table row for fast read-access.
  88
  89     Access to individual fields is via the [] operator.
  90
  91     Only read-only access is implemented.
  92
  93     '''
  94
  95     def __cinit__(self ):
  96         self.data = NULL
  97         self.fields = NULL
  98         self.index = 0
  99         self.nbytes = 0
 100         self.is_modified = 0
 101         self.nfields = 0
 102         # start counting at field offset
 103         self.offset = 0
 104
 105     def __dealloc__(self):
 106         cdef int x
 107         if self.is_modified:
 108             for x from 0 <= x < self.nfields:
 109                 if isNew( self.fields[x], self.data, self.nbytes ):
 110                     free( self.fields[x] )
 111                     self.fields[x] = NULL
 112
 113         if self.data != NULL: free(self.data)
 114         if self.fields != NULL: free( self.fields )
 115
 116     cdef take( self, char * buffer, size_t nbytes ):
 117         '''start presenting buffer.
 118
 119         Take ownership of the pointer.
 120         '''
 121         self.data = buffer
 122         self.nbytes = nbytes
 123         self.update( buffer, nbytes )
 124
 125     cdef present( self, char * buffer, size_t nbytes ):
 126         '''start presenting buffer.
 127
 128         Do not take ownership of the pointer.
 129         '''
 130         self.update( buffer, nbytes )
 131
 132     cdef copy( self, char * buffer, size_t nbytes ):
 133         '''start presenting buffer of size *nbytes*.
 134
 135         Buffer is a '\0'-terminated string without the '\n'.
 136
 137         Take a copy of buffer.
 138         '''
 139         cdef int s
 140         # +1 for '\0'
 141         s = sizeof(char) *  (nbytes + 1)
 142         self.data = <char*>malloc( s )
 143         if self.data == NULL:
 144             raise ValueError("out of memory" )
 145         self.nbytes = nbytes
 146         memcpy( <char*>self.data, buffer, s )
 147         self.update( self.data, nbytes )
 148
 149     cdef int getMaxFields( self, size_t nbytes ):
 150         '''initialize fields.'''
 151         return nbytes / 2
 152
 153     cdef update( self, char * buffer, size_t nbytes ):
 154         '''update internal data.
 155
 156         *buffer* is a \0 terminated string.
 157
 158         *nbytes* is the number of bytes in buffer (excluding
 159         the \0)
 160
 161         Update starts work in buffer, thus can be used
 162         to collect any number of fields until nbytes
 163         is exhausted.
 164
 165         If max_fields is set, the number of fields is initialized to
 166         max_fields.
 167         '''
 168         cdef char * pos
 169         cdef char * old_pos
 170         cdef int field
 171         cdef int max_fields, x
 172
 173         assert strlen(buffer) == nbytes
 174
 175         if buffer[nbytes] != 0:
 176             raise ValueError( "incomplete line at %s" % buffer )
 177
 178         #################################
 179         # remove line breaks and feeds and update number of bytes
 180         x = nbytes - 1
 181         while x > 0 and (buffer[x] == '\n' or buffer[x] == '\r'):
 182             buffer[x] = '\0'
 183             x -= 1
 184         self.nbytes = x + 1
 185
 186         #################################
 187         # clear data
 188         if self.fields != NULL: free(self.fields)
 189
 190         for field from 0 <= field < self.nfields:
 191             if isNew( self.fields[field], self.data, self.nbytes ):
 192                 free( self.fields[field] )
 193
 194         self.is_modified = self.nfields = 0
 195
 196         #################################
 197         # allocate new
 198         max_fields = self.getMaxFields( nbytes )
 199         self.fields = <char **>calloc( max_fields, sizeof(char *) )
 200         if self.fields == NULL:
 201             raise ValueError("out of memory" )
 202
 203         #################################
 204         # start filling
 205         field = 0
 206         self.fields[field] = pos = buffer
 207         field += 1
 208         old_pos = pos
 209
 210         while 1:
 211
 212             pos = <char*>memchr( pos, '\t', nbytes )
 213             if pos == NULL: break
 214             pos[0] = '\0'
 215             pos += 1
 216             self.fields[field] = pos
 217             field += 1
 218             if field > max_fields:
 219                 raise ValueError("row too large - more than %i fields" % max_fields )
 220             nbytes -= pos - old_pos
 221             if nbytes < 0: break
 222             old_pos = pos
 223
 224         self.nfields = field
 225
 226     def _getindex( self, int index ):
 227         '''return item at idx index'''
 228         cdef int i = index
 229         if i < 0: i += self.nfields
 230         if i < 0: raise IndexError( "list index out of range" )
 231         i += self.offset
 232         if i >= self.nfields:
 233             raise IndexError( "list index out of range %i >= %i" % (i, self.nfields ))
 234         return self.fields[i]
 235
 236     def __getitem__( self, key ):
 237         if type(key) == int: return self._getindex( key )
 238         # slice object
 239         start, end, step = key.indices( self.nfields )
 240         result = []
 241         for index in range( start, end, step ):
 242             result.append( self._getindex( index ) )
 243         return result
 244
 245     def _setindex( self, index, value ):
 246         '''set item at idx index.'''
 247         cdef int idx = index
 248         if idx < 0: raise IndexError( "list index out of range" )
 249         if idx >= self.nfields:
 250             raise IndexError( "list index out of range" )
 251
 252         if isNew( self.fields[idx], self.data, self.nbytes ):
 253             free( self.fields[idx] )
 254
 255         self.is_modified = 1
 256
 257         if value == None:
 258             self.fields[idx] = NULL
 259             return
 260
 261         # conversion with error checking
 262         value = _force_bytes(value)
 263         cdef char * tmp = <char*>value
 264         self.fields[idx] = <char*>malloc( (strlen( tmp ) + 1) * sizeof(char) )
 265         if self.fields[idx] == NULL:
 266             raise ValueError("out of memory" )
 267         strcpy( self.fields[idx], tmp )
 268
 269     def __setitem__(self, index, value ):
 270         '''set item at *index* to *value*'''
 271         cdef int i = index
 272         if i < 0: i += self.nfields
 273         i += self.offset
 274
 275         self._setindex( i, value )
 276
 277     def __len__(self):
 278         return self.nfields
 279
 280     def __iter__(self):
 281         self.index = 0
 282         return self
 283
 284     def __next__(self):
 285         """python version of next().
 286         """
 287         if self.index >= self.nfields:
 288             raise StopIteration
 289         cdef char * retval = self.fields[self.index]
 290         self.index += 1
 291         if retval == NULL: return None
 292         else: return retval
 293
 294     def __str__(self):
 295         '''return original data'''
 296         # copy and replace \0 bytes with \t characters
 297         if self.is_modified:
 298             # todo: treat NULL values
 299             result = []
 300             for x in xrange( 0, self.nfields ):
 301                 result.append( StrOrEmpty( self.fields[x]).decode('ascii') )
 302             return "\t".join( result )
 303         else:
 304             cpy = <char*>calloc( sizeof(char), self.nbytes+1 )
 305             if cpy == NULL:
 306                 raise ValueError("out of memory" )
 307             memcpy( cpy, self.data, self.nbytes+1)
 308             for x from 0 <= x < self.nbytes:
 309                 if cpy[x] == '\0': cpy[x] = '\t'
 310             result = cpy[:self.nbytes]
 311             free(cpy)
 312             return result.decode('ascii')
 313
 314 def toDot( v ):
 315     '''convert value to '.' if None'''
 316     if v == None: return "."
 317     else: return str(v)
 318
 319 def quote( v ):
 320     '''return a quoted attribute.'''
 321     if type(v) in types.StringTypes:
 322         return '"%s"' % v
 323     else:
 324         return str(v)
 325
 326 cdef class GTFProxy( TupleProxy ):
 327     '''Proxy class for access to GTF fields.
 328
 329     This class represents a GTF entry for fast read-access.
 330     Write-access has been added as well, though some care must
 331     be taken. If any of the string fields (contig, source, ...)
 332     are set, the new value is tied to the lifetime of the
 333     argument that was supplied.
 334
 335     The only exception is the attributes field when set from
 336     a dictionary - this field will manage its own memory.
 337     '''
 338
 339     def __cinit__(self ):
 340         # automatically calls TupleProxy.__cinit__
 341         self.hasOwnAttributes = False
 342         self._attributes = NULL
 343
 344     def __dealloc__(self):
 345         # automatically calls TupleProxy.__dealloc__
 346         if self.hasOwnAttributes:
 347             free(self._attributes)
 348
 349     cdef int getMaxFields( self, size_t nbytes ):
 350         '''return max number of fields.'''
 351         return 9
 352
 353     property contig:
 354        '''contig of feature.'''
 355        def __get__( self ): return self._getindex( 0 )
 356        def __set__( self, value ): self._setindex( 0, value )
 357
 358     property source:
 359        '''feature source.'''
 360        def __get__( self ): return self._getindex( 1 )
 361        def __set__( self, value ): self._setindex( 1, value )
 362
 363     property feature:
 364        '''feature name.'''
 365        def __get__( self ): return self._getindex( 2 )
 366        def __set__( self, value ): self._setindex( 2, value )
 367
 368     property start:
 369        '''feature start (in 0-based open/closed coordinates).'''
 370        def __get__( self ): return int( self._getindex( 3 )) - 1
 371        def __set__( self, value ): self._setindex( 3, str(value+1) )
 372
 373     property end:
 374        '''feature end (in 0-based open/closed coordinates).'''
 375        def __get__( self ): return int( self._getindex( 4 ) )
 376        def __set__( self, value ): self._setindex( 4, str(value) )
 377
 378     property score:
 379        '''feature score.'''
 380        def __get__( self ):
 381            v = self._getindex(5)
 382            if v == "" or v[0] == '.':
 383                return None
 384            else:
 385                return float(v)
 386
 387        def __set__( self, value ): self._setindex( 5, value )
 388
 389     property strand:
 390        '''feature strand.'''
 391        def __get__( self ): return self._getindex( 6 )
 392        def __set__( self, value ): self._setindex( 6, value )
 393
 394     property frame:
 395        '''feature frame.'''
 396        def __get__( self ): return self._getindex( 7 )
 397        def __set__( self, value ): self._setindex( 7, value )
 398
 399     property attributes:
 400        '''feature attributes (as a string).'''
 401        def __get__( self ):
 402            if self.hasOwnAttributes:
 403                return self._attributes
 404            else:
 405                return self._getindex( 8 )
 406        def __set__( self, value ):
 407            if self.hasOwnAttributes:
 408                free(self._attributes)
 409                self._attributes = NULL
 410                self.hasOwnAttributes = False
 411            self._setindex(8, value )
 412
 413     cdef char * getAttributes( self ):
 414        '''return pointer to attributes.'''
 415        if self.hasOwnAttributes:
 416            return self._attributes
 417        else:
 418            return self.fields[ 8 ]
 419
 420     def asDict( self ):
 421         """parse attributes - return as dict
 422         """
 423
 424         # remove comments
 425         attributes = self.attributes
 426
 427         # separate into fields
 428         fields = [ x.strip() for x in attributes.split(";")[:-1]]
 429
 430         result = {}
 431
 432         for f in fields:
 433
 434             d = [ x.strip() for x in f.split(" ")]
 435
 436             n,v = d[0], d[1]
 437             if len(d) > 2: v = d[1:]
 438
 439             if v[0] == '"' and v[-1] == '"':
 440                 v = v[1:-1]
 441             else:
 442                 ## try to convert to a value
 443                 try:
 444                     v = float( v )
 445                     v = int( v )
 446                 except ValueError:
 447                     pass
 448                 except TypeError:
 449                     pass
 450
 451             result[n] = v
 452
 453         return result
 454
 455     def fromDict( self, d ):
 456         '''set attributes from a dictionary.'''
 457         cdef char * p
 458         cdef int l
 459
 460         # clean up if this field is set twice
 461         if self.hasOwnAttributes:
 462             free(self._attributes)
 463
 464         aa = []
 465         for k,v in d.items():
 466             if type(v) in types.StringTypes:
 467                 aa.append( '%s "%s"' % (k,v) )
 468             else:
 469                 aa.append( '%s %s' % (k,str(v)) )
 470
 471         a = "; ".join( aa ) + ";"
 472         p = a
 473         l = len(a)
 474         self._attributes = <char *>calloc( l + 1, sizeof(char) )
 475         if self._attributes == NULL:
 476             raise ValueError("out of memory" )
 477         memcpy( self._attributes, p, l )
 478
 479         self.hasOwnAttributes = True
 480         self.is_modified = True
 481
 482     def __str__(self):
 483         cdef char * cpy
 484         cdef int x
 485
 486         if self.is_modified:
 487             return "\t".join(
 488                 (self.contig,
 489                  self.source,
 490                  self.feature,
 491                  str(self.start+1),
 492                  str(self.end),
 493                  toDot(self.score),
 494                  self.strand,
 495                  self.frame,
 496                  self.attributes ) )
 497         else:
 498             return TupleProxy.__str__(self)
 499
 500     def invert( self, int lcontig ):
 501         '''invert coordinates to negative strand coordinates
 502
 503         This method will only act if the feature is on the
 504         negative strand.'''
 505
 506         if self.strand[0] == '-':
 507             start = min(self.start, self.end)
 508             end = max(self.start, self.end)
 509             self.start, self.end = lcontig - end, lcontig - start
 510
 511     def keys( self ):
 512         '''return a list of attributes defined in this entry.'''
 513         r = self.attributes
 514         return [ x.strip().split(" ")[0] for x in r.split(";") if x.strip() != '' ]
 515
 516     def __getitem__(self, key):
 517         return self.__getattr__( key )
 518
 519     def __getattr__(self, item ):
 520         """Generic lookup of attribute from GFF/GTF attributes
 521         Only called if there *isn't* an attribute with this name
 522         """
 523         cdef char * start
 524         cdef char * query
 525         cdef char * cpy
 526         cdef char * end
 527         cdef int l
 528
 529         #
 530         # important to use the getAttributes function.
 531         # Using the self.attributes property to access
 532         # the attributes caused a hard-to-trace bug
 533         # in which fields in the attribute string were
 534         # set to 0.
 535         # Running through valgrind complained that
 536         # memory was accessed in the memory field
 537         # that has been released. It is not clear
 538         # why this happened and might be a cython bug
 539         # (Version 0.16). The valgrind warnings
 540         # disappeard after accessing the C data structures
 541         # directly and so did the bug.
 542         cdef char * attributes = self.getAttributes()
 543
 544         r = _force_bytes(item)
 545         query = r
 546         start = strstr( attributes, query)
 547
 548         if start == NULL:
 549             raise AttributeError("'GTFProxy' has no attribute '%s'" % item )
 550
 551         start += strlen(query) + 1
 552         # skip gaps before
 553         while start[0] == ' ': start += 1
 554
 555         if start[0] == '"':
 556             start += 1
 557             end = start
 558             while end[0] != '\0' and end[0] != '"': end += 1
 559             l = end - start
 560             result = _force_str( PyBytes_FromStringAndSize( start, l ) )
 561             return result
 562         else:
 563             return _force_str( start )
 564
 565     def setAttribute( self, name, value ):
 566         '''convenience method to set an attribute.'''
 567         r = self.asDict()
 568         r[name] = value
 569         self.fromDict( r )
 570
 571 cdef class NamedTupleProxy( TupleProxy ):
 572
 573     map_key2field = {}
 574
 575     def __setattr__(self, key, value ):
 576         '''set attribute.'''
 577         cdef int idx
 578         idx, f = self.map_key2field[key]
 579         if self.nfields < idx:
 580             raise KeyError( "field %s not set" % key )
 581         TupleProxy.__setitem__(self, idx, str(value) )
 582
 583     def __getattr__(self, key ):
 584         cdef int idx
 585         idx, f = self.map_key2field[key]
 586         if self.nfields < idx:
 587             raise KeyError( "field %s not set" % key )
 588         return f( self.fields[idx] )
 589
 590 cdef class BedProxy( NamedTupleProxy ):
 591     '''Proxy class for access to Bed fields.
 592
 593     This class represents a GTF entry for fast read-access.
 594     '''
 595     map_key2field = {
 596         'contig' : (0, bytes),
 597         'start' : (1, int),
 598         'end' : (2, int),
 599         'name' : (3, bytes),
 600         'score' : (4, float),
 601         'strand' : (5, bytes),
 602         'thickStart' : (6, int ),
 603         'thickEnd' : (7, int),
 604         'itemRGB' : (8, bytes),
 605         'blockCount': (9, int),
 606         'blockSizes': (10, bytes),
 607         'blockStarts': (11, bytes), }
 608
 609     cdef int getMaxFields( self, size_t nbytes ):
 610         '''return max number of fields.'''
 611         return 12
 612
 613     cdef update( self, char * buffer, size_t nbytes ):
 614         '''update internal data.
 615
 616         nbytes does not include the terminal '\0'.
 617         '''
 618         TupleProxy.update( self, buffer, nbytes )
 619
 620         if self.nfields < 3:
 621             raise ValueError( "bed format requires at least three columns" )
 622
 623         # determines bed format
 624         self.bedfields = self.nfields
 625
 626         # do automatic conversion
 627         self.contig = self.fields[0]
 628         self.start = atoi( self.fields[1] )
 629         self.end = atoi( self.fields[2] )
 630
 631     # __setattr__ in base class seems to take precedence
 632     # hence implement setters in __setattr__
 633     #property start:
 634     #    def __get__( self ): return self.start
 635     #property end:
 636     #    def __get__( self ): return self.end
 637
 638     def __str__(self):
 639
 640         cdef int save_fields = self.nfields
 641         # ensure fields to use correct format
 642         self.nfields = self.bedfields
 643         retval = TupleProxy.__str__( self )
 644         self.nfields = save_fields
 645         return retval
 646
 647     def __setattr__(self, key, value ):
 648         '''set attribute.'''
 649         if key == "start": self.start = value
 650         elif key == "end": self.end = value
 651
 652         cdef int idx
 653         idx, f = self.map_key2field[key]
 654         TupleProxy._setindex(self, idx, str(value) )
 655
 656 cdef class VCFProxy( NamedTupleProxy ):
 657     '''Proxy class for access to VCF fields.
 658
 659     The genotypes are accessed via index.
 660     '''
 661     map_key2field = {
 662         'contig' : (0, bytes),
 663         'pos' : (1, int),
 664         'id' : (2, bytes),
 665         'ref' : (3, bytes),
 666         'alt' : (4, bytes),
 667         'qual' : (5, bytes),
 668         'filter' : (6, bytes),
 669         'info' : (7, bytes),
 670         'format' : (8, bytes) }
 671
 672     def __cinit__(self ):
 673         # automatically calls TupleProxy.__cinit__
 674         # start indexed access at genotypes
 675         self.offset = 9
 676
 677     cdef update( self, char * buffer, size_t nbytes ):
 678         '''update internal data.
 679
 680         nbytes does not include the terminal '\0'.
 681         '''
 682         TupleProxy.update( self, buffer, nbytes )
 683
 684         self.contig = self.fields[0]
 685         # vcf counts from 1 - correct here
 686         self.pos = atoi( self.fields[1] ) - 1
 687
 688     def __len__(self):
 689         '''return number of genotype fields.'''
 690         return max(0, self.nfields - 9)
 691
 692     property pos:
 693        '''feature end (in 0-based open/closed coordinates).'''
 694        def __get__( self ):
 695            return self.pos
 696
 697     def __setattr__(self, key, value ):
 698         '''set attribute.'''
 699         if key == "pos":
 700             self.pos = value
 701             value += 1
 702
 703         cdef int idx
 704         idx, f = self.map_key2field[key]
 705         TupleProxy._setindex(self, idx, str(value) )
 706