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

Concatenate several seq within a file according to their percentage of similarities

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

这是一个 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

所以基本上这个想法是:

例如,当我比较时:

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

(如您所知,新序列的名称现在是:名称。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我想连接这2个序列如:

>Name1.1-3HSPs_-__Sp1
AAAAAACCCCCCCCCCDDDDD

(如您所见,新序列的名称现在是:名称。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_-__Sp2Name2.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_vilvetiXO009980.1_20634508-20634890_-__Agapornis_vilveti 之间的比率是:0/75 = 0

这里如你所见,一些序列没有[\d]+[-]+[\d]模式,例如YP_009186705YUUBBOX12,这些序列不必连接,它们只需要连接被添加到输出文件中。

非常感谢您的帮助。

您可以按如下方式进行。

from collections import defaultdict
with open('lines.txt','r') as fp:
    lines=fp.readlines()

dnalist = defaultdict(list)
for i,line in enumerate(lines):
    line = line.replace('\n','')
    if i%2: #'Name' in line:
        dnalist[n].append(line)
    else:
        n = line.split('-')[-1]

这为您提供了一个字典,其中的键是文件编号,值是列表中的 dna 序列。

def calc_ratio(str1,str2):
    n_skipped,n_matched,n_notmatched=0,0,0
    print(len(str1),len(str2))
    for i,ch in enumerate(str1):
        if ch=='-' or str2[i]=='-':
            n_skipped +1
        elif ch == str2[i]:
            n_matched += 1
        else:
            n_notmatched+=1
    retval = float(n_matched)/float(n_matched+n_notmatched+n_skipped)
    print(n_matched,n_notmatched,n_skipped)
    return retval

这就是你的比例;您可能想考虑序列中的字符不匹配(“-”也不匹配)的情况,这里我假设这与“-”没有不同。

连接字符串的辅助函数:这里我采用了不匹配字符的情况并放入 'X' 来标记它(如果它曾经发生过)。

def dna_concat(str1,str2):
    outstr=[]
    for i,ch in enumerate(str1):
        if ch!=str2[i]:
            if ch == '-':
                outchar = str2[i]
            elif str2[i] == '-':
                outchar = ch
            else:
                outchar = 'X'
        else:
            outchar = ch
        outstr.append(outchar)
    outstr = ''.join(outstr)
    return outstr

最后,在另一个字典中以文件编号作为键,将串联列表作为值,通过字典列表进行循环以获得串联的答案。

for filenum,dnalist in dnalist.items():
    print(dnalist)
    answers = defaultdict(list)
    for i,seq in enumerate(dnalist):
        for seq2 in dnalist[i+1:len(dnalist)]:
            ratio = calc_ratio(seq,seq2)
            print('i {} {} ration {}'.format(seq,seq2,ratio))
            if ratio<0.2:
                answers[filenum].append(dna_concat(seq,seq2))
                print(dna_concat(seq,seq2))

首先,让我们将文本文件读入 (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')]

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

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