[Biopython] Struggling with MSA...

--===============0863626562050933289==
Content-Type: multipart/alternative; boundary="0000000000000c3d9a06116d04b8"

--0000000000000c3d9a06116d04b8
Content-Type: text/plain; charset="UTF-8"

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=seed

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 = "PF08241.alignment.full.gz"
msa = "PF08241.alignment.seed.gz"

with gzip.open(msa, "rt") as handle:
align = Bio.AlignIO.read(handle, "stockholm")
ncols = align.get_alignment_length()

for col in range(ncols):
amino_acids = dict()
for s in align[:, col]:
amino_acids[s] = 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.

--0000000000000c3d9a06116d04b8
Content-Type: text/html; charset="UTF-8"
Content-Transfer-Encoding: quoted-printable

<div dir=3D"ltr">Hi,<div><br></div><div>Sorry, I can&#39;t follow the docs =
(or find the right docs).</div><div><br></div><div>I&#39;ve got the &#39;se=
ed&#39; stockholm alignment=C2=A0for this domain:</div><div><a href=3D"http=
s://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></div><div><br></div><div>and I&#39;m trying to reproduce th=
e signature it shows here:</div><div><a href=3D"https://www.ebi.ac.uk/inter=
pro/entry/pfam/PF08241/logo/">https://www.ebi.ac.uk/interpro/entry/pfam/PF0=
8241/logo/</a><br></div><div><br></div><div>I&#39;m not sure a) why the pro=
babilities 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).</div><div>=
<br></div><div>I think if I can answer b) then the answer to a) will be, &q=
uot;look at the full alignment&quot;.</div><div><br></div><div>Here is my c=
rude &#39;best guess&#39; code:</div><div><br></div><div><div style=3D"colo=
r:rgb(204,204,204);background-color:rgb(31,31,31);font-family:&quot;Droid S=
ans Mono&quot;,&quot;monospace&quot;,monospace;font-size:14px;line-height:1=
9px;white-space:pre"><div><span style=3D"color:rgb(197,134,192)">import</sp=
an> <span style=3D"color:rgb(78,201,176)">gzip</span></div><div><span style=
=3D"color:rgb(197,134,192)">import</span> <span style=3D"color:rgb(78,201,1=
76)">Bio</span>.<span style=3D"color:rgb(78,201,176)">AlignIO</span></div><=
br><div><span style=3D"color:rgb(106,153,85)"># msa =3D &quot;PF08241.align=
ment.full.gz&quot;</span></div><div><span style=3D"color:rgb(156,220,254)">=
msa</span> <span style=3D"color:rgb(212,212,212)">=3D</span> <span style=3D=
"color:rgb(206,145,120)">&quot;PF08241.alignment.seed.gz&quot;</span></div>=
<br><div><span style=3D"color:rgb(197,134,192)">with</span> <span style=3D"=
color:rgb(78,201,176)">gzip</span>.<span style=3D"color:rgb(220,220,170)">o=
pen</span>(<span style=3D"color:rgb(156,220,254)">msa</span>, <span style=
=3D"color:rgb(206,145,120)">&quot;rt&quot;</span>) <span style=3D"color:rgb=
(197,134,192)">as</span> <span style=3D"color:rgb(156,220,254)">handle</spa=
n>:</div><div>    <span style=3D"color:rgb(156,220,254)">align</span> <span=
 style=3D"color:rgb(212,212,212)">=3D</span> <span style=3D"color:rgb(78,20=
1,176)">Bio</span>.<span style=3D"color:rgb(78,201,176)">AlignIO</span>.<sp=
an style=3D"color:rgb(220,220,170)">read</span>(<span style=3D"color:rgb(15=
6,220,254)">handle</span>, <span style=3D"color:rgb(206,145,120)">&quot;sto=
ckholm&quot;</span>)</div><div>    <span style=3D"color:rgb(156,220,254)">n=
cols</span> <span style=3D"color:rgb(212,212,212)">=3D</span> <span style=
=3D"color:rgb(156,220,254)">align</span>.<span style=3D"color:rgb(220,220,1=
70)">get_alignment_length</span>()</div><br><div>    <span style=3D"color:r=
gb(197,134,192)">for</span> <span style=3D"color:rgb(156,220,254)">col</spa=
n> <span style=3D"color:rgb(197,134,192)">in</span> <span style=3D"color:rg=
b(78,201,176)">range</span>(<span style=3D"color:rgb(156,220,254)">ncols</s=
pan>):</div><div>        <span style=3D"color:rgb(156,220,254)">amino_acids=
</span> <span style=3D"color:rgb(212,212,212)">=3D</span> <span style=3D"co=
lor:rgb(78,201,176)">dict</span>()</div><div>        <span style=3D"color:r=
gb(197,134,192)">for</span> <span style=3D"color:rgb(156,220,254)">s</span>=
 <span style=3D"color:rgb(197,134,192)">in</span> <span style=3D"color:rgb(=
156,220,254)">align</span>[:, <span style=3D"color:rgb(156,220,254)">col</s=
pan>]:</div><div>            <span style=3D"color:rgb(156,220,254)">amino_a=
cids</span>[<span style=3D"color:rgb(156,220,254)">s</span>] <span style=3D=
"color:rgb(212,212,212)">=3D</span> <span style=3D"color:rgb(156,220,254)">=
amino_acids</span>.<span style=3D"color:rgb(220,220,170)">get</span>(<span =
style=3D"color:rgb(156,220,254)">s</span>, <span style=3D"color:rgb(181,206=
,168)">0</span>) <span style=3D"color:rgb(212,212,212)">+</span> <span styl=
e=3D"color:rgb(181,206,168)">1</span></div><div>        <span style=3D"colo=
r:rgb(220,220,170)">print</span>(<span style=3D"color:rgb(156,220,254)">ami=
no_acids</span>)</div><div>        <span style=3D"color:rgb(197,134,192)">f=
or</span> <span style=3D"color:rgb(156,220,254)">s</span> <span style=3D"co=
lor:rgb(197,134,192)">in</span> <span style=3D"color:rgb(156,220,254)">amin=
o_acids</span>:</div><div>            <span style=3D"color:rgb(220,220,170)=
">print</span>(<span style=3D"color:rgb(86,156,214)">f</span><span style=3D=
"color:rgb(206,145,120)">&quot;</span><span style=3D"color:rgb(86,156,214)"=
>{</span><span style=3D"color:rgb(156,220,254)">s</span><span style=3D"colo=
r:rgb(86,156,214)">}</span><span style=3D"color:rgb(206,145,120)">: </span>=
<span style=3D"color:rgb(86,156,214)">{</span><span style=3D"color:rgb(156,=
220,254)">amino_acids</span>[<span style=3D"color:rgb(156,220,254)">s</span=
>]<span style=3D"color:rgb(86,156,214)">:3d}</span><span style=3D"color:rgb=
(206,145,120)"> </span><span style=3D"color:rgb(86,156,214)">{</span><span =
style=3D"color:rgb(156,220,254)">amino_acids</span>[<span style=3D"color:rg=
b(156,220,254)">s</span>] <span style=3D"color:rgb(212,212,212)">/</span> <=
span style=3D"color:rgb(220,220,170)">len</span>(<span style=3D"color:rgb(1=
56,220,254)">align</span>)<span style=3D"color:rgb(86,156,214)">:.3f}</span=
><span style=3D"color:rgb(206,145,120)">&quot;</span>)</div><br></div></div=
><div><br></div><div><br></div><div>I have the feeling I&#39;m doin it rong=
...</div><div><br></div><div>The above is just a &#39;warm up&#39;, really =
I want to see the conservation score, base by base on a given protein in th=
e alignment (where it matches the model).</div><div><br></div><div>Many tha=
nks for any suggestions, and sorry for not being able to find the right doc=
ument to answer these questions.</div><div><br></div><div><br></div><div>kt=
hxbi,</div><div>Dan.</div></div>

--0000000000000c3d9a06116d04b8--

--===============0863626562050933289==
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

--===============0863626562050933289==--