Naive implementation of a de-Bruijn graph

Думала я как собрать риды в контиг (для задачки). Помню, на лекциях все так восторженно рассказывали про граф Де Брюйна. Но почему так делать - это оптимально не объяснили. Ну и ладно, нашла я код под названием "naive implementation of a de-Bruijn graph". Ключевой момент там в ф-ии  get_contig_forward, она очень упрощена, не рассматривается случай раздвоения вариантов (когда надо выбрать более длинный контиг), не рассматривается случай повтороу участка кода длиннее k (да и как их обработать). Но сама обертка есть...

import collections

 

def twin(km):

    return Seq.reverse_complement(km)

 

def kmers(seq,k):

    for i in xrange(len(seq)-k+1):

        yield seq[i:i+k]

 

def fw(km):

    for x in 'ACGT':

        yield km[1:]+x

 

def bw(km):

    for x in 'ACGT':

        yield x + km[:-1]

 

def build(fn,k=31,limit=1):

   

   d = collections.defaultdict(int)

   with open('tests/1', 'r') as f:
        main_string = f.readline()
        for line in f.readlines()[2:]:
            line = line.strip()

                 for km in kmers(seq,k):

                    d[km] +=1

                 seq = twin(seq)

                 for km in kmers(seq,k):

                    d[km] += 1

 

    d1 = [x for x in d if d[x] <= limit]

    for x in d1:

        del d[x]

 

    return d

 

def get_contig_forward(d,km):

    c_fw = [km]

     

    while True:

        if sum(x in d for x in fw(c_fw[-1])) != 1:

            break

         

        cand = [x for x in fw(c_fw[-1]) if x in d][0]

        if cand in c_fw:

            break

         

        if sum(x in d for x in bw(cand)) != 1:

            break

 

        c_fw.append(cand)

 

    return c_fw

 

def get_contig(d,km):

    c_fw = get_contig_forward(d,km)

    c_bw = get_contig_forward(d,twin(km))

 

    c = [twin(x) for x in c_bw[-1:0:-1]] + c_fw

    s = c[0] + ''.join(x[-1] for x in c[1:])

    return s,c

 

def all_contigs(d):

    done = set()

    r = []

    for x in d:

        if x not in done:

            s,c = get_contig(d,x)

            for y in c:

                done.add(y)

                done.add(twin(y))

            r.append(s)

    return r