如何为生物信息学查询优化python脚本

Question

我对python还很陌生，如果可能的话，我将非常感谢您的帮助。我正在比较两种密切相关的生物的基因组，并试图找出一些基本的插入和缺失。我用两种生物的序列进行了FASTA成对比对（glsearch36）。在

下面是我的python脚本的一部分，我在其中识别了一个序列（数据库）中的7个核苷酸（七肽体），它对应于另一个序列（查询）中的间隙。这是我的一个例子：

ATGCACAA-ACCTGTATG # query
ATGCAGAGGAAGAGCAAG # database
9
GAGGAAG

假设间隙位于位置9。我正在尝试改进脚本，以选择两个序列上相距20个核苷酸或更多的间隙，并且仅当周围的核苷酸也匹配时

^{pr2}$

这是我脚本的一部分，上半部分处理打开不同的文件。它还打印了一个字典，在末尾列出了每个序列的计数。在

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