在Python中连接两个fasta文件

问题描述 投票:0回答:1

我有两个数据文件(FASTA),每个文件代表一个基因,序列由物种和本地识别。我想将这些文件连接成一个例子:

psbki.fas:
>E_oleracea_Docas_de_Belm 
AACCT
ycf1b.fas: 
>E_oleracea_Docas_de_B 
GGTTC

output:
>E_oleracea_Docas_de_Belm 
AACCTGGTTC

如果你在两个文件中查看物种的名称,它们会被写成一些语法问题,然后彼此不同。另外,我还有另一个问题:某些物种不在两个文件中。

为了解决这些问题,我编写了以下代码:

ids, sequences = parse_fasta(open('psbki.fas', 'r').read().split('\n'))
ids2, sequences2 = parse_fasta(open('ycf1b.fas', 'r').read().split('\n'))

for i, j, z, h in zip(ids, sequences, sequences2, ids2):
    if i != h:
        print(">"+i + "\n"+j)
    else:
        print(">"+i + "\n"+j+z)

两个第一序列的输出是可以的。但对于其他序列,代码只打印一个文件中的文件,但它们都在两个文件中。我的代码出了什么问题?我是python的初学者

输出是:

>E_edulis_I1
AAATAGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTTTGGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAAG----ATCTTGAATTGGTATTATTCTCATTATCAGCATAAATTATCACACGTCTGGCTCTTCTTGAACGAATTTCAATATCTTCTATCGGTTTTTCCTCATTTTCTTCCTCCTGTTCTTCCAGAAGATTGGTCAATTTATATGACCATCGAGAAACCTTTTTACTGATTTCTTCTATTCCAATAGATTCATTTCTAGTTGTTTTATCATTTGGATCAATTGTCATTATATCGAATACAAATTTCAAAGATTTTGCTTGACTTTCTGAATCCATTTTTCTTTGTTCTGCCAATAAAGAACAGTTTTTCAAACAAAAATTGGGTGTGAATTCAAAAGAAAATGAAGTTAAGGAATTACCGATATAATTCAAAAATGATTTACCACCACCAAGTGAATTCTTTTGATGTTCAAATTCTCTGAAATTATTAGGAAGTAGCTCATGGATCTTATTTATCCAAAGACTTTTTATGGAATCCTCCATATAAGGGAAAAAATCATTTATGATTGTACGTAAATCAAAATCTTTTATTGCTCCACGGCATGGTCCGCTCAATAAAGGATCATATGTTTTGGTCAAGCATTTTTGTTTATTCTCATGATTGCAAAATCTAGTCTTTTTTTCGAGCATATCTAGAGCAAGAAATCCCTTTTCTTTTTTTTCTTTTTCTAGAGCTTTTATTCGACTTATTAATTCATTGCTCAAGTTGTATTTTTTTTGTTCATTGGTAAAAACCCAAAAATTATACAGGTCTCCATGGGATAATTTTTT-GTCGTGTACAAAAACATTTTTCGTTCTATCATTTCC
>E_edulis_I2
AAATAGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTTTGGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAAG----ATCTTGAATTGGTATTATTCTCATTATCAGCATAAATTATCACACGTCTGGCTCTTCTTGAACGAATTTCAATATCTTCTATCGGTTTTTCCTCAATTTCTTCCTCCTGTTCTTCCAGAAGATTGGTCAATTTATATGACCATCGAGAAACCTTTTTACTGATTTCTTCTATTCCAATAGATTCATTTCTAGTTGTTTTATCATTTGGATCAATTGTCATTATATCGAATACAAATTTCAAAGATTTTGCTTGACTTTCTGAATCCATTTTTCTTTGTTCTGCCAATAAAGAACAGTTTTTCAAACAAAAATTGGGTGTGAATTCAAAAGAAAATGAAGTTAAGGAATTACCGATATAATTCAAAAATGATTTACCACCACCAAGTGAATTCTTTTGATGTTCAAATTCTCTGAAATTATTAGGAAGTAGCTCATGGATCTTATTTATCCAAAGACTTTTTATGGAATCCTCCATATAAGGGAAAAAATCATTTATGATTGTACGTAAATCAAAATCTTTTATTGCTCCACGGCATGGTCCGCTCAATAAAGGATCATATGTTTTGGTCAAGCATTTTTGTTTATTCTCATGATTGCAAAATCTAGTCTTTTTTTCGAGCATATCTAGAGCAAGAAATCCCTTTTCTTTTTTTTCTTTTTCTAGAGCTTTTATTCGACTTATTAATTCATTGCTCAAGTTGTATTTTTTTTGTTCATTGGTAAAAACCCAAAAATTATACAGGTCTCCATGGGATAATTTTTTTGTCGTGTACAAAAACATTTTTCGTTCTATCATTTCC
>E_edulis_F7
AAATAGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCTTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTT-GGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAA-G----ATCTTG
>E_edulis_R10
AAATAGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTTTGGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAAG----ATCTTG
>E_edulis_R11
AAATAGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTTTGGGGGTGTCATGTCAAAACAAAATAGKGTATGTGGTAAAGTAAAAAATAASTAACCTATTCCCTTTTTCAAAAAAAAAAG----ATCTTG
>E_edulis_R12
AAATAGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTTTGGGGGWGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAAA----ATCTTG
>E_edulis_IFES
AAATAGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGARCAAAGACTTTATTAGGTTGCTTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTT-GGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAA-G----ATCTTG
>E_oleracea_Ilha_do_combu_1
AAATCGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTT-GGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAGATCTTATCTTG
>E_oleracea_Ilha_do_combu_2
AAATCGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTT-GGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAGATCTTATCTTG
>E_oleracea_Ilha_do_combu_3
AAATCGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTT-GGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAGATCTTATCTTG
>E_oleracea_Ilha_do_combu_5
AAATCGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTT-GGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAGATCTTATCTTG
>E_oleracea_Ilha_do_combu_10
AAATCGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTT-GGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAGATCTTATCTTG
>E_oleracea_Mangal_2
AAATCGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTT-GGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAGATCTTATCTTG
>E_oleracea_Mangal_3
AAATCGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTT-GGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAGATCTTATCTTG
>E_oleracea_Docas_de_Belm
AAATCGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTT-GGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAGATCTTATCTTG
>E_oleracea_Utinga
AAATCGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTT-GGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAGATCTTATCTTG
>E_oleracea_Canto_de_Roa_1
AAATCGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTT-GGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAGATCTTATCTTG
>E_oleracea_Canto_de_Roa_2
AAATCGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTT-GGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAGATCTTATCTTG
>E_oleracea_Canto_de_Roa_3
AAATCGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTT-GGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAGATCTTATCTTG
>E_oleracea_IFES
AAATCGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTT-GGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAGATCTTATCTTG

换句话说,我想连接两个文件中的基因而不是打印,而不是连接只出现在一个文件中的物种。我不知道如何解决物种写错的问题。

编辑1:我已经改变了代码,使用Levenshtein比率来解决某些物种名称中的写入错误,但输出是相同的。

新代码是:

import Levenshtein as lev
Str1 = str(ids)
Str2 = str(ids2)
Ratio = lev.ratio(Str1.lower(),Str2.lower())

for i, j, z, h in zip(ids, sequences, sequences2, ids2):
    if lev.ratio(i,h) > 0.70 and i in h:
         print(">"+i + "\n"+j+z)
    else:
        print(">"+i + "\n"+j)

编辑2

    Input File1: gene 1
>E_edulis_I1
AAATAGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTTTGGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAAG----ATCTTG
>E_edulis_I2
AAATAGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTTTGGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAAG----ATCTTG
>E_edulis_F7
AAATAGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCTTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTT-GGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAA-G----ATCTTG

Input File 2: gene 2
>E_edulis_I1
AAATAGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTTTGGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAAG----ATCTTG
>E_ed_I2
AAATAGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTTTGGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAAG----ATCTTG

My desired output:
 >E_edulis_I1
    AAATAGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTTTGGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAAG----ATCTTGAAATAGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTTTGGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAAG----ATCTTG
    >E_edulis_I2
    AAATAGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTTTGGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAAG----ATCTTGAAATAGAAATTCTTGTATATTGAATAACCGCGGCGATGAATTTTGATCAACTTATTTCCTCGTTCTGACCTTACAGTGAGCAAAGACTTTATTAGGTTGCCTACAATACCTAATTATTCATATGACAAGAAATTTTTGATAACGAAGGAATCAAAATCTTATTCCAAAGAAATTCGTGAAAATGACTTTCTTTTCAAAAAACACTTCATTTTTTTTGGGGGTGTCATGTCAAAACAAAATAGTGTATGTGGTAAAGTAAAAAATAAGTAACCTATTCCCTTTTTCAAAAAAAAAAG----ATCTTG

附:在第二个文件中,我有相同的物种E_edulis_I2,名称不完整 - > E_ed_I2。我希望脚本能够识别并将序列与第一个序列连接起来(文件1 = E_edulis_I2)。另一个问题是E_edulis_F7硬币只出现在文件1中,所以我不想在我的输出中使用thar specie。

python for-loop if-statement concatenation fasta
1个回答
0
投票

这比我想象的要复杂一点。问题是,例如,"E_edulis_I1""E_edulis_I2"更接近"E_ed_I"。我认为最好的解决方案是比较两个文件之间的每对名称,并确定它们是否足够相似,称为匹配。然后,一旦您拥有这组匹配的名称及其相似性级别,您就可以按照最相似和最不相似的顺序进行搜索,并在结合时将结果添加到组合的FASTA文件中。因此,因为两个文件之间的"E_edulis_I1"完全匹配,所以将首先放入组合的FASTA中。然后,当我们找到匹配的名称"E_edulis_I1""E_edulis_I2"时,我们会看到我们已经使用过"E_edulis_I1",所以这对不能匹配。

这对我来说似乎有点脆弱。您只需要注意相似度函数和相似度阈值是什么。您可以添加的一件事是在匹配完成时打印出名称而没有相似度1.这样,您可以快速扫描这些(希望不会太多)并确定是否有任何名称匹配已经过了。

无论如何,这是代码。它至少可以用你给出的例子。

from pathlib import Path
from typing import Dict
import Levenshtein as lev

def fasta_to_dict(path: str) -> Dict[str, str]:
    """Read in a fasta file and return a dict mapping names to sequences."""
    fasta = Path(path).read_text().split("\n")
    return {fasta[i]: fasta[i + 1] for i in range(0, len(fasta) - 1, 2)}

# put in the paths to your fasta files here
fasta1 = fasta_to_dict("f1.fasta")
fasta2 = fasta_to_dict("f2.fasta")

def get_similarity(name1: str, name2: str) -> float:
    """Return a number describing how similar the names are."""
    # this could be any string comparison function you want
    return lev.ratio(name1, name2)

similarity_threshold = 0.7 # must be at least this similar to be called a match

def get_correct_name(name1: str, name2: str) -> str:
    """Determine which of 2 close-match names is the correct one."""
    # this can also be whatever function makes sense in your application
    # I'm just using the longer name as the "correct" one
    return name1 if len(name1) > len(name2) else  name2

possible_matches = []
for name1 in fasta1:
    for name2 in fasta2:
        similarity = get_similarity(name1, name2)
        if similarity > similarity_threshold:
            possible_matches.append((name1, name2, similarity))

# sort by most similar matches first
possible_matches = sorted(possible_matches, key=lambda match: -match[-1])

combined_fasta = {}
used_names1 = set()
used_names2 = set()
for name1, name2, _ in possible_matches:
    if name1 in used_names1 or name2 in used_names2:
        continue
    correct_name = get_correct_name(name1, name2)
    combined_fasta[correct_name] = fasta1[name1] + fasta2[name2]
    used_names1.add(name1)
    used_names2.add(name2)

# put the path to your output fasta file here
with open("combined.fasta", "w") as f:
    for name, seq in combined_fasta.items():
        f.write(name + "\n")
        f.write(seq + "\n")
© www.soinside.com 2019 - 2024. All rights reserved.