python - 如何为生物信息学查询优化 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 nt 以尝试减少噪音来解决我之前的一个问题，这改善了我的结果。上面编辑的脚本。

这是我的结果的一个例子，最后我有一个字典，它计算每个序列的出现次数。

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”替换为您想比较的字符串。