根据相似性的百分比将一个文件中的几个seq连接起来

>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

>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

2条回答

网友

1楼 · 编辑于 2024-04-20 11:27:05

你可以这样做。你知道吗

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

网友

2楼 · 编辑于 2024-04-20 11:27:05

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

然后我们定义helper函数，每个函数用于检查concat，然后concat用于seq，merge用于name（使用nest helper获取HSPs计数）：

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\u检查，并为下一次比较转发新的配对，必要时添加到最终列表中：

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

请注意，我只检查了文本文件中连续序列的串联，您必须重新使用我在不同脚本中编写的方法来说明组合。我让你自己决定。你知道吗

希望这有帮助。：）

相关问题更多 >

编程相关推荐

热门问题

热门文章