Неожиданный результат при рандомизированном поиске мотивов в цепочках ДНК

У меня есть следующие t=5 цепочки ДНК:

DNA = '''CGCCCCTCTCGGGGGTGTTCAGTAAACGGCCA
GGGCGAGGTATGTGTAAGTGCCAAGGTGCCAG
TAGTACCGAGACCGAAAGAAGTATACAGGCGT
TAGATCAAGTTTCAGGTGCACGTCGGTGAACC
AATCCACCAGCTCCACGTGCAATGTTGGCCTA'''
k = 8
t = 5

Я пытаюсь найти лучшие мотивы длины k=8 из коллекции строк, используя преобразование Лапласа для случайной выборки фрагментов длины k из каждой из t строк.

Мои вспомогательные функции следующие:

def window(s, k):
    for i in range(1 + len(s) - k):
        yield s[i:i+k]

def HammingDistance(seq1, seq2):
    if len(seq1) != len(seq2):
        raise ValueError('Undefined for sequences of unequal length.')
    return sum(ch1 != ch2 for ch1, ch2 in zip(seq1, seq2))

def score(motifs):
    score = 0
    for i in range(len(motifs[0])):
        motif = ''.join([motifs[j][i] for j in range(len(motifs))])
        score += min([HammingDistance(motif, homogeneous*len(motif)) for homogeneous in 'ACGT'])
    return score

def profile(motifs):
    prof = []
    for i in range(len(motifs[0])):
        col = ''.join([motifs[j][i] for j in range(len(motifs))])
        prof.append([float(col.count(nuc))/float(len(col)) for nuc in 'ACGT'])
    return prof

def profile_most_probable_kmer(dna, k, prof):
    dna = dna.splitlines()
    nuc_loc = {nucleotide:index for index,nucleotide in enumerate('ACGT')}
    motif_matrix = []
    max_prob = [-1, None]
    for i in range(len(dna)):
        motif_matrix.append(max_prob)
    for i in range(len(dna)):
        for chunk in window(dna[i],K):
            current_prob = 1
            for j, nuc in enumerate(chunk):
                current_prob*=prof[j][nuc_loc[nuc]]
            if current_prob>motif_matrix[i][0]:
                motif_matrix[i] = [current_prob,chunk]
    return list(list(zip(*motif_matrix))[1])

def profile_with_pseudocounts(motifs):
    prof = []
    for i in range(len(motifs[0])):
        col = ''.join([motifs[j][i] for j in range(len(motifs))])
        prof.append([float(col.count(nuc)+1)/float(len(col)+4) for nuc in 'ACGT'])
    return prof

from random import randint

def SampleMotifs(Dna,k,t):
    Dna = Dna.splitlines()
    BestMotifs = []
    for line in Dna:
        position = randint(0,len(line)-k)
        BestMotifs.append(line[position:position+k])
    return BestMotifs

def motifs_from_profile(profile, dna, k):
    return [profile_most_probable_kmer(seq,k,profile) for seq in dna]


def randomized_motif_search(dna,k,t):
    from random import randint
    dna = dna.splitlines()
    rand_ints = [randint(0,len(dna[0])-k)) for a in range(len(dna))]
    motifs = [dna[i][r:r+k] for i,r in enumerate(rand_ints)]
    best_score = [score(motifs), motifs]
    while True:
        current_profile = profile_with_pseudocounts(motifs)
        motifs = motifs_from_profile(current_profile, dna, k)
        current_score = score(motifs)
        if current_score < best_score[0]:
            best_score = [current_score, motifs]
        else:
            return best_score[1]


def Laplace(dna,k,t):
    i = 0
    LastMotifs = randomized_motif_search(dna,k,t)
    while i < 1000:
        try:
            BestMotifs = randomized_motif_search(dna,k,t)
            if score(BestMotifs)<score(LastMotifs):
                LastMotifs = BestMotifs
        except:
            pass
        i+=1
    print(*LastMotifs)

Результат, который я должен получить для этого:

TCTCGGGG
CCAAGGTG
TACAGGCG
TTCAGGTG
TCCACGTG

Я получаю разные результаты каждый раз, чего я ожидаю от метода со случайным элементом, но он должен сходиться, когда я повторяю 1000 раз и обновляю свои лучшие мотивы только в том случае, если оценка ниже. Тот факт, что мне пришлось добавить обработчик ошибок в лаплас, поскольку я получаю индекс при вызове randomized_motif_search(dna,k,t), говорит мне, что это может быть источником проблемы. Я провел последние два дня, просматривая код, чтобы убедиться, что все в правильной форме, но факт, что я получаю неправильный ответ или следующие ошибки:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-811-ee735449bb8e> in <module>
----> 1 Laplace(DNA,K,T)

<ipython-input-810-31301d79eb95> in Laplace(dna, k, t)
      3     LastMotifs = randomized_motif_search(dna,k,t)
      4     while i < 2000:
----> 5         BestMotifs = randomized_motif_search(dna,k,t)
      6         if score(BestMotifs)<score(LastMotifs):
      7             LastMotifs = BestMotifs

<ipython-input-809-43600d882734> in randomized_motif_search(dna, k, t)
      8     while True:
      9         current_profile = profile_with_pseudocounts(motifs)
---> 10         motifs = motifs_from_profile(current_profile, dna, k)
     11         current_score = score(motifs)
     12         if current_score < best_score[0]:

<ipython-input-408-7c866045d839> in motifs_from_profile(profile, dna, k)
      1 def motifs_from_profile(profile, dna, k):
----> 2     return [profile_most_probable_kmer(seq,k,profile) for seq in dna]

<ipython-input-408-7c866045d839> in <listcomp>(.0)
      1 def motifs_from_profile(profile, dna, k):
----> 2     return [profile_most_probable_kmer(seq,k,profile) for seq in dna]

<ipython-input-795-56f83ba5ee2b> in profile_most_probable_kmer(dna, k, prof)
     25             current_prob = 1
     26             for j, nuc in enumerate(chunk):
---> 27                 current_prob*=prof[j][nuc_loc[nuc]]
     28             if current_prob>motif_matrix[i][0]:
     29                 motif_matrix[i] = [current_prob,chunk]

IndexError: list index out of range

Это более чем досадно. Реальная помощь будет принята с благодарностью.

РЕДАКТИРОВАТЬ: проблема заключалась в моем индексировании случайных целых чисел, которые отбирали случайные мотивы из строк ДНК, и что моя функция motifs_from_profile возвращает список списков, а не просто список, от которого зависит код. Я обновил функции ниже: Хотя эти исправления разрешили проблемы, которые вызывали ошибки в коде, теперь я получаю вывод каждый раз, когда запускаю функцию Laplace, но результат не такой, как я ожидал, даже когда я кормлю правильный ответ на первой итерации. Я постараюсь изо всех сил отладить, что происходит в функции подсчета очков, и, как мне кажется, просмотреть литературу. Может быть, поможет более неопределенный вклад сообщества, но кто знает?

обновленный randomized_motif_search:

def randomized_motif_search(dna,k,t):
    from random import randint
    from itertools import chain
    dna = dna.splitlines()
#     Randomly generate k-mers from each sequence in the dna list.
    rand_ints = [randint(0,len(dna[0])-(k)) for a in range(len(dna))]
    motifs = [dna[i][r:r+k] for i,r in enumerate(rand_ints)]
    best_score = [score(motifs), motifs]
    while True:
        current_profile = profile_with_pseudocounts(motifs)
        mfp = motifs_from_profile(current_profile,dna,k)
        motifs = []
        for i in range(len(mfp)):
            motifs.append(mfp[i][0])
        current_score = score(motifs)
        if current_score < best_score[0]:
            best_score = [current_score, motifs]
        else:
            return best_score[1]

и новый Laplace:

def Laplace(dna,k,t):
    i = 0

    LastMotifs = randomized_motif_search(dna,k,t)
    while i < 1000:
        BestMotifs = randomized_motif_search(dna,k,t)
        if score(BestMotifs)<score(LastMotifs):
            LastMotifs = BestMotifs
        i+=1
    print(*LastMotifs)

РЕДАКТИРОВАТЬ: Дорогой журнал, я испортил метод оценки и выяснил, как вернуть оценку мотива, но все еще не прихожу к правильным ответам на вопрос. Я официально застрял, вот обновленный код функции score с правильной индексацией:

def score(motifs):
    score = 0
    for i in range(len(motifs[0])):
        motif = ''.join([Motifs[j][i] for j in range(len(Motifs))])
        score+=min([HammingDistance(motif,homogenous*len(motif)) for homogenous in 'ACGT'])
    return score

python random search bioinformatics dna-sequence
person slap-a-da-bias    schedule 20.02.2020    source источник
comment
Тот факт, что мне пришлось поместить обработчик ошибок в лаплас. Выявление всех исключений и их отбрасывание, безусловно, является смелым способом обработки ошибок (см .: stackoverflow.com/questions/54948548/, stackoverflow.com/questions/4990718/). В чем именно заключалась ошибка в Laplace?   -  person AMC    schedule 20.02.2020
comment
@AMC Это индекс списка вне допустимого диапазона, который исходит из вызова profile_most_probable_kmer в строке: _1 _... Он только иногда вызывает ошибку, но, может быть, это как-то связано с генерируемыми мной rand_ints?   -  person slap-a-da-bias    schedule 20.02.2020
comment
Это только иногда вызывает ошибку, так что, может быть, это как-то связано с генерируемыми мной rand_ints? Да, это, безусловно, возможно.   -  person AMC    schedule 20.02.2020
comment
Вы получаете разные результаты при каждом запуске (что кажется естественным, поскольку вы randint).   -  person CristiFati    schedule 28.02.2020

Ответы (1)


arrow_upward
3
arrow_downward

Я понял! Вся индексация возникла из вызова splitlines(), который у меня был в начале многих моих вспомогательных функций. Я также не отладил полностью свою score() функцию, поэтому моя оценка происходила не так, как меня просили сделать это в литературе.

все вспомогательные функции следующие:

def window(s, k):
    for i in range(1 + len(s) - k):
        yield s[i:i+k]

def HammingDistance(seq1, seq2):
    if len(seq1) != len(seq2):
        raise ValueError('Undefined for sequences of unequal length.')
    return sum(ch1 != ch2 for ch1, ch2 in zip(seq1, seq2))

def score(motifs):
    score = 0
    for i in range(len(motifs[0])):
        motif = ''.join(motifs[j][i] for j in range(len(motifs)))
        score += min([HammingDistance(motif,homogenous*len(motif)) for homogenous in 'ACGT'])
    return score

def profile(motifs):
    prof = []
    for i in range(len(motifs[0])):
        col = ''.join([motifs[j][i] for j in range(len(motifs))])
        prof.append([float(col.count(nuc))/float(len(col)) for nuc in 'ACGT'])
    return prof

def profile_most_probable_kmer(dna, k, prof):
    nuc_loc = {nucleotide:index for index,nucleotide in enumerate('ACGT')}
    motif_matrix = []
    max_prob = [-1, None]
    for i in range(len(dna)):
        motif_matrix.append(max_prob)
    for i in range(len(dna)):
        for chunk in window(dna[i],k):
            current_prob = 1
            for j, nuc in enumerate(chunk):
                current_prob*=prof[j][nuc_loc[nuc]]
            if current_prob>motif_matrix[i][0]:
                motif_matrix[i] = [current_prob,chunk]
    return list(list(zip(*motif_matrix))[1])

def profile_with_pseudocounts(motifs):
    prof = []
    for i in range(len(motifs[0])):
        col = ''.join([motifs[j][i] for j in range(len(motifs))])
        prof.append([float(col.count(nuc)+1)/float(len(col)+4) for nuc in 'ACGT'])
    return prof

from random import randint

def SampleMotifs(Dna,k,t):
    Dna = Dna.splitlines()
    BestMotifs = []
    for line in Dna:
        position = randint(0,len(line)-k)
        BestMotifs.append(line[position:position+k])
    return BestMotifs

def randomized_motif_search(dna,k,t):
    from random import randint
    from itertools import chain
    dna = dna.splitlines()
#     Randomly generate k-mers from each sequence in the dna list.
    rand_ints = [randint(0,len(dna[0])-(k)) for a in range(len(dna))]
    motifs = [dna[i][r:r+k] for i,r in enumerate(rand_ints)]
    best_score = [score(motifs), motifs]
    while True:
        current_profile = profile_with_pseudocounts(motifs)
        motifs = profile_most_probable_kmer(dna,k,current_profile)
#         motifs = []
#         for i in range(len(mfp)):
#             motifs.append(mfp[i][0])
        current_score = score(motifs)
        if current_score < best_score[0]:
            best_score = [current_score, motifs]
        else:
            return best_score[1]

def Laplace(dna,k,t):
    import random
    random.seed(0)
    i = 0

    LastMotifs = randomized_motif_search(dna,k,t)
    while i < 1000:
        BestMotifs = randomized_motif_search(dna,k,t)
        if score(BestMotifs)<score(LastMotifs):
            LastMotifs = BestMotifs
        i+=1
    print('\n'.join(LastMotifs))
person slap-a-da-bias    schedule 23.02.2020