tests/tabix_test.py

   1 #!/usr/bin/env python
   2 '''unit testing code for pysam.
   3
   4 Execute in the :file:`tests` directory as it requires the Makefile
   5 and data files located there.
   6 '''
   7
   8 import sys, os, shutil, gzip
   9 import pysam
  10 import unittest
  11 import itertools
  12 import subprocess
  13
  14 def checkBinaryEqual( filename1, filename2 ):
  15     '''return true if the two files are binary equal.'''
  16     if os.path.getsize( filename1 ) !=  os.path.getsize( filename2 ):
  17         return False
  18
  19     infile1 = open(filename1, "rb")
  20     infile2 = open(filename2, "rb")
  21
  22     def chariter( infile ):
  23         while 1:
  24             c = infile.read(1)
  25             if c == "": break
  26             yield c
  27
  28     found = False
  29     for c1,c2 in itertools.izip( chariter( infile1), chariter( infile2) ):
  30         if c1 != c2: break
  31     else:
  32         found = True
  33
  34     infile1.close()
  35     infile2.close()
  36     return found
  37
  38 class TestIndexing(unittest.TestCase):
  39     filename = "example.gtf.gz"
  40     filename_idx = "example.gtf.gz.tbi"
  41
  42     def setUp( self ):
  43
  44         self.tmpfilename = "tmp_%i.gtf.gz" % id(self)
  45         shutil.copyfile( self.filename, self.tmpfilename )
  46
  47     def testIndexPreset( self ):
  48         '''test indexing via preset.'''
  49
  50         pysam.tabix_index( self.tmpfilename, preset = "gff" )
  51         checkBinaryEqual( self.tmpfilename + ".tbi", self.filename_idx )
  52
  53     def tearDown( self ):
  54         os.unlink( self.tmpfilename )
  55         os.unlink( self.tmpfilename + ".tbi" )
  56
  57 class TestCompression(unittest.TestCase):
  58     filename = "example.gtf.gz"
  59     filename_idx = "example.gtf.gz.tbi"
  60
  61     def setUp( self ):
  62
  63         self.tmpfilename = "tmp_%i.gtf" % id(self)
  64         infile = gzip.open( self.filename, "r")
  65         outfile = open( self.tmpfilename, "w" )
  66         outfile.write( "".join(infile.readlines()) )
  67         outfile.close()
  68         infile.close()
  69
  70     def testIndexPreset( self ):
  71         '''test indexing via preset.'''
  72
  73         pysam.tabix_index( self.tmpfilename, preset = "gff" )
  74         checkBinaryEqual( self.tmpfilename + ".gz", self.filename )
  75         checkBinaryEqual( self.tmpfilename + ".gz.tbi", self.filename_idx )
  76
  77     def tearDown( self ):
  78         os.unlink( self.tmpfilename + ".gz" )
  79         os.unlink( self.tmpfilename + ".gz.tbi" )
  80
  81 class TestIteration( unittest.TestCase ):
  82
  83     filename = "example.gtf.gz"
  84
  85     def setUp( self ):
  86
  87         self.tabix = pysam.Tabixfile( self.filename )
  88         lines = [ x for x in gzip.open(self.filename).readlines() if not x.startswith("#") ]
  89         # creates index of contig, start, end, adds content without newline.
  90         self.compare = [
  91             (x[0][0], int(x[0][3]), int(x[0][4]), x[1])
  92             for x in [ (y.split("\t"), y[:-1]) for y in lines ] ]
  93
  94     def getSubset( self, contig = None, start = None, end = None):
  95
  96         if contig == None:
  97             # all lines
  98             subset = [ x[3] for x in self.compare ]
  99         else:
 100             if start != None and end == None:
 101                 # until end of contig
 102                 subset = [ x[3] for x in self.compare if x[0] == contig and x[2] > start ]
 103             elif start == None and end != None:
 104                 # from start of contig
 105                 subset = [ x[3] for x in self.compare if x[0] == contig and x[1] <= end ]
 106             elif start == None and end == None:
 107                 subset = [ x[3] for x in self.compare if x[0] == contig ]
 108             else:
 109                 # all within interval
 110                 subset = [ x[3] for x in self.compare if x[0] == contig and \
 111                                min( x[2], end) - max(x[1], start) > 0 ]
 112
 113         return subset
 114
 115     def checkPairwise( self, result, ref ):
 116
 117         result.sort()
 118         ref.sort()
 119
 120         a = set(result)
 121         b = set(ref)
 122
 123         self.assertEqual( len(result), len(ref),
 124                           "unexpected number of results: %i, expected %i, differences are %s: %s" \
 125                               % (len(result), len(ref),
 126                                  a.difference(b),
 127                                  b.difference(a) ))
 128
 129         for x, d in enumerate( zip( result, ref )):
 130
 131             self.assertEqual( d[0], d[1],
 132                               "unexpected results in pair %i: '%s', expected '%s'" % \
 133                                   (x,
 134                                    d[0],
 135                                    d[1]) )
 136
 137
 138     def testAll( self ):
 139         result = list(self.tabix.fetch())
 140         ref = self.getSubset( )
 141         self.checkPairwise( result, ref )
 142
 143     def testPerContig( self ):
 144         for contig in ("chr1", "chr2", "chr1", "chr2" ):
 145             result = list(self.tabix.fetch( contig ))
 146             ref = self.getSubset( contig )
 147             self.checkPairwise( result, ref )
 148
 149     def testPerContigToEnd( self ):
 150
 151         end = None
 152         for contig in ("chr1", "chr2", "chr1", "chr2" ):
 153             for start in range( 0, 200000, 1000):
 154                 result = list(self.tabix.fetch( contig, start, end ))
 155                 ref = self.getSubset( contig, start, end )
 156                 self.checkPairwise( result, ref )
 157
 158     def testPerContigFromStart( self ):
 159
 160         start = None
 161         for contig in ("chr1", "chr2", "chr1", "chr2" ):
 162             for end in range( 0, 200000, 1000):
 163                 result = list(self.tabix.fetch( contig, start, end ))
 164                 ref = self.getSubset( contig, start, end )
 165                 self.checkPairwise( result, ref )
 166
 167     def testPerContig( self ):
 168
 169         start, end  = None, None
 170         for contig in ("chr1", "chr2", "chr1", "chr2" ):
 171             result = list(self.tabix.fetch( contig, start, end ))
 172             ref = self.getSubset( contig, start, end )
 173             self.checkPairwise( result, ref )
 174
 175     def testPerInterval( self ):
 176
 177         start, end  = None, None
 178         for contig in ("chr1", "chr2", "chr1", "chr2" ):
 179             for start in range( 0, 200000, 2000):
 180                 for end in range( start, start + 2000, 500):
 181                     result = list(self.tabix.fetch( contig, start, end ))
 182                     ref = self.getSubset( contig, start, end )
 183                     self.checkPairwise( result, ref )
 184
 185
 186     def testInvalidIntervals( self ):
 187
 188         self.assertRaises( ValueError, self.tabix.fetch, "chr1", 0, -10)
 189         self.assertRaises( ValueError, self.tabix.fetch, "chr1", -10, 200)
 190         self.assertRaises( ValueError, self.tabix.fetch, "chr1", 200, 0)
 191         self.assertRaises( ValueError, self.tabix.fetch, "chr1", -10, -20)
 192         self.assertRaises( ValueError, self.tabix.fetch, "chrUn" )
 193
 194     def testGetContigs( self ):
 195         self.assertEqual( sorted(self.tabix.contigs), ["chr1", "chr2"] )
 196         # check that contigs is read-only
 197         self.assertRaises( AttributeError, setattr, self.tabix, "contigs", ["chr1", "chr2"] )
 198
 199     def testHeader( self ):
 200         ref = []
 201         for x in gzip.open( self.filename ):
 202             if not x.startswith("#"): break
 203             ref.append( x[:-1] )
 204         header = list( self.tabix.header )
 205         self.assertEqual( ref, header )
 206
 207 class TestParser( unittest.TestCase ):
 208
 209     filename = "example.gtf.gz"
 210
 211     def setUp( self ):
 212
 213         self.tabix = pysam.Tabixfile( self.filename )
 214         self.compare = [ x[:-1].split("\t") for x in gzip.open( self.filename, "r") if not x.startswith("#") ]
 215
 216     def testRead( self ):
 217
 218         for x, r in enumerate(self.tabix.fetch( parser = pysam.asTuple() )):
 219             self.assertEqual( self.compare[x], list(r) )
 220             self.assertEqual( len(self.compare[x]), len(r) )
 221
 222             for c in range(0,len(r)):
 223                 self.assertEqual( self.compare[x][c], r[c] )
 224
 225     def testWrite( self ):
 226
 227         for x, r in enumerate(self.tabix.fetch( parser = pysam.asTuple() )):
 228             self.assertEqual( self.compare[x], list(r) )
 229             c = list(r)
 230             for y in range(len(r)):
 231                 r[y] = "test_%05i" % y
 232                 c[y] = "test_%05i" % y
 233             self.assertEqual( c, list(r) )
 234             self.assertEqual( "\t".join( c ), str(r) )
 235             # check second assignment
 236             for y in range(len(r)):
 237                 r[y] = "test_%05i" % y
 238             self.assertEqual( c, list(r) )
 239             self.assertEqual( "\t".join( c ), str(r) )
 240
 241     def testUnset( self ):
 242         for x, r in enumerate(self.tabix.fetch( parser = pysam.asTuple() )):
 243             self.assertEqual( self.compare[x], list(r) )
 244             c = list(r)
 245             e = list(r)
 246             for y in range(len(r)):
 247                 r[y] = c[y] = None
 248                 e[y] = ""
 249                 self.assertEqual( c, list(r) )
 250                 self.assertEqual( "\t".join(e), str(r) )
 251
 252 class TestGTF( TestParser ):
 253     def testRead( self ):
 254
 255         for x, r in enumerate(self.tabix.fetch( parser = pysam.asGTF() )):
 256             self.assertEqual( "\t".join( self.compare[x]), str(r) )
 257
 258 class TestBed( unittest.TestCase ):
 259     filename = "example.bed.gz"
 260
 261     def setUp( self ):
 262
 263         self.tabix = pysam.Tabixfile( self.filename )
 264         self.compare = [ x[:-1].split("\t") for x in gzip.open( self.filename, "r") if not x.startswith("#") ]
 265
 266     def testRead( self ):
 267
 268         for x, r in enumerate(self.tabix.fetch( parser = pysam.asBed() )):
 269             c = self.compare[x]
 270             self.assertEqual( "\t".join( c ), str(r) )
 271             self.assertEqual( list(c), list(r) )
 272             self.assertEqual( c[0], r.contig)
 273             self.assertEqual( int(c[1]), r.start)
 274             self.assertEqual( int(c[2]), r.end)
 275
 276     def testWrite( self ):
 277
 278         for x, r in enumerate(self.tabix.fetch( parser = pysam.asBed() )):
 279             c = self.compare[x]
 280             self.assertEqual( "\t".join( c ), str(r) )
 281             self.assertEqual( list(c), list(r) )
 282
 283             r.contig = "test"
 284             self.assertEqual( "test", r.contig)
 285             self.assertEqual( "test", r[0])
 286
 287             r.start += 1
 288             self.assertEqual( int(c[1]) + 1, r.start )
 289             self.assertEqual( str(int(c[1]) + 1), r[1] )
 290
 291             r.end += 1
 292             self.assertEqual( int(c[2]) + 1, r.end )
 293             self.assertEqual( str(int(c[2]) + 1), r[2] )
 294
 295 class TestVCF( TestParser ):
 296
 297     filename = "example.vcf40.gz"
 298     columns = ("contig", "pos", "id",
 299                "ref", "alt", "qual",
 300                "filter", "info", "format" )
 301
 302     def testRead( self ):
 303
 304         ncolumns = len(self.columns)
 305
 306         for x, r in enumerate(self.tabix.fetch( parser = pysam.asVCF() )):
 307             c = self.compare[x]
 308             for y, field in enumerate( self.columns ):
 309                 if field == "pos":
 310                     self.assertEqual( int(c[y])-1, getattr( r, field ) )
 311                     self.assertEqual( int(c[y])-1, r.pos )
 312                 else:
 313                     self.assertEqual( c[y], getattr( r, field ),
 314                                       "mismatch in field %s: %s != %s" %\
 315                                           ( field,c[y], getattr( r, field ) ) )
 316             self.assertEqual( len(c), len( r ) + ncolumns )
 317
 318             for y in range(len(c) - ncolumns):
 319                 self.assertEqual( c[ncolumns+y], r[y] )
 320
 321     def testWrite( self ):
 322
 323         ncolumns = len(self.columns)
 324
 325         for x, r in enumerate(self.tabix.fetch( parser = pysam.asVCF() )):
 326             c = self.compare[x]
 327             for y, field in enumerate( self.columns ):
 328                 if field == "pos":
 329                     r.pos += 1
 330                     self.assertEqual( int(c[y]), getattr( r, field ) )
 331                     self.assertEqual( int(c[y]), r.pos )
 332                 else:
 333                     setattr( r, field, "test_%i" % y)
 334                     c[y] = "test_%i" % y
 335                     self.assertEqual( c[y], getattr( r, field ),
 336                                       "mismatch in field %s: %s != %s" %\
 337                                           ( field,c[y], getattr( r, field ) ) )
 338
 339             self.assertEqual( len(c), len( r ) + ncolumns )
 340
 341             for y in range(len(c) - ncolumns):
 342                 c[ncolumns+y] = "test_%i" % y
 343                 r[y] = "test_%i" % y
 344                 self.assertEqual( c[ncolumns+y], r[y] )
 345
 346 if __name__ == "__main__":
 347
 348     unittest.main()
 349
 350