[Biopython] Struggling with MSA...

[Dan Bolser]

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

Methyltransf_11 PF08241.15 CL0063 63 173 63 172 1 95 51.30 1.5e-13 3.1e-17
Methyltransferase domain

---------------------------------------------------------------------------------------------------------------------------------

MODEL
 LdvGcGtGrlaealakrg.arvvgvDlskemlklakekaseeglkvefvvadaeklpfednsfDlvvssevlhhv...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-QLRWETMDVRKLDFPSASFDVVLEKGTLDALlagErdpwtvssegvhTVDQVLSEVSRVLVPGGRFI


So I guess it's possible to work backwards from there...



On Thu, 15 Feb 2024 at 16:02, Peter Cock <p.j.a.cock at googlemail.com> 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 <= to match, or perhaps they use a more
> sophisticated criteria.
>
> Peter
>
> On Thu, Feb 15, 2024 at 3:56 PM Dan Bolser <dan.bolser at outsee.co.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 <dan.bolser at outsee.co.uk>
> wrote:
> >>
> >> Hi Peter,
> >>
> >> On Thu, 15 Feb 2024 at 15:28, Peter Cock <p.j.a.cock at googlemail.com>
> wrote:
> >>>
> >>> Hello Dan,
> >>>
> >>> It that is a probability at the end, you have the wrong denominator -
> >>> should be len(align) == 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 PM Dan Bolser <dan.bolser at 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=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.
> >>> > _______________________________________________
> >>> > Biopython mailing list  -  Biopython at biopython.org
> >>> > https://mailman.open-bio.org/mailman/listinfo/biopython
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.open-bio.org/pipermail/biopython/attachments/20240215/9f8844fb/attachment.htm>