Re: [Biopython] Comparing two already aligned sequences quickly to count gaps/matches/mismatches

--===============5038451685463138993==
Content-Type: multipart/alternative; boundary="00000000000044e756062d7707b8"

--00000000000044e756062d7707b8
Content-Type: text/plain; charset="UTF-8"
Content-Transfer-Encoding: quoted-printable

This is with a FASTA input MSA (all sequences the same length with gap
characters, and some pairs will have common gaps).

I have been loading this incrementally so only two sequences were in RAM at
any point - and then switched to multiple threads (so at any point with K
threads only 2*K sequences were loaded).

The test case is only about 100MB on disk, 100 sequences each of 1 million
base pairs, and takes about a minute or two (multi-threaded laptop).

I don't think this will be memory constrained, so working on the entire MSA
with a single thread is fine (if faster).

Thank you,

Peter

On Thu, Feb 6, 2025 at 10:36=E2=80=AFAM Michiel de Hoon <[email protected]=
m> wrote:

> > 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 <
> [email protected]> wrote:
>
>
> Hello all, and Michiel in particular,
>
> I am wondering if any of the pairwise alignment code in Bio.Align (much o=
f
> 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
> MSA file), for which I am doing something like this:
>
> ```python
> def count_matches_etc(query_seq, subject_seq):
>     assert len(query_seq) =3D=3D len(subject_seq), "Should be same length=
"
>     matches =3D non_gap_mismatches =3D either_gapped =3D both_gapped =3D =
0
>     for q, s in zip(query_seq, subject_seq, strict=3DTrue):
>         if q =3D=3D "-" and s =3D=3D "-":
>             both_gapped +=3D 1
>         elif q =3D=3D "-" or s =3D=3D "-":
>             either_gapped +=3D 1
>         elif q =3D=3D s:
>             matches +=3D 1
>         else:
>             non_gap_mismatches +=3D 1
>     assert matches + non_gap_mismatches + either_gapped + both_gapped =3D=
=3D
> len(query_seq)
>     return matches, non_gap_mismatches, either_gapped, both_gapped
>
>
> # Test case
> assert (9, 1, 2, 1) =3D=3D count_matches_etc("ACGTAC-TAC-GT", "AGGT-CGTAC=
-GT")
> ```
>
> 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.
>
> Thank you,
>
> Peter
> _______________________________________________
> Biopython mailing list  -  [email protected]
> https://mailman.open-bio.org/mailman/listinfo/biopython
>

--00000000000044e756062d7707b8
Content-Type: text/html; charset="UTF-8"
Content-Transfer-Encoding: quoted-printable

<div dir=3D"ltr"><div>This is with a FASTA input MSA (all sequences the sam=
e length with gap characters, and some pairs will have common gaps).</div><=
div><br></div><div>I have been loading this incrementally so only two seque=
nces were in RAM at any point - and then switched to multiple threads (so a=
t any point with K threads only 2*K sequences were loaded).</div><div><br><=
/div><div>The test case is only about 100MB on disk, 100 sequences each of =
1 million base pairs, and takes about a minute or two (multi-threaded lapto=
p).<br></div><div><br></div><div><div>I don&#39;t think this will be memory=
 constrained, so working on the entire MSA with a single thread is fine (if=
 faster).<br></div><div><br></div><div>Thank you,</div><div><br></div><div>=
Peter<br></div></div></div><br><div class=3D"gmail_quote gmail_quote_contai=
ner"><div dir=3D"ltr" class=3D"gmail_attr">On Thu, Feb 6, 2025 at 10:36=E2=
=80=AFAM Michiel de Hoon &lt;<a href=3D"mailto:[email protected]">mjldeho=
[email protected]</a>&gt; wrote:<br></div><blockquote class=3D"gmail_quote" styl=
e=3D"margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);paddin=
g-left:1ex"><div><div style=3D"font-family:Helvetica Neue,Helvetica,Arial,s=
ans-serif;font-size:10px"><div></div>
        <div dir=3D"ltr">&gt; <span>I have lots of pairs of pre-aligned seq=
uences (imported from an external MSA file),</span></div><div dir=3D"ltr"><=
span><br></span></div><div dir=3D"ltr"><span>In which format is your MSA fi=
le?</span></div><div dir=3D"ltr"><span><br></span></div><div dir=3D"ltr"><s=
pan>-Michiel<br></span></div><div><br></div>
       =20
        </div><div id=3D"m_8475683859613579565yahoo_quoted_9753999089">
            <div style=3D"font-family:&quot;Helvetica Neue&quot;,Helvetica,=
Arial,sans-serif;font-size:13px;color:rgb(38,40,42)">
               =20
                <div>
                        On Thursday, January 30, 2025 at 11:59:33 PM GMT+9,=
 Peter Cock &lt;<a href=3D"mailto:[email protected]" target=3D"_bla=
nk">[email protected]</a>&gt; wrote:
                    </div>
                    <div><br></div>
                    <div><br></div>
               =20
               =20
                <div><div id=3D"m_8475683859613579565yiv8603226725"><div di=
r=3D"ltr">Hello all, and Michiel in particular,<br><br>I am wondering if an=
y of the pairwise alignment code in Bio.Align (much of which is written in =
C for speed) could help with this use case?:<br><br>I have lots of pairs of=
 pre-aligned sequences (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):<br>=C2=A0 =C2=A0 assert len(query_seq) =3D=3D len(subje=
ct_seq), &quot;Should be same length&quot;<br>=C2=A0 =C2=A0 matches =3D non=
_gap_mismatches =3D either_gapped =3D both_gapped =3D 0<br>=C2=A0 =C2=A0 fo=
r q, s in zip(query_seq, subject_seq, strict=3DTrue):<br>=C2=A0 =C2=A0 =C2=
=A0 =C2=A0 if q =3D=3D &quot;-&quot; and s =3D=3D &quot;-&quot;:<br>=C2=A0 =
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 both_gapped +=3D 1<br>=C2=A0 =C2=A0 =C2=
=A0 =C2=A0 elif q =3D=3D &quot;-&quot; or s =3D=3D &quot;-&quot;:<br>=C2=A0=
 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 either_gapped +=3D 1<br>=C2=A0 =C2=A0 =
=C2=A0 =C2=A0 elif q =3D=3D s:<br>=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0=
 matches +=3D 1<br>=C2=A0 =C2=A0 =C2=A0 =C2=A0 else:<br>=C2=A0 =C2=A0 =C2=
=A0 =C2=A0 =C2=A0 =C2=A0 non_gap_mismatches +=3D 1<br>=C2=A0 =C2=A0 assert =
matches + non_gap_mismatches + either_gapped + both_gapped =3D=3D len(query=
_seq)<br>=C2=A0 =C2=A0 return matches, non_gap_mismatches, either_gapped, b=
oth_gapped<br><br><br># Test case<br>assert (9, 1, 2, 1) =3D=3D count_match=
es_etc(&quot;ACGTAC-TAC-GT&quot;, &quot;AGGT-CGTAC-GT&quot;)<br>```<br><br>=
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.<br><div><br></div><div>Than=
k you,</div><div><br></div>Peter</div>
</div>_______________________________________________<br>Biopython mailing =
list=C2=A0 -=C2=A0 <a href=3D"mailto:[email protected]" target=3D"_bl=
ank">[email protected]</a><br><a href=3D"https://mailman.open-bio.org=
/mailman/listinfo/biopython" target=3D"_blank">https://mailman.open-bio.org=
/mailman/listinfo/biopython</a><br></div>
            </div>
        </div></div></blockquote></div>

--00000000000044e756062d7707b8--

--===============5038451685463138993==
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

--===============5038451685463138993==--