3

我对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

对此的任何帮助将不胜感激！

4

1 回答 1

1

我想我理解你想要做什么，但正如@alko 所说 - 你的代码中的注释肯定会有很大帮助。

至于在间隙周围找到完全匹配，您可以运行字符串比较：

类似于以下内容：

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

您需要将“查询”和“数据库”替换为您所称的要比较的字符串。

于 2013-11-06T15:43:13.427 回答