Re: [Biopython] Struggling with MSA...
Dan Bolser <[email protected]> Fri, 16 Feb 2024 10:59:32 +0000
| Newsgroups | gmane.comp.python.bio.general |
|---|---|
| Message-ID | <CANs1yPLHoL27Jp3+WFza=X3wqN+sASabpAFMnCRMHtwL2QA5mA@mail.gmail.com> |
--===============3842091600841923411== Content-Type: multipart/alternative; boundary="000000000000234ef306117da295" --000000000000234ef306117da295 Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: quoted-printable I don't suppose there is a package for reading in the HMM? =F0=9F=98=85 On Thu, 15 Feb 2024 at 16:18, Dan Bolser <[email protected]> wrote: > Recreating the profile was just a sanity check. Really I want to look at > conservation across the protein of interest. > > On Thu, 15 Feb 2024 at 16:13, Peter Cock <[email protected]> > wrote: > >> You may have more joy starting from the HMM for the model? >> >> https://www.ebi.ac.uk/interpro/wwwapi//entry/pfam/PF08241?annotation=3Dh= mm >> >> Peter >> >> On Thu, Feb 15, 2024 at 4:06=E2=80=AFPM Dan Bolser <[email protected]= .uk> >> wrote: >> > >> > Searching HMMER gives me (cryptic?) details of the model match: >> > >> > Family-Id Family-Accession Clan Start End Ali-Start Ali-End Model-Star= t >> Model-End Bit-Score Ind.-E-value Cond.-E-value Description >> > >> =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D >> > >> > Methyltransf_11 PF08241.15 CL0063 63 173 63 172 1 95 51.30 1.5e-13 >> 3.1e-17 Methyltransferase domain >> > >> > >> ------------------------------------------------------------------------= --------------------------------------------------------- >> > >> > MODEL >> LdvGcGtGrlaealakrg.arvvgvDlskemlklakekaseeglkvefvvadaeklpfednsfDlvvssevl= hhv...e............dpekalkeiaRvLkpgGllv >> > MATCH L +GcG+ +l+ +l g +v+ vD+s ++++ +++ + ++++ +d++kl >> f+++sfD+v+ + +l+ + e ++++l+e+ RvL pgG+++ >> > PPL >> 569************7779*************77766666666.69**************************= ***8662667778888888******************97 >> > SEQ >> LVLGCGNSALSYELFLGGfPNVTSVDYSSVVVAAMQARHAHVP-QLRWETMDVRKLDFPSASFDVVLEKGTL= DALlagErdpwtvssegvhTVDQVLSEVSRVLVPGGRFI >> > >> > >> > So I guess it's possible to work backwards from there... >> > >> > >> > >> > On Thu, 15 Feb 2024 at 16:02, Peter Cock <[email protected]> >> wrote: >> >> >> >> Sorry yes, I used this as a threshold for if a column was gappy, and >> >> likely to be in the model or not: >> >> >> >> amino_acids.get("-", 0) < len(align)*0.5 >> >> >> >> Might need to use <=3D to match, or perhaps they use a more >> >> sophisticated criteria. >> >> >> >> Peter >> >> >> >> On Thu, Feb 15, 2024 at 3:56=E2=80=AFPM Dan Bolser <dan.bolser@outsee= .co.uk> >> wrote: >> >> > >> >> > Sorry, by 'half gaps' I thought you meant a specific thing, like, >> "look at that half gap!". I guess you mean, columns where half of the >> values in the column are gaps. >> >> > >> >> > Looking directly at the alignment, I guess there is no way to work >> out precisely which columns are gaps or not in the model? >> >> > >> >> > >> >> > >> >> > >> >> > >> >> > >> >> > On Thu, 15 Feb 2024 at 15:47, Dan Bolser <[email protected]> >> wrote: >> >> >> >> >> >> Hi Peter, >> >> >> >> >> >> On Thu, 15 Feb 2024 at 15:28, Peter Cock <[email protected]= m> >> wrote: >> >> >>> >> >> >>> Hello Dan, >> >> >>> >> >> >>> It that is a probability at the end, you have the wrong >> denominator - >> >> >>> should be len(align) =3D=3D sum(amino_acids.values() >> >> >> >> >> >> >> >> >> Makes no difference, they are the same (including -s). The profile >> reports -s as 'Insert Probability'. >> >> >> >> >> >> >> >> >>> Note this seed alignment has 164 columns, yet the model page says >> the >> >> >>> model has 95 columns. If I count the number of columns which are >> half >> >> >>> gaps that's pretty close... >> >> >> >> >> >> >> >> >> I don't understand what that means or implies... Is something wron= g >> somewhere? >> >> >> >> >> >> >> >> >> Cheers, >> >> >> Dan. >> >> >> >> >> >> >> >> >>> Peter >> >> >>> >> >> >>> On Thu, Feb 15, 2024 at 3:09=E2=80=AFPM Dan Bolser < >> [email protected]> wrote: >> >> >>> > >> >> >>> > Hi, >> >> >>> > >> >> >>> > Sorry, I can't follow the docs (or find the right docs). >> >> >>> > >> >> >>> > I've got the 'seed' stockholm alignment for this domain: >> >> >>> > >> https://www.ebi.ac.uk/interpro/entry/pfam/PF08241/entry_alignments/?type= =3Dseed >> >> >>> > >> >> >>> > and I'm trying to reproduce the signature it shows here: >> >> >>> > https://www.ebi.ac.uk/interpro/entry/pfam/PF08241/logo/ >> >> >>> > >> >> >>> > I'm not sure a) why the probabilities differ in the profile >> relative to the seed alignment, or b) how to filter columns in the >> alignment by those that have a match in the model (see columns 4-6 in th= e >> alignment, which are gaps in the model). >> >> >>> > >> >> >>> > I think if I can answer b) then the answer to a) will be, "look >> at the full alignment". >> >> >>> > >> >> >>> > Here is my crude 'best guess' code: >> >> >>> > >> >> >>> > import gzip >> >> >>> > import Bio.AlignIO >> >> >>> > >> >> >>> > # msa =3D "PF08241.alignment.full.gz" >> >> >>> > msa =3D "PF08241.alignment.seed.gz" >> >> >>> > >> >> >>> > with gzip.open(msa, "rt") as handle: >> >> >>> > align =3D Bio.AlignIO.read(handle, "stockholm") >> >> >>> > ncols =3D align.get_alignment_length() >> >> >>> > >> >> >>> > for col in range(ncols): >> >> >>> > amino_acids =3D dict() >> >> >>> > for s in align[:, col]: >> >> >>> > amino_acids[s] =3D amino_acids.get(s, 0) + 1 >> >> >>> > print(amino_acids) >> >> >>> > for s in amino_acids: >> >> >>> > print(f"{s}: {amino_acids[s]:3d} {amino_acids[s] / >> len(align):.3f}") >> >> >>> > >> >> >>> > >> >> >>> > >> >> >>> > I have the feeling I'm doin it rong... >> >> >>> > >> >> >>> > The above is just a 'warm up', really I want to see the >> conservation score, base by base on a given protein in the alignment (wh= ere >> it matches the model). >> >> >>> > >> >> >>> > Many thanks for any suggestions, and sorry for not being able t= o >> find the right document to answer these questions. >> >> >>> > >> >> >>> > >> >> >>> > kthxbi, >> >> >>> > Dan. >> >> >>> > _______________________________________________ >> >> >>> > Biopython mailing list - [email protected] >> >> >>> > https://mailman.open-bio.org/mailman/listinfo/biopython >> > --000000000000234ef306117da295 Content-Type: text/html; charset="UTF-8" Content-Transfer-Encoding: quoted-printable <div dir=3D"ltr">I don't suppose there is a package for reading in the = HMM?=C2=A0=F0=9F=98=85<div><br></div></div><br><div class=3D"gmail_quote"><= div dir=3D"ltr" class=3D"gmail_attr">On Thu, 15 Feb 2024 at 16:18, Dan Bols= er <<a href=3D"mailto:[email protected]">[email protected]</= a>> wrote:<br></div><blockquote class=3D"gmail_quote" style=3D"margin:0p= x 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><d= iv dir=3D"ltr">Recreating the profile was just a sanity check. Really I wan= t to look at conservation across the protein of interest.</div><br><div cla= ss=3D"gmail_quote"><div dir=3D"ltr" class=3D"gmail_attr">On Thu, 15 Feb 202= 4 at 16:13, Peter Cock <<a href=3D"mailto:[email protected]" tar= get=3D"_blank">[email protected]</a>> wrote:<br></div><blockquot= e class=3D"gmail_quote" style=3D"margin:0px 0px 0px 0.8ex;border-left:1px s= olid rgb(204,204,204);padding-left:1ex">You may have more joy starting from= the HMM for the model?<br> <br> <a href=3D"https://www.ebi.ac.uk/interpro/wwwapi//entry/pfam/PF08241?annota= tion=3Dhmm" rel=3D"noreferrer" target=3D"_blank">https://www.ebi.ac.uk/inte= rpro/wwwapi//entry/pfam/PF08241?annotation=3Dhmm</a><br> <br> Peter<br> <br> On Thu, Feb 15, 2024 at 4:06=E2=80=AFPM Dan Bolser <<a href=3D"mailto:da= [email protected]" target=3D"_blank">[email protected]</a>> wr= ote:<br> ><br> > Searching HMMER gives me (cryptic?) details of the model match:<br> ><br> > Family-Id Family-Accession Clan Start End Ali-Start Ali-End Model-Star= t Model-End Bit-Score Ind.-E-value Cond.-E-value Description<br> > =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D<br> ><br> > Methyltransf_11 PF08241.15 CL0063 63 173 63 172 1 95 51.30 1.5e-13 3.1= e-17 Methyltransferase domain<br> ><br> > ----------------------------------------------------------------------= -----------------------------------------------------------<br> ><br> > MODEL=C2=A0 LdvGcGtGrlaealakrg.arvvgvDlskemlklakekaseeglkvefvvadaeklpf= ednsfDlvvssevlhhv...e............dpekalkeiaRvLkpgGllv<br> > MATCH=C2=A0 L +GcG+ +l+ +l=C2=A0 =C2=A0g=C2=A0 +v+ vD+s=C2=A0 ++++ +++= =C2=A0 =C2=A0 + ++++=C2=A0 +d++kl f+++sfD+v+ + +l+ +=C2=A0 =C2=A0e=C2=A0 = =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0++++l+e+ RvL pgG+++<br> > PPL=C2=A0 =C2=A0 569************7779*************77766666666.69*******= **********************8662667778888888******************97<br> > SEQ=C2=A0 =C2=A0 LVLGCGNSALSYELFLGGfPNVTSVDYSSVVVAAMQARHAHVP-QLRWETMDV= RKLDFPSASFDVVLEKGTLDALlagErdpwtvssegvhTVDQVLSEVSRVLVPGGRFI<br> ><br> ><br> > So I guess it's possible to work backwards from there...<br> ><br> ><br> ><br> > On Thu, 15 Feb 2024 at 16:02, Peter Cock <<a href=3D"mailto:p.j.a.c= [email protected]" target=3D"_blank">[email protected]</a>> wro= te:<br> >><br> >> Sorry yes, I used this as a threshold for if a column was gappy, a= nd<br> >> likely to be in the model or not:<br> >><br> >> amino_acids.get("-", 0) < len(align)*0.5<br> >><br> >> Might need to use <=3D to match, or perhaps they use a more<br> >> sophisticated criteria.<br> >><br> >> Peter<br> >><br> >> On Thu, Feb 15, 2024 at 3:56=E2=80=AFPM Dan Bolser <<a href=3D"= mailto:[email protected]" target=3D"_blank">[email protected]</= a>> wrote:<br> >> ><br> >> > Sorry, by 'half gaps' I thought you meant a specific = thing, like, "look at that half gap!". I guess you mean, columns = where half of the values in the column are gaps.<br> >> ><br> >> > Looking directly at the alignment, I guess there is no way to= work out precisely which columns are gaps or not in the model?<br> >> ><br> >> ><br> >> ><br> >> ><br> >> ><br> >> ><br> >> > On Thu, 15 Feb 2024 at 15:47, Dan Bolser <<a href=3D"mailt= o:[email protected]" target=3D"_blank">[email protected]</a>>= ; wrote:<br> >> >><br> >> >> Hi Peter,<br> >> >><br> >> >> On Thu, 15 Feb 2024 at 15:28, Peter Cock <<a href=3D"m= ailto:[email protected]" target=3D"_blank">[email protected]= m</a>> wrote:<br> >> >>><br> >> >>> Hello Dan,<br> >> >>><br> >> >>> It that is a probability at the end, you have the wro= ng denominator -<br> >> >>> should be len(align) =3D=3D sum(amino_acids.values()<= br> >> >><br> >> >><br> >> >> Makes no difference, they are the same (including -s). Th= e profile reports -s as 'Insert Probability'.<br> >> >><br> >> >><br> >> >>> Note this seed alignment has 164 columns, yet the mod= el page says the<br> >> >>> model has 95 columns. If I count the number of column= s which are half<br> >> >>> gaps that's pretty close...<br> >> >><br> >> >><br> >> >> I don't understand what that means or implies... Is s= omething wrong somewhere?<br> >> >><br> >> >><br> >> >> Cheers,<br> >> >> Dan.<br> >> >><br> >> >><br> >> >>> Peter<br> >> >>><br> >> >>> On Thu, Feb 15, 2024 at 3:09=E2=80=AFPM Dan Bolser &l= t;<a href=3D"mailto:[email protected]" target=3D"_blank">dan.bolser@o= utsee.co.uk</a>> wrote:<br> >> >>> ><br> >> >>> > Hi,<br> >> >>> ><br> >> >>> > Sorry, I can't follow the docs (or find the = right docs).<br> >> >>> ><br> >> >>> > I've got the 'seed' stockholm alignm= ent for this domain:<br> >> >>> > <a href=3D"https://www.ebi.ac.uk/interpro/entry/= pfam/PF08241/entry_alignments/?type=3Dseed" rel=3D"noreferrer" target=3D"_b= lank">https://www.ebi.ac.uk/interpro/entry/pfam/PF08241/entry_alignments/?t= ype=3Dseed</a><br> >> >>> ><br> >> >>> > and I'm trying to reproduce the signature it= shows here:<br> >> >>> > <a href=3D"https://www.ebi.ac.uk/interpro/entry/= pfam/PF08241/logo/" rel=3D"noreferrer" target=3D"_blank">https://www.ebi.ac= .uk/interpro/entry/pfam/PF08241/logo/</a><br> >> >>> ><br> >> >>> > I'm not sure a) why the probabilities differ= in the profile relative to the seed alignment, or b) how to filter columns= in the alignment by those that have a match in the model (see columns 4-6 = in the alignment, which are gaps in the model).<br> >> >>> ><br> >> >>> > I think if I can answer b) then the answer to a)= will be, "look at the full alignment".<br> >> >>> ><br> >> >>> > Here is my crude 'best guess' code:<br> >> >>> ><br> >> >>> > import gzip<br> >> >>> > import Bio.AlignIO<br> >> >>> ><br> >> >>> > # msa =3D "PF08241.alignment.full.gz"<= br> >> >>> > msa =3D "PF08241.alignment.seed.gz"<br= > >> >>> ><br> >> >>> > with gzip.open(msa, "rt") as handle:<b= r> >> >>> > align =3D Bio.AlignIO.read(handle, "stockho= lm")<br> >> >>> > ncols =3D align.get_alignment_length()<br> >> >>> ><br> >> >>> > for col in range(ncols):<br> >> >>> > amino_acids =3D dict()<br> >> >>> > for s in align[:, col]:<br> >> >>> > amino_acids[s] =3D amino_acids.get(s, 0) + 1<br> >> >>> > print(amino_acids)<br> >> >>> > for s in amino_acids:<br> >> >>> > print(f"{s}: {amino_acids[s]:3d} {amino_aci= ds[s] / len(align):.3f}")<br> >> >>> ><br> >> >>> ><br> >> >>> ><br> >> >>> > I have the feeling I'm doin it rong...<br> >> >>> ><br> >> >>> > The above is just a 'warm up', really I = want to see the conservation score, base by base on a given protein in the = alignment (where it matches the model).<br> >> >>> ><br> >> >>> > Many thanks for any suggestions, and sorry for n= ot being able to find the right document to answer these questions.<br> >> >>> ><br> >> >>> ><br> >> >>> > kthxbi,<br> >> >>> > Dan.<br> >> >>> > _______________________________________________<= br> >> >>> > Biopython mailing list=C2=A0 -=C2=A0 <a href=3D"= mailto:[email protected]" target=3D"_blank">[email protected]</= a><br> >> >>> > <a href=3D"https://mailman.open-bio.org/mailman/= listinfo/biopython" rel=3D"noreferrer" target=3D"_blank">https://mailman.op= en-bio.org/mailman/listinfo/biopython</a><br> </blockquote></div> </blockquote></div> --000000000000234ef306117da295-- --===============3842091600841923411== 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 --===============3842091600841923411==--