[Biopython] Struggling with MSA...

[Dan Bolser]

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/3e6ffa91/attachment-0001.htm>