Re: [Biopython] Struggling with MSA...
Michiel de Hoon <[email protected]> Fri, 16 Feb 2024 23:42:01 +0000 (UTC)
| Newsgroups | gmane.comp.python.bio.general |
|---|---|
| Message-ID | <[email protected]> |
--===============7844582813674361105==
Content-Type: multipart/alternative;
boundary="----=_Part_3582539_1018162062.1708126921125"
------=_Part_3582539_1018162062.1708126921125
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: quoted-printable
You could use the new parser in Bio.Align:
>>> import gzip>>> from Bio import Align>>> from Bio.motifs import Motif
>>> from Bio.Data import IUPACData
>>> stream =3D gzip.open("PF08241.alignment.seed.gz", "rt")>>> alignment =
=3D Align.read(stream, "Stockholm")>>> indices =3D [index for index, c in e=
numerate(alignment.column_annotations['consensus sequence']) if c!=3D"."]>>=
> alignment =3D alignment[:, indices]>>> m.consensus
Seq('LDVGCGTGLLTRALARLGARVTGVDLSPEMLELARERAPRAVVGDAEDLPFPDN...LII')
>>> m.weblogo("PF08241.png")
-Michiel
On Friday, February 16, 2024 at 07:59:14 PM GMT+9, Dan Bolser <dan.bols=
[email protected]> wrote: =20
=20
I don't suppose there is a package for reading in the HMM?=C2=A0=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 co=
nservation 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=3Dhmm
Peter
On Thu, Feb 15, 2024 at 4:06=E2=80=AFPM Dan Bolser <[email protected]=
> wrote:
>
> Searching HMMER gives me (cryptic?) details of the model match:
>
> Family-Id Family-Accession Clan Start End Ali-Start Ali-End Model-Start M=
odel-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-1=
7 Methyltransferase domain
>
> -------------------------------------------------------------------------=
--------------------------------------------------------
>
> MODEL=C2=A0 LdvGcGtGrlaealakrg.arvvgvDlskemlklakekaseeglkvefvvadaeklpfedn=
sfDlvvssevlhhv...e............dpekalkeiaRvLkpgGllv
> 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+++
> PPL=C2=A0 =C2=A0 569************7779*************77766666666.69**********=
*******************8662667778888888******************97
> SEQ=C2=A0 =C2=A0 LVLGCGNSALSYELFLGGfPNVTSVDYSSVVVAAMQARHAHVP-QLRWETMDVRKL=
DFPSASFDVVLEKGTLDALlagErdpwtvssegvhTVDQVLSEVSRVLVPGGRFI
>
>
> So I guess it's possible to work backwards from there...
>
>
>
> On Thu, 15 Feb 2024 at 16:02, Peter Cock <[email protected]> wrot=
e:
>>
>> 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, "loo=
k 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]> wro=
te:
>> >>
>> >> 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 re=
ports -s as 'Insert Probability'.
>> >>
>> >>
>> >>> Note this seed alignment has 164 columns, yet the model page says th=
e
>> >>> model has 95 columns. If I count the number of columns which are hal=
f
>> >>> gaps that's pretty close...
>> >>
>> >>
>> >> I don't understand what that means or implies... Is something wrong s=
omewhere?
>> >>
>> >>
>> >> Cheers,
>> >> Dan.
>> >>
>> >>
>> >>> Peter
>> >>>
>> >>> On Thu, Feb 15, 2024 at 3:09=E2=80=AFPM Dan Bolser <dan.bolser@outse=
e.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 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).
>> >>> >
>> >>> > 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 conservati=
on 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 f=
ind the right document to answer these questions.
>> >>> >
>> >>> >
>> >>> > kthxbi,
>> >>> > Dan.
>> >>> > _______________________________________________
>> >>> > Biopython mailing list=C2=A0 -=C2=A0 [email protected]
>> >>> > https://mailman.open-bio.org/mailman/listinfo/biopython
_______________________________________________
Biopython mailing list=C2=A0 -=C2=A0 [email protected]
https://mailman.open-bio.org/mailman/listinfo/biopython
=20
------=_Part_3582539_1018162062.1708126921125
Content-Type: text/html; charset=UTF-8
Content-Transfer-Encoding: quoted-printable
<html><head></head><body><div class=3D"ydpce4c934eyahoo-style-wrap" style=
=3D"font-family:Helvetica Neue, Helvetica, Arial, sans-serif;font-size:10px=
;"><div></div>
<div dir=3D"ltr" data-setdir=3D"false">You could use the new parser=
in Bio.Align:</div><div dir=3D"ltr" data-setdir=3D"false"><br></div><div d=
ir=3D"ltr" data-setdir=3D"false"><div dir=3D"ltr" data-setdir=3D"false"><sp=
an><span><span><span><span>>>> </span></span></span></span></span>=
import gzip</div><div dir=3D"ltr" data-setdir=3D"false"><span><span><span><=
span><span>>>> </span></span></span></span></span>from Bio import =
Align</div><div dir=3D"ltr" data-setdir=3D"false"><div dir=3D"ltr" data-set=
dir=3D"false"><span><span><span><span><span>>>> </span></span></sp=
an></span></span>from Bio.motifs import Motif<br><div dir=3D"ltr" data-setd=
ir=3D"false"><span><span><span><span><span>>>> </span></span></spa=
n></span></span>from Bio.Data import IUPACData<br></div></div></div><div di=
r=3D"ltr" data-setdir=3D"false"><span><span><span><span><span><span>>>=
;> </span></span></span></span></span>stream =3D gzip.open("PF08241.alig=
nment.seed.gz", "rt")</span></div><div dir=3D"ltr" data-setdir=3D"false"><s=
pan><span><span><span><span><span><span>>>> </span></span></span><=
/span></span>alignment =3D Align.read(stream, "Stockholm")</span></span></d=
iv><div dir=3D"ltr" data-setdir=3D"false"><span><span><span><span><span><sp=
an><span><span>>>> </span></span></span></span></span>indices =3D =
[index for index, c in enumerate(alignment.column_annotations['consensus se=
quence']) if c!=3D"."]</span></span></span></div><div dir=3D"ltr" data-setd=
ir=3D"false"><span><span><span><span>>>> alignment =3D alignment[:=
, indices]</span></span></span></span></div><div dir=3D"ltr" data-setdir=3D=
"false"><span><span><div dir=3D"ltr" data-setdir=3D"false">>>> m.c=
onsensus<br>Seq('LDVGCGTGLLTRALARLGARVTGVDLSPEMLELARERAPRAVVGDAEDLPFPDN...L=
II')<br><div>>>> m.weblogo("PF08241.png")<br><br></div></div><div =
dir=3D"ltr" data-setdir=3D"false">-Michiel<br></div></span></span></div><br=
></div><div><br></div>
=20
</div><div id=3D"yahoo_quoted_8854620181" class=3D"yahoo_quoted">
<div style=3D"font-family:'Helvetica Neue', Helvetica, Arial, s=
ans-serif;font-size:13px;color:#26282a;">
=20
<div>
On Friday, February 16, 2024 at 07:59:14 PM GMT+9, Dan =
Bolser <[email protected]> wrote:
</div>
<div><br></div>
<div><br></div>
<div><div id=3D"yiv4136689180"><div><div dir=3D"ltr">I don'=
t suppose there is a package for reading in the HMM? =F0=9F=98=85<div>=
<br clear=3D"none"></div></div><br clear=3D"none"><div id=3D"yiv4136689180y=
qt23012" class=3D"yiv4136689180yqt8247348347"><div class=3D"yiv4136689180gm=
ail_quote"><div dir=3D"ltr" class=3D"yiv4136689180gmail_attr">On Thu, 15 Fe=
b 2024 at 16:18, Dan Bolser <<a rel=3D"nofollow noopener noreferrer" sha=
pe=3D"rect" ymailto=3D"mailto:[email protected]" target=3D"_blank" hr=
ef=3D"mailto:[email protected]">[email protected]</a>> wrote=
:<br clear=3D"none"></div><blockquote style=3D"margin:0px 0px 0px 0.8ex;bor=
der-left:1px solid rgb(204,204,204);padding-left:1ex;" class=3D"yiv41366891=
80gmail_quote"><div dir=3D"ltr">Recreating the profile was just a sanity ch=
eck. Really I want to look at conservation across the protein of interest.<=
/div><br clear=3D"none"><div class=3D"yiv4136689180gmail_quote"><div dir=3D=
"ltr" class=3D"yiv4136689180gmail_attr">On Thu, 15 Feb 2024 at 16:13, Peter=
Cock <<a rel=3D"nofollow noopener noreferrer" shape=3D"rect" ymailto=3D=
"mailto:[email protected]" target=3D"_blank" href=3D"mailto:p.j.a.c=
[email protected]">[email protected]</a>> wrote:<br clear=3D"no=
ne"></div><blockquote style=3D"margin:0px 0px 0px 0.8ex;border-left:1px sol=
id rgb(204,204,204);padding-left:1ex;" class=3D"yiv4136689180gmail_quote">Y=
ou may have more joy starting from the HMM for the model?<br clear=3D"none"=
>
<br clear=3D"none">
<a rel=3D"nofollow noopener noreferrer" shape=3D"rect" target=3D"_blank" hr=
ef=3D"https://www.ebi.ac.uk/interpro/wwwapi//entry/pfam/PF08241?annotation=
=3Dhmm">https://www.ebi.ac.uk/interpro/wwwapi//entry/pfam/PF08241?annotatio=
n=3Dhmm</a><br clear=3D"none">
<br clear=3D"none">
Peter<br clear=3D"none">
<br clear=3D"none">
On Thu, Feb 15, 2024 at 4:06=E2=80=AFPM Dan Bolser <<a rel=3D"nofollow n=
oopener noreferrer" shape=3D"rect" ymailto=3D"mailto:[email protected]=
k" target=3D"_blank" href=3D"mailto:[email protected]">dan.bolser@out=
see.co.uk</a>> wrote:<br clear=3D"none">
><br clear=3D"none">
> Searching HMMER gives me (cryptic?) details of the model match:<br cle=
ar=3D"none">
><br clear=3D"none">
> 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 clear=3D"no=
ne">
> =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=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 clear=3D"none">
><br clear=3D"none">
> Methyltransf_11 PF08241.15 CL0063 63 173 63 172 1 95 51.30 1.5e-13 3.1=
e-17 Methyltransferase domain<br clear=3D"none">
><br clear=3D"none">
> ----------------------------------------------------------------------=
-----------------------------------------------------------<br clear=3D"non=
e">
><br clear=3D"none">
> MODEL LdvGcGtGrlaealakrg.arvvgvDlskemlklakekaseeglkvefvvadaeklpf=
ednsfDlvvssevlhhv...e............dpekalkeiaRvLkpgGllv<br clear=3D"none">
> MATCH L +GcG+ +l+ +l g +v+ vD+s ++++ +++=
+ ++++ +d++kl f+++sfD+v+ + +l+ + e &n=
bsp; ++++l+e+ RvL pgG+++<br clear=3D"none=
">
> PPL 569************7779*************77766666666.69*******=
**********************8662667778888888******************97<br clear=3D"none=
">
> SEQ LVLGCGNSALSYELFLGGfPNVTSVDYSSVVVAAMQARHAHVP-QLRWETMDV=
RKLDFPSASFDVVLEKGTLDALlagErdpwtvssegvhTVDQVLSEVSRVLVPGGRFI<br clear=3D"none=
">
><br clear=3D"none">
><br clear=3D"none">
> So I guess it's possible to work backwards from there...<br clear=3D"n=
one">
><br clear=3D"none">
><br clear=3D"none">
><br clear=3D"none">
> On Thu, 15 Feb 2024 at 16:02, Peter Cock <<a rel=3D"nofollow noopen=
er noreferrer" shape=3D"rect" ymailto=3D"mailto:[email protected]" =
target=3D"_blank" href=3D"mailto:[email protected]">p.j.a.cock@goog=
lemail.com</a>> wrote:<br clear=3D"none">
>><br clear=3D"none">
>> Sorry yes, I used this as a threshold for if a column was gappy, a=
nd<br clear=3D"none">
>> likely to be in the model or not:<br clear=3D"none">
>><br clear=3D"none">
>> amino_acids.get("-", 0) < len(align)*0.5<br clear=3D"none">
>><br clear=3D"none">
>> Might need to use <=3D to match, or perhaps they use a more<br =
clear=3D"none">
>> sophisticated criteria.<br clear=3D"none">
>><br clear=3D"none">
>> Peter<br clear=3D"none">
>><br clear=3D"none">
>> On Thu, Feb 15, 2024 at 3:56=E2=80=AFPM Dan Bolser <<a rel=3D"n=
ofollow noopener noreferrer" shape=3D"rect" ymailto=3D"mailto:dan.bolser@ou=
tsee.co.uk" target=3D"_blank" href=3D"mailto:[email protected]">dan.b=
[email protected]</a>> wrote:<br clear=3D"none">
>> ><br clear=3D"none">
>> > Sorry, by 'half gaps' I thought you meant a specific thing, l=
ike, "look at that half gap!". I guess you mean, columns where half of the =
values in the column are gaps.<br clear=3D"none">
>> ><br clear=3D"none">
>> > 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 clear=3D=
"none">
>> ><br clear=3D"none">
>> ><br clear=3D"none">
>> ><br clear=3D"none">
>> ><br clear=3D"none">
>> ><br clear=3D"none">
>> ><br clear=3D"none">
>> > On Thu, 15 Feb 2024 at 15:47, Dan Bolser <<a rel=3D"nofoll=
ow noopener noreferrer" shape=3D"rect" ymailto=3D"mailto:dan.bolser@outsee.=
co.uk" target=3D"_blank" href=3D"mailto:[email protected]">dan.bolser=
@outsee.co.uk</a>> wrote:<br clear=3D"none">
>> >><br clear=3D"none">
>> >> Hi Peter,<br clear=3D"none">
>> >><br clear=3D"none">
>> >> On Thu, 15 Feb 2024 at 15:28, Peter Cock <<a rel=3D"no=
follow noopener noreferrer" shape=3D"rect" ymailto=3D"mailto:p.j.a.cock@goo=
glemail.com" target=3D"_blank" href=3D"mailto:[email protected]">p.=
[email protected]</a>> wrote:<br clear=3D"none">
>> >>><br clear=3D"none">
>> >>> Hello Dan,<br clear=3D"none">
>> >>><br clear=3D"none">
>> >>> It that is a probability at the end, you have the wro=
ng denominator -<br clear=3D"none">
>> >>> should be len(align) =3D=3D sum(amino_acids.values()<=
br clear=3D"none">
>> >><br clear=3D"none">
>> >><br clear=3D"none">
>> >> Makes no difference, they are the same (including -s). Th=
e profile reports -s as 'Insert Probability'.<br clear=3D"none">
>> >><br clear=3D"none">
>> >><br clear=3D"none">
>> >>> Note this seed alignment has 164 columns, yet the mod=
el page says the<br clear=3D"none">
>> >>> model has 95 columns. If I count the number of column=
s which are half<br clear=3D"none">
>> >>> gaps that's pretty close...<br clear=3D"none">
>> >><br clear=3D"none">
>> >><br clear=3D"none">
>> >> I don't understand what that means or implies... Is somet=
hing wrong somewhere?<br clear=3D"none">
>> >><br clear=3D"none">
>> >><br clear=3D"none">
>> >> Cheers,<br clear=3D"none">
>> >> Dan.<br clear=3D"none">
>> >><br clear=3D"none">
>> >><br clear=3D"none">
>> >>> Peter<br clear=3D"none">
>> >>><br clear=3D"none">
>> >>> On Thu, Feb 15, 2024 at 3:09=E2=80=AFPM Dan Bolser &l=
t;<a rel=3D"nofollow noopener noreferrer" shape=3D"rect" ymailto=3D"mailto:=
[email protected]" target=3D"_blank" href=3D"mailto:dan.bolser@outsee=
.co.uk">[email protected]</a>> wrote:<br clear=3D"none">
>> >>> ><br clear=3D"none">
>> >>> > Hi,<br clear=3D"none">
>> >>> ><br clear=3D"none">
>> >>> > Sorry, I can't follow the docs (or find the righ=
t docs).<br clear=3D"none">
>> >>> ><br clear=3D"none">
>> >>> > I've got the 'seed' stockholm alignment for this=
domain:<br clear=3D"none">
>> >>> > <a rel=3D"nofollow noopener noreferrer" shape=3D=
"rect" target=3D"_blank" href=3D"https://www.ebi.ac.uk/interpro/entry/pfam/=
PF08241/entry_alignments/?type=3Dseed">https://www.ebi.ac.uk/interpro/entry=
/pfam/PF08241/entry_alignments/?type=3Dseed</a><br clear=3D"none">
>> >>> ><br clear=3D"none">
>> >>> > and I'm trying to reproduce the signature it sho=
ws here:<br clear=3D"none">
>> >>> > <a rel=3D"nofollow noopener noreferrer" shape=3D=
"rect" target=3D"_blank" href=3D"https://www.ebi.ac.uk/interpro/entry/pfam/=
PF08241/logo/">https://www.ebi.ac.uk/interpro/entry/pfam/PF08241/logo/</a><=
br clear=3D"none">
>> >>> ><br clear=3D"none">
>> >>> > 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 t=
he alignment, which are gaps in the model).<br clear=3D"none">
>> >>> ><br clear=3D"none">
>> >>> > I think if I can answer b) then the answer to a)=
will be, "look at the full alignment".<br clear=3D"none">
>> >>> ><br clear=3D"none">
>> >>> > Here is my crude 'best guess' code:<br clear=3D"=
none">
>> >>> ><br clear=3D"none">
>> >>> > import gzip<br clear=3D"none">
>> >>> > import Bio.AlignIO<br clear=3D"none">
>> >>> ><br clear=3D"none">
>> >>> > # msa =3D "PF08241.alignment.full.gz"<br clear=
=3D"none">
>> >>> > msa =3D "PF08241.alignment.seed.gz"<br clear=3D"=
none">
>> >>> ><br clear=3D"none">
>> >>> > with gzip.open(msa, "rt") as handle:<br clear=3D=
"none">
>> >>> > align =3D Bio.AlignIO.read(handle, "stockholm")<=
br clear=3D"none">
>> >>> > ncols =3D align.get_alignment_length()<br clear=
=3D"none">
>> >>> ><br clear=3D"none">
>> >>> > for col in range(ncols):<br clear=3D"none">
>> >>> > amino_acids =3D dict()<br clear=3D"none">
>> >>> > for s in align[:, col]:<br clear=3D"none">
>> >>> > amino_acids[s] =3D amino_acids.get(s, 0) + 1<br =
clear=3D"none">
>> >>> > print(amino_acids)<br clear=3D"none">
>> >>> > for s in amino_acids:<br clear=3D"none">
>> >>> > print(f"{s}: {amino_acids[s]:3d} {amino_acids[s]=
/ len(align):.3f}")<br clear=3D"none">
>> >>> ><br clear=3D"none">
>> >>> ><br clear=3D"none">
>> >>> ><br clear=3D"none">
>> >>> > I have the feeling I'm doin it rong...<br clear=
=3D"none">
>> >>> ><br clear=3D"none">
>> >>> > The above is just a 'warm up', really I want to =
see the conservation score, base by base on a given protein in the alignmen=
t (where it matches the model).<br clear=3D"none">
>> >>> ><br clear=3D"none">
>> >>> > Many thanks for any suggestions, and sorry for n=
ot being able to find the right document to answer these questions.<br clea=
r=3D"none">
>> >>> ><br clear=3D"none">
>> >>> ><br clear=3D"none">
>> >>> > kthxbi,<br clear=3D"none">
>> >>> > Dan.<br clear=3D"none">
>> >>> > _______________________________________________<=
br clear=3D"none">
>> >>> > Biopython mailing list - <a rel=3D"n=
ofollow noopener noreferrer" shape=3D"rect" ymailto=3D"mailto:Biopython@bio=
python.org" target=3D"_blank" href=3D"mailto:[email protected]">Biopy=
[email protected]</a><br clear=3D"none">
>> >>> > <a rel=3D"nofollow noopener noreferrer" shape=3D=
"rect" target=3D"_blank" href=3D"https://mailman.open-bio.org/mailman/listi=
nfo/biopython">https://mailman.open-bio.org/mailman/listinfo/biopython</a><=
br clear=3D"none">
</blockquote></div>
</blockquote></div></div>
</div></div><div class=3D"yqt8247348347" id=3D"yqt57419">__________________=
_____________________________<br clear=3D"none">Biopython mailing list =
; - <a shape=3D"rect" ymailto=3D"mailto:[email protected]" href=
=3D"mailto:[email protected]">[email protected]</a><br clear=3D=
"none"><a shape=3D"rect" href=3D"https://mailman.open-bio.org/mailman/listi=
nfo/biopython" target=3D"_blank">https://mailman.open-bio.org/mailman/listi=
nfo/biopython</a><br clear=3D"none"></div></div>
</div>
</div></body></html>
------=_Part_3582539_1018162062.1708126921125--
--===============7844582813674361105==
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
--===============7844582813674361105==--