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 回答