Re: [Biopython] Struggling with MSA...
Dan Bolser <[email protected]> Thu, 15 Feb 2024 15:47:33 +0000
| Newsgroups | gmane.comp.python.bio.general |
|---|---|
| Message-ID | <CANs1yPKGqPKQxK+6ab2wNH-HsSQkEkoa5pm12cFsqVZrGeC7ow@mail.gmail.com> |
--===============4371104272496301092== Content-Type: multipart/alternative; boundary="0000000000007679e806116d8a7b" --0000000000007679e806116d8a7b Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: quoted-printable 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 <[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 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 a= re > 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 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 > --0000000000007679e806116d8a7b Content-Type: text/html; charset="UTF-8" Content-Transfer-Encoding: quoted-printable <div dir=3D"ltr"><div dir=3D"ltr">Hi Peter,</div><br><div class=3D"gmail_qu= ote"><div dir=3D"ltr" class=3D"gmail_attr">On Thu, 15 Feb 2024 at 15:28, Pe= ter Cock <<a href=3D"mailto:[email protected]">p.j.a.cock@google= mail.com</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-le= ft: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> --0000000000007679e806116d8a7b-- --===============4371104272496301092== 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 --===============4371104272496301092==--