Re: [Biopython] Comparing two already aligned sequences quickly to count gaps/matches/mismatches
Peter Cock <[email protected]> Thu, 6 Feb 2025 12:56:36 +0000
| Newsgroups | gmane.comp.python.bio.general |
|---|---|
| Message-ID | <CAKVJ-_5eOmsoYJWN0eEDSrcC-qci5h1DnUJeCq+rgCFr6DbMkA@mail.gmail.com> |
--===============4688907330968466005== Content-Type: multipart/alternative; boundary="0000000000006ff1d9062d78c6d9" --0000000000006ff1d9062d78c6d9 Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: quoted-printable Sadly I don't think that answers my query. 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. However, I am specifically needing the number of matches and mismatches. That means I can't use the faster mode you're suggesting (which ignores the non-gap characters in the sequences) :( Moreover, I don't want the global counts, I want the pairwise counts. I rather 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 > characters, and some pairs will have common gaps). > > Then you can use the `"fasta"` parser in `Bio.Align`, which parses the > aligned sequences in C, so it should be fast. > Then call the new `.counts` method in Biopython 1.86dev with > ignore_sequences=3DTrue to calculate the number of insertions and deletio= ns > 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 < > [email protected]> wrote: > > > 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 wit= h > K threads only 2*K sequences were loaded). > > The test case is only about 100MB on disk, 100 sequences each of 1 millio= n > 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 <mjldehoon@yahoo.= com> > 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 > > --0000000000006ff1d9062d78c6d9 Content-Type: text/html; charset="UTF-8" Content-Transfer-Encoding: quoted-printable <div dir=3D"ltr"><div>Sadly I don't think that answers my query.</div><= div><br></div><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></div><div>However,= I am specifically needing the number of matches and mismatches. That means= I can't use the faster mode you're suggesting (which ignores the n= on-gap characters in the sequences) :(</div><div><br></div><div>Moreover, I= don't want the global counts, I want the pairwise counts. I rather sus= pect based on the Biopython 1.85 code that I do just have to do the inner l= oop over pairs of bases:</div><div><br></div><div><a href=3D"https://github= .com/biopython/biopython/blob/biopython-185/Bio/Align/__init__.py#L3576">ht= tps://github.com/biopython/biopython/blob/biopython-185/Bio/Align/__init__.= py#L3576</a></div><div><br></div><div>Thanks though - it looks like I wasn&= #39;t overlooking something "off the shelf" here.</div><div><br><= /div><div>Peter</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 12:29=E2= =80=AFPM Michiel de Hoon <<a href=3D"mailto:[email protected]">mjldeho= [email protected]</a>> 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">> <span><div>This is with a FASTA input MSA (al= l sequences the same length with gap characters, and some pairs will have c= ommon gaps).</div></span></div><div dir=3D"ltr"><br></div><div dir=3D"ltr">= Then you can use the `"fasta"` parser in `Bio.Align`, which parse= s the aligned sequences in C, so it should be fast.</div><div dir=3D"ltr">T= hen call the new `.counts`=C2=A0 method in Biopython 1.86dev with ignore_se= quences=3DTrue to calculate the number of insertions and deletions quickly.= </div><div dir=3D"ltr">This won't give you the number of gap-against-ga= p alignments, but those are not meaningful anyway.</div><div dir=3D"ltr"><b= r></div><div dir=3D"ltr">-Michiel<br></div><div><br></div> =20 </div><div id=3D"m_7082227319688521718yahoo_quoted_8899767467"> <div style=3D"font-family:"Helvetica Neue",Helvetica,= Arial,sans-serif;font-size:13px;color:rgb(38,40,42)"> =20 <div> On Thursday, February 6, 2025 at 07:51:48 PM GMT+9,= Peter Cock <<a href=3D"mailto:[email protected]" target=3D"_bla= nk">[email protected]</a>> wrote: </div> <div><br></div> <div><br></div> =20 =20 <div><div id=3D"m_7082227319688521718yiv7863050864"><div><d= iv dir=3D"ltr"><div>This is with a FASTA input MSA (all sequences the same = length with gap characters, and some pairs will have common gaps).</div><di= v><br clear=3D"none"></div><div>I have been loading this incrementally so o= nly 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).</d= iv><div><br clear=3D"none"></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 laptop).<br clear=3D"none"></div><div><br clear=3D"n= one"></div><div><div>I don't think this will be memory constrained, so = working on the entire MSA with a single thread is fine (if faster).<br clea= r=3D"none"></div><div><br clear=3D"none"></div><div>Thank you,</div><div><b= r clear=3D"none"></div><div>Peter<br clear=3D"none"></div></div></div><br c= lear=3D"none"><div id=3D"m_7082227319688521718yiv7863050864yqt97629"><div><= div dir=3D"ltr">On Thu, Feb 6, 2025 at 10:36=E2=80=AFAM Michiel de Hoon <= ;<a rel=3D"nofollow noopener noreferrer" shape=3D"rect" href=3D"mailto:mjld= [email protected]" target=3D"_blank">[email protected]</a>> wrote:<br cl= ear=3D"none"></div><blockquote style=3D"margin:0px 0px 0px 0.8ex;border-lef= t:1px solid rgb(204,204,204);padding-left:1ex"><div><div style=3D"font-fami= ly: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"m_7082227319688521718yiv7863050864m_84756838596135= 79565yahoo_quoted_9753999089"> <div style=3D"font-family: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 <<a rel=3D"nofollow noopener noreferrer" shape=3D"rect" href= =3D"mailto:[email protected]" target=3D"_blank">p.j.a.cock@googlema= il.com</a>> wrote: </div> <div><br clear=3D"none"></div> <div><br clear=3D"none"></div> =20 =20 <div><div id=3D"m_7082227319688521718yiv7863050864m_8475683= 859613579565yiv8603226725"><div dir=3D"ltr">Hello all, and Michiel in parti= cular,<br clear=3D"none"><br clear=3D"none">I am wondering if any of the pa= irwise alignment code in Bio.Align (much of which is written in C for speed= ) could help with this use case?:<br clear=3D"none"><br clear=3D"none">I ha= ve lots of pairs of pre-aligned sequences (imported from an external MSA fi= le), for which I am doing something like this:<br clear=3D"none"><br clear= =3D"none">```python<br clear=3D"none">def count_matches_etc(query_seq, subj= ect_seq):<br clear=3D"none">=C2=A0 =C2=A0 assert len(query_seq) =3D=3D len(= subject_seq), "Should be same length"<br clear=3D"none">=C2=A0 = =C2=A0 matches =3D non_gap_mismatches =3D either_gapped =3D both_gapped =3D= 0<br clear=3D"none">=C2=A0 =C2=A0 for q, s in zip(query_seq, subject_seq, = strict=3DTrue):<br clear=3D"none">=C2=A0 =C2=A0 =C2=A0 =C2=A0 if q =3D=3D &= quot;-" and s =3D=3D "-":<br clear=3D"none">=C2=A0 =C2=A0 = =C2=A0 =C2=A0 =C2=A0 =C2=A0 both_gapped +=3D 1<br clear=3D"none">=C2=A0 =C2= =A0 =C2=A0 =C2=A0 elif q =3D=3D "-" or s =3D=3D "-":<br= clear=3D"none">=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 either_gapped += =3D 1<br clear=3D"none">=C2=A0 =C2=A0 =C2=A0 =C2=A0 elif q =3D=3D s:<br cle= ar=3D"none">=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 matches +=3D 1<br cle= ar=3D"none">=C2=A0 =C2=A0 =C2=A0 =C2=A0 else:<br clear=3D"none">=C2=A0 =C2= =A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 non_gap_mismatches +=3D 1<br clear=3D"none"= >=C2=A0 =C2=A0 assert matches + non_gap_mismatches + either_gapped + both_g= apped =3D=3D len(query_seq)<br clear=3D"none">=C2=A0 =C2=A0 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 using this with s= equences of a million base pairs but few gaps), however I have written this= example with clarity foremost in mind.<br clear=3D"none"><div><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=C2=A0 -=C2=A0 <a rel=3D"nofollow noopener noreferrer" s= hape=3D"rect" href=3D"mailto:[email protected]" target=3D"_blank">Bio= [email protected]</a><br clear=3D"none"><a rel=3D"nofollow noopener nore= ferrer" shape=3D"rect" href=3D"https://mailman.open-bio.org/mailman/listinf= o/biopython" target=3D"_blank">https://mailman.open-bio.org/mailman/listinf= o/biopython</a><br clear=3D"none"></div> </div> </div></div></blockquote></div></div> </div></div></div> </div> </div></div></blockquote></div> --0000000000006ff1d9062d78c6d9-- --===============4688907330968466005== 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 --===============4688907330968466005==--