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 __getitem__( self, key ):
 145
 146         cdef int i = key
 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" )
 152         return self.fields[i]
 153
 154     def _setindex( self, index, value ):
 155         '''set item at idx index.'''
 156         cdef int idx = index
 157         if idx < 0: raise IndexError( "list index out of range" )
 158         if idx >= self.nfields:
 159             raise IndexError( "list index out of range" )
 160
 161         if isNew( self.fields[idx], self.data, self.nbytes ):
 162             free( self.fields[idx] )
 163
 164         self.is_modified = 1
 165
 166         if value == None:
 167             self.fields[idx] = NULL
 168             return
 169
 170         # conversion with error checking
 171         cdef char * tmp = PyString_AsString( value )
 172         self.fields[idx] = <char*>malloc( (strlen( tmp ) + 1) * sizeof(char) )
 173         if self.fields[idx] == NULL:
 174             raise ValueError("out of memory" )
 175         strcpy( self.fields[idx], tmp )
 176
 177     def __setitem__(self, index, value ):
 178         '''set item at *index* to *value*'''
 179         cdef int i = index
 180         if i < 0: i += self.nfields
 181         i += self.offset
 182
 183         self._setindex( i, value )
 184
 185     def __len__(self):
 186         return self.nfields
 187
 188     def __iter__(self):
 189         self.index = 0
 190         return self
 191
 192     def __next__(self):
 193         """python version of next().
 194         """
 195         if self.index >= self.nfields:
 196             raise StopIteration
 197         cdef char * retval = self.fields[self.index]
 198         self.index += 1
 199         if retval == NULL: return None
 200         else: return retval
 201
 202     def __str__(self):
 203         '''return original data'''
 204         # copy and replace \0 bytes with \t characters
 205         if self.is_modified:
 206             # todo: treat NULL values
 207             return "\t".join( [StrOrEmpty( self.fields[x]) for x in xrange(0, self.nfields ) ] )
 208         else:
 209             cpy = <char*>calloc( sizeof(char), self.nbytes+1 )
 210             if cpy == NULL:
 211                 raise ValueError("out of memory" )
 212             memcpy( cpy, self.data, self.nbytes+1)
 213             for x from 0 <= x < self.nbytes:
 214                 if cpy[x] == '\0': cpy[x] = '\t'
 215             result = PyString_FromStringAndSize(cpy, self.nbytes)
 216             free(cpy)
 217             return result
 218
 219 def toDot( v ):
 220     '''convert value to '.' if None'''
 221     if v == None: return "."
 222     else: return str(v)
 223
 224 def quote( v ):
 225     '''return a quoted attribute.'''
 226     if type(v) in types.StringTypes:
 227         return '"%s"' % v
 228     else:
 229         return str(v)
 230
 231 cdef class GTFProxy( TupleProxy ):
 232     '''Proxy class for access to GTF fields.
 233
 234     This class represents a GTF entry for fast read-access.
 235     Write-access has been added as well, though some care must
 236     be taken. If any of the string fields (contig, source, ...)
 237     are set, the new value is tied to the lifetime of the
 238     argument that was supplied.
 239
 240     The only exception is the attributes field when set from
 241     a dictionary - this field will manage its own memory.
 242
 243     '''
 244
 245     def __cinit__(self ):
 246         # automatically calls TupleProxy.__cinit__
 247         self.hasOwnAttributes = False
 248
 249     def __dealloc__(self):
 250         # automatically calls TupleProxy.__dealloc__
 251         if self.hasOwnAttributes:
 252             free(self.attributes)
 253
 254     cdef int getMaxFields( self, size_t nbytes ):
 255         '''return max number of fields.'''
 256         return 9
 257
 258     cdef update( self, char * buffer, size_t nbytes ):
 259         '''update internal data.
 260
 261         nbytes does not include the terminal '\0'.
 262         '''
 263         cdef int end
 264         cdef char * cstart, * cend, * cscore
 265         self.contig = buffer
 266         cdef char * pos
 267
 268         if buffer[nbytes] != 0:
 269             raise ValueError( "incomplete line at %s" % buffer )
 270
 271         self.source = pos = nextItem( buffer )
 272         self.feature = pos = nextItem( pos )
 273         cstart = pos = nextItem( pos )
 274         cend = pos = nextItem( pos )
 275         self.score = pos = nextItem( pos )
 276         self.strand = pos = nextItem( pos )
 277         self.frame = pos = nextItem( pos )
 278         self.attributes = pos = nextItem( pos )
 279
 280         self.start = atoi( cstart ) - 1
 281         self.end = atoi( cend )
 282
 283     property contig:
 284        '''contig of feature.'''
 285        def __get__( self ): return self.contig
 286        def __set__( self, value ):
 287            self.is_modified = True
 288            self.contig = value
 289
 290     property feature:
 291        '''feature name.'''
 292        def __get__( self ): return self.feature
 293        def __set__( self, value ):
 294            self.is_modified = True
 295            self.feature = value
 296
 297     property source:
 298        '''feature source.'''
 299        def __get__( self ): return self.source
 300        def __set__( self, value ):
 301            self.is_modified = True
 302            self.source = value
 303
 304     property start:
 305        '''feature start (in 0-based open/closed coordinates).'''
 306        def __get__( self ): return self.start
 307        def __set__( self, value ):
 308            self.is_modified = True
 309            self.start = value
 310
 311     property end:
 312        '''feature end (in 0-based open/closed coordinates).'''
 313        def __get__( self ): return self.end
 314        def __set__( self, value ):
 315            self.is_modified = True
 316            self.end = value
 317
 318     property score:
 319        '''feature score.'''
 320        def __get__( self ):
 321            if self.score[0] == '.' and self.score[1] == '\0' :
 322                return None
 323            else:
 324                return atof(self.score)
 325        def __set__( self, value ):
 326            self.is_modified = True
 327            self.score = value
 328
 329     property strand:
 330        '''feature strand.'''
 331        def __get__( self ): return self.strand
 332        def __set__( self, value ):
 333            self.is_modified = True
 334            self.strand = value
 335
 336     property frame:
 337        '''feature frame.'''
 338        def __get__( self ): return self.frame
 339        def __set__( self, value ):
 340            self.is_modified = True
 341            self.frame = value
 342
 343     property attributes:
 344        '''feature attributes (as a string).'''
 345        def __get__( self ): return self.attributes
 346        def __set__( self, value ):
 347            self.is_modified = True
 348            self.attributes = value
 349
 350     def asDict( self ):
 351         """parse attributes - return as dict
 352         """
 353
 354         # remove comments
 355         attributes = self.attributes
 356
 357         # separate into fields
 358         fields = [ x.strip() for x in attributes.split(";")[:-1]]
 359
 360         result = {}
 361
 362         for f in fields:
 363
 364             d = [ x.strip() for x in f.split(" ")]
 365
 366             n,v = d[0], d[1]
 367             if len(d) > 2: v = d[1:]
 368
 369             if v[0] == '"' and v[-1] == '"':
 370                 v = v[1:-1]
 371             else:
 372                 ## try to convert to a value
 373                 try:
 374                     v = float( v )
 375                     v = int( v )
 376                 except ValueError:
 377                     pass
 378                 except TypeError:
 379                     pass
 380
 381             result[n] = v
 382
 383         return result
 384
 385     def fromDict( self, d ):
 386         '''set attributes from a dictionary.'''
 387         cdef char * p
 388         cdef int l
 389
 390         # clean up if this field is set twice
 391         if self.hasOwnAttributes:
 392             free(self.attributes)
 393
 394         aa = []
 395         for k,v in d.items():
 396             if type(v) == types.StringType:
 397                 aa.append( '%s "%s"' % (k,v) )
 398             else:
 399                 aa.append( '%s %s' % (k,str(v)) )
 400
 401         a = "; ".join( aa ) + ";"
 402         p = a
 403         l = len(a)
 404         self.attributes = <char *>calloc( l + 1, sizeof(char) )
 405         if self.attributes == NULL:
 406             raise ValueError("out of memory" )
 407         memcpy( self.attributes, p, l )
 408
 409         self.hasOwnAttributes = True
 410         self.is_modified = True
 411
 412     def __str__(self):
 413         cdef char * cpy
 414         cdef int x
 415
 416         if self.is_modified:
 417             return "\t".join(
 418                 (self.contig,
 419                  self.source,
 420                  self.feature,
 421                  str(self.start+1),
 422                  str(self.end),
 423                  toDot(self.score),
 424                  self.strand,
 425                  self.frame,
 426                  self.attributes ) )
 427         else:
 428             return TupleProxy.__str__(self)
 429
 430     def invert( self, int lcontig ):
 431         '''invert coordinates to negative strand coordinates
 432
 433         This method will only act if the feature is on the
 434         negative strand.'''
 435
 436         if self.strand[0] == '-':
 437             start = min(self.start, self.end)
 438             end = max(self.start, self.end)
 439             self.start, self.end = lcontig - end, lcontig - start
 440
 441     def keys( self ):
 442         '''return a list of attributes defined in this entry.'''
 443         r = self.attributes
 444         return [ x.strip().split(" ")[0] for x in r.split(";") if x.strip() != '' ]
 445
 446     def __getitem__(self, item):
 447         return self.__getattr__( item )
 448
 449     def __getattr__(self, item ):
 450         """Generic lookup of attribute from GFF/GTF attributes
 451         Only called if there *isn't* an attribute with this name
 452         """
 453         cdef char * start
 454         cdef char * query
 455         cdef char * cpy
 456         cdef char * end
 457         cdef int l
 458         query = item
 459
 460         start = strstr( self.attributes, query)
 461         if start == NULL:
 462             raise AttributeError("'GTFProxy' has no attribute '%s'" % item )
 463
 464         start += strlen(query) + 1
 465         # skip gaps before
 466         while start[0] == " ": start += 1
 467         if start[0] == '"':
 468             start += 1
 469             end = start
 470             while end[0] != '\0' and end[0] != '"': end += 1
 471             l = end - start + 1
 472             cpy = <char*>calloc( l, sizeof(char ) )
 473             if cpy == NULL: raise ValueError("out of memory" )
 474             memcpy( cpy, start, l )
 475             cpy[l-1] = '\0'
 476             result = cpy
 477             free(cpy)
 478             return result
 479         else:
 480             return start
 481
 482     def setAttribute( self, name, value ):
 483         '''convenience method to set an attribute.'''
 484         r = self.asDict()
 485         r[name] = value
 486         self.fromDict( r )
 487
 488 cdef class NamedTupleProxy( TupleProxy ):
 489
 490     map_key2field = {}
 491
 492     def __setattr__(self, key, value ):
 493         '''set attribute.'''
 494         cdef int idx
 495         idx, f = self.map_key2field[key]
 496         if self.nfields < idx:
 497             raise KeyError( "field %s not set" % key )
 498         TupleProxy.__setitem__(self, idx, str(value) )
 499
 500     def __getattr__(self, key ):
 501         cdef int idx
 502         idx, f = self.map_key2field[key]
 503         if self.nfields < idx:
 504             raise KeyError( "field %s not set" % key )
 505         return f( self.fields[idx] )
 506
 507 cdef class BedProxy( NamedTupleProxy ):
 508     '''Proxy class for access to Bed fields.
 509
 510     This class represents a GTF entry for fast read-access.
 511     '''
 512     map_key2field = {
 513         'contig' : (0, str),
 514         'start' : (1, int),
 515         'end' : (2, int),
 516         'name' : (3, str),
 517         'score' : (4, float),
 518         'strand' : (5, str),
 519         'thickStart' : (6,int ),
 520         'thickEnd' : (7,int),
 521         'itemRGB' : (8,str),
 522         'blockCount': (9,int),
 523         'blockSizes': (10,str),
 524         'blockStarts': (11,str), }
 525
 526     cdef int getMaxFields( self, size_t nbytes ):
 527         '''return max number of fields.'''
 528         return 12
 529
 530     cdef update( self, char * buffer, size_t nbytes ):
 531         '''update internal data.
 532
 533         nbytes does not include the terminal '\0'.
 534         '''
 535         TupleProxy.update( self, buffer, nbytes )
 536
 537         if self.nfields < 3:
 538             raise ValueError( "bed format requires at least three columns" )
 539
 540         # determines bed format
 541         self.bedfields = self.nfields
 542
 543         # do automatic conversion
 544         self.contig = self.fields[0]
 545         self.start = atoi( self.fields[1] )
 546         self.end = atoi( self.fields[2] )
 547
 548     # __setattr__ in base class seems to take precedence
 549     # hence implement setters in __setattr__
 550     #property start:
 551     #    def __get__( self ): return self.start
 552     #property end:
 553     #    def __get__( self ): return self.end
 554
 555     def __str__(self):
 556
 557         cdef int save_fields = self.nfields
 558         # ensure fields to use correct format
 559         self.nfields = self.bedfields
 560         retval = TupleProxy.__str__( self )
 561         self.nfields = save_fields
 562         return retval
 563
 564     def __setattr__(self, key, value ):
 565         '''set attribute.'''
 566         if key == "start": self.start = value
 567         elif key == "end": self.end = value
 568
 569         cdef int idx
 570         idx, f = self.map_key2field[key]
 571         TupleProxy._setindex(self, idx, str(value) )
 572
 573 cdef class VCFProxy( NamedTupleProxy ):
 574     '''Proxy class for access to VCF fields.
 575
 576     The genotypes are accessed via index.
 577     '''
 578     map_key2field = {
 579         'contig' : (0, str),
 580         'pos' : (1, int),
 581         'id' : (2, str),
 582         'ref' : (3, str),
 583         'alt' : (4, str),
 584         'qual' : (5, str),
 585         'filter' : (6,str),
 586         'info' : (7,str),
 587         'format' : (8,str) }
 588
 589     def __cinit__(self ):
 590         # automatically calls TupleProxy.__cinit__
 591         # start indexed access at genotypes
 592         self.offset = 9
 593
 594     cdef update( self, char * buffer, size_t nbytes ):
 595         '''update internal data.
 596
 597         nbytes does not include the terminal '\0'.
 598         '''
 599         TupleProxy.update( self, buffer, nbytes )
 600
 601         self.contig = self.fields[0]
 602         # vcf counts from 1 - correct here
 603         self.pos = atoi( self.fields[1] ) - 1
 604
 605     def __len__(self):
 606         return max(0, self.nfields - 9)
 607
 608     def __setattr__(self, key, value ):
 609         '''set attribute.'''
 610         if key == "pos":
 611             self.pos = value
 612             value += 1
 613
 614         cdef int idx
 615         idx, f = self.map_key2field[key]
 616         TupleProxy._setindex(self, idx, str(value) )
 617