tests/blast_test.py

   1 from itertools import *
   2 import re
   3 import unittest
   4 import glob
   5 import os
   6 from testlib import testutil, SkipTest, PygrTestProgram
   7 from pygr import worldbase
   8 from pygr import sequence, cnestedlist, seqdb, blast, logger, parse_blast
   9 from pygr.nlmsa_utils import CoordsGroupStart,CoordsGroupEnd
  10 from pygr import translationDB
  11
  12 def check_results(results, correct, formatter, delta=0.01,
  13                   reformatCorrect=False, reformatResults=True):
  14     if reformatResults:
  15         results = reformat_results(results, formatter)
  16
  17     if reformatCorrect: # reformat these data too
  18         correct = reformat_results(correct, formatter)
  19     else:
  20         correct.sort()
  21
  22     # this is to help troubleshooting the mismatches if there are any
  23     mismatch = [ (a, b) for a, b in zip(correct, results) if
  24                  testutil.approximate_cmp([a], [b], delta)]
  25     if mismatch:
  26         logger.warn('blast mismatches found')
  27         for m in mismatch:
  28             logger.warn('%s != %s' % m)
  29
  30     # this is the actual test
  31     assert testutil.approximate_cmp(correct, results, delta) == 0
  32
  33 def check_results_relaxed_blastp(results, correct, formatter, delta=0.01,
  34                                  reformatCorrect=False, allowedLengthDiff=0,
  35                                  identityMin=0.6, reformatResults=True):
  36     if reformatResults:
  37         results = reformat_results(results, formatter)
  38
  39     if reformatCorrect: # reformat these data too
  40         correct = reformat_results(correct, formatter)
  41     else:
  42         correct.sort()
  43
  44     # Length of output
  45     assert abs(len(results) - len(correct)) <= allowedLengthDiff
  46
  47     # Format check
  48     key_re = re.compile('^[A-Z]{3}[A-Z0-9]?_[A-Z]{2,5}$')
  49     for result in results:
  50         assert key_re.search(result[0])
  51         assert key_re.search(result[1])
  52         assert (0. < result[2] and result[2] <= 1.)
  53
  54     # High-identity comparison
  55     results_high = []
  56     correct_high = []
  57     for result in results:
  58         if result[2] > identityMin:
  59             results_high.append(result)
  60     for result in correct:
  61         if result[2] > identityMin:
  62             correct_high.append(result)
  63     assert testutil.approximate_cmp(correct_high, results_high, delta) == 0
  64
  65 def check_results_relaxed_blastx(results, correct, formatter, delta=0.01,
  66                                  reformatCorrect=False, allowedLengthDiff=0,
  67                                  identityMin=0.6):
  68     results = reformat_results(results, formatter)
  69
  70     if reformatCorrect: # reformat these data too
  71         correct = reformat_results(correct, formatter)
  72     else:
  73         correct.sort()
  74
  75     # Length of output
  76     assert abs(len(results) - len(correct)) <= allowedLengthDiff
  77
  78     # Format check
  79     for result in results:
  80         assert 3 * result[0] == result[2]
  81         assert (0. < result[3] and result[3] <= 1.)
  82
  83     # High-identity comparison
  84     results_high = []
  85     correct_high = []
  86     for result in results:
  87         if result[3] > identityMin:
  88             results_high.append(result)
  89     for result in correct:
  90         if result[3] > identityMin:
  91             correct_high.append(result)
  92     assert testutil.approximate_cmp(correct_high, results_high, delta) == 0
  93
  94 def reformat_results(results, formatter):
  95     reffed = []
  96     for result in results:
  97         for t in result.edges(mergeMost=True):
  98             reffed.append(formatter(t))
  99     reffed.sort()
 100     return reffed
 101
 102 def pair_identity_tuple(t):
 103     'standard formatter for blast matches'
 104     return (t[0].id, t[1].id, t[2].pIdentity())
 105
 106 class BlastBase(unittest.TestCase):
 107     def setUp(self):
 108         hbb1_mouse = testutil.datafile('hbb1_mouse.fa')
 109         hbb1_mouse_rc = testutil.datafile('hbb1_mouse_rc.fa')
 110         sp_hbb1 = testutil.datafile('sp_hbb1')
 111         gapping = testutil.datafile('gapping.fa')
 112
 113         self.dna = seqdb.SequenceFileDB(hbb1_mouse)
 114         self.dna_rc = seqdb.SequenceFileDB(hbb1_mouse_rc)
 115         self.prot = seqdb.SequenceFileDB(sp_hbb1)
 116         self.gapping = seqdb.SequenceFileDB(gapping)
 117
 118     def tearDown(self):
 119         'do the RIGHT thing... close resources that have been opened!'
 120         self.dna.close()
 121         self.dna_rc.close()
 122         self.prot.close()
 123         self.gapping.close()
 124
 125
 126 _multiblast_results = None
 127
 128 class Blast_Test(BlastBase):
 129     """
 130     Test basic BLAST stuff (using blastp).
 131     """
 132     def test_blastp(self):
 133         "Testing blastp"
 134         if not testutil.blast_enabled():
 135             raise SkipTest, "no BLAST installed"
 136
 137         blastmap = blast.BlastMapping(self.prot, verbose=False)
 138         results = blastmap[self.prot['HBB1_XENLA']]
 139
 140         check_results_relaxed_blastp([results], blastp_correct_results,
 141                                      pair_identity_tuple,
 142                                      allowedLengthDiff=2)
 143
 144     def test_repr(self):
 145         blastmap = blast.BlastMapping(self.prot, verbose=False)
 146         assert '<BlastMapping' in repr(blastmap)
 147
 148     def test_no_query(self):
 149         blastmap = blast.BlastMapping(self.dna, verbose=False)
 150         try:
 151             blastmap()
 152             assert 0, "should fail before this"
 153         except ValueError:
 154             pass
 155
 156     def test_both_seq_and_db(self):
 157         "Testing db arg present"
 158         blastmap = blast.BlastMapping(self.prot, verbose=False)
 159         seq = self.prot['HBB1_XENLA']
 160
 161         try:
 162             blastmap(seq=seq, queryDB=self.prot)
 163             assert 0, "should fail before this"
 164         except ValueError:
 165             pass
 166
 167     def test_multiblast(self):
 168         "testing multi sequence blast"
 169         results = self.get_multiblast_results()
 170         check_results_relaxed_blastp(results, correct_multiblast_results,
 171                                      None, reformatResults=False,
 172                                      allowedLengthDiff=10)
 173
 174     def get_multiblast_results(self):
 175         """return saved results or generate them if needed;
 176         results are saved so we only do this time-consuming operation once"""
 177         global _multiblast_results
 178
 179         if not testutil.blast_enabled():
 180             raise SkipTest, "no BLAST installed"
 181
 182         if not _multiblast_results:
 183             logger.info("running expensive multiblast")
 184             blastmap = blast.BlastMapping(self.prot, verbose=False)
 185             al = cnestedlist.NLMSA('blasthits', 'memory', pairwiseMode=True,
 186                                    bidirectional=False)
 187
 188             blastmap(al=al, queryDB=self.prot) # all vs all
 189
 190             al.build() # construct the alignment indexes
 191             results = [al[seq] for seq in self.prot.values()]
 192             _multiblast_results = reformat_results(results,
 193                                                    pair_identity_tuple)
 194
 195         return _multiblast_results
 196
 197     def test_multiblast_single(self):
 198         "Test multi-sequence BLAST results, for BLASTs run one by one."
 199         if not testutil.blast_enabled():
 200             raise SkipTest, "no BLAST installed"
 201
 202         blastmap = blast.BlastMapping(self.prot, verbose=False)
 203         al = cnestedlist.NLMSA('blasthits', 'memory', pairwiseMode=True,
 204                                bidirectional=False)
 205
 206         for seq in self.prot.values():
 207             blastmap(seq, al) # all vs all, one by one
 208
 209         al.build() # construct the alignment indexes
 210         results = [al[seq] for seq in self.prot.values()]
 211         results_multi = self.get_multiblast_results()
 212         # Strict check must work here even on live BLAST output
 213         check_results(results, results_multi, pair_identity_tuple)
 214
 215     def test_multiblast_long(self):
 216         "testing multi sequence blast with long db"
 217         if not testutil.blast_enabled():
 218             raise SkipTest, "no BLAST installed"
 219
 220         longerFile = testutil.datafile('sp_all_hbb')
 221
 222         sp_all_hbb = seqdb.SequenceFileDB(longerFile)
 223         blastmap = blast.BlastMapping(self.prot, verbose=False)
 224         al = cnestedlist.NLMSA('blasthits', 'memory', pairwiseMode=True,
 225                                bidirectional=False)
 226         blastmap(None, al, queryDB=sp_all_hbb) # all vs all
 227         al.build() # construct the alignment indexes
 228
 229     def test_maskEnd(self):
 230         """
 231         This tests against a minor bug in cnestedlist where maskEnd
 232         is used to clip the end to the mask region.
 233         """
 234         if not testutil.blast_enabled():
 235             raise SkipTest, "no BLAST installed"
 236
 237         db = self.gapping
 238         blastmap = blast.BlastMapping(db)
 239         ungapped = db['ungapped']
 240         gapped = db['gapped']
 241         results = blastmap[gapped]
 242
 243         results[ungapped]
 244
 245     def test_no_bidirectional(self):
 246         if not testutil.blast_enabled():
 247             raise SkipTest, "no BLAST installed"
 248
 249         db = self.gapping
 250         gapped = db['gapped']
 251         ungapped = db['ungapped']
 252
 253         blastmap = blast.BlastMapping(db)
 254         al = blastmap(queryDB=db)
 255         slice = al[gapped]
 256
 257         found_once = False
 258         for src, dest, edge in al[gapped].edges():
 259             if src == gapped[0:40] and dest == ungapped[0:40]:
 260                 assert not found_once, \
 261                        "BLAST results should not be bidirectional"
 262                 found_once = True
 263
 264         assert found_once, "should have found this match exactly once!"
 265
 266     def test_formatdb_fail(self):
 267         db = self.gapping
 268         try:
 269             blastmap = blast.BlastMapping(db, filepath='foobarbaz.fa',
 270                                           blastReady=True,
 271                                           showFormatdbMessages=False)
 272             assert 0, "should not reach this point"
 273         except IOError:                 # should fail with 'cannot build'
 274             pass
 275
 276         remnants = glob.glob('foobarbaz.fa.n??')
 277         for filename in remnants:
 278             os.unlink(filename)
 279
 280     def test_seq_without_db(self):
 281         "Check that sequences without associated DBs work as query strings"
 282         if not testutil.blast_enabled():
 283             raise SkipTest, "no BLAST installed"
 284         blastmap = blast.BlastMapping(self.prot, verbose=False)
 285
 286         seq = self.prot['HBB1_XENLA']
 287         seq_no_db = sequence.Sequence(str(seq), 'HBB1_XENLA_no_db')
 288         slice = blastmap(seq=seq_no_db)[seq_no_db]
 289         assert len(slice)
 290
 291 class Blastx_Test(BlastBase):
 292     def test_blastx(self):
 293         "Testing blastx"
 294         if not testutil.blast_enabled():
 295             raise SkipTest, "no BLAST installed"
 296
 297         blastmap = blast.BlastxMapping(self.prot, verbose=False)
 298
 299         correct = [(143, 143, 429, 0.53146853146853146),
 300                    (143, 145, 429, 0.28275862068965518),
 301                    (143, 145, 429, 0.28965517241379313),
 302                    (143, 145, 429, 0.29655172413793102),
 303                    (143, 145, 429, 0.30344827586206896),
 304                    (144, 144, 432, 0.4513888888888889),
 305                    (144, 144, 432, 0.4513888888888889),
 306                    (145, 145, 435, 0.45517241379310347),
 307                    (145, 145, 435, 0.51034482758620692),
 308                    (146, 142, 438, 0.35616438356164382),
 309                    (146, 146, 438, 0.4589041095890411),
 310                    (146, 146, 438, 0.46575342465753422),
 311                    (146, 146, 438, 0.4726027397260274),
 312                    (146, 146, 438, 0.4726027397260274),
 313                    (146, 146, 438, 0.4863013698630137),
 314                    (146, 146, 438, 0.59589041095890416),
 315                    (146, 146, 438, 0.62328767123287676),
 316                    (146, 146, 438, 0.66438356164383561),
 317                    (146, 146, 438, 0.74657534246575341),
 318                    (146, 146, 438, 0.91095890410958902),
 319                    (146, 146, 438, 0.97945205479452058)]
 320
 321         results = blastmap[self.dna['gi|171854975|dbj|AB364477.1|']]
 322         check_results_relaxed_blastx(results, correct,
 323                       lambda t:(len(t[0]), len(t[1]), len(t[0].sequence),
 324                                 t[2].pIdentity()), allowedLengthDiff=2)
 325
 326     def test_blastx_rc(self):
 327         "Testing blastx with negative frames"
 328         if not testutil.blast_enabled():
 329             raise SkipTest, "no BLAST installed"
 330
 331         blastmap = blast.BlastxMapping(self.prot, verbose=False)
 332
 333         correct = [(143, 143, 429, 0.53146853146853146),
 334                    (143, 145, 429, 0.28275862068965518),
 335                    (143, 145, 429, 0.28965517241379313),
 336                    (143, 145, 429, 0.29655172413793102),
 337                    (143, 145, 429, 0.30344827586206896),
 338                    (144, 144, 432, 0.4513888888888889),
 339                    (144, 144, 432, 0.4513888888888889),
 340                    (145, 145, 435, 0.45517241379310347),
 341                    (145, 145, 435, 0.51034482758620692),
 342                    (146, 142, 438, 0.35616438356164382),
 343                    (146, 146, 438, 0.4589041095890411),
 344                    (146, 146, 438, 0.46575342465753422),
 345                    (146, 146, 438, 0.4726027397260274),
 346                    (146, 146, 438, 0.4726027397260274),
 347                    (146, 146, 438, 0.4863013698630137),
 348                    (146, 146, 438, 0.59589041095890416),
 349                    (146, 146, 438, 0.62328767123287676),
 350                    (146, 146, 438, 0.66438356164383561),
 351                    (146, 146, 438, 0.74657534246575341),
 352                    (146, 146, 438, 0.91095890410958902),
 353                    (146, 146, 438, 0.97945205479452058)]
 354
 355         results = blastmap[self.dna_rc['hbb1_mouse_RC']]
 356         check_results_relaxed_blastx(results, correct,
 357                       lambda t:(len(t[0]), len(t[1]), len(t[0].sequence),
 358                                 t[2].pIdentity()), allowedLengthDiff=2)
 359
 360         results = blastmap[self.dna_rc['hbb1_mouse_RC_2']]
 361         check_results_relaxed_blastx(results, correct,
 362                       lambda t:(len(t[0]), len(t[1]), len(t[0].sequence),
 363                                 t[2].pIdentity()), allowedLengthDiff=2)
 364
 365         results = blastmap[self.dna_rc['hbb1_mouse_RC_3']]
 366         check_results_relaxed_blastx(results, correct,
 367                       lambda t:(len(t[0]), len(t[1]), len(t[0].sequence),
 368                                 t[2].pIdentity()), allowedLengthDiff=2)
 369
 370     def test_repr(self):
 371         blastmap = blast.BlastxMapping(self.prot, verbose=False)
 372         assert '<BlastxMapping' in repr(blastmap)
 373
 374     def test_blastx_no_blastp(self):
 375         if not testutil.blast_enabled():
 376             raise SkipTest, "no BLAST installed"
 377
 378         blastmap = blast.BlastxMapping(self.prot, verbose=False)
 379
 380         try:
 381             results = blastmap(self.prot['HBB1_MOUSE'])
 382             raise AssertionError('failed to trap blastp in BlastxMapping')
 383         except ValueError:
 384             pass
 385
 386     def test_no_query(self):
 387         blastmap = blast.BlastxMapping(self.prot, verbose=False)
 388         try:
 389             blastmap()
 390             assert 0, "should fail before this"
 391         except ValueError:
 392             pass
 393
 394     def test_both_seq_and_db(self):
 395         "Testing blastp"
 396         blastmap = blast.BlastxMapping(self.prot, verbose=False)
 397         seq = self.prot['HBB1_XENLA']
 398
 399         try:
 400             blastmap(seq=seq, queryDB=self.prot)
 401             assert 0, "should fail before this"
 402         except ValueError:
 403             pass
 404
 405     def test_translation_db_in_results_of_db_search(self):
 406         """
 407         Test that the NLMSA in a BlastxMapping properly picks up the
 408         translationDB from the query sequence dict.
 409         """
 410         if not testutil.blast_enabled():
 411             raise SkipTest, "no BLAST installed"
 412         blastmap = blast.BlastxMapping(self.prot, verbose=False)
 413         results = blastmap(queryDB=self.dna)
 414
 415         tdb = translationDB.get_translation_db(self.dna)
 416         assert tdb.annodb in results.seqDict.dicts
 417
 418     def test_translation_db_in_results_of_seq_search(self):
 419         """
 420         Test that the NLMSA in a BlastxMapping properly picks up the
 421         translationDB from a single input sequence.
 422         """
 423         if not testutil.blast_enabled():
 424             raise SkipTest, "no BLAST installed"
 425         blastmap = blast.BlastxMapping(self.prot, verbose=False)
 426
 427         query_seq = self.dna['gi|171854975|dbj|AB364477.1|']
 428         results = blastmap(seq=query_seq)
 429
 430         tdb = translationDB.get_translation_db(self.dna)
 431         assert tdb.annodb in results.seqDict.dicts
 432
 433     def test_translated_seqs_in_results(self):
 434         """
 435         Only NLMSASlices for the query sequence should show up in
 436         BlastxMapping.__getitem__, right?
 437         """
 438         if not testutil.blast_enabled():
 439             raise SkipTest, "no BLAST installed"
 440         blastmap = blast.BlastxMapping(self.prot, verbose=False)
 441
 442         query_seq = self.dna['gi|171854975|dbj|AB364477.1|']
 443         results = blastmap[query_seq]
 444
 445         tdb = translationDB.get_translation_db(self.dna)
 446         annodb = tdb.annodb
 447
 448         for slice in results:
 449             assert slice.seq.id in annodb, '%s not in annodb!' % slice.seq.id
 450
 451     def test_non_consumable_results(self):
 452         if not testutil.blast_enabled():
 453             raise SkipTest, "no BLAST installed"
 454         blastmap = blast.BlastxMapping(self.prot, verbose=False)
 455
 456         query_seq = self.dna['gi|171854975|dbj|AB364477.1|']
 457         results = blastmap[query_seq]
 458
 459         x = list(results)
 460         y = list(results)
 461
 462         assert len(x), x
 463         assert x == y, "BlastxMapping.__getitem__ should return list"
 464
 465 class Tblastn_Test(BlastBase):
 466     def test_tblastn(self):
 467         "tblastn test"
 468         if not testutil.blast_enabled():
 469             raise SkipTest, "no BLAST installed"
 470
 471         blastmap = blast.BlastMapping(self.dna, verbose=False)
 472         correct = [(144, 144, 432, 0.451)]
 473
 474         result = blastmap[self.prot['HBB1_XENLA']]
 475         check_results_relaxed_blastx([result], correct,
 476                       lambda t:(len(t[1]), len(t[0]), len(t[1].sequence),
 477                                 t[2].pIdentity()))
 478
 479     def test_tblastn_no_blastx(self):
 480         if not testutil.blast_enabled():
 481             raise SkipTest, "no BLAST installed"
 482
 483         blastmap = blast.BlastMapping(self.prot)
 484         try:
 485             results = blastmap[self.dna['gi|171854975|dbj|AB364477.1|']]
 486             raise AssertionError('failed to trap blastx in BlastMapping')
 487         except ValueError:
 488             pass
 489
 490     def test_megablast(self):
 491         '''test megablast'''
 492         if not testutil.blast_enabled():
 493             raise SkipTest, "no BLAST installed"
 494
 495         blastmap = blast.MegablastMapping(self.dna, verbose=False)
 496         # must use copy of sequence to get "self matches" from NLMSA...
 497         query = seqdb.Sequence(str(self.dna['gi|171854975|dbj|AB364477.1|']),
 498                                'foo')
 499         try:
 500             result = blastmap[query]
 501         except OSError: # silently ignore missing RepeatMasker, megablast
 502             return
 503         found = [(len(t[0]), len(t[1])) for t in result.edges()]
 504         assert found == [(444, 444)]
 505
 506     def test_megablast_repr(self):
 507         blastmap = blast.MegablastMapping(self.dna, verbose=False)
 508         assert '<MegablastMapping' in repr(blastmap)
 509
 510     def test_bad_subject(self):
 511         "Test bad subjects"
 512
 513         correctCoords = ((12, 63, 99508, 99661),
 514                          (65, 96, 99661, 99754),
 515                          (96, 108, 99778, 99814),
 516                          (108, 181, 99826, 100045))
 517
 518         fp = file(testutil.datafile('bad_tblastn.txt'))
 519         try:
 520             p = parse_blast.BlastHitParser()
 521             it = iter(correctCoords)
 522             for ival in p.parse_file(fp):
 523                 if not isinstance(ival, (CoordsGroupStart, CoordsGroupEnd)):
 524                     assert (ival.src_start,ival.src_end,
 525                             ival.dest_start,ival.dest_end) \
 526                         == it.next()
 527         finally:
 528             fp.close()
 529
 530 class BlastParsers_Test(BlastBase):
 531     def test_blastp_parser(self):
 532         "Testing blastp parser"
 533         blastp_output = open(testutil.datafile('blastp_output.txt'), 'r')
 534
 535         seq_dict = { 'HBB1_XENLA' : self.prot['HBB1_XENLA'] }
 536         prot_index = blast.BlastIDIndex(self.prot)
 537         try:
 538             alignment = blast.read_blast_alignment(blastp_output, seq_dict,
 539                                                    prot_index)
 540             results = alignment[self.prot['HBB1_XENLA']]
 541         finally:
 542             blastp_output.close()
 543
 544         check_results([results], blastp_correct_results, pair_identity_tuple)
 545
 546     def test_multiblast_parser(self):
 547         "Testing multiblast parser"
 548         multiblast_output = open(testutil.datafile('multiblast_output.txt'),
 549                                  'r')
 550
 551         try:
 552             al = cnestedlist.NLMSA('blasthits', 'memory', pairwiseMode=True,
 553                                    bidirectional=False)
 554             al = blast.read_blast_alignment(multiblast_output, self.prot,
 555                                             blast.BlastIDIndex(self.prot), al)
 556         finally:
 557             multiblast_output.close()
 558         al.build()
 559         results = [al[seq] for seq in self.prot.values()]
 560
 561         check_results(results, correct_multiblast_results,
 562                       pair_identity_tuple)
 563
 564     def test_multiblast_parser_long(self):
 565         "Testing multiblast parser with long input"
 566         longerFile = testutil.datafile('sp_all_hbb')
 567         sp_all_hbb = seqdb.SequenceFileDB(longerFile)
 568
 569         filename = testutil.datafile('multiblast_long_output.txt')
 570         multiblast_output = open(filename, 'r')
 571         try:
 572             al = cnestedlist.NLMSA('blasthits', 'memory', pairwiseMode=True,
 573                                    bidirectional=False)
 574             al = blast.read_blast_alignment(multiblast_output, sp_all_hbb,
 575                                             self.prot, al)
 576         finally:
 577             multiblast_output.close()
 578         al.build()
 579
 580         results = []
 581         for seq in sp_all_hbb.values():
 582             try:
 583                 results.append(al[seq])
 584             except KeyError:
 585                 pass
 586         correctfile = file(testutil.datafile('multiblast_long_correct.txt'),
 587                            'r')
 588         try:
 589             correct = []
 590             for line in correctfile:
 591                 t = line.split()
 592                 correct.append((t[0], t[1], float(t[2])))
 593         finally:
 594             correctfile.close()
 595         check_results(results, correct, pair_identity_tuple)
 596
 597     def test_blastx_parser(self):
 598         "Testing blastx parser"
 599         blastx_output = open(testutil.datafile('blastx_output.txt'), 'r')
 600         seq_dict =  { 'gi|171854975|dbj|AB364477.1|' :
 601                       self.dna['gi|171854975|dbj|AB364477.1|'] }
 602         try:
 603             results = blast.read_blast_alignment(blastx_output,
 604                                                  seq_dict,
 605                                                  blast.BlastIDIndex(self.prot),
 606                                                  translateSrc=True)
 607         finally:
 608             blastx_output.close()
 609         correct = [(143, 143, 429, 0.53146853146853146),
 610                    (143, 145, 429, 0.28275862068965518),
 611                    (143, 145, 429, 0.28965517241379313),
 612                    (143, 145, 429, 0.29655172413793102),
 613                    (143, 145, 429, 0.30344827586206896),
 614                    (144, 144, 432, 0.4513888888888889),
 615                    (144, 144, 432, 0.4513888888888889),
 616                    (145, 145, 435, 0.45517241379310347),
 617                    (145, 145, 435, 0.51034482758620692),
 618                    (146, 142, 438, 0.35616438356164382),
 619                    (146, 146, 438, 0.4589041095890411),
 620                    (146, 146, 438, 0.46575342465753422),
 621                    (146, 146, 438, 0.4726027397260274),
 622                    (146, 146, 438, 0.4726027397260274),
 623                    (146, 146, 438, 0.4863013698630137),
 624                    (146, 146, 438, 0.59589041095890416),
 625                    (146, 146, 438, 0.62328767123287676),
 626                    (146, 146, 438, 0.66438356164383561),
 627                    (146, 146, 438, 0.74657534246575341),
 628                    (146, 146, 438, 0.91095890410958902),
 629                    (146, 146, 438, 0.97945205479452058)]
 630
 631         check_results([results], correct,
 632                       lambda t:(len(t[0]), len(t[1]), len(t[0].sequence),
 633                                 t[2].pIdentity()))
 634
 635     def test_tblastn_parser(self):
 636         "Testing tblastn parser"
 637         seq_dict = { 'HBB1_XENLA' : self.prot['HBB1_XENLA'] }
 638         dna_db = blast.BlastIDIndex(self.dna)
 639         tblastn_output = open(testutil.datafile('tblastn_output.txt'), 'r')
 640         try:
 641             al = blast.read_blast_alignment(tblastn_output, seq_dict,
 642                                             dna_db, translateDest=True)
 643             result = al[self.prot['HBB1_XENLA']]
 644         finally:
 645             tblastn_output.close()
 646         src, dest, edge = iter(result.edges()).next()
 647
 648         self.assertEqual(str(src),
 649             'LTAHDRQLINSTWGKLCAKTIGQEALGRLLWTYPWTQRYFSSFGNLNSADAVFHNEAVAAHGEK'
 650             'VVTSIGEAIKHMDDIKGYYAQLSKYHSETLHVDPLNFKRFGGCLSIALARHFHEEYTPELHAAY'
 651             'EHLFDAIADALGKGYH')
 652         self.assertEqual(str(dest),
 653             'LTDAEKAAVSGLWGKVNSDEVGGEALGRLLVVYPWTQRYFDSFGDLSSASAIMGNAKVKAHGKK'
 654             'VITAFNEGLNHLDSLKGTFASLSELHCDKLHVDPENFRLLGNMIVIVLGHHLGKDFTPAAQAAF'
 655             'QKVMAGVATALAHKYH')
 656         self.assertEqual(str(dest.sequence),
 657             'CTGACTGATGCTGAGAAGGCTGCTGTCTCTGGCCTGTGGGGAAAGGTGAACTCCGATGAAGTTG'
 658             'GTGGTGAGGCCCTGGGCAGGCTGCTGGTTGTCTACCCTTGGACCCAGAGGTACTTTGATAGCTT'
 659             'TGGAGACCTATCCTCTGCCTCTGCTATCATGGGTAATGCCAAAGTGAAGGCCCATGGCAAGAAA'
 660             'GTGATAACTGCCTTTAACGAGGGCCTGAATCACTTGGACAGCCTCAAGGGCACCTTTGCCAGCC'
 661             'TCAGTGAGCTCCACTGTGACAAGCTCCATGTGGATCCTGAGAACTTCAGGCTCCTGGGCAATAT'
 662             'GATCGTGATTGTGCTGGGCCACCACCTGGGCAAGGATTTCACCCCCGCTGCACAGGCTGCCTTC'
 663             'CAGAAGGTGATGGCTGGAGTGGCCACTGCCCTGGCTCACAAGTACCAC')
 664
 665         self.assertAlmostEqual(edge.pIdentity(), 0.451, 3)
 666
 667
 668 # not used currently
 669 def all_vs_all_blast_save():
 670     """
 671     Creates the blast files used during testing.
 672     Must be called before running the tests
 673     """
 674
 675     tempdir = testutil.TempDir('blast-test')
 676     testutil.change_pygrdatapath(tempdir.path)
 677
 678     sp_hbb1 = testutil.datafile('sp_hbb1')
 679     all_vs_all = testutil.tempdatafile('all_vs_all')
 680
 681     sp = seqdb.BlastDB(sp_hbb1)
 682     msa = cnestedlist.NLMSA(all_vs_all ,mode='w', pairwiseMode=True, bidirectional=False)
 683
 684     # get strong homologs, save alignment in msa for every sequence
 685     reader = islice(sp.iteritems(), None)
 686     for id, s in reader:
 687         sp.blast(s, msa, expmax=1e-10, verbose=False)
 688
 689     # done constructing the alignment, so build the alignment db indexes
 690     msa.build(saveSeqDict=True)
 691
 692     db = msa.seqDict.dicts.keys()[0]
 693     working, result = {}, {}
 694     for k in db.values():
 695         edges = msa[k].edges(minAlignSize=12, pIdentityMin=0.5)
 696         for t in edges:
 697             assert len(t[0]) >= 12
 698         tmpdict = dict(map(lambda x:(x, None), [(str(t[0]), str(t[1]), t[2].pIdentity(trapOverflow=False)) for t in edges]))
 699         result[repr(k)] = tmpdict.keys()
 700         result[repr(k)].sort()
 701
 702     # save it into worldbase
 703     data = testutil.TestData()
 704     data.__doc__ = 'sp_allvall'
 705     data.result = result
 706     worldbase.Bio.Blast = data
 707     worldbase.commit()
 708
 709     #return msa
 710
 711 ###
 712
 713 blastp_correct_results = \
 714         [('HBB1_XENLA', 'HBB0_PAGBO', 0.44055944055944057),
 715          ('HBB1_XENLA', 'HBB1_ANAMI', 0.45323741007194246),
 716          ('HBB1_XENLA', 'HBB1_CYGMA', 0.46715328467153283),
 717          ('HBB1_XENLA', 'HBB1_IGUIG', 0.48951048951048953),
 718          ('HBB1_XENLA', 'HBB1_MOUSE', 0.44444444444444442),
 719          ('HBB1_XENLA', 'HBB1_ONCMY', 0.39436619718309857),
 720          ('HBB1_XENLA', 'HBB1_PAGBO', 0.44055944055944057),
 721          ('HBB1_XENLA', 'HBB1_RAT', 0.45833333333333331),
 722          ('HBB1_XENLA', 'HBB1_SPHPU', 0.4825174825174825),
 723          ('HBB1_XENLA', 'HBB1_TAPTE', 0.47222222222222221),
 724          ('HBB1_XENLA', 'HBB1_TORMA', 0.33793103448275863),
 725          ('HBB1_XENLA', 'HBB1_TRICR', 0.49305555555555558),
 726          ('HBB1_XENLA', 'HBB1_UROHA', 0.3776223776223776),
 727          ('HBB1_XENLA', 'HBB1_VAREX', 0.5174825174825175),
 728          ('HBB1_XENLA', 'HBB1_XENBO', 0.96551724137931039),
 729          ('HBB1_XENLA', 'HBB1_XENLA', 1.0),
 730          ('HBB1_XENLA', 'HBB1_XENTR', 0.75),
 731          ('HBB1_XENLA', 'MYG_DIDMA', 0.32075471698113206),
 732          ('HBB1_XENLA', 'MYG_ELEMA', 0.26415094339622641),
 733          ('HBB1_XENLA', 'MYG_ERIEU', 0.32075471698113206),
 734          ('HBB1_XENLA', 'MYG_ESCGI', 0.28282828282828282),
 735          ('HBB1_XENLA', 'MYG_GALCR', 0.32075471698113206)]
 736
 737
 738
 739
 740
 741 correct_multiblast_results = \
 742         [('HBB0_PAGBO', 'HBB0_PAGBO', 1.0),
 743          ('HBB0_PAGBO', 'HBB1_ANAMI', 0.66896551724137931),
 744          ('HBB0_PAGBO', 'HBB1_CYGMA', 0.68493150684931503),
 745          ('HBB0_PAGBO', 'HBB1_IGUIG', 0.4863013698630137),
 746          ('HBB0_PAGBO', 'HBB1_MOUSE', 0.45205479452054792),
 747          ('HBB0_PAGBO', 'HBB1_ONCMY', 0.55172413793103448),
 748          ('HBB0_PAGBO', 'HBB1_PAGBO', 0.69178082191780821),
 749          ('HBB0_PAGBO', 'HBB1_RAT', 0.4589041095890411),
 750          ('HBB0_PAGBO', 'HBB1_SPHPU', 0.4589041095890411),
 751          ('HBB0_PAGBO', 'HBB1_TAPTE', 0.4863013698630137),
 752          ('HBB0_PAGBO', 'HBB1_TORMA', 0.31506849315068491),
 753          ('HBB0_PAGBO', 'HBB1_TRICR', 0.4375),
 754          ('HBB0_PAGBO', 'HBB1_UROHA', 0.4041095890410959),
 755          ('HBB0_PAGBO', 'HBB1_VAREX', 0.49315068493150682),
 756          ('HBB0_PAGBO', 'HBB1_XENBO', 0.43356643356643354),
 757          ('HBB0_PAGBO', 'HBB1_XENLA', 0.44055944055944057),
 758          ('HBB0_PAGBO', 'HBB1_XENTR', 0.4726027397260274),
 759          ('HBB0_PAGBO', 'MYG_DIDMA', 0.22222222222222221),
 760          ('HBB0_PAGBO', 'MYG_ELEMA', 0.20833333333333334),
 761          ('HBB0_PAGBO', 'MYG_ERIEU', 0.21527777777777779),
 762          ('HBB0_PAGBO', 'MYG_ESCGI', 0.25),
 763          ('HBB0_PAGBO', 'MYG_GALCR', 0.24305555555555555),
 764          ('HBB1_ANAMI', 'HBB0_PAGBO', 0.66896551724137931),
 765          ('HBB1_ANAMI', 'HBB1_ANAMI', 1.0),
 766          ('HBB1_ANAMI', 'HBB1_CYGMA', 0.75862068965517238),
 767          ('HBB1_ANAMI', 'HBB1_IGUIG', 0.47586206896551725),
 768          ('HBB1_ANAMI', 'HBB1_MOUSE', 0.45517241379310347),
 769          ('HBB1_ANAMI', 'HBB1_ONCMY', 0.59310344827586203),
 770          ('HBB1_ANAMI', 'HBB1_PAGBO', 0.75862068965517238),
 771          ('HBB1_ANAMI', 'HBB1_RAT', 0.48965517241379308),
 772          ('HBB1_ANAMI', 'HBB1_SPHPU', 0.46206896551724136),
 773          ('HBB1_ANAMI', 'HBB1_TAPTE', 0.48965517241379308),
 774          ('HBB1_ANAMI', 'HBB1_TORMA', 0.32413793103448274),
 775          ('HBB1_ANAMI', 'HBB1_TRICR', 0.41258741258741261),
 776          ('HBB1_ANAMI', 'HBB1_UROHA', 0.38620689655172413),
 777          ('HBB1_ANAMI', 'HBB1_VAREX', 0.48275862068965519),
 778          ('HBB1_ANAMI', 'HBB1_XENBO', 0.4460431654676259),
 779          ('HBB1_ANAMI', 'HBB1_XENLA', 0.45323741007194246),
 780          ('HBB1_ANAMI', 'HBB1_XENTR', 0.4689655172413793),
 781          ('HBB1_CYGMA', 'HBB0_PAGBO', 0.68493150684931503),
 782          ('HBB1_CYGMA', 'HBB1_ANAMI', 0.75862068965517238),
 783          ('HBB1_CYGMA', 'HBB1_CYGMA', 1.0),
 784          ('HBB1_CYGMA', 'HBB1_IGUIG', 0.5),
 785          ('HBB1_CYGMA', 'HBB1_MOUSE', 0.47945205479452052),
 786          ('HBB1_CYGMA', 'HBB1_ONCMY', 0.53103448275862064),
 787          ('HBB1_CYGMA', 'HBB1_PAGBO', 0.86986301369863017),
 788          ('HBB1_CYGMA', 'HBB1_RAT', 0.50684931506849318),
 789          ('HBB1_CYGMA', 'HBB1_SPHPU', 0.47945205479452052),
 790          ('HBB1_CYGMA', 'HBB1_TAPTE', 0.4726027397260274),
 791          ('HBB1_CYGMA', 'HBB1_TORMA', 0.33561643835616439),
 792          ('HBB1_CYGMA', 'HBB1_TRICR', 0.4375),
 793          ('HBB1_CYGMA', 'HBB1_UROHA', 0.36986301369863012),
 794          ('HBB1_CYGMA', 'HBB1_VAREX', 0.4863013698630137),
 795          ('HBB1_CYGMA', 'HBB1_XENBO', 0.45985401459854014),
 796          ('HBB1_CYGMA', 'HBB1_XENLA', 0.46715328467153283),
 797          ('HBB1_CYGMA', 'HBB1_XENTR', 0.47945205479452052),
 798          ('HBB1_CYGMA', 'MYG_ESCGI', 0.2361111111111111),
 799          ('HBB1_IGUIG', 'HBB0_PAGBO', 0.4863013698630137),
 800          ('HBB1_IGUIG', 'HBB1_ANAMI', 0.47586206896551725),
 801          ('HBB1_IGUIG', 'HBB1_CYGMA', 0.5),
 802          ('HBB1_IGUIG', 'HBB1_IGUIG', 1.0),
 803          ('HBB1_IGUIG', 'HBB1_MOUSE', 0.63013698630136983),
 804          ('HBB1_IGUIG', 'HBB1_ONCMY', 0.51034482758620692),
 805          ('HBB1_IGUIG', 'HBB1_PAGBO', 0.4863013698630137),
 806          ('HBB1_IGUIG', 'HBB1_RAT', 0.61643835616438358),
 807          ('HBB1_IGUIG', 'HBB1_SPHPU', 0.71232876712328763),
 808          ('HBB1_IGUIG', 'HBB1_TAPTE', 0.64383561643835618),
 809          ('HBB1_IGUIG', 'HBB1_TORMA', 0.36301369863013699),
 810          ('HBB1_IGUIG', 'HBB1_TRICR', 0.47916666666666669),
 811          ('HBB1_IGUIG', 'HBB1_UROHA', 0.64383561643835618),
 812          ('HBB1_IGUIG', 'HBB1_VAREX', 0.77397260273972601),
 813          ('HBB1_IGUIG', 'HBB1_XENBO', 0.4825174825174825),
 814          ('HBB1_IGUIG', 'HBB1_XENLA', 0.48951048951048953),
 815          ('HBB1_IGUIG', 'HBB1_XENTR', 0.49315068493150682),
 816          ('HBB1_IGUIG', 'MYG_DIDMA', 0.25179856115107913),
 817          ('HBB1_IGUIG', 'MYG_ERIEU', 0.28368794326241137),
 818          ('HBB1_IGUIG', 'MYG_ESCGI', 0.27659574468085107),
 819          ('HBB1_IGUIG', 'MYG_GALCR', 0.28368794326241137),
 820          ('HBB1_MOUSE', 'HBB0_PAGBO', 0.45205479452054792),
 821          ('HBB1_MOUSE', 'HBB1_ANAMI', 0.45517241379310347),
 822          ('HBB1_MOUSE', 'HBB1_CYGMA', 0.47945205479452052),
 823          ('HBB1_MOUSE', 'HBB1_IGUIG', 0.63013698630136983),
 824          ('HBB1_MOUSE', 'HBB1_MOUSE', 1.0),
 825          ('HBB1_MOUSE', 'HBB1_ONCMY', 0.50344827586206897),
 826          ('HBB1_MOUSE', 'HBB1_PAGBO', 0.4726027397260274),
 827          ('HBB1_MOUSE', 'HBB1_RAT', 0.9178082191780822),
 828          ('HBB1_MOUSE', 'HBB1_SPHPU', 0.65753424657534243),
 829          ('HBB1_MOUSE', 'HBB1_TAPTE', 0.76027397260273977),
 830          ('HBB1_MOUSE', 'HBB1_TORMA', 0.35616438356164382),
 831          ('HBB1_MOUSE', 'HBB1_TRICR', 0.52083333333333337),
 832          ('HBB1_MOUSE', 'HBB1_UROHA', 0.47945205479452052),
 833          ('HBB1_MOUSE', 'HBB1_VAREX', 0.6095890410958904),
 834          ('HBB1_MOUSE', 'HBB1_XENBO', 0.44444444444444442),
 835          ('HBB1_MOUSE', 'HBB1_XENLA', 0.44444444444444442),
 836          ('HBB1_MOUSE', 'HBB1_XENTR', 0.4589041095890411),
 837          ('HBB1_MOUSE', 'MYG_DIDMA', 0.29655172413793102),
 838          ('HBB1_MOUSE', 'MYG_ELEMA', 0.27586206896551724),
 839          ('HBB1_MOUSE', 'MYG_ERIEU', 0.30344827586206896),
 840          ('HBB1_MOUSE', 'MYG_ESCGI', 0.28965517241379313),
 841          ('HBB1_MOUSE', 'MYG_GALCR', 0.28275862068965518),
 842          ('HBB1_ONCMY', 'HBB0_PAGBO', 0.55172413793103448),
 843          ('HBB1_ONCMY', 'HBB1_ANAMI', 0.59310344827586203),
 844          ('HBB1_ONCMY', 'HBB1_CYGMA', 0.53103448275862064),
 845          ('HBB1_ONCMY', 'HBB1_IGUIG', 0.51034482758620692),
 846          ('HBB1_ONCMY', 'HBB1_MOUSE', 0.50344827586206897),
 847          ('HBB1_ONCMY', 'HBB1_ONCMY', 1.0),
 848          ('HBB1_ONCMY', 'HBB1_PAGBO', 0.56551724137931036),
 849          ('HBB1_ONCMY', 'HBB1_RAT', 0.50344827586206897),
 850          ('HBB1_ONCMY', 'HBB1_SPHPU', 0.46206896551724136),
 851          ('HBB1_ONCMY', 'HBB1_TAPTE', 0.50344827586206897),
 852          ('HBB1_ONCMY', 'HBB1_TORMA', 0.33793103448275863),
 853          ('HBB1_ONCMY', 'HBB1_TRICR', 0.41258741258741261),
 854          ('HBB1_ONCMY', 'HBB1_UROHA', 0.44827586206896552),
 855          ('HBB1_ONCMY', 'HBB1_VAREX', 0.48965517241379308),
 856          ('HBB1_ONCMY', 'HBB1_XENBO', 0.40140845070422537),
 857          ('HBB1_ONCMY', 'HBB1_XENLA', 0.39436619718309857),
 858          ('HBB1_ONCMY', 'HBB1_XENTR', 0.39310344827586208),
 859          ('HBB1_ONCMY', 'MYG_DIDMA', 0.25694444444444442),
 860          ('HBB1_ONCMY', 'MYG_ERIEU', 0.2361111111111111),
 861          ('HBB1_ONCMY', 'MYG_ESCGI', 0.25),
 862          ('HBB1_ONCMY', 'MYG_GALCR', 0.24305555555555555),
 863          ('HBB1_PAGBO', 'HBB0_PAGBO', 0.69178082191780821),
 864          ('HBB1_PAGBO', 'HBB1_ANAMI', 0.75862068965517238),
 865          ('HBB1_PAGBO', 'HBB1_CYGMA', 0.86986301369863017),
 866          ('HBB1_PAGBO', 'HBB1_IGUIG', 0.4863013698630137),
 867          ('HBB1_PAGBO', 'HBB1_MOUSE', 0.4726027397260274),
 868          ('HBB1_PAGBO', 'HBB1_ONCMY', 0.56551724137931036),
 869          ('HBB1_PAGBO', 'HBB1_PAGBO', 1.0),
 870          ('HBB1_PAGBO', 'HBB1_RAT', 0.4863013698630137),
 871          ('HBB1_PAGBO', 'HBB1_SPHPU', 0.4726027397260274),
 872          ('HBB1_PAGBO', 'HBB1_TAPTE', 0.46575342465753422),
 873          ('HBB1_PAGBO', 'HBB1_TORMA', 0.34931506849315069),
 874          ('HBB1_PAGBO', 'HBB1_TRICR', 0.4375),
 875          ('HBB1_PAGBO', 'HBB1_UROHA', 0.35616438356164382),
 876          ('HBB1_PAGBO', 'HBB1_VAREX', 0.4726027397260274),
 877          ('HBB1_PAGBO', 'HBB1_XENBO', 0.43356643356643354),
 878          ('HBB1_PAGBO', 'HBB1_XENLA', 0.44055944055944057),
 879          ('HBB1_PAGBO', 'HBB1_XENTR', 0.47945205479452052),
 880          ('HBB1_RAT', 'HBB0_PAGBO', 0.4589041095890411),
 881          ('HBB1_RAT', 'HBB1_ANAMI', 0.48965517241379308),
 882          ('HBB1_RAT', 'HBB1_CYGMA', 0.50684931506849318),
 883          ('HBB1_RAT', 'HBB1_IGUIG', 0.61643835616438358),
 884          ('HBB1_RAT', 'HBB1_MOUSE', 0.9178082191780822),
 885          ('HBB1_RAT', 'HBB1_ONCMY', 0.50344827586206897),
 886          ('HBB1_RAT', 'HBB1_PAGBO', 0.4863013698630137),
 887          ('HBB1_RAT', 'HBB1_RAT', 1.0),
 888          ('HBB1_RAT', 'HBB1_SPHPU', 0.66438356164383561),
 889          ('HBB1_RAT', 'HBB1_TAPTE', 0.76712328767123283),
 890          ('HBB1_RAT', 'HBB1_TORMA', 0.35616438356164382),
 891          ('HBB1_RAT', 'HBB1_TRICR', 0.52777777777777779),
 892          ('HBB1_RAT', 'HBB1_UROHA', 0.5),
 893          ('HBB1_RAT', 'HBB1_VAREX', 0.62328767123287676),
 894          ('HBB1_RAT', 'HBB1_XENBO', 0.45833333333333331),
 895          ('HBB1_RAT', 'HBB1_XENLA', 0.45833333333333331),
 896          ('HBB1_RAT', 'HBB1_XENTR', 0.45205479452054792),
 897          ('HBB1_RAT', 'MYG_DIDMA', 0.29655172413793102),
 898          ('HBB1_RAT', 'MYG_ELEMA', 0.28275862068965518),
 899          ('HBB1_RAT', 'MYG_ERIEU', 0.29655172413793102),
 900          ('HBB1_RAT', 'MYG_ESCGI', 0.28275862068965518),
 901          ('HBB1_RAT', 'MYG_GALCR', 0.27586206896551724),
 902          ('HBB1_SPHPU', 'HBB0_PAGBO', 0.4589041095890411),
 903          ('HBB1_SPHPU', 'HBB1_ANAMI', 0.46206896551724136),
 904          ('HBB1_SPHPU', 'HBB1_CYGMA', 0.47945205479452052),
 905          ('HBB1_SPHPU', 'HBB1_IGUIG', 0.71232876712328763),
 906          ('HBB1_SPHPU', 'HBB1_MOUSE', 0.65753424657534243),
 907          ('HBB1_SPHPU', 'HBB1_ONCMY', 0.46206896551724136),
 908          ('HBB1_SPHPU', 'HBB1_PAGBO', 0.4726027397260274),
 909          ('HBB1_SPHPU', 'HBB1_RAT', 0.66438356164383561),
 910          ('HBB1_SPHPU', 'HBB1_SPHPU', 1.0),
 911          ('HBB1_SPHPU', 'HBB1_TAPTE', 0.63698630136986301),
 912          ('HBB1_SPHPU', 'HBB1_TORMA', 0.38356164383561642),
 913          ('HBB1_SPHPU', 'HBB1_TRICR', 0.47916666666666669),
 914          ('HBB1_SPHPU', 'HBB1_UROHA', 0.54109589041095896),
 915          ('HBB1_SPHPU', 'HBB1_VAREX', 0.69178082191780821),
 916          ('HBB1_SPHPU', 'HBB1_XENBO', 0.48951048951048953),
 917          ('HBB1_SPHPU', 'HBB1_XENLA', 0.4825174825174825),
 918          ('HBB1_SPHPU', 'HBB1_XENTR', 0.4726027397260274),
 919          ('HBB1_TAPTE', 'HBB0_PAGBO', 0.4863013698630137),
 920          ('HBB1_TAPTE', 'HBB1_ANAMI', 0.48965517241379308),
 921          ('HBB1_TAPTE', 'HBB1_CYGMA', 0.4726027397260274),
 922          ('HBB1_TAPTE', 'HBB1_IGUIG', 0.64383561643835618),
 923          ('HBB1_TAPTE', 'HBB1_MOUSE', 0.76027397260273977),
 924          ('HBB1_TAPTE', 'HBB1_ONCMY', 0.50344827586206897),
 925          ('HBB1_TAPTE', 'HBB1_PAGBO', 0.46575342465753422),
 926          ('HBB1_TAPTE', 'HBB1_RAT', 0.76712328767123283),
 927          ('HBB1_TAPTE', 'HBB1_SPHPU', 0.63698630136986301),
 928          ('HBB1_TAPTE', 'HBB1_TAPTE', 1.0),
 929          ('HBB1_TAPTE', 'HBB1_TORMA', 0.34931506849315069),
 930          ('HBB1_TAPTE', 'HBB1_TRICR', 0.4861111111111111),
 931          ('HBB1_TAPTE', 'HBB1_UROHA', 0.51369863013698636),
 932          ('HBB1_TAPTE', 'HBB1_VAREX', 0.62328767123287676),
 933          ('HBB1_TAPTE', 'HBB1_XENBO', 0.4861111111111111),
 934          ('HBB1_TAPTE', 'HBB1_XENLA', 0.47222222222222221),
 935          ('HBB1_TAPTE', 'HBB1_XENTR', 0.45205479452054792),
 936          ('HBB1_TAPTE', 'MYG_DIDMA', 0.26277372262773724),
 937          ('HBB1_TAPTE', 'MYG_ERIEU', 0.27007299270072993),
 938          ('HBB1_TAPTE', 'MYG_ESCGI', 0.30344827586206896),
 939          ('HBB1_TAPTE', 'MYG_GALCR', 0.27007299270072993),
 940          ('HBB1_TORMA', 'HBB0_PAGBO', 0.31506849315068491),
 941          ('HBB1_TORMA', 'HBB1_ANAMI', 0.32413793103448274),
 942          ('HBB1_TORMA', 'HBB1_CYGMA', 0.33561643835616439),
 943          ('HBB1_TORMA', 'HBB1_IGUIG', 0.36301369863013699),
 944          ('HBB1_TORMA', 'HBB1_MOUSE', 0.35616438356164382),
 945          ('HBB1_TORMA', 'HBB1_ONCMY', 0.33793103448275863),
 946          ('HBB1_TORMA', 'HBB1_PAGBO', 0.34931506849315069),
 947          ('HBB1_TORMA', 'HBB1_RAT', 0.35616438356164382),
 948          ('HBB1_TORMA', 'HBB1_SPHPU', 0.38356164383561642),
 949          ('HBB1_TORMA', 'HBB1_TAPTE', 0.34931506849315069),
 950          ('HBB1_TORMA', 'HBB1_TORMA', 1.0),
 951          ('HBB1_TORMA', 'HBB1_TRICR', 0.31724137931034485),
 952          ('HBB1_TORMA', 'HBB1_UROHA', 0.29452054794520549),
 953          ('HBB1_TORMA', 'HBB1_VAREX', 0.35616438356164382),
 954          ('HBB1_TORMA', 'HBB1_XENBO', 0.34482758620689657),
 955          ('HBB1_TORMA', 'HBB1_XENLA', 0.33793103448275863),
 956          ('HBB1_TORMA', 'HBB1_XENTR', 0.33561643835616439),
 957          ('HBB1_TORMA', 'MYG_ESCGI', 0.25675675675675674),
 958          ('HBB1_TRICR', 'HBB0_PAGBO', 0.4375),
 959          ('HBB1_TRICR', 'HBB1_ANAMI', 0.41258741258741261),
 960          ('HBB1_TRICR', 'HBB1_CYGMA', 0.4375),
 961          ('HBB1_TRICR', 'HBB1_IGUIG', 0.47916666666666669),
 962          ('HBB1_TRICR', 'HBB1_MOUSE', 0.52083333333333337),
 963          ('HBB1_TRICR', 'HBB1_ONCMY', 0.41258741258741261),
 964          ('HBB1_TRICR', 'HBB1_PAGBO', 0.4375),
 965          ('HBB1_TRICR', 'HBB1_RAT', 0.52777777777777779),
 966          ('HBB1_TRICR', 'HBB1_SPHPU', 0.47916666666666669),
 967          ('HBB1_TRICR', 'HBB1_TAPTE', 0.4861111111111111),
 968          ('HBB1_TRICR', 'HBB1_TORMA', 0.31724137931034485),
 969          ('HBB1_TRICR', 'HBB1_TRICR', 1.0),
 970          ('HBB1_TRICR', 'HBB1_UROHA', 0.3611111111111111),
 971          ('HBB1_TRICR', 'HBB1_VAREX', 0.4513888888888889),
 972          ('HBB1_TRICR', 'HBB1_XENBO', 0.4861111111111111),
 973          ('HBB1_TRICR', 'HBB1_XENLA', 0.49305555555555558),
 974          ('HBB1_TRICR', 'HBB1_XENTR', 0.49305555555555558),
 975          ('HBB1_UROHA', 'HBB0_PAGBO', 0.4041095890410959),
 976          ('HBB1_UROHA', 'HBB1_ANAMI', 0.42857142857142855),
 977          ('HBB1_UROHA', 'HBB1_CYGMA', 0.36986301369863012),
 978          ('HBB1_UROHA', 'HBB1_IGUIG', 0.64383561643835618),
 979          ('HBB1_UROHA', 'HBB1_MOUSE', 0.51666666666666672),
 980          ('HBB1_UROHA', 'HBB1_ONCMY', 0.50420168067226889),
 981          ('HBB1_UROHA', 'HBB1_PAGBO', 0.38333333333333336),
 982          ('HBB1_UROHA', 'HBB1_RAT', 0.54166666666666663),
 983          ('HBB1_UROHA', 'HBB1_SPHPU', 0.54109589041095896),
 984          ('HBB1_UROHA', 'HBB1_TAPTE', 0.55833333333333335),
 985          ('HBB1_UROHA', 'HBB1_TORMA', 0.31034482758620691),
 986          ('HBB1_UROHA', 'HBB1_TRICR', 0.39316239316239315),
 987          ('HBB1_UROHA', 'HBB1_UROHA', 1.0),
 988          ('HBB1_UROHA', 'HBB1_VAREX', 0.59589041095890416),
 989          ('HBB1_UROHA', 'HBB1_XENBO', 0.42608695652173911),
 990          ('HBB1_UROHA', 'HBB1_XENLA', 0.41739130434782606),
 991          ('HBB1_UROHA', 'HBB1_XENTR', 0.40000000000000002),
 992          ('HBB1_UROHA', 'MYG_ERIEU', 0.27927927927927926),
 993          ('HBB1_VAREX', 'HBB0_PAGBO', 0.49315068493150682),
 994          ('HBB1_VAREX', 'HBB1_ANAMI', 0.48275862068965519),
 995          ('HBB1_VAREX', 'HBB1_CYGMA', 0.4863013698630137),
 996          ('HBB1_VAREX', 'HBB1_IGUIG', 0.77397260273972601),
 997          ('HBB1_VAREX', 'HBB1_MOUSE', 0.6095890410958904),
 998          ('HBB1_VAREX', 'HBB1_ONCMY', 0.48965517241379308),
 999          ('HBB1_VAREX', 'HBB1_PAGBO', 0.4726027397260274),
1000          ('HBB1_VAREX', 'HBB1_RAT', 0.62328767123287676),
1001          ('HBB1_VAREX', 'HBB1_SPHPU', 0.69178082191780821),
1002          ('HBB1_VAREX', 'HBB1_TAPTE', 0.62328767123287676),
1003          ('HBB1_VAREX', 'HBB1_TORMA', 0.35616438356164382),
1004          ('HBB1_VAREX', 'HBB1_TRICR', 0.4513888888888889),
1005          ('HBB1_VAREX', 'HBB1_UROHA', 0.59589041095890416),
1006          ('HBB1_VAREX', 'HBB1_VAREX', 1.0),
1007          ('HBB1_VAREX', 'HBB1_XENBO', 0.51048951048951052),
1008          ('HBB1_VAREX', 'HBB1_XENLA', 0.5174825174825175),
1009          ('HBB1_VAREX', 'HBB1_XENTR', 0.4726027397260274),
1010          ('HBB1_VAREX', 'MYG_DIDMA', 0.25531914893617019),
1011          ('HBB1_VAREX', 'MYG_ERIEU', 0.25531914893617019),
1012          ('HBB1_VAREX', 'MYG_ESCGI', 0.24822695035460993),
1013          ('HBB1_VAREX', 'MYG_GALCR', 0.24822695035460993),
1014          ('HBB1_XENBO', 'HBB0_PAGBO', 0.43356643356643354),
1015          ('HBB1_XENBO', 'HBB1_ANAMI', 0.4460431654676259),
1016          ('HBB1_XENBO', 'HBB1_CYGMA', 0.45985401459854014),
1017          ('HBB1_XENBO', 'HBB1_IGUIG', 0.4825174825174825),
1018          ('HBB1_XENBO', 'HBB1_MOUSE', 0.44444444444444442),
1019          ('HBB1_XENBO', 'HBB1_ONCMY', 0.40140845070422537),
1020          ('HBB1_XENBO', 'HBB1_PAGBO', 0.43356643356643354),
1021          ('HBB1_XENBO', 'HBB1_RAT', 0.45833333333333331),
1022          ('HBB1_XENBO', 'HBB1_SPHPU', 0.48951048951048953),
1023          ('HBB1_XENBO', 'HBB1_TAPTE', 0.4861111111111111),
1024          ('HBB1_XENBO', 'HBB1_TORMA', 0.34482758620689657),
1025          ('HBB1_XENBO', 'HBB1_TRICR', 0.4861111111111111),
1026          ('HBB1_XENBO', 'HBB1_UROHA', 0.38461538461538464),
1027          ('HBB1_XENBO', 'HBB1_VAREX', 0.51048951048951052),
1028          ('HBB1_XENBO', 'HBB1_XENBO', 1.0),
1029          ('HBB1_XENBO', 'HBB1_XENLA', 0.96551724137931039),
1030          ('HBB1_XENBO', 'HBB1_XENTR', 0.76388888888888884),
1031          ('HBB1_XENBO', 'MYG_DIDMA', 0.32075471698113206),
1032          ('HBB1_XENBO', 'MYG_ELEMA', 0.27358490566037735),
1033          ('HBB1_XENBO', 'MYG_ERIEU', 0.32075471698113206),
1034          ('HBB1_XENBO', 'MYG_GALCR', 0.32075471698113206),
1035          ('HBB1_XENLA', 'HBB0_PAGBO', 0.44055944055944057),
1036          ('HBB1_XENLA', 'HBB1_ANAMI', 0.45323741007194246),
1037          ('HBB1_XENLA', 'HBB1_CYGMA', 0.46715328467153283),
1038          ('HBB1_XENLA', 'HBB1_IGUIG', 0.48951048951048953),
1039          ('HBB1_XENLA', 'HBB1_MOUSE', 0.44444444444444442),
1040          ('HBB1_XENLA', 'HBB1_ONCMY', 0.39436619718309857),
1041          ('HBB1_XENLA', 'HBB1_PAGBO', 0.44055944055944057),
1042          ('HBB1_XENLA', 'HBB1_RAT', 0.45833333333333331),
1043          ('HBB1_XENLA', 'HBB1_SPHPU', 0.4825174825174825),
1044          ('HBB1_XENLA', 'HBB1_TAPTE', 0.47222222222222221),
1045          ('HBB1_XENLA', 'HBB1_TORMA', 0.33793103448275863),
1046          ('HBB1_XENLA', 'HBB1_TRICR', 0.49305555555555558),
1047          ('HBB1_XENLA', 'HBB1_UROHA', 0.3776223776223776),
1048          ('HBB1_XENLA', 'HBB1_VAREX', 0.5174825174825175),
1049          ('HBB1_XENLA', 'HBB1_XENBO', 0.96551724137931039),
1050          ('HBB1_XENLA', 'HBB1_XENLA', 1.0),
1051          ('HBB1_XENLA', 'HBB1_XENTR', 0.75),
1052          ('HBB1_XENLA', 'MYG_DIDMA', 0.32075471698113206),
1053          ('HBB1_XENLA', 'MYG_ELEMA', 0.26415094339622641),
1054          ('HBB1_XENLA', 'MYG_ERIEU', 0.32075471698113206),
1055          ('HBB1_XENLA', 'MYG_ESCGI', 0.28282828282828282),
1056          ('HBB1_XENLA', 'MYG_GALCR', 0.32075471698113206),
1057          ('HBB1_XENTR', 'HBB0_PAGBO', 0.4726027397260274),
1058          ('HBB1_XENTR', 'HBB1_ANAMI', 0.4689655172413793),
1059          ('HBB1_XENTR', 'HBB1_CYGMA', 0.47945205479452052),
1060          ('HBB1_XENTR', 'HBB1_IGUIG', 0.49315068493150682),
1061          ('HBB1_XENTR', 'HBB1_MOUSE', 0.4589041095890411),
1062          ('HBB1_XENTR', 'HBB1_ONCMY', 0.39310344827586208),
1063          ('HBB1_XENTR', 'HBB1_PAGBO', 0.47945205479452052),
1064          ('HBB1_XENTR', 'HBB1_RAT', 0.45205479452054792),
1065          ('HBB1_XENTR', 'HBB1_SPHPU', 0.4726027397260274),
1066          ('HBB1_XENTR', 'HBB1_TAPTE', 0.45205479452054792),
1067          ('HBB1_XENTR', 'HBB1_TORMA', 0.33561643835616439),
1068          ('HBB1_XENTR', 'HBB1_TRICR', 0.49305555555555558),
1069          ('HBB1_XENTR', 'HBB1_UROHA', 0.35616438356164382),
1070          ('HBB1_XENTR', 'HBB1_VAREX', 0.4726027397260274),
1071          ('HBB1_XENTR', 'HBB1_XENBO', 0.76388888888888884),
1072          ('HBB1_XENTR', 'HBB1_XENLA', 0.75),
1073          ('HBB1_XENTR', 'HBB1_XENTR', 1.0),
1074          ('HBB1_XENTR', 'MYG_DIDMA', 0.2857142857142857),
1075          ('HBB1_XENTR', 'MYG_ERIEU', 0.27067669172932329),
1076          ('HBB1_XENTR', 'MYG_ESCGI', 0.27272727272727271),
1077          ('HBB1_XENTR', 'MYG_GALCR', 0.2781954887218045),
1078          ('MYG_DIDMA', 'HBB0_PAGBO', 0.22222222222222221),
1079          ('MYG_DIDMA', 'HBB1_IGUIG', 0.25179856115107913),
1080          ('MYG_DIDMA', 'HBB1_MOUSE', 0.29655172413793102),
1081          ('MYG_DIDMA', 'HBB1_ONCMY', 0.25694444444444442),
1082          ('MYG_DIDMA', 'HBB1_RAT', 0.29655172413793102),
1083          ('MYG_DIDMA', 'HBB1_TAPTE', 0.26277372262773724),
1084          ('MYG_DIDMA', 'HBB1_VAREX', 0.25531914893617019),
1085          ('MYG_DIDMA', 'HBB1_XENBO', 0.32075471698113206),
1086          ('MYG_DIDMA', 'HBB1_XENLA', 0.32075471698113206),
1087          ('MYG_DIDMA', 'HBB1_XENTR', 0.2857142857142857),
1088          ('MYG_DIDMA', 'MYG_DIDMA', 1.0),
1089          ('MYG_DIDMA', 'MYG_ELEMA', 0.81045751633986929),
1090          ('MYG_DIDMA', 'MYG_ERIEU', 0.87581699346405228),
1091          ('MYG_DIDMA', 'MYG_ESCGI', 0.83552631578947367),
1092          ('MYG_DIDMA', 'MYG_GALCR', 0.83006535947712423),
1093          ('MYG_ELEMA', 'HBB0_PAGBO', 0.20833333333333334),
1094          ('MYG_ELEMA', 'HBB1_MOUSE', 0.27586206896551724),
1095          ('MYG_ELEMA', 'HBB1_RAT', 0.28275862068965518),
1096          ('MYG_ELEMA', 'HBB1_XENBO', 0.27358490566037735),
1097          ('MYG_ELEMA', 'HBB1_XENLA', 0.26415094339622641),
1098          ('MYG_ELEMA', 'MYG_DIDMA', 0.81045751633986929),
1099          ('MYG_ELEMA', 'MYG_ELEMA', 1.0),
1100          ('MYG_ELEMA', 'MYG_ERIEU', 0.82352941176470584),
1101          ('MYG_ELEMA', 'MYG_ESCGI', 0.83552631578947367),
1102          ('MYG_ELEMA', 'MYG_GALCR', 0.84313725490196079),
1103          ('MYG_ERIEU', 'HBB0_PAGBO', 0.21527777777777779),
1104          ('MYG_ERIEU', 'HBB1_IGUIG', 0.28368794326241137),
1105          ('MYG_ERIEU', 'HBB1_MOUSE', 0.30344827586206896),
1106          ('MYG_ERIEU', 'HBB1_ONCMY', 0.2361111111111111),
1107          ('MYG_ERIEU', 'HBB1_RAT', 0.29655172413793102),
1108          ('MYG_ERIEU', 'HBB1_TAPTE', 0.27007299270072993),
1109          ('MYG_ERIEU', 'HBB1_UROHA', 0.27927927927927926),
1110          ('MYG_ERIEU', 'HBB1_VAREX', 0.25531914893617019),
1111          ('MYG_ERIEU', 'HBB1_XENBO', 0.32075471698113206),
1112          ('MYG_ERIEU', 'HBB1_XENLA', 0.32075471698113206),
1113          ('MYG_ERIEU', 'HBB1_XENTR', 0.27067669172932329),
1114          ('MYG_ERIEU', 'MYG_DIDMA', 0.87581699346405228),
1115          ('MYG_ERIEU', 'MYG_ELEMA', 0.82352941176470584),
1116          ('MYG_ERIEU', 'MYG_ERIEU', 1.0),
1117          ('MYG_ERIEU', 'MYG_ESCGI', 0.83552631578947367),
1118          ('MYG_ERIEU', 'MYG_GALCR', 0.85620915032679734),
1119          ('MYG_ESCGI', 'HBB0_PAGBO', 0.25),
1120          ('MYG_ESCGI', 'HBB1_CYGMA', 0.2361111111111111),
1121          ('MYG_ESCGI', 'HBB1_IGUIG', 0.27659574468085107),
1122          ('MYG_ESCGI', 'HBB1_MOUSE', 0.28965517241379313),
1123          ('MYG_ESCGI', 'HBB1_ONCMY', 0.25),
1124          ('MYG_ESCGI', 'HBB1_RAT', 0.28275862068965518),
1125          ('MYG_ESCGI', 'HBB1_TAPTE', 0.3611111111111111),
1126          ('MYG_ESCGI', 'HBB1_TORMA', 0.25675675675675674),
1127          ('MYG_ESCGI', 'HBB1_VAREX', 0.24822695035460993),
1128          ('MYG_ESCGI', 'HBB1_XENLA', 0.28282828282828282),
1129          ('MYG_ESCGI', 'HBB1_XENTR', 0.27272727272727271),
1130          ('MYG_ESCGI', 'MYG_DIDMA', 0.83552631578947367),
1131          ('MYG_ESCGI', 'MYG_ELEMA', 0.83552631578947367),
1132          ('MYG_ESCGI', 'MYG_ERIEU', 0.83552631578947367),
1133          ('MYG_ESCGI', 'MYG_ESCGI', 1.0),
1134          ('MYG_ESCGI', 'MYG_GALCR', 0.84210526315789469),
1135          ('MYG_GALCR', 'HBB0_PAGBO', 0.24305555555555555),
1136          ('MYG_GALCR', 'HBB1_IGUIG', 0.28368794326241137),
1137          ('MYG_GALCR', 'HBB1_MOUSE', 0.28275862068965518),
1138          ('MYG_GALCR', 'HBB1_ONCMY', 0.24305555555555555),
1139          ('MYG_GALCR', 'HBB1_RAT', 0.27586206896551724),
1140          ('MYG_GALCR', 'HBB1_TAPTE', 0.27007299270072993),
1141          ('MYG_GALCR', 'HBB1_VAREX', 0.24822695035460993),
1142          ('MYG_GALCR', 'HBB1_XENBO', 0.32075471698113206),
1143          ('MYG_GALCR', 'HBB1_XENLA', 0.32075471698113206),
1144          ('MYG_GALCR', 'HBB1_XENTR', 0.2781954887218045),
1145          ('MYG_GALCR', 'MYG_DIDMA', 0.83006535947712423),
1146          ('MYG_GALCR', 'MYG_ELEMA', 0.84313725490196079),
1147          ('MYG_GALCR', 'MYG_ERIEU', 0.85620915032679734),
1148          ('MYG_GALCR', 'MYG_ESCGI', 0.84210526315789469),
1149          ('MYG_GALCR', 'MYG_GALCR', 1.0),
1150          ('PRCA_ANASP', 'PRCA_ANASP', 1.0),
1151          ('PRCA_ANASP', 'PRCA_ANAVA', 0.97222222222222221),
1152          ('PRCA_ANAVA', 'PRCA_ANASP', 0.97222222222222221),
1153          ('PRCA_ANAVA', 'PRCA_ANAVA', 1.0)]
1154
1155 ###
1156
1157 if __name__ == '__main__':
1158     PygrTestProgram(verbosity=2)
1159