pysam/TabProxies.pyx

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