tests/pysam_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 pysam
   9 import unittest
  10 import os, re
  11 import itertools
  12 import subprocess
  13 import shutil
  14
  15
  16 def checkBinaryEqual( filename1, filename2 ):
  17     '''return true if the two files are binary equal.'''
  18     if os.path.getsize( filename1 ) !=  os.path.getsize( filename2 ):
  19         return False
  20
  21     infile1 = open(filename1, "rb")
  22     infile2 = open(filename2, "rb")
  23
  24     def chariter( infile ):
  25         while 1:
  26             c = infile.read(1)
  27             if c == "": break
  28             yield c
  29
  30     found = False
  31     for c1,c2 in itertools.izip( chariter( infile1), chariter( infile2) ):
  32         if c1 != c2: break
  33     else:
  34         found = True
  35
  36     infile1.close()
  37     infile2.close()
  38     return found
  39
  40 def runSamtools( cmd ):
  41     '''run a samtools command'''
  42
  43     try:
  44         retcode = subprocess.call(cmd, shell=True)
  45         if retcode < 0:
  46             print >>sys.stderr, "Child was terminated by signal", -retcode
  47     except OSError, e:
  48         print >>sys.stderr, "Execution failed:", e
  49
  50 def getSamtoolsVersion():
  51     '''return samtools version'''
  52
  53     pipe = subprocess.Popen("samtools", shell=True, stderr=subprocess.PIPE).stderr
  54     lines = "".join(pipe.readlines())
  55     return re.search( "Version:\s+(\S+)", lines).groups()[0]
  56
  57 class BinaryTest(unittest.TestCase):
  58     '''test samtools command line commands and compare
  59     against pysam commands.
  60
  61     Tests fail, if the output is not binary identical.
  62     '''
  63
  64     first_time = True
  65
  66     # a list of commands to test
  67     mCommands = \
  68         { "faidx" : \
  69         (
  70             ("ex1.fa.fai", "samtools faidx ex1.fa"),
  71             ("pysam_ex1.fa.fai", (pysam.faidx, "ex1.fa") ),
  72             ),
  73           "import" :
  74               (
  75                 ("ex1.bam", "samtools import ex1.fa.fai ex1.sam.gz ex1.bam" ),
  76                 ("pysam_ex1.bam", (pysam.samimport, "ex1.fa.fai ex1.sam.gz pysam_ex1.bam") ),
  77                 ),
  78           "index":
  79               (
  80                 ("ex1.bam.bai", "samtools index ex1.bam" ),
  81                 ("pysam_ex1.bam.bai", (pysam.index, "pysam_ex1.bam" ) ),
  82                 ),
  83           "pileup1" :
  84           (
  85                 ("ex1.pileup", "samtools pileup -cf ex1.fa ex1.bam > ex1.pileup" ),
  86                 ("pysam_ex1.pileup", (pysam.pileup, "-c -f ex1.fa ex1.bam" ) )
  87                 ),
  88           "pileup2" :
  89           (
  90                 ("ex1.glf", "samtools pileup -gf ex1.fa ex1.bam > ex1.glf" ),
  91                 ("pysam_ex1.glf", (pysam.pileup, "-g -f ex1.fa ex1.bam" ) )
  92                 ),
  93           "glfview" :
  94         (
  95                 ("ex1.glfview", "samtools glfview ex1.glf > ex1.glfview"),
  96                 ("pysam_ex1.glfview", (pysam.glfview, "ex1.glf" ) ),
  97                 ),
  98           "view" :
  99         (
 100                 ("ex1.view", "samtools view ex1.bam > ex1.view"),
 101                 ("pysam_ex1.view", (pysam.view, "ex1.bam" ) ),
 102                 ),
 103         }
 104
 105     # some tests depend on others. The order specifies in which order
 106     # the samtools commands are executed.
 107     mOrder = ('faidx', 'import', 'index', 'pileup1', 'pileup2', 'glfview', 'view' )
 108
 109     def setUp( self ):
 110         '''setup tests.
 111
 112         For setup, all commands will be run before the first test is
 113         executed. Individual tests will then just compare the output
 114         files.
 115         '''
 116         if BinaryTest.first_time:
 117             # copy the source
 118             shutil.copy( "ex1.fa", "pysam_ex1.fa" )
 119
 120             for label in self.mOrder:
 121                 command = self.mCommands[label]
 122                 samtools_target, samtools_command = command[0]
 123                 pysam_target, pysam_command = command[1]
 124                 runSamtools( samtools_command )
 125                 pysam_method, pysam_options = pysam_command
 126                 output = pysam_method( *pysam_options.split(" "), raw=True)
 127                 if ">" in samtools_command:
 128                     outfile = open( pysam_target, "w" )
 129                     for line in output: outfile.write( line )
 130                     outfile.close()
 131
 132             BinaryTest.first_time = False
 133
 134         samtools_version = getSamtoolsVersion()
 135         if samtools_version != pysam.__samtools_version__:
 136             raise ValueError("versions of pysam/samtools and samtools differ: %s != %s" % \
 137                                  (pysam.__samtools_version__,
 138                                   samtools_version ))
 139
 140     def checkCommand( self, command ):
 141
 142         if command:
 143             samtools_target, pysam_target = self.mCommands[command][0][0], self.mCommands[command][1][0]
 144             self.assertTrue( checkBinaryEqual( samtools_target, pysam_target ),
 145                              "%s failed: files %s and %s are not the same" % (command, samtools_target, pysam_target) )
 146
 147     def testImport( self ):
 148         self.checkCommand( "import" )
 149
 150     def testIndex( self ):
 151         self.checkCommand( "index" )
 152
 153     def testPileup1( self ):
 154         self.checkCommand( "pileup1" )
 155
 156     def testPileup2( self ):
 157         self.checkCommand( "pileup2" )
 158
 159     def testGLFView( self ):
 160         self.checkCommand( "glfview" )
 161
 162     def testView( self ):
 163         self.checkCommand( "view" )
 164
 165     def testEmptyIndex( self ):
 166         self.assertRaises( pysam.SamtoolsError, pysam.index, "exdoesntexist.bam" )
 167
 168     def __del__(self):
 169         return
 170         for label, command in self.mCommands.iteritems():
 171             samtools_target, samtools_command = command[0]
 172             pysam_target, pysam_command = command[1]
 173             if os.path.exists( samtools_target): os.remove( samtools_target )
 174             if os.path.exists( pysam_target): os.remove( pysam_target )
 175         if os.path.exists( "pysam_ex1.fa" ): os.remove( "pysam_ex1.fa" )
 176
 177 class IOTest(unittest.TestCase):
 178     '''check if reading samfile and writing a samfile are consistent.'''
 179
 180     def checkEcho( self, input_filename, reference_filename,
 181                    output_filename,
 182                    input_mode, output_mode, use_template = True):
 183         '''iterate through *input_filename* writing to *output_filename* and
 184         comparing the output to *reference_filename*.
 185
 186         The files are opened according to the *input_mode* and *output_mode*.
 187
 188         If *use_template* is set, the header is copied from infile using the
 189         template mechanism, otherwise target names and lengths are passed explicitely.
 190         '''
 191
 192         infile = pysam.Samfile( input_filename, input_mode )
 193         if use_template:
 194             outfile = pysam.Samfile( output_filename, output_mode, template = infile )
 195         else:
 196             outfile = pysam.Samfile( output_filename, output_mode,
 197                                      referencenames = infile.references,
 198                                      referencelengths = infile.lengths )
 199
 200         iter = infile.fetch()
 201         for x in iter: outfile.write( x )
 202         infile.close()
 203         outfile.close()
 204
 205         self.assertTrue( checkBinaryEqual( reference_filename, output_filename),
 206                          "files %s and %s are not the same" % (reference_filename, output_filename) )
 207
 208     def testReadWriteBam( self ):
 209
 210         input_filename = "ex1.bam"
 211         output_filename = "pysam_ex1.bam"
 212         reference_filename = "ex1.bam"
 213
 214         self.checkEcho( input_filename, reference_filename, output_filename,
 215                         "rb", "wb" )
 216
 217     def testReadWriteBamWithTargetNames( self ):
 218
 219         input_filename = "ex1.bam"
 220         output_filename = "pysam_ex1.bam"
 221         reference_filename = "ex1.bam"
 222
 223         self.checkEcho( input_filename, reference_filename, output_filename,
 224                         "rb", "wb", use_template = False )
 225
 226     def testReadWriteSamWithHeader( self ):
 227
 228         input_filename = "ex2.sam"
 229         output_filename = "pysam_ex2.sam"
 230         reference_filename = "ex2.sam"
 231
 232         self.checkEcho( input_filename, reference_filename, output_filename,
 233                         "r", "wh" )
 234
 235     def testReadWriteSamWithoutHeader( self ):
 236
 237         input_filename = "ex2.sam"
 238         output_filename = "pysam_ex2.sam"
 239         reference_filename = "ex1.sam"
 240
 241         self.checkEcho( input_filename, reference_filename, output_filename,
 242                         "r", "w" )
 243
 244     def testFetchFromClosedFile( self ):
 245
 246         samfile = pysam.Samfile( "ex1.bam", "rb" )
 247         samfile.close()
 248         self.assertRaises( ValueError, samfile.fetch, 'chr1', 100, 120)
 249
 250     def testPileupFromClosedFile( self ):
 251
 252         samfile = pysam.Samfile( "ex1.bam", "rb" )
 253         samfile.close()
 254         self.assertRaises( ValueError, samfile.pileup, 'chr1', 100, 120)
 255
 256     def testBinaryReadFromSamfile( self ):
 257         pass
 258         # needs to re-activated, see issue 19
 259         #samfile = pysam.Samfile( "ex1.bam", "r" )
 260         #samfile.fetch().next()
 261
 262     def testReadingFromFileWithoutIndex( self ):
 263         '''read from bam file without index.'''
 264
 265         assert not os.path.exists( "ex2.bam.bai" )
 266         samfile = pysam.Samfile( "ex2.bam", "rb" )
 267         self.assertRaises( ValueError, samfile.fetch )
 268         self.assertEqual( len(list( samfile.fetch(until_eof = True) )), 3270 )
 269
 270     def testReadingFromFileWithWrongMode( self ):
 271
 272         assert not os.path.exists( "ex2.bam.bai" )
 273         samfile = pysam.Samfile( "ex2.bam", "r" )
 274         self.assertRaises( ValueError, samfile.fetch )
 275
 276 class TestIteratorRow(unittest.TestCase):
 277
 278     def setUp(self):
 279         self.samfile=pysam.Samfile( "ex1.bam","rb" )
 280
 281     def checkRange( self, rnge ):
 282         '''compare results from iterator with those from samtools.'''
 283         ps = list(self.samfile.fetch(region=rnge))
 284         sa = list(pysam.view( "ex1.bam", rnge , raw = True) )
 285         self.assertEqual( len(ps), len(sa), "unequal number of results for range %s: %i != %i" % (rnge, len(ps), len(sa) ))
 286         # check if the same reads are returned and in the same order
 287         for line, pair in enumerate( zip( ps, sa ) ):
 288             data = pair[1].split("\t")
 289             self.assertEqual( pair[0].qname, data[0], "read id mismatch in line %i: %s != %s" % (line, pair[0].rname, data[0]) )
 290
 291     def testIteratePerContig(self):
 292         '''check random access per contig'''
 293         for contig in self.samfile.references:
 294             self.checkRange( contig )
 295
 296     def testIterateRanges(self):
 297         '''check random access per range'''
 298         for contig, length in zip(self.samfile.references, self.samfile.lengths):
 299             for start in range( 1, length, 90):
 300                 self.checkRange( "%s:%i-%i" % (contig, start, start + 90) ) # this includes empty ranges
 301
 302     def tearDown(self):
 303         self.samfile.close()
 304
 305
 306 class TestIteratorRowAll(unittest.TestCase):
 307
 308     def setUp(self):
 309         self.samfile=pysam.Samfile( "ex1.bam","rb" )
 310
 311     def testIterate(self):
 312         '''compare results from iterator with those from samtools.'''
 313         ps = list(self.samfile.fetch())
 314         sa = list(pysam.view( "ex1.bam", raw = True) )
 315         self.assertEqual( len(ps), len(sa), "unequal number of results: %i != %i" % (len(ps), len(sa) ))
 316         # check if the same reads are returned
 317         for line, pair in enumerate( zip( ps, sa ) ):
 318             data = pair[1].split("\t")
 319             self.assertEqual( pair[0].qname, data[0], "read id mismatch in line %i: %s != %s" % (line, pair[0].rname, data[0]) )
 320
 321     def tearDown(self):
 322         self.samfile.close()
 323
 324 class TestIteratorColumn(unittest.TestCase):
 325     '''test iterator column against contents of ex3.bam.'''
 326
 327     # note that samfile contains 1-based coordinates
 328     # 1D means deletion with respect to reference sequence
 329     #
 330     mCoverages = { 'chr1' : [ 0 ] * 20 + [1] * 36 + [0] * (100 - 20 -35 ),
 331                    'chr2' : [ 0 ] * 20 + [1] * 35 + [0] * (100 - 20 -35 ),
 332                    }
 333
 334     def setUp(self):
 335         self.samfile=pysam.Samfile( "ex4.bam","rb" )
 336
 337     def checkRange( self, rnge ):
 338         '''compare results from iterator with those from samtools.'''
 339         # check if the same reads are returned and in the same order
 340         for column in self.samfile.pileup(region=rnge):
 341             thiscov = len(column.pileups)
 342             refcov = self.mCoverages[self.samfile.getrname(column.tid)][column.pos]
 343             self.assertEqual( thiscov, refcov, "wrong coverage at pos %s:%i %i should be %i" % (self.samfile.getrname(column.tid), column.pos, thiscov, refcov))
 344
 345     def testIterateAll(self):
 346         '''check random access per contig'''
 347         self.checkRange( None )
 348
 349     def testIteratePerContig(self):
 350         '''check random access per contig'''
 351         for contig in self.samfile.references:
 352             self.checkRange( contig )
 353
 354     def testIterateRanges(self):
 355         '''check random access per range'''
 356         for contig, length in zip(self.samfile.references, self.samfile.lengths):
 357             for start in range( 1, length, 90):
 358                 self.checkRange( "%s:%i-%i" % (contig, start, start + 90) ) # this includes empty ranges
 359
 360     def testInverse( self ):
 361         '''test the inverse, is point-wise pileup accurate.'''
 362         for contig, refseq in self.mCoverages.items():
 363             refcolumns = sum(refseq)
 364             for pos, refcov in enumerate( refseq ):
 365                 columns = list(self.samfile.pileup( contig, pos, pos+1) )
 366                 if refcov == 0:
 367                     # if no read, no coverage
 368                     self.assertEqual( len(columns), refcov, "wrong number of pileup columns returned for position %s:%i, %i should be %i" %(contig,pos,len(columns), refcov) )
 369                 elif refcov == 1:
 370                     # one read, all columns of the read are returned
 371                     self.assertEqual( len(columns), refcolumns, "pileup incomplete at position %i: got %i, expected %i " %\
 372                                           (pos, len(columns), refcolumns))
 373
 374     def tearDown(self):
 375         self.samfile.close()
 376
 377 class TestAlignedReadFromBam(unittest.TestCase):
 378
 379     def setUp(self):
 380         self.samfile=pysam.Samfile( "ex3.bam","rb" )
 381         self.reads=list(self.samfile.fetch())
 382
 383     def testARqname(self):
 384         self.assertEqual( self.reads[0].qname, "read_28833_29006_6945", "read name mismatch in read 1: %s != %s" % (self.reads[0].qname, "read_28833_29006_6945") )
 385         self.assertEqual( self.reads[1].qname, "read_28701_28881_323b", "read name mismatch in read 2: %s != %s" % (self.reads[1].qname, "read_28701_28881_323b") )
 386
 387     def testARflag(self):
 388         self.assertEqual( self.reads[0].flag, 99, "flag mismatch in read 1: %s != %s" % (self.reads[0].flag, 99) )
 389         self.assertEqual( self.reads[1].flag, 147, "flag mismatch in read 2: %s != %s" % (self.reads[1].flag, 147) )
 390
 391     def testARrname(self):
 392         self.assertEqual( self.reads[0].rname, 0, "chromosome/target id mismatch in read 1: %s != %s" % (self.reads[0].rname, 0) )
 393         self.assertEqual( self.reads[1].rname, 1, "chromosome/target id mismatch in read 2: %s != %s" % (self.reads[1].rname, 1) )
 394
 395     def testARpos(self):
 396         self.assertEqual( self.reads[0].pos, 33-1, "mapping position mismatch in read 1: %s != %s" % (self.reads[0].pos, 33-1) )
 397         self.assertEqual( self.reads[1].pos, 88-1, "mapping position mismatch in read 2: %s != %s" % (self.reads[1].pos, 88-1) )
 398
 399     def testARmapq(self):
 400         self.assertEqual( self.reads[0].mapq, 20, "mapping quality mismatch in read 1: %s != %s" % (self.reads[0].mapq, 20) )
 401         self.assertEqual( self.reads[1].mapq, 30, "mapping quality mismatch in read 2: %s != %s" % (self.reads[1].mapq, 30) )
 402
 403     def testARcigar(self):
 404         self.assertEqual( self.reads[0].cigar, [(0, 10), (2, 1), (0, 25)], "read name length mismatch in read 1: %s != %s" % (self.reads[0].cigar, [(0, 10), (2, 1), (0, 25)]) )
 405         self.assertEqual( self.reads[1].cigar, [(0, 35)], "read name length mismatch in read 2: %s != %s" % (self.reads[1].cigar, [(0, 35)]) )
 406
 407     def testARmrnm(self):
 408         self.assertEqual( self.reads[0].mrnm, 0, "mate reference sequence name mismatch in read 1: %s != %s" % (self.reads[0].mrnm, 0) )
 409         self.assertEqual( self.reads[1].mrnm, 1, "mate reference sequence name mismatch in read 2: %s != %s" % (self.reads[1].mrnm, 1) )
 410
 411     def testARmpos(self):
 412         self.assertEqual( self.reads[0].mpos, 200-1, "mate mapping position mismatch in read 1: %s != %s" % (self.reads[0].mpos, 200-1) )
 413         self.assertEqual( self.reads[1].mpos, 500-1, "mate mapping position mismatch in read 2: %s != %s" % (self.reads[1].mpos, 500-1) )
 414
 415     def testARisize(self):
 416         self.assertEqual( self.reads[0].isize, 167, "insert size mismatch in read 1: %s != %s" % (self.reads[0].isize, 167) )
 417         self.assertEqual( self.reads[1].isize, 412, "insert size mismatch in read 2: %s != %s" % (self.reads[1].isize, 412) )
 418
 419     def testARseq(self):
 420         self.assertEqual( self.reads[0].seq, "AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG", "sequence mismatch in read 1: %s != %s" % (self.reads[0].seq, "AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG") )
 421         self.assertEqual( self.reads[1].seq, "ACCTATATCTTGGCCTTGGCCGATGCGGCCTTGCA", "sequence size mismatch in read 2: %s != %s" % (self.reads[1].seq, "ACCTATATCTTGGCCTTGGCCGATGCGGCCTTGCA") )
 422         self.assertEqual( self.reads[3].seq, "AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG", "sequence mismatch in read 4: %s != %s" % (self.reads[3].seq, "AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG") )
 423
 424     def testARqual(self):
 425         self.assertEqual( self.reads[0].qual, "<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<", "quality string mismatch in read 1: %s != %s" % (self.reads[0].qual, "<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<") )
 426         self.assertEqual( self.reads[1].qual, "<<<<<;<<<<7;:<<<6;<<<<<<<<<<<<7<<<<", "quality string mismatch in read 2: %s != %s" % (self.reads[1].qual, "<<<<<;<<<<7;:<<<6;<<<<<<<<<<<<7<<<<") )
 427         self.assertEqual( self.reads[3].qual, "<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<", "quality string mismatch in read 3: %s != %s" % (self.reads[3].qual, "<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<") )
 428
 429     def testARquery(self):
 430         self.assertEqual( self.reads[0].query, "AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG", "query mismatch in read 1: %s != %s" % (self.reads[0].query, "AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG") )
 431         self.assertEqual( self.reads[1].query, "ACCTATATCTTGGCCTTGGCCGATGCGGCCTTGCA", "query size mismatch in read 2: %s != %s" % (self.reads[1].query, "ACCTATATCTTGGCCTTGGCCGATGCGGCCTTGCA") )
 432         self.assertEqual( self.reads[3].query, "TAGCTAGCTACCTATATCTTGGTCTT", "query mismatch in read 4: %s != %s" % (self.reads[3].query, "TAGCTAGCTACCTATATCTTGGTCTT") )
 433
 434     def testARqqual(self):
 435         self.assertEqual( self.reads[0].qqual, "<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<", "qquality string mismatch in read 1: %s != %s" % (self.reads[0].qqual, "<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<") )
 436         self.assertEqual( self.reads[1].qqual, "<<<<<;<<<<7;:<<<6;<<<<<<<<<<<<7<<<<", "qquality string mismatch in read 2: %s != %s" % (self.reads[1].qqual, "<<<<<;<<<<7;:<<<6;<<<<<<<<<<<<7<<<<") )
 437         self.assertEqual( self.reads[3].qqual, "<<<<<<<<<<<<<<<<<:<9/,&,22", "qquality string mismatch in read 3: %s != %s" % (self.reads[3].qqual, "<<<<<<<<<<<<<<<<<:<9/,&,22") )
 438
 439     def testPresentOptionalFields(self):
 440         self.assertEqual( self.reads[0].opt('NM'), 1, "optional field mismatch in read 1, NM: %s != %s" % (self.reads[0].opt('NM'), 1) )
 441         self.assertEqual( self.reads[0].opt('RG'), 'L1', "optional field mismatch in read 1, RG: %s != %s" % (self.reads[0].opt('RG'), 'L1') )
 442         self.assertEqual( self.reads[1].opt('RG'), 'L2', "optional field mismatch in read 2, RG: %s != %s" % (self.reads[1].opt('RG'), 'L2') )
 443         self.assertEqual( self.reads[1].opt('MF'), 18, "optional field mismatch in read 2, MF: %s != %s" % (self.reads[1].opt('MF'), 18) )
 444
 445     def testPairedBools(self):
 446         self.assertEqual( self.reads[0].is_paired, True, "is paired mismatch in read 1: %s != %s" % (self.reads[0].is_paired, True) )
 447         self.assertEqual( self.reads[1].is_paired, True, "is paired mismatch in read 2: %s != %s" % (self.reads[1].is_paired, True) )
 448         self.assertEqual( self.reads[0].is_proper_pair, True, "is proper pair mismatch in read 1: %s != %s" % (self.reads[0].is_proper_pair, True) )
 449         self.assertEqual( self.reads[1].is_proper_pair, True, "is proper pair mismatch in read 2: %s != %s" % (self.reads[1].is_proper_pair, True) )
 450
 451     def testTags( self ):
 452         self.assertEqual( self.reads[0].tags,
 453                           [('NM', 1), ('RG', 'L1'),
 454                            ('PG', 'P1'), ('XT', 'U')] )
 455         self.assertEqual( self.reads[1].tags,
 456                           [('MF', 18), ('RG', 'L2'),
 457                            ('PG', 'P2'),('XT', 'R') ] )
 458
 459     def testOpt( self ):
 460         self.assertEqual( self.reads[0].opt("XT"), "U" )
 461         self.assertEqual( self.reads[1].opt("XT"), "R" )
 462
 463     def testMissingOpt( self ):
 464         self.assertRaises( KeyError, self.reads[0].opt, "XP" )
 465
 466     def testEmptyOpt( self ):
 467         self.assertRaises( KeyError, self.reads[2].opt, "XT" )
 468
 469     def tearDown(self):
 470         self.samfile.close()
 471
 472 class TestAlignedReadFromSam(TestAlignedReadFromBam):
 473
 474     def setUp(self):
 475         self.samfile=pysam.Samfile( "ex3.sam","r" )
 476         self.reads=list(self.samfile.fetch())
 477
 478 # needs to be implemented
 479 # class TestAlignedReadFromSamWithoutHeader(TestAlignedReadFromBam):
 480 #
 481 #     def setUp(self):
 482 #         self.samfile=pysam.Samfile( "ex7.sam","r" )
 483 #         self.reads=list(self.samfile.fetch())
 484
 485 class TestHeaderSam(unittest.TestCase):
 486
 487     header = {'SQ': [{'LN': 1575, 'SN': 'chr1'},
 488                      {'LN': 1584, 'SN': 'chr2'}],
 489               'RG': [{'LB': 'SC_1', 'ID': 'L1', 'SM': 'NA12891', 'PU': 'SC_1_10', "CN":"name:with:colon"},
 490                      {'LB': 'SC_2', 'ID': 'L2', 'SM': 'NA12891', 'PU': 'SC_2_12', "CN":"name:with:colon"}],
 491               'PG': [{'ID': 'P1', 'VN': '1.0'}, {'ID': 'P2', 'VN': '1.1'}],
 492               'HD': {'VN': '1.0'},
 493               'CO' : [ 'this is a comment', 'this is another comment'],
 494               }
 495
 496     def compareHeaders( self, a, b ):
 497         '''compare two headers a and b.'''
 498         for ak,av in a.iteritems():
 499             self.assertTrue( ak in b, "key '%s' not in '%s' " % (ak,b) )
 500             self.assertEqual( av, b[ak] )
 501
 502     def setUp(self):
 503         self.samfile=pysam.Samfile( "ex3.sam","r" )
 504
 505     def testHeaders(self):
 506         self.compareHeaders( self.header, self.samfile.header )
 507         self.compareHeaders( self.samfile.header, self.header )
 508
 509     def tearDown(self):
 510         self.samfile.close()
 511
 512 class TestHeaderBam(TestHeaderSam):
 513
 514     def setUp(self):
 515         self.samfile=pysam.Samfile( "ex3.bam","rb" )
 516
 517 class TestUnmappedReads(unittest.TestCase):
 518
 519     def testSAM(self):
 520         samfile=pysam.Samfile( "ex5.sam","r" )
 521         self.assertEqual( len(list(samfile.fetch( until_eof = True))), 2 )
 522         samfile.close()
 523
 524     def testBAM(self):
 525         samfile=pysam.Samfile( "ex5.bam","rb" )
 526         self.assertEqual( len(list(samfile.fetch( until_eof = True))), 2 )
 527         samfile.close()
 528
 529 class TestPileupObjects(unittest.TestCase):
 530
 531     def setUp(self):
 532         self.samfile=pysam.Samfile( "ex1.bam","rb" )
 533
 534     def testPileupColumn(self):
 535         for pcolumn1 in self.samfile.pileup( region="chr1:105" ):
 536             if pcolumn1.pos == 104:
 537                 self.assertEqual( pcolumn1.tid, 0, "chromosome/target id mismatch in position 1: %s != %s" % (pcolumn1.tid, 0) )
 538                 self.assertEqual( pcolumn1.pos, 105-1, "position mismatch in position 1: %s != %s" % (pcolumn1.pos, 105-1) )
 539                 self.assertEqual( pcolumn1.n, 2, "# reads mismatch in position 1: %s != %s" % (pcolumn1.n, 2) )
 540         for pcolumn2 in self.samfile.pileup( region="chr2:1480" ):
 541             if pcolumn2.pos == 1479:
 542                 self.assertEqual( pcolumn2.tid, 1, "chromosome/target id mismatch in position 1: %s != %s" % (pcolumn2.tid, 1) )
 543                 self.assertEqual( pcolumn2.pos, 1480-1, "position mismatch in position 1: %s != %s" % (pcolumn2.pos, 1480-1) )
 544                 self.assertEqual( pcolumn2.n, 12, "# reads mismatch in position 1: %s != %s" % (pcolumn2.n, 12) )
 545
 546     def testPileupRead(self):
 547         for pcolumn1 in self.samfile.pileup( region="chr1:105" ):
 548             if pcolumn1.pos == 104:
 549                 self.assertEqual( len(pcolumn1.pileups), 2, "# reads aligned to column mismatch in position 1: %s != %s" % (len(pcolumn1.pileups), 2) )
 550 #                self.assertEqual( pcolumn1.pileups[0]  # need to test additional properties here
 551
 552     def tearDown(self):
 553         self.samfile.close()
 554
 555 class TestContextManager(unittest.TestCase):
 556
 557     def testManager( self ):
 558         with pysam.Samfile('ex1.bam', 'rb') as samfile:
 559             samfile.fetch()
 560         self.assertEqual( samfile._isOpen(), False )
 561
 562 class TestExceptions(unittest.TestCase):
 563
 564     def setUp(self):
 565         self.samfile=pysam.Samfile( "ex1.bam","rb" )
 566
 567     def testMissingFile(self):
 568
 569         self.assertRaises( IOError, pysam.Samfile, "exdoesntexist.bam", "rb" )
 570         self.assertRaises( IOError, pysam.Samfile, "exdoesntexist.sam", "r" )
 571         self.assertRaises( IOError, pysam.Samfile, "exdoesntexist.bam", "r" )
 572         self.assertRaises( IOError, pysam.Samfile, "exdoesntexist.sam", "rb" )
 573
 574     def testBadContig(self):
 575         self.assertRaises( ValueError, self.samfile.fetch, "chr88" )
 576
 577     def testMeaninglessCrap(self):
 578         self.assertRaises( ValueError, self.samfile.fetch, "skljf" )
 579
 580     def testBackwardsOrderNewFormat(self):
 581         self.assertRaises( ValueError, self.samfile.fetch, 'chr1', 100, 10 )
 582
 583     def testBackwardsOrderOldFormat(self):
 584         self.assertRaises( ValueError, self.samfile.fetch, region="chr1:100-10")
 585
 586     def testOutOfRangeNegativeNewFormat(self):
 587         self.assertRaises( ValueError, self.samfile.fetch, "chr1", 5, -10 )
 588         self.assertRaises( ValueError, self.samfile.fetch, "chr1", 5, 0 )
 589         self.assertRaises( ValueError, self.samfile.fetch, "chr1", -5, -10 )
 590
 591     def testOutOfRangeNegativeOldFormat(self):
 592         self.assertRaises( ValueError, self.samfile.fetch, region="chr1:-5-10" )
 593         self.assertRaises( ValueError, self.samfile.fetch, region="chr1:-5-0" )
 594         self.assertRaises( ValueError, self.samfile.fetch, region="chr1:-5--10" )
 595
 596     def testOutOfRangNewFormat(self):
 597         self.assertRaises( ValueError, self.samfile.fetch, "chr1", 9999999999, 99999999999 )
 598
 599     def testOutOfRangeLargeNewFormat(self):
 600         self.assertRaises( ValueError, self.samfile.fetch, "chr1", 9999999999999999999999999999999, 9999999999999999999999999999999999999999 )
 601
 602     def testOutOfRangeLargeOldFormat(self):
 603         self.assertRaises( ValueError, self.samfile.fetch, "chr1:99999999999999999-999999999999999999" )
 604
 605     def tearDown(self):
 606         self.samfile.close()
 607
 608 class TestFastaFile(unittest.TestCase):
 609
 610     mSequences = { 'chr1' :
 611                        "CACTAGTGGCTCATTGTAAATGTGTGGTTTAACTCGTCCATGGCCCAGCATTAGGGAGCTGTGGACCCTGCAGCCTGGCTGTGGGGGCCGCAGTGGCTGAGGGGTGCAGAGCCGAGTCACGGGGTTGCCAGCACAGGGGCTTAACCTCTGGTGACTGCCAGAGCTGCTGGCAAGCTAGAGTCCCATTTGGAGCCCCTCTAAGCCGTTCTATTTGTAATGAAAACTATATTTATGCTATTCAGTTCTAAATATAGAAATTGAAACAGCTGTGTTTAGTGCCTTTGTTCAACCCCCTTGCAACAACCTTGAGAACCCCAGGGAATTTGTCAATGTCAGGGAAGGAGCATTTTGTCAGTTACCAAATGTGTTTATTACCAGAGGGATGGAGGGAAGAGGGACGCTGAAGAACTTTGATGCCCTCTTCTTCCAAAGATGAAACGCGTAACTGCGCTCTCATTCACTCCAGCTCCCTGTCACCCAATGGACCTGTGATATCTGGATTCTGGGAAATTCTTCATCCTGGACCCTGAGAGATTCTGCAGCCCAGCTCCAGATTGCTTGTGGTCTGACAGGCTGCAACTGTGAGCCATCACAATGAACAACAGGAAGAAAAGGTCTTTCAAAAGGTGATGTGTGTTCTCATCAACCTCATACACACACATGGTTTAGGGGTATAATACCTCTACATGGCTGATTATGAAAACAATGTTCCCCAGATACCATCCCTGTCTTACTTCCAGCTCCCCAGAGGGAAAGCTTTCAACGCTTCTAGCCATTTCTTTTGGCATTTGCCTTCAGACCCTACACGAATGCGTCTCTACCACAGGGGGCTGCGCGGTTTCCCATCATGAAGCACTGAACTTCCACGTCTCATCTAGGGGAACAGGGAGGTGCACTAATGCGCTCCACGCCCAAGCCCTTCTCACAGTTTCTGCCCCCAGCATGGTTGTACTGGGCAATACATGAGATTATTAGGAAATGCTTTACTGTCATAACTATGAAGAGACTATTGCCAGATGAACCACACATTAATACTATGTTTCTTATCTGCACATTACTACCCTGCAATTAATATAATTGTGTCCATGTACACACGCTGTCCTATGTACTTATCATGACTCTATCCCAAATTCCCAATTACGTCCTATCTTCTTCTTAGGGAAGAACAGCTTAGGTATCAATTTGGTGTTCTGTGTAAAGTCTCAGGGAGCCGTCCGTGTCCTCCCATCTGGCCTCGTCCACACTGGTTCTCTTGAAAGCTTGGGCTGTAATGATGCCCCTTGGCCATCACCCAGTCCCTGCCCCATCTCTTGTAATCTCTCTCCTTTTTGCTGCATCCCTGTCTTCCTCTGTCTTGATTTACTTGTTGTTGGTTTTCTGTTTCTTTGTTTGATTTGGTGGAAGACATAATCCCACGCTTCCTATGGAAAGGTTGTTGGGAGATTTTTAATGATTCCTCAATGTTAAAATGTCTATTTTTGTCTTGACACCCAACTAATATTTGTCTGAGCAAAACAGTCTAGATGAGAGAGAACTTCCCTGGAGGTCTGATGGCGTTTCTCCCTCGTCTTCTTA",
 612                    'chr2' :
 613                        "TTCAAATGAACTTCTGTAATTGAAAAATTCATTTAAGAAATTACAAAATATAGTTGAAAGCTCTAACAATAGACTAAACCAAGCAGAAGAAAGAGGTTCAGAACTTGAAGACAAGTCTCTTATGAATTAACCCAGTCAGACAAAAATAAAGAAAAAAATTTTAAAAATGAACAGAGCTTTCAAGAAGTATGAGATTATGTAAAGTAACTGAACCTATGAGTCACAGGTATTCCTGAGGAAAAAGAAAAAGTGAGAAGTTTGGAAAAACTATTTGAGGAAGTAATTGGGGAAAACCTCTTTAGTCTTGCTAGAGATTTAGACATCTAAATGAAAGAGGCTCAAAGAATGCCAGGAAGATACATTGCAAGACAGACTTCATCAAGATATGTAGTCATCAGACTATCTAAAGTCAACATGAAGGAAAAAAATTCTAAAATCAGCAAGAGAAAAGCATACAGTCATCTATAAAGGAAATCCCATCAGAATAACAATGGGCTTCTCAGCAGAAACCTTACAAGCCAGAAGAGATTGGATCTAATTTTTGGACTTCTTAAAGAAAAAAAAACCTGTCAAACACGAATGTTATGCCCTGCTAAACTAAGCATCATAAATGAAGGGGAAATAAAGTCAAGTCTTTCCTGACAAGCAAATGCTAAGATAATTCATCATCACTAAACCAGTCCTATAAGAAATGCTCAAAAGAATTGTAAAAGTCAAAATTAAAGTTCAATACTCACCATCATAAATACACACAAAAGTACAAAACTCACAGGTTTTATAAAACAATTGAGACTACAGAGCAACTAGGTAAAAAATTAACATTACAACAGGAACAAAACCTCATATATCAATATTAACTTTGAATAAAAAGGGATTAAATTCCCCCACTTAAGAGATATAGATTGGCAGAACAGATTTAAAAACATGAACTAACTATATGCTGTTTACAAGAAACTCATTAATAAAGACATGAGTTCAGGTAAAGGGGTGGAAAAAGATGTTCTACGCAAACAGAAACCAAATGAGAGAAGGAGTAGCTATACTTATATCAGATAAAGCACACTTTAAATCAACAACAGTAAAATAAAACAAAGGAGGTCATCATACAATGATAAAAAGATCAATTCAGCAAGAAGATATAACCATCCTACTAAATACATATGCACCTAACACAAGACTACCCAGATTCATAAAACAAATACTACTAGACCTAAGAGGGATGAGAAATTACCTAATTGGTACAATGTACAATATTCTGATGATGGTTACACTAAAAGCCCATACTTTACTGCTACTCAATATATCCATGTAACAAATCTGCGCTTGTACTTCTAAATCTATAAAAAAATTAAAATTTAACAAAAGTAAATAAAACACATAGCTAAAACTAAAAAAGCAAAAACAAAAACTATGCTAAGTATTGGTAAAGATGTGGGGAAAAAAGTAAACTCTCAAATATTGCTAGTGGGAGTATAAATTGTTTTCCACTTTGGAAAACAATTTGGTAATTTCGTTTTTTTTTTTTTCTTTTCTCTTTTTTTTTTTTTTTTTTTTGCATGCCAGAAAAAAATATTTACAGTAACT",
 614                    }
 615
 616     def setUp(self):
 617         self.file=pysam.Fastafile( "ex1.fa" )
 618
 619     def testFetch(self):
 620         for id, seq in self.mSequences.items():
 621             self.assertEqual( seq, self.file.fetch( id ) )
 622             for x in range( 0, len(seq), 10):
 623                 self.assertEqual( seq[x:x+10], self.file.fetch( id, x, x+10) )
 624                 # test x:end
 625                 self.assertEqual( seq[x:], self.file.fetch( id, x) )
 626                 # test 0:x
 627                 self.assertEqual( seq[:x], self.file.fetch( id, None, x) )
 628
 629
 630         # unknown sequence returns ""
 631         self.assertEqual( "", self.file.fetch("chr12") )
 632
 633     def testFetchErrors( self ):
 634         self.assertRaises( ValueError, self.file.fetch )
 635         self.assertRaises( ValueError, self.file.fetch, "chr1", -1, 10 )
 636         self.assertRaises( ValueError, self.file.fetch, "chr1", 20, 10 )
 637
 638     def testLength( self ):
 639         self.assertEqual( len(self.file), 2 )
 640
 641     def tearDown(self):
 642         self.file.close()
 643
 644 class TestAlignedRead(unittest.TestCase):
 645     '''tests to check if aligned read can be constructed
 646     and manipulated.
 647     '''
 648
 649     def checkFieldEqual( self, read1, read2, exclude = []):
 650         '''check if two reads are equal by comparing each field.'''
 651
 652         for x in ("qname", "seq", "flag",
 653                   "rname", "pos", "mapq", "cigar",
 654                   "mrnm", "mpos", "isize", "qual",
 655                   "is_paired", "is_proper_pair",
 656                   "is_unmapped", "mate_is_unmapped",
 657                   "is_reverse", "mate_is_reverse",
 658                   "is_read1", "is_read2",
 659                   "is_secondary", "is_qcfail",
 660                   "is_duplicate", "bin"):
 661             if x in exclude: continue
 662             self.assertEqual( getattr(read1, x), getattr(read2,x), "attribute mismatch for %s: %s != %s" %
 663                               (x, getattr(read1, x), getattr(read2,x)))
 664
 665     def testEmpty( self ):
 666         a = pysam.AlignedRead()
 667         self.assertEqual( a.qname, None )
 668         self.assertEqual( a.seq, None )
 669         self.assertEqual( a.qual, None )
 670         self.assertEqual( a.flag, 0 )
 671         self.assertEqual( a.rname, 0 )
 672         self.assertEqual( a.mapq, 0 )
 673         self.assertEqual( a.cigar, None )
 674         self.assertEqual( a.tags, None )
 675         self.assertEqual( a.mrnm, 0 )
 676         self.assertEqual( a.mpos, 0 )
 677         self.assertEqual( a.isize, 0 )
 678
 679     def buildRead( self ):
 680         '''build an example read.'''
 681
 682         a = pysam.AlignedRead()
 683         a.qname = "read_12345"
 684         a.seq="ACGT" * 3
 685         a.flag = 0
 686         a.rname = 0
 687         a.pos = 33
 688         a.mapq = 20
 689         a.cigar = ( (0,10), (2,1), (0,25) )
 690         a.mrnm = 0
 691         a.mpos=200
 692         a.isize=167
 693         a.qual="1234" * 3
 694
 695         return a
 696
 697     def testUpdate( self ):
 698         '''check if updating fields affects other variable length data
 699         '''
 700         a = self.buildRead()
 701         b = self.buildRead()
 702
 703         # check qname
 704         b.qname = "read_123"
 705         self.checkFieldEqual( a, b, "qname" )
 706         b.qname = "read_12345678"
 707         self.checkFieldEqual( a, b, "qname" )
 708         b.qname = "read_12345"
 709         self.checkFieldEqual( a, b)
 710
 711         # check cigar
 712         b.cigar = ( (0,10), )
 713         self.checkFieldEqual( a, b, "cigar" )
 714         b.cigar = ( (0,10), (2,1), (0,25), (2,1), (0,25) )
 715         self.checkFieldEqual( a, b, "cigar" )
 716         b.cigar = ( (0,10), (2,1), (0,25) )
 717         self.checkFieldEqual( a, b)
 718
 719         # check seq
 720         b.seq = "ACGT"
 721         self.checkFieldEqual( a, b, ("seq", "qual") )
 722         b.seq = "ACGT" * 10
 723         self.checkFieldEqual( a, b, ("seq", "qual") )
 724         b.seq = "ACGT" * 3
 725         self.checkFieldEqual( a, b, ("qual",))
 726
 727         # reset qual
 728         b = self.buildRead()
 729
 730         # check flags:
 731         for x in (
 732             "is_paired", "is_proper_pair",
 733             "is_unmapped", "mate_is_unmapped",
 734             "is_reverse", "mate_is_reverse",
 735             "is_read1", "is_read2",
 736             "is_secondary", "is_qcfail",
 737             "is_duplicate"):
 738             setattr( b, x, True )
 739             self.assertEqual( getattr(b, x), True )
 740             self.checkFieldEqual( a, b, ("flag", x,) )
 741             setattr( b, x, False )
 742             self.assertEqual( getattr(b, x), False )
 743             self.checkFieldEqual( a, b )
 744
 745     def testLargeRead( self ):
 746         '''build an example read.'''
 747
 748         a = pysam.AlignedRead()
 749         a.qname = "read_12345"
 750         a.seq="ACGT" * 200
 751         a.flag = 0
 752         a.rname = 0
 753         a.pos = 33
 754         a.mapq = 20
 755         a.cigar = ( (0,10), (2,1), (0,25) )
 756         a.mrnm = 0
 757         a.mpos=200
 758         a.isize=167
 759         a.qual="1234" * 200
 760
 761         return a
 762
 763 class TestDeNovoConstruction(unittest.TestCase):
 764     '''check BAM/SAM file construction using ex3.sam
 765
 766     (note these are +1 coordinates):
 767
 768     read_28833_29006_6945       99      chr1    33      20      10M1D25M        =       200     167     AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG     <<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<     NM:i:1  RG:Z:L1
 769     read_28701_28881_323b       147     chr2    88      30      35M     =       500     412     ACCTATATCTTGGCCTTGGCCGATGCGGCCTTGCA     <<<<<;<<<<7;:<<<6;<<<<<<<<<<<<7<<<<     MF:i:18 RG:Z:L2
 770     '''
 771
 772     header = { 'HD': {'VN': '1.0'},
 773                'SQ': [{'LN': 1575, 'SN': 'chr1'},
 774                       {'LN': 1584, 'SN': 'chr2'}], }
 775
 776     bamfile = "ex6.bam"
 777     samfile = "ex6.sam"
 778
 779     def checkFieldEqual( self, read1, read2, exclude = []):
 780         '''check if two reads are equal by comparing each field.'''
 781
 782         for x in ("qname", "seq", "flag",
 783                   "rname", "pos", "mapq", "cigar",
 784                   "mrnm", "mpos", "isize", "qual",
 785                   "bin",
 786                   "is_paired", "is_proper_pair",
 787                   "is_unmapped", "mate_is_unmapped",
 788                   "is_reverse", "mate_is_reverse",
 789                   "is_read1", "is_read2",
 790                   "is_secondary", "is_qcfail",
 791                   "is_duplicate"):
 792             if x in exclude: continue
 793             self.assertEqual( getattr(read1, x), getattr(read2,x), "attribute mismatch for %s: %s != %s" %
 794                               (x, getattr(read1, x), getattr(read2,x)))
 795
 796     def setUp( self ):
 797
 798
 799         a = pysam.AlignedRead()
 800         a.qname = "read_28833_29006_6945"
 801         a.seq="AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG"
 802         a.flag = 99
 803         a.rname = 0
 804         a.pos = 32
 805         a.mapq = 20
 806         a.cigar = ( (0,10), (2,1), (0,25) )
 807         a.mrnm = 0
 808         a.mpos=199
 809         a.isize=167
 810         a.qual="<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<"
 811         a.tags = ( ("NM", 1),
 812                    ("RG", "L1") )
 813
 814         b = pysam.AlignedRead()
 815         b.qname = "read_28701_28881_323b"
 816         b.seq="ACCTATATCTTGGCCTTGGCCGATGCGGCCTTGCA"
 817         b.flag = 147
 818         b.rname = 1
 819         b.pos = 87
 820         b.mapq = 30
 821         b.cigar = ( (0,35), )
 822         b.mrnm = 1
 823         b.mpos=499
 824         b.isize=412
 825         b.qual="<<<<<;<<<<7;:<<<6;<<<<<<<<<<<<7<<<<"
 826         b.tags = ( ("MF", 18),
 827                    ("RG", "L2") )
 828
 829         self.reads = (a,b)
 830
 831     def testSAMWholeFile( self ):
 832
 833         tmpfilename = "tmp_%i.sam" % id(self)
 834
 835         outfile = pysam.Samfile( tmpfilename, "wh", header = self.header )
 836
 837         for x in self.reads: outfile.write( x )
 838         outfile.close()
 839
 840         self.assertTrue( checkBinaryEqual( tmpfilename, self.samfile ),
 841                          "mismatch when construction SAM file, see %s %s" % (tmpfilename, self.samfile))
 842
 843         os.unlink( tmpfilename )
 844
 845     def testBAMPerRead( self ):
 846         '''check if individual reads are binary equal.'''
 847         infile = pysam.Samfile( self.bamfile, "rb")
 848
 849         others = list(infile)
 850         for denovo, other in zip( others, self.reads):
 851             self.checkFieldEqual( other, denovo )
 852             self.assertEqual( other.compare( denovo ), 0 )
 853
 854     def testSAMPerRead( self ):
 855         '''check if individual reads are binary equal.'''
 856         infile = pysam.Samfile( self.samfile, "r")
 857
 858         others = list(infile)
 859         for denovo, other in zip( others, self.reads):
 860             self.checkFieldEqual( other, denovo )
 861             self.assertEqual( other.compare( denovo), 0 )
 862
 863     def testBAMWholeFile( self ):
 864
 865         tmpfilename = "tmp_%i.bam" % id(self)
 866
 867         outfile = pysam.Samfile( tmpfilename, "wb", header = self.header )
 868
 869         for x in self.reads: outfile.write( x )
 870         outfile.close()
 871
 872         self.assertTrue( checkBinaryEqual( tmpfilename, self.bamfile ),
 873                          "mismatch when construction BAM file, see %s %s" % (tmpfilename, self.bamfile))
 874
 875         os.unlink( tmpfilename )
 876
 877 class TestDoubleFetch(unittest.TestCase):
 878     '''check if two iterators on the same bamfile are independent.'''
 879
 880     def testDoubleFetch( self ):
 881
 882         samfile1 = pysam.Samfile('ex1.bam', 'rb')
 883
 884         for a,b in zip(samfile1.fetch(), samfile1.fetch()):
 885             self.assertEqual( a.compare( b ), 0 )
 886
 887     def testDoubleFetchWithRegion( self ):
 888
 889         samfile1 = pysam.Samfile('ex1.bam', 'rb')
 890         chr, start, stop = 'chr1', 200, 3000000
 891         self.assertTrue(len(list(samfile1.fetch ( chr, start, stop))) > 0) #just making sure the test has something to catch
 892
 893         for a,b in zip(samfile1.fetch( chr, start, stop), samfile1.fetch( chr, start, stop)):
 894             self.assertEqual( a.compare( b ), 0 )
 895
 896     def testDoubleFetchUntilEOF( self ):
 897
 898         samfile1 = pysam.Samfile('ex1.bam', 'rb')
 899
 900         for a,b in zip(samfile1.fetch( until_eof = True),
 901                        samfile1.fetch( until_eof = True )):
 902             self.assertEqual( a.compare( b), 0 )
 903
 904 class TestRemoteFileFTP(unittest.TestCase):
 905     '''test remote access.
 906
 907     '''
 908
 909     # Need to find an ftp server without password on standard
 910     # port.
 911
 912     url = "ftp://ftp.sanger.ac.uk/pub/rd/humanSequences/CV.bam"
 913     region = "1:1-1000"
 914
 915     def testFTPView( self ):
 916         result = pysam.view( self.url, self.region )
 917         self.assertEqual( len(result), 36 )
 918
 919     def testFTPFetch( self ):
 920         samfile = pysam.Samfile(self.url, "rb")
 921         result = list(samfile.fetch( region = self.region ))
 922         self.assertEqual( len(result), 36 )
 923
 924 class TestRemoteFileHTTP( unittest.TestCase):
 925
 926     url = "http://genserv.anat.ox.ac.uk/downloads/pysam/test/ex1.bam"
 927     region = "chr1:1-1000"
 928     local = "ex1.bam"
 929
 930     def testView( self ):
 931         self.assertRaises( pysam.SamtoolsError, pysam.view, self.url, self.region )
 932
 933     def testFetch( self ):
 934         samfile = pysam.Samfile(self.url, "rb")
 935         result = list(samfile.fetch( region = self.region ))
 936         samfile_local = pysam.Samfile(self.local, "rb")
 937         ref = list(samfile_local.fetch( region = self.region ))
 938
 939         self.assertEqual( len(ref), len(result) )
 940         for x, y in zip(result, ref):
 941             self.assertEqual( x.compare( y ), 0 )
 942
 943     def testFetchAll( self ):
 944         samfile = pysam.Samfile(self.url, "rb")
 945         result = list(samfile.fetch())
 946         samfile_local = pysam.Samfile(self.local, "rb")
 947         ref = list(samfile_local.fetch() )
 948
 949         self.assertEqual( len(ref), len(result) )
 950         for x, y in zip(result, ref):
 951             self.assertEqual( x.compare( y ), 0 )
 952
 953
 954 # TODOS
 955 # 1. finish testing all properties within pileup objects
 956 # 2. check exceptions and bad input problems (missing files, optional fields that aren't present, etc...)
 957
 958 if __name__ == "__main__":
 959     # build data files
 960     print "building data files"
 961     subprocess.call( "make", shell=True)
 962     print "starting tests"
 963     unittest.main()