Как усовершенствовать скрипт на Python для запроса биоинформатики

Question

Я довольно новичок в Python, и я был бы благодарен за некоторую помощь, если это возможно. Я сравниваю геномы двух близкородственных организмов [E_C & E_F] и пытаюсь определить некоторые основные вставки и делеции. Я провел парное выравнивание FASTA (glsearch36), используя последовательности из обоих организмов.

Ниже приведен раздел моего скрипта на python, где я смог идентифицировать 7 нуклеотидов (гептамер) в одной последовательности (базе данных), что соответствует пробелу в другой последовательности (запрос). Это пример того, что у меня есть:

ATGCACAA-ACCTGTATG # query
ATGCAGAGGAAGAGCAAG # database
9
GAGGAAG

Предположим, что разрыв находится в положении 9. Я пытаюсь уточнить сценарий, чтобы выбрать разрывы, которые расположены на расстоянии 20 или более нуклеотидов в обеих последовательностях, и только в том случае, если окружающие нуклеотиды также совпадают.

ATGCACAAGTAAGGTTACCG-ACCTGTATGTGAACTCAACA
                 ||| |||
GTGCTCGGGTCACCTTACCGGACCGCCCAGGGCGGCCCAAG
21
CCGGACC

Это раздел моего скрипта, верхняя часть которого посвящена открытию разных файлов. он также печатает словарь с количеством каждой последовательности в конце.

list_of_positions = []

for match in re.finditer(r'(?=(%s))' % re.escape("-"), dict_seqs[E_C]): 
    list_of_positions.append(match.start())

set_of_positions = set(list_of_positions)
for position in list_of_positions:
    list_no_indels = []
    for number in range(position-20, position) :
        list_no_indels.append(number)
    for number in range(position+1, position+21) :
        list_no_indels.append(number)
    set_no_indels = set(list_no_indels)
    if len(set_no_indels.intersection(set_of_positions))> 0 : continue

    if len(set_no_indels.intersection(set_of_positions_EF))> 0 : continue


    print position 
    #print match.start()

    print dict_seqs[E_F][position -3:position +3]

    key = dict_seqs[E_F][position -3: position +3]

    if nt_dict.has_key(key):
        nt_dict[key] += 1 
    else:
        nt_dict[key] = 1


print nt_dict

По сути, я пытался отредактировать результаты парных выравниваний, чтобы попытаться идентифицировать нуклеотиды напротив пробелов в последовательности запроса и базы данных, чтобы провести некоторый базовый анализ вставки / удаления.

Я смог решить одну из моих предыдущих проблем, увеличив расстояние между пробелами "-" до 20 нт, пытаясь уменьшить шум, это улучшило мои результаты. Скрипт отредактирован выше.

Это пример моих результатов, и в конце у меня есть словарь, в котором подсчитаны случаи каждой последовательности.

ATGCACAA-ACCTGTATG # query
ATGCAGAGGAAGAGCAAG # database
9 (position on the sequence)
GAGGAA (hexamer)


ATGCACAAGACCTGTATG # query
ATGCAGAG-AAGAGCAAG # database
9 (position)
CAAGAC (hexamer)

Тем не менее, я все еще пытаюсь исправить сценарий, где я получаю нуклеотиды вокруг пробела, чтобы точно соответствовать, как это, где | просто чтобы показать соответствующие NT на каждой последовательности:

GGTTACCG-ACCTGTATGTGAACTCAACA # query
     ||| ||
CCTTACCGGACCGCCCAGGGCGGCCCAAG # database

9
ACCGAC

Любая помощь с этим будет с благодарностью!

Answer 1

Я думаю, что понимаю, что вы пытаетесь сделать, но, как сказал @alko, комментарии в вашем коде определенно помогут.

Что касается нахождения точного соответствия вокруг пропуска, вы можете выполнить сравнение строк:

Что-то вроде:

if query[position -3: position] == database[position -3: position] and query[position +1: position +3] == database[position +1: position +3]:
   # Do something

Вам нужно будет заменить "query" и "database" на то, что вы назвали своими строками, которые вы хотите сравнить.