Re: [Biopython] Comparing two already aligned sequences quickly to count gaps/matches/mismatches
Michiel de Hoon <[email protected]> Thu, 6 Feb 2025 10:36:01 +0000 (UTC)
| Newsgroups | gmane.comp.python.bio.general |
|---|---|
| Message-ID | <[email protected]> |
--===============5283190597617378020==
Content-Type: multipart/alternative;
boundary="----=_Part_9218259_733540587.1738838161137"
------=_Part_9218259_733540587.1738838161137
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: quoted-printable
> I have lots of pairs of pre-aligned sequences (imported from an external=
MSA file),
In which format is your MSA file?
-Michiel
On Thursday, January 30, 2025 at 11:59:33 PM GMT+9, Peter Cock <p.j.a.c=
[email protected]> wrote: =20
=20
Hello all, and Michiel in particular,
I am wondering if any of the pairwise alignment code in Bio.Align (much of =
which is written in C for speed) could help with this use case?:
I have lots of pairs of pre-aligned sequences (imported from an external MS=
A file), for which I am doing something like this:
```python
def count_matches_etc(query_seq, subject_seq):
=C2=A0 =C2=A0 assert len(query_seq) =3D=3D len(subject_seq), "Should be sam=
e length"
=C2=A0 =C2=A0 matches =3D non_gap_mismatches =3D either_gapped =3D both_gap=
ped =3D 0
=C2=A0 =C2=A0 for q, s in zip(query_seq, subject_seq, strict=3DTrue):
=C2=A0 =C2=A0 =C2=A0 =C2=A0 if q =3D=3D "-" and s =3D=3D "-":
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 both_gapped +=3D 1
=C2=A0 =C2=A0 =C2=A0 =C2=A0 elif q =3D=3D "-" or s =3D=3D "-":
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 either_gapped +=3D 1
=C2=A0 =C2=A0 =C2=A0 =C2=A0 elif q =3D=3D s:
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 matches +=3D 1
=C2=A0 =C2=A0 =C2=A0 =C2=A0 else:
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 non_gap_mismatches +=3D 1
=C2=A0 =C2=A0 assert matches + non_gap_mismatches + either_gapped + both_ga=
pped =3D=3D len(query_seq)
=C2=A0 =C2=A0 return matches, non_gap_mismatches, either_gapped, both_gappe=
d
# Test case
assert (9, 1, 2, 1) =3D=3D count_matches_etc("ACGTAC-TAC-GT", "AGGT-CGTAC-G=
T")
```
Sticking with Python that could be optimized (e.g. I am currently using thi=
s with sequences of a million base pairs but few gaps), however I have writ=
ten this example with clarity foremost in mind.
Thank you,
Peter_______________________________________________
Biopython mailing list=C2=A0 -=C2=A0 [email protected]
https://mailman.open-bio.org/mailman/listinfo/biopython
=20
------=_Part_9218259_733540587.1738838161137
Content-Type: text/html; charset=UTF-8
Content-Transfer-Encoding: quoted-printable
<html><head></head><body><div class=3D"ydpcc8aacfayahoo-style-wrap" style=
=3D"font-family:Helvetica Neue, Helvetica, Arial, sans-serif;font-size:10px=
;"><div></div>
<div dir=3D"ltr" data-setdir=3D"false">> <span>I have lots of pa=
irs of pre-aligned sequences (imported from an external MSA file),</span></=
div><div dir=3D"ltr" data-setdir=3D"false"><span><br></span></div><div dir=
=3D"ltr" data-setdir=3D"false"><span>In which format is your MSA file?</spa=
n></div><div dir=3D"ltr" data-setdir=3D"false"><span><br></span></div><div =
dir=3D"ltr" data-setdir=3D"false"><span>-Michiel<br></span></div><div><br><=
/div>
=20
</div><div id=3D"yahoo_quoted_9753999089" class=3D"yahoo_quoted">
<div style=3D"font-family:'Helvetica Neue', Helvetica, Arial, s=
ans-serif;font-size:13px;color:#26282a;">
=20
<div>
On Thursday, January 30, 2025 at 11:59:33 PM GMT+9,=
Peter Cock <[email protected]> wrote:
</div>
<div><br></div>
<div><br></div>
=20
=20
<div><div id=3D"yiv8603226725"><div dir=3D"ltr">Hello all, =
and Michiel in particular,<br><br>I am wondering if any of the pairwise ali=
gnment code in Bio.Align (much of which is written in C for speed) could he=
lp with this use case?:<br><br>I have lots of pairs of pre-aligned sequence=
s (imported from an external MSA file), for which I am doing something like=
this:<br><br>```python<br>def count_matches_etc(query_seq, subject_seq):<b=
r> assert len(query_seq) =3D=3D len(subject_seq), "Should be s=
ame length"<br> matches =3D non_gap_mismatches =3D either_gapp=
ed =3D both_gapped =3D 0<br> for q, s in zip(query_seq, subjec=
t_seq, strict=3DTrue):<br> if q =3D=3D "-" and s=
=3D=3D "-":<br> both_gapped +=3D =
1<br> elif q =3D=3D "-" or s =3D=3D "-":<br>&nbs=
p; either_gapped +=3D 1<br> =
elif q =3D=3D s:<br>  =
; matches +=3D 1<br> else:<br> &nbs=
p; non_gap_mismatches +=3D 1<br> assert m=
atches + non_gap_mismatches + either_gapped + both_gapped =3D=3D len(query_=
seq)<br> return matches, non_gap_mismatches, either_gapped, bo=
th_gapped<br><br><br># Test case<br>assert (9, 1, 2, 1) =3D=3D count_matche=
s_etc("ACGTAC-TAC-GT", "AGGT-CGTAC-GT")<br>```<br><br>Sticking with Python =
that could be optimized (e.g. I am currently using this with sequences of a=
million base pairs but few gaps), however I have written this example with=
clarity foremost in mind.<br><div><br></div><div>Thank you,</div><div><br>=
</div>Peter</div>
</div>_______________________________________________<br>Biopython mailing =
list - <a ymailto=3D"mailto:[email protected]" href=3D"ma=
ilto:[email protected]">[email protected]</a><br><a href=3D"htt=
ps://mailman.open-bio.org/mailman/listinfo/biopython" target=3D"_blank">htt=
ps://mailman.open-bio.org/mailman/listinfo/biopython</a><br></div>
</div>
</div></body></html>
------=_Part_9218259_733540587.1738838161137--
--===============5283190597617378020==
Content-Type: text/plain; charset="us-ascii"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Content-Disposition: inline
_______________________________________________
Biopython mailing list - [email protected]
https://mailman.open-bio.org/mailman/listinfo/biopython
--===============5283190597617378020==--