Python seqdb.SequenceFileDB类代码示例

 def test_basic_construction(self):
     db = SequenceFileDB(self.dbfile)
     try:
         assert str(db.get('seq1')).startswith('atggtgtca')
         assert str(db.get('seq2')).startswith('GTGTTGAA')
     finally:
         db.close()

    def test_inverse_add_behavior(self):
        dnaseq = testutil.datafile('dnaseq.fasta')
        seqdb = SequenceFileDB(dnaseq)
        try:
            seq = seqdb['seq1']

            name = (~self.db)[seq]
        finally:
            seqdb.close() # only need to close if exception occurs

 def test_funny_key2(self):
     "check handling of ID containing multiple separators"
     dnaseq = testutil.datafile('funnyseq.fasta')
     seqdb = SequenceFileDB(dnaseq)     # contains 'seq1', 'seq2'
     try:
         pudb = PrefixUnionDict({'prefix': seqdb})
         seq = pudb['prefix.seq.2.even.longer']
     finally:
         seqdb.close()

    def test_cache(self):
        "Sequence slice cache mechanics."

        dnaseq = testutil.datafile('dnaseq.fasta')
        db = SequenceFileDB(dnaseq)

        try:
            # create cache components
            cacheDict = {}
            cacheHint = db.cacheHint

            # get seq1
            seq1 = db['seq1']

            # _cache is only created on first cache attempt
            assert not hasattr(db, '_cache')

            # build an 'owner' object
            class AnonymousOwner(object):
                pass
            owner = AnonymousOwner()

            # save seq1 in cache
            cacheDict['seq1'] = (seq1.start, seq1.stop)
            cacheHint(cacheDict, owner)
            del cacheDict                   # 'owner' now holds reference

            # peek into _cache and assert that only the ival coordinates
            # are stored
            v = db._cache.values()[0]
            assert len(v['seq1']) == 2
            del v

            # force a cache access & check that now we've stored actual string
            ival = str(seq1[5:10])
            v = db._cache.values()[0]
            # ...check that we've stored actual string
            assert len(v['seq1']) == 3

            # again force cache access, this time to the stored sequence string
            ival = str(seq1[5:10])

            # now, eliminate all references to the cache proxy dict
            del owner

            # trash unused objects - not strictly necessary, because there are
            # no islands of circular references & so all objects are already
            # deallocated, but that's implementation dependent.
            gc.collect()

            # ok, cached values should now be gone.
            v = db._cache.values()
            assert len(v) == 0
        finally:
            db.close()

 def test_headerfile_create_conflict(self):
     "test non-empty prefixDict with a passed in PUD header file: conflict"
     subdb = SequenceFileDB(self.dbfile)
     try:
         header = testutil.datafile('prefixUnionDict-1.txt')
         try:
             db = PrefixUnionDict(filename=header, prefixDict={ 'foo' : subdb })
             assert 0, "should not get here"
         except TypeError:
             pass
     finally:
         subdb.close()

    def test_nlmsaslice_cache(self):
        "NLMSASlice sequence caching & removal"

        # set up sequences
        dnaseq = testutil.datafile('dnaseq.fasta')

        db = SequenceFileDB(dnaseq, autoGC=-1) # use pure WeakValueDict...
        try:
            gc.collect()
            assert len(db._weakValueDict)==0, '_weakValueDict should be empty'
            seq1, seq2 = db['seq1'], db['seq2']
            assert len(db._weakValueDict)==2, \
                    '_weakValueDict should have 2 seqs'

            # build referencing NLMSA
            mymap = NLMSA('test', 'memory', db, pairwiseMode=True)
            mymap += seq1
            mymap[seq1] += seq2
            mymap.build()

            # check: no cache
            assert not hasattr(db, '_cache'), 'should be no cache yet'

            seq1, seq2 = db['seq1'], db['seq2'] # re-retrieve
            # now retrieve a NLMSASlice, forcing entry of seq into cache
            ival = seq1[5:10]
            x = mymap[ival]

            assert len(db._cache.values()) != 0

            n1 = len(db._cache)
            assert n1 == 1, "should be exactly one cache entry, not %d" % \
                    (n1, )

            # ok, now trash referencing arguments & make sure of cleanup
            del x
            gc.collect()

            assert len(db._cache.values()) == 0


            n2 = len(db._cache)
            assert n2 == 0, '%d objects remain; cache memory leak!' % n2
            # FAIL because of __dealloc__ error in cnestedlist.NLMSASlice.

            # Drop our references, the cache should empty.
            del mymap, ival, seq1, seq2
            gc.collect()
            # check that db._weakValueDict cache is empty
            assert len(db._weakValueDict)==0, '_weakValueDict should be empty'
        finally:
            db.close()

    def test_iadd_db_twice(self):
        dnaseq = testutil.datafile('dnaseq.fasta')
        seqdb = SequenceFileDB(dnaseq)
        try:
            new_seq = seqdb['seq1']

            self.db += new_seq
            name1 = (~self.db)[new_seq]

            self.db += new_seq              # should do nothing...
            name2 = (~self.db)[new_seq]
            assert name1 == name2           # ...leaving seq with same name.
        finally:
            seqdb.close()

    def test_no_db_info(self):
        dnaseq = testutil.datafile('dnaseq.fasta')
        seqdb = SequenceFileDB(dnaseq)
        try:
            new_seq = seqdb['seq1']

            assert getattr(seqdb, '_persistent_id', None) is None
            del seqdb.filepath

            self.db += new_seq
            name = (~self.db)[new_seq]
            assert name == 'noname0.seq1'
        finally:
            seqdb.close()

    def test_build_seqLenDict_with_reader(self):
        "Test that building things works properly when specifying a reader."

        class InfoBag(object):

            def __init__(self, **kw):
                self.__dict__.update(kw)

        # first, load the db & save the sequence info in a list
        l = []
        db = SequenceFileDB(self.dbfile)
        try:
            for k, v in db.items():
                info = InfoBag(id=k, length=len(v), sequence=str(v))
                l.append(info)
        finally:
            # now, erase the existing files, and recreate the db.
            db.close()
        self.trash_intermediate_files()

        # create a fake reader with access to the saved info
        def my_fake_reader(fp, filename, info_list=l):
            return info_list

        # now try creating with the fake reader
        db = SequenceFileDB(self.dbfile, reader=my_fake_reader)

        # did it work?
        try:
            assert str(db.get('seq1')).startswith('atggtgtca')
            assert str(db.get('seq2')).startswith('GTGTTGAA')
        finally:
            db.close()

    def test_inverse_noadd_behavior(self):
        # compare with test_inverse_add_behavior...
        db = SeqPrefixUnionDict(addAll=False)
        dnaseq = testutil.datafile('dnaseq.fasta')
        seqdb = SequenceFileDB(dnaseq)
        try:
            seq = seqdb['seq1']

            try:
                name = (~db)[seq]
                assert 0, "should not get here"
            except KeyError:
                pass
        finally:
            seqdb.close()

    def test_headerfile_write_fail(self):
        subdb = SequenceFileDB(self.dbfile)
        try:
            del subdb.filepath  # remove 'filepath' attribute for test
            db = PrefixUnionDict({'prefix': subdb})

            assert len(db) == 2
            assert 'prefix.seq1' in db

            output = testutil.tempdatafile('prefixUnionDict-write-fail.txt')
            try:
                db.writeHeaderFile(output)
            except AttributeError:
                pass
        finally:
            subdb.close() # closes both db and subdb

    def test_build_seqLenDict_with_bad_reader(self):
        "Test that building things fails properly with a bad reader."

        class InfoBag(object):

            def __init__(self, **kw):
                self.__dict__.update(kw)

        # first, load the db & save the sequence info in a list
        l = []
        db = SequenceFileDB(self.dbfile)
        try:
            for k, v in db.items():
                info = InfoBag(id=k, length=0, sequence=str(v))
                l.append(info)
        finally:
            # now, erase the existing files, and recreate the db.
            db.close()
        self.trash_intermediate_files()

        # create a fake reader with access to the saved info
        def my_fake_reader(fp, filename, info_list=l):
            return info_list

        # now try creating with the fake reader
        try:
            db = SequenceFileDB(self.dbfile, reader=my_fake_reader)
            try:
                assert 0, "should not reach here; db construction should fail!"
            finally:
                db.close()
        except ValueError:
            pass                        # ValueError is expected

    def test_basic_iadd(self):
        dnaseq = testutil.datafile('dnaseq.fasta')
        seqdb = SequenceFileDB(dnaseq)
        try:
            new_seq = seqdb['seq1']

            self.db += new_seq

            assert new_seq in self.db
            name = (~self.db)[new_seq]
            assert name == 'dnaseq.seq1', name

            ###

            seqdb2 = SequenceFileDB(dnaseq)
            try:
                # Munge the filepath for testing.
                seqdb2.filepath = 'foo'
                new_seq2 = seqdb2['seq1']

                self.db += new_seq2
                name2 = (~self.db)[new_seq2]
                assert name2 == 'foo.seq1', name2
            finally:
                seqdb2.close()
        finally:
            seqdb.close()

def main(argv=None):
    parser = make_parser()
    args = parser.parse_args(argv)
    genome = SequenceFileDB(args.genome)


    pwm = [IUPAC_SCORES[l] for l in args.consensus]
    pwm.extend([REQUIRED_SCORES[l] for l in args.required_3p_seq])
    pwm = motility.PWM(pwm)

    # find all matches
    with open(args.outfile, 'w') as outfile:
        for chrom in genome.keys():
            chromseq = str(genome[chrom])
            print "searching ", chrom, "of length", len(chromseq)
            if len(chromseq) < len(pwm):
                print 'chromosome/fragment', chrom, 'is too short'
                continue
            matches = pwm.find(chromseq, -args.mismatches)
            for start, stop, strand, seq in matches:
                score = pwm.calc_score(seq)
                outfile.write('\t'.join(
                    [chrom, str(start), str(stop), seq, str(score),
                     '+' if strand == 1 else '-']) + '\n')

    def test_iadd_duplicate_seqdb(self):
        dnaseq = testutil.datafile('dnaseq.fasta')
        seqdb = SequenceFileDB(dnaseq)
        try:
            seqdb2 = SequenceFileDB(dnaseq)
            try:
                new_seq = seqdb['seq1']
                new_seq2 = seqdb2['seq1']

                self.db += new_seq
                try:
                    self.db += new_seq2
                    assert 0, "should never reach this point"
                except ValueError:
                    pass
            finally:
                seqdb2.close()
        finally:
            seqdb.close()