1
$\begingroup$

I'm trying to programmatically replace a set of amino acid residues on an MSA with a "Z" from a list of unaligned positions. Any ideas on how I could do this?

Input: a list of unaligned positions for each sequence + an MSA
Output: an MSA with those unaligned positions on the aligned MSA replaced with a Z for each sequence.

I have a multiple sequence alignment i'm working with and a list of sites on that alignment that are glycosylated. I was trying to run a script that would take the input of the positions as a list of characters along with the MSA to replace the Nth character with a "Z" to easily visualize the sites. So essentially if I feed it the following sequence

M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K

and positions

2, 5, 10, 11, 14

it should spit back

M----ZFLZFLCLZZFAZAQ-------------D-GKCG----TLSNKSPS--------------------K
Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ Hi @DannyMorales and welcome. This is likely fairly easy to do, but the question is difficult to understand. Any clarity would be helpful. $\endgroup$
    – M__ ♦
    Oct 23, 2020 at 21:26
  • $\begingroup$ Hey Michael! Thanks for the reply. I have a multiple sequence alignment i'm working with and a list of sites on that alignment that are glycosylated. I was trying to run a script that would take the input of the positions as a list of characters along with the MSA to replace the Nth character with a "Z" to easily visualize the sites. So essentially if I feed it the following sequence "M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K" and positions "2, 5, 10, 11, 14" it should spit back "M----ZFLZFLCLZZFAZAQ-------------D-GKCG----TLSNKSPS--------------------K" $\endgroup$
    – Danny Morales
    Oct 23, 2020 at 22:08
  • 2
    $\begingroup$ Understand, what you need to do is update your question with the information above and then provide the code you have so far. The obvious choice here is biopython AlignIO $\endgroup$
    – M__ ♦
    Oct 23, 2020 at 22:28
  • 2
    $\begingroup$ Could you please show your code, so that is clear what you tried so far? $\endgroup$
    – Mr_Z
    Oct 24, 2020 at 13:48

2 Answers 2

Reset to default
1
$\begingroup$

Another alternative, without the need for re:

seq = 'M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K'
pos = [2, 5, 10, 11, 14]

new_seq = []
count = 0
for s in seq:
    if s == "-":
        new_seq.append(s)
    else:
        count += 1
        if count in pos:
            new_seq.append('Z')
        else:
            new_seq.append(s)

print(''.join(new_seq))

Result:

M----ZFLZFLCLZZFAZAQ-------------D-GKCG----TLSNKSPS--------------------K
Improve this answer
$\endgroup$
10
  • 1
    $\begingroup$ @M__, to be honest I didn't know re was particularly slow in Python. In any case, I think that if efficiency is an issue, it would also be important to define the final size of the new_seq list beforehand (e.g. [None]*len(seq)), and then replace the Nones as one loops over seq. $\endgroup$
    – gaspanic
    Jul 28, 2022 at 13:40
  • 1
    $\begingroup$ Hi @gaspanic its very famous for it. It should be compiled via re.compile but its still a problem. The inbuilt regex - split replace that sort of thing - is usually used instead not for functionality (it ain't as good at all) but performance. Predefining the list size - I will look into that in detail, thanks and do some tests and timings. I'm well aware for C and Cython it is important. If I'm doing re I parallelise it. $\endgroup$
    – M__ ♦
    Jul 28, 2022 at 13:57
  • 1
    $\begingroup$ ... okay I tried predefining the list but didn't noticeably improve in my code. deque is another way in. Would make a good question. $\endgroup$
    – M__ ♦
    Jul 28, 2022 at 14:23
  • 1
    $\begingroup$ Yes, looks like you're right. See this thread: stackoverflow.com/a/63600055/8861903 . I guess I must have inherited that idea from other languages (or I'm just plain wrong). $\endgroup$
    – gaspanic
    Jul 28, 2022 at 22:45
  • 1
    $\begingroup$ @M__ I actually tried extending the sequence 100000x, but the difference was still minimal. I even removed the if statements to see if that was the bottleneck, but still not much difference. It appears that the explanation lies in how Python allocates memory to lists :). $\endgroup$
    – gaspanic
    Jul 28, 2022 at 23:52
1
$\begingroup$

The code is

import re

pos = [2, 5, 10, 11, 14]
newseq = []
prot = 'M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K'
count = 0
for n,s in enumerate(prot):    
    if str(s) == '-':
        count += 1
        newseq.append('-')
    else:
        if (n + 1 - count) in pos:
            s = re.sub(r'[A-Za-z]', 'Z', s)
        newseq.append(s)
print(''.join(newseq))

output

M----ZFLZFLCLZZFAZAQ-------------D-GKCG----TLSNKSPS--------------------K

Precisely as requested, apologies for the delay.

Happy to explain why it works - its just a couple of Python tricks and a method. It must be run under Python 3 (unlikely to work in Python 2) and probably >3.6.

BTW I normally avoid re for complicated reasons but its easy, hence its here.

As @gaspanic has pointed out re is not only unnecessary it will slow the code (re is very inefficient). It is in fact one of the 'legendary inefficiencies' of Python. Thus removing it ...

# variable declarations omitted
newseq = list()
count = 0
for n,s in enumerate(prot):    
    if str(s) == '-':
        count += 1
        newseq.append('-')
    else:
        if (n + 1 - count) in pos:
            s = 'Z'
        newseq.append(s)
print(''.join(newseq))

Same output.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.