Re: [Biopython] Struggling with MSA...
Dan Bolser <[email protected]> Thu, 15 Feb 2024 16:18:49 +0000
| Newsgroups | gmane.comp.python.bio.general |
|---|---|
| Message-ID | <CANs1yPLBAVMt_-7RkW8dKkzSrLZaNBLm9kWA7q4bqZnc4ZUpNA@mail.gmail.com> |
--===============5800428202233143121== Content-Type: multipart/alternative; boundary="000000000000474af006116dfaa6" --000000000000474af006116dfaa6 Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: quoted-printable 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=3Dhm= m > > 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-Start > 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.arvvgvDlskemlklakekaseeglkvefvvadaeklpfednsfDlvvssevlh= hv...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-QLRWETMDVRKLDFPSASFDVVLEKGTLD= ALlagErdpwtvssegvhTVDQVLSEVSRVLVPGGRFI > > > > > > 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]= > > wrote: > >> >>> > >> >>> Hello Dan, > >> >>> > >> >>> It that is a probability at the end, you have the wrong denominato= r > - > >> >>> 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 wrong > somewhere? > >> >> > >> >> > >> >> Cheers, > >> >> Dan. > >> >> > >> >> > >> >>> Peter > >> >>> > >> >>> On Thu, Feb 15, 2024 at 3:09=E2=80=AFPM Dan Bolser <dan.bolser@out= see.co.uk> > 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 the > 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 (whe= re > it matches the model). > >> >>> > > >> >>> > Many thanks for any suggestions, and sorry for not being able to > find the right document to answer these questions. > >> >>> > > >> >>> > > >> >>> > kthxbi, > >> >>> > Dan. > >> >>> > _______________________________________________ > >> >>> > Biopython mailing list - [email protected] > >> >>> > https://mailman.open-bio.org/mailman/listinfo/biopython > --000000000000474af006116dfaa6 Content-Type: text/html; charset="UTF-8" Content-Transfer-Encoding: quoted-printable <div dir=3D"ltr">Recreating the profile was just a sanity check. Really I w= ant to look at conservation across the protein of interest.</div><br><div c= lass=3D"gmail_quote"><div dir=3D"ltr" class=3D"gmail_attr">On Thu, 15 Feb 2= 024 at 16:13, Peter Cock <<a href=3D"mailto:[email protected]">p= [email protected]</a>> wrote:<br></div><blockquote class=3D"gmail= _quote" style=3D"margin:0px 0px 0px 0.8ex;border-left:1px solid 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> --000000000000474af006116dfaa6-- --===============5800428202233143121== 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 --===============5800428202233143121==--