3

我有一个多 FASTA 文件，其中包含来自下一代测序的 10 000 多个 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 文件，而不仅仅是一个代码序列。另外：如何根据需要有效地存储结果。

4

1 回答 1

2

第一步：生成要对齐的对

您可能只想将序列相互对齐一次，例如，如果您有序列 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从这个元组创建对象，以便能够将其保存到文件中，将分数添加为 adescription并保留nameand 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并写入 yieldedSeqRecords到一个打开的句柄：

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")
于 2019-08-20T07:14:42.257 回答