Re: [Biopython] Comparing two already aligned sequences quickly to count gaps/matches/mismatches
Michiel de Hoon <[email protected]> Thu, 6 Feb 2025 23:30:42 +0000 (UTC)
| Newsgroups | gmane.comp.python.bio.general |
|---|---|
| Message-ID | <[email protected]> |
--===============6029306226636094164==
Content-Type: multipart/alternative;
boundary="----=_Part_9558251_1314731363.1738884642904"
------=_Part_9558251_1314731363.1738884642904
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: quoted-printable
It may be faster to create a numpy array with the sequence (as bytes), and=
then check for equality. Then the loop is done in C.
Best,-Michiel
On Thursday, February 6, 2025 at 09:56:50 PM GMT+9, Peter Cock <p.j.a.c=
[email protected]> wrote: =20
=20
Sadly I don't think that answers my query.
I don't actually need the number of gap-pairs (which can be inferred from t=
he other numbers anyway as they sum to the number of columns in the MSA). S=
o that's fine.
However, I am specifically needing the number of matches and mismatches. Th=
at means I can't use the faster mode you're suggesting (which ignores the n=
on-gap characters in the sequences) :(
Moreover, I don't want the global counts, I want the pairwise counts. I rat=
her suspect based on the Biopython 1.85 code that I do just have to do the =
inner loop over pairs of bases:
https://github.com/biopython/biopython/blob/biopython-185/Bio/Align/__init_=
_.py#L3576
Thanks though - it looks like I wasn't overlooking something "off the shelf=
" here.
Peter
On Thu, Feb 6, 2025 at 12:29=E2=80=AFPM Michiel de Hoon <[email protected]=
m> wrote:
> This is with a FASTA input MSA (all sequences the same length with gap c=
haracters, and some pairs will have common gaps).
Then you can use the `"fasta"` parser in `Bio.Align`, which parses the alig=
ned sequences in C, so it should be fast.Then call the new `.counts`=C2=A0 =
method in Biopython 1.86dev with ignore_sequences=3DTrue to calculate the n=
umber of insertions and deletions quickly.This won't give you the number of=
gap-against-gap alignments, but those are not meaningful anyway.
-Michiel
On Thursday, February 6, 2025 at 07:51:48 PM GMT+9, Peter Cock <p.j.a.c=
[email protected]> wrote: =20
=20
This is with a FASTA input MSA (all sequences the same length with gap cha=
racters, 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 <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
=20
=20
------=_Part_9558251_1314731363.1738884642904
Content-Type: text/html; charset=UTF-8
Content-Transfer-Encoding: quoted-printable
<html><head></head><body><div class=3D"ydp3abe364yahoo-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">It may be faster to create a=
numpy array with the sequence (as bytes), and then check for equality. The=
n the loop is done in C.</div><div dir=3D"ltr" data-setdir=3D"false"><br></=
div><div dir=3D"ltr" data-setdir=3D"false">Best,</div><div dir=3D"ltr" data=
-setdir=3D"false">-Michiel</div><div><br></div>
=20
</div><div id=3D"yahoo_quoted_9269853774" 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, February 6, 2025 at 09:56:50 PM GMT+9,=
Peter Cock <[email protected]> wrote:
</div>
<div><br></div>
<div><br></div>
=20
=20
<div><div id=3D"yiv4997287671"><div><div dir=3D"ltr"><div>S=
adly I don't think that answers my query.</div><div><br clear=3D"none"></di=
v><div>I don't actually need the number of gap-pairs (which can be inferred=
from the other numbers anyway as they sum to the number of columns in the =
MSA). So that's fine.</div><div><br clear=3D"none"></div><div>However, I am=
specifically needing the number of matches and mismatches. That means I ca=
n't use the faster mode you're suggesting (which ignores the non-gap charac=
ters in the sequences) :(</div><div><br clear=3D"none"></div><div>Moreover,=
I don't want the global counts, I want the pairwise counts. I rather suspe=
ct based on the Biopython 1.85 code that I do just have to do the inner loo=
p over pairs of bases:</div><div><br clear=3D"none"></div><div><a rel=3D"no=
follow noopener noreferrer" shape=3D"rect" target=3D"_blank" href=3D"https:=
//github.com/biopython/biopython/blob/biopython-185/Bio/Align/__init__.py#L=
3576">https://github.com/biopython/biopython/blob/biopython-185/Bio/Align/_=
_init__.py#L3576</a></div><div><br clear=3D"none"></div><div>Thanks though =
- it looks like I wasn't overlooking something "off the shelf" here.</div><=
div><br clear=3D"none"></div><div>Peter</div></div><br clear=3D"none"><div =
id=3D"yiv4997287671yqt21426" class=3D"yiv4997287671yqt6878398598"><div clas=
s=3D"yiv4997287671gmail_quote yiv4997287671gmail_quote_container"><div dir=
=3D"ltr" class=3D"yiv4997287671gmail_attr">On Thu, Feb 6, 2025 at 12:29=E2=
=80=AFPM Michiel de Hoon <<a rel=3D"nofollow noopener noreferrer" shape=
=3D"rect" ymailto=3D"mailto:[email protected]" target=3D"_blank" href=3D"=
mailto:[email protected]">[email protected]</a>> wrote:<br clear=3D"=
none"></div><blockquote style=3D"margin:0px 0px 0px 0.8ex;border-left:1px s=
olid rgb(204,204,204);padding-left:1ex;" class=3D"yiv4997287671gmail_quote"=
><div><div style=3D"font-family:Helvetica Neue, Helvetica, Arial, sans-seri=
f;font-size:10px;"><div></div>
<div dir=3D"ltr">> <span></span><div>This is with a FASTA input =
MSA (all sequences the same length with gap characters, and some pairs will=
have common gaps).</div></div><div dir=3D"ltr"><br clear=3D"none"></div><d=
iv dir=3D"ltr">Then you can use the `"fasta"` parser in `Bio.Align`, which =
parses the aligned sequences in C, so it should be fast.</div><div dir=3D"l=
tr">Then call the new `.counts` method in Biopython 1.86dev with igno=
re_sequences=3DTrue to calculate the number of insertions and deletions qui=
ckly.</div><div dir=3D"ltr">This won't give you the number of gap-against-g=
ap alignments, but those are not meaningful anyway.</div><div dir=3D"ltr"><=
br clear=3D"none"></div><div dir=3D"ltr">-Michiel<br clear=3D"none"></div><=
div><br clear=3D"none"></div>
=20
</div><div id=3D"yiv4997287671m_7082227319688521718yahoo_quoted_889=
9767467">
<div style=3D"font-family:Helvetica, Arial, sans-serif;font-siz=
e:13px;color:rgb(38,40,42);">
=20
<div>
On Thursday, February 6, 2025 at 07:51:48 PM GMT+9,=
Peter Cock <<a rel=3D"nofollow noopener noreferrer" shape=3D"rect" ymai=
lto=3D"mailto:[email protected]" target=3D"_blank" href=3D"mailto:p=
[email protected]">[email protected]</a>> wrote:
</div>
<div><br clear=3D"none"></div>
<div><br clear=3D"none"></div>
=20
=20
<div><div id=3D"yiv4997287671m_7082227319688521718yiv786305=
0864"><div><div dir=3D"ltr"><div>This is with a FASTA input MSA (all sequen=
ces the same length with gap characters, and some pairs will have common ga=
ps).</div><div><br clear=3D"none"></div><div>I have been loading this incre=
mentally 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 wer=
e loaded).</div><div><br clear=3D"none"></div><div>The test case is only ab=
out 100MB on disk, 100 sequences each of 1 million base pairs, and takes ab=
out a minute or two (multi-threaded laptop).<br clear=3D"none"></div><div><=
br clear=3D"none"></div><div><div>I don't think this will be memory constra=
ined, so working on the entire MSA with a single thread is fine (if faster)=
.<br clear=3D"none"></div><div><br clear=3D"none"></div><div>Thank you,</di=
v><div><br clear=3D"none"></div><div>Peter<br clear=3D"none"></div></div></=
div><br clear=3D"none"><div id=3D"yiv4997287671m_7082227319688521718yiv7863=
050864yqt97629"><div><div dir=3D"ltr">On Thu, Feb 6, 2025 at 10:36=E2=80=AF=
AM Michiel de Hoon <<a rel=3D"nofollow noopener noreferrer" shape=3D"rec=
t" ymailto=3D"mailto:[email protected]" target=3D"_blank" href=3D"mailto:=
[email protected]">[email protected]</a>> wrote:<br clear=3D"none"><=
/div><blockquote style=3D"margin:0px 0px 0px 0.8ex;border-left:1px solid rg=
b(204,204,204);padding-left:1ex;"><div><div style=3D"font-family:Helvetica =
Neue, Helvetica, Arial, sans-serif;font-size:10px;"><div></div>
<div dir=3D"ltr">> <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 clear=3D"none"></span></div><div dir=3D"ltr"><span>In which format=
is your MSA file?</span></div><div dir=3D"ltr"><span><br clear=3D"none"></=
span></div><div dir=3D"ltr"><span>-Michiel<br clear=3D"none"></span></div><=
div><br clear=3D"none"></div>
=20
</div><div id=3D"yiv4997287671m_7082227319688521718yiv7863050864m_8=
475683859613579565yahoo_quoted_9753999089">
<div style=3D"font-family:Helvetica, Arial, sans-serif;font-siz=
e:13px;color:rgb(38,40,42);">
=20
<div>
On Thursday, January 30, 2025 at 11:59:33 PM GMT+9,=
Peter Cock <<a rel=3D"nofollow noopener noreferrer" shape=3D"rect" ymai=
lto=3D"mailto:[email protected]" target=3D"_blank" href=3D"mailto:p=
[email protected]">[email protected]</a>> wrote:
</div>
<div><br clear=3D"none"></div>
<div><br clear=3D"none"></div>
=20
=20
<div><div id=3D"yiv4997287671m_7082227319688521718yiv786305=
0864m_8475683859613579565yiv8603226725"><div dir=3D"ltr">Hello all, and Mic=
hiel in particular,<br clear=3D"none"><br clear=3D"none">I am wondering if =
any of the pairwise alignment code in Bio.Align (much of which is written i=
n C for speed) could help with this use case?:<br clear=3D"none"><br clear=
=3D"none">I have lots of pairs of pre-aligned sequences (imported from an e=
xternal MSA file), for which I am doing something like this:<br clear=3D"no=
ne"><br clear=3D"none">```python<br clear=3D"none">def count_matches_etc(qu=
ery_seq, subject_seq):<br clear=3D"none"> assert len(query_seq=
) =3D=3D len(subject_seq), "Should be same length"<br clear=3D"none"> =
matches =3D non_gap_mismatches =3D either_gapped =3D both_gapped =
=3D 0<br clear=3D"none"> for q, s in zip(query_seq, subject_se=
q, strict=3DTrue):<br clear=3D"none"> if q =3D=
=3D "-" and s =3D=3D "-":<br clear=3D"none"> &nb=
sp; both_gapped +=3D 1<br clear=3D"none"> =
elif q =3D=3D "-" or s =3D=3D "-":<br clear=3D"none"> =
either_gapped +=3D 1<br clear=3D"none"> &=
nbsp; elif q =3D=3D s:<br clear=3D"none"> =
matches +=3D 1<br clear=3D"none"> =
else:<br clear=3D"none"> non_gap_=
mismatches +=3D 1<br clear=3D"none"> assert matches + non_gap_=
mismatches + either_gapped + both_gapped =3D=3D len(query_seq)<br clear=3D"=
none"> return matches, non_gap_mismatches, either_gapped, both=
_gapped<br clear=3D"none"><br clear=3D"none"><br clear=3D"none"># Test case=
<br clear=3D"none">assert (9, 1, 2, 1) =3D=3D count_matches_etc("ACGTAC-TAC=
-GT", "AGGT-CGTAC-GT")<br clear=3D"none">```<br clear=3D"none"><br clear=3D=
"none">Sticking with Python that could be optimized (e.g. I am currently us=
ing this with sequences of a million base pairs but few gaps), however I ha=
ve written this example with clarity foremost in mind.<br clear=3D"none"><d=
iv><br clear=3D"none"></div><div>Thank you,</div><div><br clear=3D"none"></=
div>Peter</div>
</div>_______________________________________________<br clear=3D"none">Bio=
python mailing list - <a rel=3D"nofollow noopener noreferrer" s=
hape=3D"rect" ymailto=3D"mailto:[email protected]" target=3D"_blank" =
href=3D"mailto:[email protected]">[email protected]</a><br clea=
r=3D"none"><a rel=3D"nofollow noopener noreferrer" shape=3D"rect" target=3D=
"_blank" href=3D"https://mailman.open-bio.org/mailman/listinfo/biopython">h=
ttps://mailman.open-bio.org/mailman/listinfo/biopython</a><br clear=3D"none=
"></div>
</div>
</div></div></blockquote></div></div>
</div></div></div>
</div>
</div></div></blockquote></div></div>
</div></div></div>
</div>
</div></body></html>
------=_Part_9558251_1314731363.1738884642904--
--===============6029306226636094164==
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
--===============6029306226636094164==--