python - 根据相似度百分比连接文件中的多个序列

Question

你好，我需要你的帮助来完成一项复杂的任务。

这是一个file1.txt :

>Name1.1_1-40_-__Sp1
AAAAAACC-------------
>Name1.1_67-90_-__Sp1
------CCCCCCCCC------
>Name1.1_90-32_-__Sp1
--------------CCDDDDD
>Name2.1_20-89_-__Sp2
AAAAAACCCCCCCCCCC----
>Name2.1_78-200_-__Sp2
-------CCCCCCCCCCDDDD

想法是创建一个名为 file1.txt_Hsp 的新文件比如:

>Name1.1-3HSPs-__Sp1
AAAAAACCCCCCCCCCDDDDD
>Name3.1_-__Sp2
AAAAAACCCCCCCCCCC----
>Name4.1_-__Sp2
-------CCCCCCCCCCCCCC

所以基本上这个想法是:

• 比较每个序列from the same SpN <--(这里很重要，只有相同的 SpN 名称)彼此在 file1.txt 中.例如，我将不得不比较:
• Name1.1_1-40_-__Sp1 vs Name1.1_67-90_-__Sp1
• Name1.1_1-40_-__Sp1 vs Name1.1_90-32_-__Sp1
• Name1.1_67-90_-__Sp1 vs Name1.1_90-32_-__Sp1
• Name2.1_20-89_-__Sp2 vs Name2.1_78-200_-__Sp2

例如，当我比较时:

Name1.1_1-40_-__Sp1 vs Name1.1_67-90_-__Sp1我得到:

>Name1.1_1-40_-__Sp1
AAAAAACC-------------
>Name1.1_67-90_-__Sp1
------CCCCCCCCC------

如果ratio between number of letter matching with another letter / nb letter matching with a (-)，我想连接两个序列是 < 0.20`。

例如这里有21 characters , 以及与另一个字母匹配的字母数 = 2 ( C and C )。以及与 - 匹配的字母数, 是 13 ( AAAAAA+CCCCCCC )

所以

ratio = 2/15 : 0.1538462

如果这个ratio < 0.20然后我想连接这两个序列，例如:

>Name1.1-2HSPs_-__Sp1
AAAAAACCCCCCCCC------

(如您所知，新序列的名称现在是:Name.1-2HSPs_-__Sp1，其中 2 表示有 2 个序列连接在一起)所以我们删除 XHSPS 的数字部分，其中 X 是数字的序列连接。并获取 file1.txt_Hsp :

>Name1.1-2HSPs_-__Sp1
AAAAAACCCCCCCCC------
>Name1.1_90-32_-__Sp1
--------------CCDDDDD
>Name2.1_20-89_-__Sp2
AAAAAACCCCCCCCCCC----
>Name2.1_78-200_-__Sp2
-------CCCCCCCCCCDDDD

然后我用 Name1.1-2HSPs_-__Sp1 vs Name1.1_90-32_-__Sp1 再做一次

>Name1.1-2HSPs_-__Sp1
AAAAAACCCCCCCCC------
>Name1.1_90-32-__Sp1
--------------CCDDDDD

Where ratio = 1/20 = 0.05 

那么因为 ratio is < 0.20我想连接这两个序列，例如:

>Name1.1-3HSPs_-__Sp1
AAAAAACCCCCCCCCCDDDDD

(如您所见，新序列的名称现在是:Name.1-3HSPs_-__Sp1，其中 3 表示有 3 个序列连接在一起)

file1.txt_Hsp:

>Name1.1-3HSPs_-__Sp1
AAAAAACCCCCCCCCCDDDDD
>Name2.1_20-89_-__Sp2
AAAAAACCCCCCCCCCC----
>Name2.1_78-200_-__Sp2
-------CCCCCCCCCCDDDD

然后我用 Name2.1_20-89_-__Sp2 再做一次对比Name2.1_78-200_-__Sp2

>Name2.1_20-89_-__Sp2
AAAAAACCCCCCCCCCC----
>Name2.1_78-200_-__Sp2
-------CCCCCCCCCCDDDD

Where ratio = 10/11 = 0.9090909

那么因为 ratio is > 0.20我什么都不做，得到了最后的 file1.txt_Hsp :

>Name1.1-3HSPs_-__Sp1
AAAAAACCCCCCCCCCDDDDD
>Name2.1_20-89_-__Sp2
AAAAAACCCCCCCCCCC----
>Name2.1_78-200_-__Sp2
-------CCCCCCCCCCDDDD

这是我需要的最终结果。

---

一个最简单的例子是:

>Name1.1_10-60_-__Seq1
AAA------
>Name1.1_70-120_-__Seq1
--AAAAAAA
>Name2.1_12-78_-__Seq2
--AAAAAAA

比率是1/8 = 0.125因为只有 1 个字母匹配而 8 因为 8 个字母与 (-) 匹配

因为 ratio < 0.20我将两个序列 Seq1 连接到:

>Name1.1_2HSPs_-__Seq1
AAAAAAAAA

新文件应该是:

>Name1.1_2HSPs_-__Seq1
AAAAAAAAA
>Name2.1_-__Seq2
--AAAAAAA

** 这是我真实数据的一个例子**

>YP_009186705
MMSCQSWMMKYFTKVCNRSNLALPFDQSVNPVSFSMISSHDVMLKLDDEIFYKSLNQSNL
ALPFDQSVNPVSFSMISSHDLIA
>XO009980.1_26784332-20639090_-__Agapornis_vilveti
------------------------------------------------------LNQSNL
ALPFDQSVNPVSFSMISSHDLIA
>CM009917.1_20634332-20634508_-__Neodiprion_lecontei
---CDSWMIKFFARISQMC---IKIHSKYEEVSFFLFQSK--KKKIADSHFFRSLNQDTA
-------LNTVSY----------
>XO009980.1_20634508-20634890_-__Agapornis_vilveti
MMSCQSWMMKYFTKVCNRSNLALPFDQSVNPVSFSMISSHDVMLKL--------------
-----------------------
>YUUBBOX12
MMSCQSWMMKYFTKVCNRSNLALPFDQSVNPVSFSMISSHDVMLKLDDEIFYKSLNQSNL
ALPFDQSVNPVSFSMISSHDLIA

我应该得到:

>YP_009186705
MMSCQSWMMKYFTKVCNRSNLALPFDQSVNPVSFSMISSHDVMLKLDDEIFYKSLNQSNL
ALPFDQSVNPVSFSMISSHDLIA
>XO009980.1_2HSPs_-__Agapornis_vilveti
MMSCQSWMMKYFTKVCNRSNLALPFDQSVNPVSFSMISSHDVMLKLLNQSNL
ALPFDQSVNPVSFSMISSHDLIA
>CM009917.1_20634332-20634508_-__Neodiprion_lecontei
---CDSWMIKFFARISQMC---IKIHSKYEEVSFFLFQSK--KKKIADSHFFRSLNQDTA
-------LNTVSY----------
>YUUBBOX12
MMSCQSWMMKYFTKVCNRSNLALPFDQSVNPVSFSMISSHDVMLKLDDEIFYKSLNQSNL
ALPFDQSVNPVSFSMISSHDLIA

XO009980.1_26784332-20639090_-__Agapornis_vilveti 之间的比率和 XO009980.1_20634508-20634890_-__Agapornis_vilveti是:0/75 = 0

如你所见，一些序列没有[\d]+[-]+[\d] YP_009186705 等模式或 YUUBBOX12 ，这些不必连接，只需将它们添加到输出文件中即可。

非常感谢您的帮助。

Answer 1

首先，让我们将文本文件读入(name, seq)的元组中:

with open('seq.txt', 'r+') as f:
    lines = f.readlines()
seq_map = []
for i in range(0, len(lines), 2):
    seq_map.append((lines[i].strip('\n'), lines[i+1].strip('\n')))

#[('>Name1.1_10-60_-__Seq1', 'AAA------'),
# ('>Name1.1_70-120_-__Seq1', '--AAAAAAA'),
# ('>Name2.1_12-78_-__Seq2', '--AAAAAAA')]
#
# or
#
# [('>Name1.1_1-40_-__Sp1', 'AAAAAACC-------------'),
#  ('>Name1.1_67-90_-__Sp1', '------CCCCCCCCC------'),
#  ('>Name1.1_90-32_-__Sp1', '--------------CCDDDDD'),
#  ('>Name2.1_20-89_-__Sp2', 'AAAAAACCCCCCCCCCC----'),
#  ('>Name2.1_78-200_-__Sp2', '-------CCCCCCCCCCDDDD')]

然后我们定义辅助函数，每个辅助函数用于检查 concat，然后 concat 用于 seq，合并用于名称(使用嵌套辅助函数来获取 HSP 计数):

import re

def count_num(x):
    num = re.findall(r'[\d]+?(?=HSPs)', x)
    count = int(num[0]) if num and 'HSPs' in x else 1
    return count

def concat_name(nx, ny):
    count, new_name = 0, []
    count += count_num(nx)
    count += count_num(ny)
    for ind, x in enumerate(nx.split('_')):
        if ind == 1:
            new_name.append('{}HSPs'.format(count))
        else:
            new_name.append(x)
    new_name = '_'.join([x for x in new_name])
    return new_name

def concat_seq(x, y):
    mash, new_seq = zip(x, y), ''
    for i in mash:
        if i.count('-') > 1:
            new_seq += '-'
        else:
            new_seq += i[0] if i[1] == '-' else i[1]
    return new_seq

def check_concat(x, y):
    mash, sim, dissim = zip(x, y), 0 ,0
    for i in mash:
        if i[0] == i[1] and '-' not in i:
            sim += 1
        if '-' in i and i.count('-') == 1:
            dissim += 1
    return False if not dissim or float(sim)/float(dissim) >= 0.2 else True

然后我们将编写一个脚本来按顺序运行元组，检查 spn 匹配，然后是 concat_checks，并为下一次比较采用新的配对，在必要时添加到最终列表:

tmp_seq_map = seq_map[:]
final_seq = []

for ind in range(1, len(seq_map)):
    end = True if ind == len(seq_map)-1 else False
    pair_a = tmp_seq_map[ind-1]
    pair_b = tmp_seq_map[ind]

    name_a = pair_a[0][:]
    name_b = pair_b[0][:]

    if name_a.split('__')[1] == name_b.split('__')[1]:

        if check_concat(pair_a[1], pair_b[1]):

            new_name = concat_name(pair_a[0], pair_b[0])
            new_seq = concat_seq(pair_a[1], pair_b[1])
            tmp_seq_map[ind] = (((new_name, new_seq)))

            if end:
                final_seq.append(tmp_seq_map[ind])
                end = False
        else:
            final_seq.append(pair_a)
    else:
        final_seq.append(pair_a)
    if end:
        final_seq.append(pair_b)
print(final_seq)

#[('>Name1.1_2HSPs_-__Seq1', 'AAAAAAAAA'),
# ('>Name2.1_12-78_-__Seq2', '--AAAAAAA')]
#
# or
#
#[('>Name1.1_3HSPs_-__Sp1', 'AAAAAACCCCCCCCCCDDDDD'),
# ('>Name2.1_20-89_-__Sp2', 'AAAAAACCCCCCCCCCC----'),
#  ('>Name2.1_78-200_-__Sp2', '-------CCCCCCCCCCDDDD')]

请注意，我只检查了文本文件中连续序列的串联，您将不得不重新使用我在不同脚本中编写的方法来计算组合.我将其留给您自行决定。

希望这对您有所帮助。 :)