Re: [Biopython] Struggling with MSA...
Dan Bolser <[email protected]> Thu, 15 Feb 2024 15:56:58 +0000
| Newsgroups | gmane.comp.python.bio.general |
|---|---|
| Message-ID | <CANs1yPLyQs0cd8G_eRNY2+OM-ry52rbB3Y5Wn7ebXCD6a+6P3w@mail.gmail.com> |
--===============6991427981698599123== Content-Type: multipart/alternative; boundary="00000000000017859c06116dac96" --00000000000017859c06116dac96 Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: quoted-printable 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 report= s > -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 <[email protected]= .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 t= o >> the seed alignment, or b) how to filter columns in the alignment by thos= e >> 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 (where it matche= s >> 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 >> > --00000000000017859c06116dac96 Content-Type: text/html; charset="UTF-8" Content-Transfer-Encoding: quoted-printable <div dir=3D"ltr">Sorry, by 'half gaps' I thought you meant a specif= ic thing, like, "look at that half gap!". I guess you mean, colum= ns where half of the values in the column are gaps.<div><br></div><div>Look= ing directly at the alignment, I guess there is no way to work out precisel= y=C2=A0which columns are gaps or not in the model?</div><div><br></div><div= ><br></div><div><br></div><div><br></div><div><br></div></div><br><div clas= s=3D"gmail_quote"><div dir=3D"ltr" class=3D"gmail_attr">On Thu, 15 Feb 2024= at 15:47, Dan Bolser <<a href=3D"mailto:[email protected]">dan.bo= [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);p= adding-left:1ex"><div dir=3D"ltr"><div dir=3D"ltr">Hi Peter,</div><br><div = class=3D"gmail_quote"><div dir=3D"ltr" class=3D"gmail_attr">On Thu, 15 Feb = 2024 at 15:28, Peter Cock <<a href=3D"mailto:[email protected]" = target=3D"_blank">[email protected]</a>> wrote:<br></div><blockq= uote class=3D"gmail_quote" style=3D"margin:0px 0px 0px 0.8ex;border-left:1p= x solid rgb(204,204,204);padding-left:1ex">Hello Dan,<br> <br> It that is a probability at the end, you have the wrong denominator -<br> should be len(align) =3D=3D sum(amino_acids.values()<br></blockquote><div><= br></div><div>Makes no difference, they are the same (including -s). The pr= ofile reports -s as 'Insert Probability'.</div><div><br></div><div>= <br></div><blockquote class=3D"gmail_quote" style=3D"margin:0px 0px 0px 0.8= ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Note this seed = alignment has 164 columns, yet the model page says the<br> model has 95 columns. If I count the number of columns which are half<br> gaps that's pretty close...<br></blockquote><div><br></div><div>I don&#= 39;t understand what that means or implies... Is something=C2=A0wrong somew= here?</div><div><br></div><div><br></div><div>Cheers,</div><div>Dan.</div><= div><br></div><div><br></div><blockquote class=3D"gmail_quote" style=3D"mar= gin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1= ex"> Peter<br> <br> On Thu, Feb 15, 2024 at 3:09=E2=80=AFPM Dan Bolser <<a href=3D"mailto:da= [email protected]" target=3D"_blank">[email protected]</a>> wr= ote:<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 alignment for this domain:<b= r> > <a href=3D"https://www.ebi.ac.uk/interpro/entry/pfam/PF08241/entry_ali= gnments/?type=3Dseed" rel=3D"noreferrer" target=3D"_blank">https://www.ebi.= ac.uk/interpro/entry/pfam/PF08241/entry_alignments/?type=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/" re= l=3D"noreferrer" target=3D"_blank">https://www.ebi.ac.uk/interpro/entry/pfa= m/PF08241/logo/</a><br> ><br> > I'm not sure a) why the probabilities differ in the profile relati= ve to the seed alignment, or b) how to filter columns in the alignment by t= hose that have a match in the model (see columns 4-6 in the alignment, whic= h are gaps in the model).<br> ><br> > I think if I can answer b) then the answer to a) will be, "look a= t 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] / 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 conser= vation score, base by base on a given protein in the alignment (where it ma= tches 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:Biopython@biopy= thon.org" target=3D"_blank">[email protected]</a><br> > <a href=3D"https://mailman.open-bio.org/mailman/listinfo/biopython" re= l=3D"noreferrer" target=3D"_blank">https://mailman.open-bio.org/mailman/lis= tinfo/biopython</a><br> </blockquote></div></div> </blockquote></div> --00000000000017859c06116dac96-- --===============6991427981698599123== 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 --===============6991427981698599123==--