Extracting gapped sequences from Bio.Align.PairwiseAligner

Question

[biomarco]

Hi everyone, I noticed that the Bio.pairwise2 module has been deprecated, so I wanted to switch to Bio.Align.PairwiseAligner in my scripts, as suggested by the warning message.

I used to do the following to calculate sequence identity and query coverage:

aln1, aln2, _, start_aln, end_aln = pairwise2.align.localds(x_seq, y_seq, blosum_matrix, -11, -1)[0]
identities = len([i for i,k in zip(aln1[start_aln:end_aln], aln2[start_aln:end_aln]) if i == k and i != '-'])
perc_ident = round(100 * (identities / (end_aln - start_aln)), 1)
query_cov = round(100 * ((end_aln - start_aln) / len(aln1)), 1)

With PairwiseAligner it looks like you can just print the aligned sequences, and I could not find a way to access the query and target gapped sequences (aln1 and aln2 variables above), so there is no way for me to calculate the sequence identity.

Can somebody help?

[Mensur Dlakic]

I think this point of view is different among generations: I never update a program or an operating system unless there is a very good reason to do so. To be clear, in my book having a newer version is not a good reason to update.

So while you are waiting for someone to answer your question, you can always go back to the BioPython version that worked for you.

[biomarco]

I get your point, but if I can get rid of warning messages without silencing or ignoring them, I'm very happy to do so, especially if a module I'm using has been deprecated (I assume they had a good reason to do so). Anyway the current BioPython version is still working fine, I just get a warning message when I import pairwise2 suggesting to use something else and to let the developers know if you really still need that module. That's why I'm wondering whether it is possible to do the same I used to do with PairwiseAligner instead of pairwise2.

[Ram]

Do not delete posts that have received feedback. Redacting content is extremely bad etiquette and if done again, will get your user account suspended.

I've restored the content from Google cache.

[Joe]

FWIW, there are good reasons beyond it simply being a new version to update here. Pairwise in particular has been overhauled quite a bit in more recent versions of Biopython.

Answer 1

[biomarco]

Sorry for deleting my question, I was not sure anymore it was a question worthy to be asked. Won't do that again! I will also post the solution I found, for anyone stumbling across the same problem:

from Bio.Align import substitution_matrices, PairwiseAligner

aligner_local = PairwiseAligner(
    mode = 'local',
    open_gap_score = -11,
    extend_gap_score = -1,
    substitution_matrix = substitution_matrices.load("BLOSUM62")
    )

def get_pairwise_identity(query: str, subject: str, aligner: PairwiseAligner):
    # NOTE using local aligner is recommended!
    alignment = aligner.align(query, subject)[0]  # looks just at the first alignment
    aln1, aln2 = alignment  ###### <--- THIS LINE ANSWERS MY INITIAL QUESTION! ######
    identities = len([i for i, k in zip(aln1, aln2) if i == k and i != '-'])
    identity = 100 * identities / len(aln1)
    query_coverage = 100 * (alignment.coordinates[0][-1] - alignment.coordinates[0][0]) / len(query)
    subject_coverage = 100 * (alignment.coordinates[1][-1] - alignment.coordinates[1][0]) / len(subject)
    return round(identity, 2), round(query_coverage, 2), round(subject_coverage, 2)

identity, q_cov, s_cov = get_pairwise_identity(seq1, seq2, aligner_local)

[Ram]

Thank you for posting the solution - that is the right way to go!

[Ram]

You can accept your own answer to resolve the post.