python - 多 FASTA 文件序列的成对比对

Question

我有一个多 FASTA 文件，其中包含超过 10000 个由下一代测序产生的 fasta 序列，我想将每个序列与文件内的每个序列进行成对比对，并将所有结果存储在同一个新文件中，以便之后进行聚类分析。 FASTA 序列的示例和我使用 python 执行成对序列比对的代码如下所示。

FASTA序列

>m180921_230442_42149_c101464342550000001823297908121882_s1_X0/538/ccs
AGAAGCCACTCATCCATCCAGGCAGGAAGACTCTTAGGATCCTGACTTTCTCCTGGTCCCCACATCCCCT
AAACCGAGGAAGGGGTCCAGCAGGGTCCGAGTCCCTGAAGCAAGGATTCTCCGTGGTCGTGTCCCCACAG

请忽略第一行，因为它包含序列的描述摘要。

我的代码

    from Bio import pairwise2
    from Bio.pairwise2 import format_alignment

    X = "ACGGGT"
    Y = "ACG"

    #match score = 2, mismatch score = -1, gap opening = -5, gap extension = -2
    alignments = pairwise2.align.globalms(X, Y, 2, -1, -5, -2)

    for a in alignments:
        print(format_alignment(*a))

问题

我想知道如何修改它以遍历整个多 FASTA 文件而不仅仅是一个代码序列。 另外:如何根据需要有效地存储结果。

Answer 1

第一步:生成要对齐的对

您可能只想将序列相互比对一次，例如如果你有序列 1、2 和 3，你只想对齐 1vs2、1vs3 和 2vs3(即所有组合)，丢弃 2vs1 和 3vs2 以及自对齐。这将为您节省一些运行时间:

from itertools import combinations
from Bio import SeqIO

combinations(SeqIO.parse(infile, "fasta"), 2)

第 2 步:对齐生成的对

函数 pairwise2.align.globalms 返回一个 (seqA, seqB, score, begin, end) 的元组。我们需要从这个元组创建 SeqRecord 对象，以便能够将其保存到文件中，将分数添加为 description 并保留 name 和 id:

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord

def align(pair):
    """Yield two Seqs aligned."""
    al = pairwise2.align.globalms(pair[0].seq, pair[1].seq, 2, -1, -5, -2)[0]

    yield SeqRecord(Seq(al[0]), id=pair[0].id, name=pair[0].name,
                    description=f"Score={al[2]}")
    yield SeqRecord(Seq(al[1]), id=pair[1].id, name=pair[1].name,
                    description=f"Score={al[2]}")

第 3 步:将它们缝合在一起

可以看出上面的函数都是生成器。 Biopython writer 干净地处理生成的序列，所以你可以只要求在第一个函数中生成对，将它发送到 align 并将产生的 SeqRecords 写入一个打开的句柄:

input = "your_fasta.fas"

with open("output.fas", "w") as output:
    for pair in combinations(SeqIO.parse(infile, "fasta"), 2):
        SeqIO.write(align(pair), output, "fasta")