报告与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)

但我只想报告每一次阅读的最佳排列。如何从a2的分数中选择这种对齐方式。我想选择a2值最高的对齐方式。

Answer 1

您可以根据a2对对齐进行排序

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)

为了找到它，我不得不在文档中做一些调查，但这给出了最佳的分数。