This is a test version of Biostars. For the public version, visit https://www.biostars.org.
trying to upgrade from Bio.pairwise2.align to Bio.Align

Hi folks,

I'm trying to upgrade some old code to use Bio.Align that currently uses Bio.pairwise2.align.

I'm working through the tests that exist in the code and most are failing. One case that fails expects these two alignments to return equal (i.e., ==) results:

from Bio.pairwise2 import align as bio_align
fasta_dict = {"query": "SCASSWVTVTTVSS", "target": "QVQLCXXXWGQGVTVSS"}
exp_fasta_dict = {"query": "---SCASSWVTVTTVSS", "target": "QVQLCXXXWGQGVTVSS"}

Here are the commands I am running for comparison :

nt_ms = 1
nt_mmp = -1
nt_gop = -3
nt_gep = -1

align1a = bio_align.localms(
    fasta_dict["query"],
    fasta_dict["target"],
    nt_ms,
    nt_mmp,
    nt_gop,
    nt_gep,
    one_alignment_only=True,
)[0]

align1b = bio_align.localms(
    exp_fasta_dict["query"],
    exp_fasta_dict["target"],
    nt_ms,
    nt_mmp,
    nt_gop,
    nt_gep,
    one_alignment_only=True,
)[0]

align1b==align1a #both seqAs returned have a leading "---" and are "---SCASSWVTVTTVSS"

Then I switch to Bio.Align:

from Bio import Align
aligner = Align.PairwiseAligner()

# Set the match score, mismatch score, gap open penalty, and gap extension penalty
aligner.match_score = nt_ms
aligner.mismatch_score = nt_mmp
aligner.open_gap_score = nt_gop
aligner.extend_gap_score = nt_gep
aligner.end_open_gap_score = nt_gop
aligner.end_extend_gap_score = nt_gep
aligner.mode = "local"

# Perform the alignment
aligns2 = aligner.align(fasta_dict["query"], fasta_dict["target"])
align2a = next(aligns2)  # returns best match first

aligns2 = aligner.align(exp_fasta_dict["query"], exp_fasta_dict["target"])
align2b = next(aligns2)  # returns best match first

But here align2a is:

print(align2a)
target           10 TVSS 14
                  0 ||||  4
query            13 TVSS 17

and align2b had different coodinates

print(align2b)
target           13 TVSS 17
                  0 ||||  4
query            13 TVSS 17

Is there a way to get Bio.Align to report the same results as `Bio.pairwise2.align?

pairwise2 bio.align python biopython

0 answers

No answers yet.

Log in to answer this question.