报告与pairwise2的最佳对齐

Question

我有一个fastq的读取文件，比如“reads.fastq”。我想将sequnces与保存为fasta文件ref.faa的字符串对齐。我正在使用以下代码

reads_array = []

for x in Bio.SeqIO.parse("reads.fastq","fastq"):
    reads_array.append(x)

for x in Bio.SeqIO.parse("ref.faa","fasta"):
    refseq = x

result = open("alignments_G10_final","w")

aligned_reads = []

for x in reads_array:
    alignments =pairwise2.align.globalms(str(refseq.seq).upper(),str(x.seq),2,-1,-5,-0.05)
    for a in alignments:
        result.write(format_alignment(*a))
        aligned_reads.append(x)

但我想报告每次读取的最佳对齐方式。如何从[2]中的分数中选择这种对齐方式。我想选择具有最高值a [2]的对齐方式

Answer 1

您可以根据[2]对对齐进行排序：

for x in reads_array:
    alignments =pairwise2.align.globalms(str(refseq.seq).upper(),str(x.seq),2,-1,-5,-0.05)
    sorted_alignments = sorted(alignments, key=operator.itemgetter(2))

    result.write(format_alignment(*sorted_alignments[0]))
    aligned_reads.append(x)

Answer 2

我知道这是一个老问题，但对于仍在寻找正确答案的人，请在你的对齐方法中添加one_alignment_only=True参数：

alignments =pairwise2.align.globalms(str(refseq.seq).upper(),
    str(x.seq),
    2,-1,-5,-0.05,
    one_alignment_only=True)

我不得不在文档中进行一些挖掘以找到它，但这会给出最佳分数。