python - 根据相似性百分比连接文件中的多个 seq

标签 python concatenation string-concatenation

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

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

这里我想连接两个序列 if 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 ，这些不必串联，只需将它们添加到输出文件中即可。

非常感谢您的帮助。

最佳答案

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

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

希望这有帮助。 :)

关于python - 根据相似性百分比连接文件中的多个 seq，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/59114817/

上一篇：python - 使用seaborn绘制与外部数据源的标准差

下一篇：python - 在 Pandas 中创建类似相关矩阵的数据框

相关文章：

python - 无法在不转换为 ascii 的情况下拆分 unicode 字符串 - python 2.7

python - __future__ 模块在 exec 中不生效

python - 通过 Pandas 连接(许多)CSV 文件

c# - 将两个锯齿状列表合并为一个

javascript - 字符串数组的逐元素连接

python - 字符串连接 Python 2.5(使用 cStringIO)

python - 为什么 Apache Beam python 中的 GroupByKey 之后的 FlatMap 这么慢？

python - 在 for-in 循环之间添加暂停

java - 为什么 ".concat(String)"比 "+"快这么多？

sql - 通过 xml 路径连接字符串

©2024 IT工具网  联系我们