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

Question

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

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

所以基本上这个想法是：

• 在file1.txt中将每个序列from the same SpN<--(此处仅具有相同的SpN名称非常重要(相互比较。 例如，我将不得不比较：
• 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')]

请注意，我已经检查了文本文件中仅连续序列的串联，并且您必须重用我用不同脚本编写的方法来解释组合。我让你们自行决定。

希望这有帮助。 :)

Answer 2

您可以按如下方式执行此操作。

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