Re: [Biopython] Struggling with MSA...
Dan Bolser <[email protected]> Thu, 15 Feb 2024 16:07:07 +0000
| Newsgroups | gmane.comp.python.bio.general |
|---|---|
| Message-ID | <CANs1yPJgoEZPEaFocWkm94pkMC_yXzO4C9xscK3t_4Ec1_0sgg@mail.gmail.com> |
--===============7035859266241881958== Content-Type: multipart/alternative; boundary="00000000000063b61506116dd028" --00000000000063b61506116dd028 Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: quoted-printable 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.arvvgvDlskemlklakekaseeglkvefvvadaeklpfednsfDlvvssevlhh= v...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-QLRWETMDVRKLDFPSASFDVVLEKGTLDA= LlagErdpwtvssegvhTVDQVLSEVSRVLVPGGRFI 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 <[email protected].= 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 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 wrong > somewhere? > >> > >> > >> Cheers, > >> Dan. > >> > >> > >>> Peter > >>> > >>> On Thu, Feb 15, 2024 at 3:09=E2=80=AFPM Dan Bolser <dan.bolser@outsee= .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 relativ= e > 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 conservatio= n > score, base by base on a given protein in the alignment (where 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 > --00000000000063b61506116dd028 Content-Type: text/html; charset="UTF-8" Content-Transfer-Encoding: quoted-printable <div dir=3D"ltr">Searching HMMER gives me (cryptic?) details of the model m= atch:<div><br></div><div>Family-Id Family-Accession Clan Start End Ali-Star= t Ali-End Model-Start Model-End Bit-Score Ind.-E-value Cond.-E-value Descri= ption<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.1e-17 Methyltransferase domain<br><br>-----------= ---------------------------------------------------------------------------= -------------------------------------------<br><br>MODEL =C2=A0LdvGcGtGrlae= alakrg.arvvgvDlskemlklakekaseeglkvefvvadaeklpfednsfDlvvssevlhhv...e........= ....dpekalkeiaRvLkpgGllv<br>MATCH =C2=A0L +GcG+ +l+ +l =C2=A0 g =C2=A0+v+ v= D+s =C2=A0++++ +++ =C2=A0 =C2=A0+ ++++ =C2=A0+d++kl f+++sfD+v+ + +l+ + =C2= =A0 e =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 ++++l+e+ RvL pgG+++<br>PPL = =C2=A0 =C2=A0569************7779*************77766666666.69****************= *************8662667778888888******************97<br>SEQ =C2=A0 =C2=A0LVLGC= GNSALSYELFLGGfPNVTSVDYSSVVVAAMQARHAHVP-QLRWETMDVRKLDFPSASFDVVLEKGTLDALlagEr= dpwtvssegvhTVDQVLSEVSRVLVPGGRFI<br></div><div><br></div><div><br></div><div= >So I guess it's possible to work backwards from there...</div><div><br= ></div><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:02, Peter Cock <<a href= =3D"mailto:[email protected]">[email protected]</a>> wro= te:<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">Sorry yes, I= used this as a threshold for if a column was gappy, and<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:da= [email protected]" target=3D"_blank">[email protected]</a>> wr= ote:<br> ><br> > Sorry, by 'half gaps' I thought you meant a specific thing, li= ke, "look at that half gap!". I guess you mean, columns where hal= f 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"mailto:dan.bol= [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"mailto:p.j= [email protected]" target=3D"_blank">[email protected]</a>>= wrote:<br> >>><br> >>> Hello Dan,<br> >>><br> >>> It that is a probability at the end, you have the wrong denomi= nator -<br> >>> should be len(align) =3D=3D sum(amino_acids.values()<br> >><br> >><br> >> Makes no difference, they are the same (including -s). The profile= reports -s as 'Insert Probability'.<br> >><br> >><br> >>> Note this seed alignment has 164 columns, yet the model page s= ays the<br> >>> model has 95 columns. If I count the number of columns which a= re half<br> >>> gaps that's pretty close...<br> >><br> >><br> >> I don't understand what that means or implies... Is something = 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 <<a href= =3D"mailto:[email protected]" target=3D"_blank">[email protected].= uk</a>> wrote:<br> >>> ><br> >>> > Hi,<br> >>> ><br> >>> > Sorry, I can't follow the docs (or find the right doc= s).<br> >>> ><br> >>> > I've got the 'seed' stockholm alignment for t= his domain:<br> >>> > <a href=3D"https://www.ebi.ac.uk/interpro/entry/pfam/PF08= 241/entry_alignments/?type=3Dseed" rel=3D"noreferrer" target=3D"_blank">htt= ps://www.ebi.ac.uk/interpro/entry/pfam/PF08241/entry_alignments/?type=3Dsee= d</a><br> >>> ><br> >>> > and I'm trying to reproduce the signature it shows he= re:<br> >>> > <a href=3D"https://www.ebi.ac.uk/interpro/entry/pfam/PF08= 241/logo/" rel=3D"noreferrer" target=3D"_blank">https://www.ebi.ac.uk/inter= pro/entry/pfam/PF08241/logo/</a><br> >>> ><br> >>> > I'm not sure a) why the probabilities differ in the p= rofile relative to the seed alignment, or b) how to filter columns in the a= lignment by those that have a match in the model (see columns 4-6 in the al= ignment, 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:<br> >>> > align =3D Bio.AlignIO.read(handle, "stockholm")= <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_acids[s] / l= en(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 s= ee 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 not 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:Bi= [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.open-bio.or= g/mailman/listinfo/biopython</a><br> </blockquote></div> --00000000000063b61506116dd028-- --===============7035859266241881958== 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 --===============7035859266241881958==--