<p>The code is</p>
<pre><code>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))
</code></pre>
<p>output</p>
<blockquote>
<p>M----ZFLZFLCLZZFAZAQ-------------D-GKCG----TLSNKSPS--------------------K</p>
</blockquote>
<p>Precisely as requested, apologies for the delay.</p>
<p>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 &gt;3.6.</p>
<p>BTW I normally avoid <code>re</code> for complicated reasons but its easy, hence its here.</p>
<hr />
<p>As @gaspanic has pointed out <code>re</code> is not only unnecessary it will slow the code (<code>re</code> is very inefficient). It is in fact one of the 'legendary inefficiencies' of Python. Thus removing it ...</p>
<pre><code># 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))
</code></pre>
<p>Same output.</p>

CC BY-SA 4.0
2022-07-27T20:03:23.030Z
2022-07-28T02:16:22.137Z
2022-07-28T02:16:22.137Z
M__
<p>Another alternative, without the need for <code>re</code>:</p>
<pre><code>seq = 'M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K'
pos = [2, 5, 10, 11, 14]

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

print(''.join(new_seq))
</code></pre>
<p>Result:</p>
<pre><code>M----ZFLZFLCLZZFAZAQ-------------D-GKCG----TLSNKSPS--------------------K
</code></pre>

CC BY-SA 4.0
2022-07-28T02:11:23.647Z
Thanks @gaspanic. You are right `re` is very slow and removing is of course the better solution altogether. I came to coding via Perl so I tend to over regex and `re` is one the downsides of Python.
M__
CC BY-SA 4.0
2022-07-28T13:40:37.480Z
@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 ```None```s as one loops over ```seq```.
gaspanic
CC BY-SA 4.0
2022-07-28T13:57:45.857Z
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.
M__
CC BY-SA 4.0
2022-07-28T14:23:18.267Z
... okay I tried predefining the list but didn't noticeably improve in my code. deque is another way in. Would make a good question.
M__
CC BY-SA 4.0
2022-07-28T22:45:33.403Z
Yes, looks like you're right. See this thread: https://stackoverflow.com/a/63600055/8861903 . I guess I must have inherited that idea from other languages (or I'm just plain wrong).
gaspanic
CC BY-SA 4.0
2022-07-28T23:46:33.490Z
Thanks @gaspanic those timings were better than mine. I'll take another look later. It could be the data set size, this 'data set' is very small and there could be an overhead.
M__
CC BY-SA 4.0
2022-07-28T23:52:25.677Z
@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 :).
gaspanic
CC BY-SA 4.0
2022-08-03T19:39:36.307Z
Hi @gaspanic as promised I've placed a question here: https://bioinformatics.stackexchange.com/questions/19506/improving-list-speed-glycosylation-example
M__
CC BY-SA 4.0
2022-08-04T11:59:17.377Z
Great, thx for taking the time! :)
gaspanic
CC BY-SA 4.0
2022-08-05T02:33:41.993Z
Hi @gaspanic I've placed a 'final answer' in the thread https://bioinformatics.stackexchange.com/questions/19506/improving-list-speed-glycosylation-example . This combines an `iterator` with a `list comprehension` ... Many thanks for your input.
M__
CC BY-SA 4.0
2022-07-28T00:55:26.400Z
2022-07-28T01:21:56.273Z
2022-07-28T01:21:56.273Z
gaspanic
<p>I'm trying to programmatically replace a set of amino acid residues on an MSA with a &quot;Z&quot; from a list of unaligned positions. Any ideas on how I could do this?</p>
<p>Input: a list of unaligned positions for each sequence + an MSA<br />
Output: an MSA with those unaligned positions on the aligned MSA replaced with a Z for each sequence.</p>
<p>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 &quot;Z&quot; to easily visualize the sites. So essentially if I feed it the following sequence</p>
<pre><code>M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K
</code></pre>
<p>and positions</p>
<pre><code>2, 5, 10, 11, 14
</code></pre>
<p>it should spit back</p>
<pre><code>M----ZFLZFLCLZZFAZAQ-------------D-GKCG----TLSNKSPS--------------------K
</code></pre>

CC BY-SA 4.0
2020-10-23T21:26:43.007Z
Hi @DannyMorales and welcome. This is likely fairly easy to do, but the question is difficult to understand. Any clarity would be helpful.
M__
CC BY-SA 4.0
2020-10-23T22:08:01.593Z
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"
Danny Morales
CC BY-SA 4.0
2020-10-23T22:28:41.497Z
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
M__
CC BY-SA 4.0
2020-10-24T13:48:06.557Z
Could you please show your code, so that is clear what you tried so far?
Mr_Z
<fasta><python><proteins>
Python/Biopython - Replace amino acid residue on MSA with "z" from a list of unaligned positions
# Python/Biopython - Replace amino acid residue on MSA with "z" from a list of unaligned positions

(Asked by: username_0 on 2020-10-23)

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

**Tags:** fasta, python, proteins
> **Comments:**
> **username_1**: Hi @DannyMorales and welcome. This is likely fairly easy to do, but the question is difficult to understand. Any clarity would be helpful.  
> **username_0**: 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"  
> **username_1**: 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  
> **username_2**: Could you please show your code, so that is clear what you tried so far?  

---

## 2 answers

### Answer by username_1 on 2022-07-27. Score: 1

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.


---

### Answer by username_3 on 2022-07-28. Score: 1

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

> **Comments:**
> **username_1**: Thanks @gaspanic. You are right `re` is very slow and removing is of course the better solution altogether. I came to coding via Perl so I tend to over regex and `re` is one the downsides of Python.  
> **username_3**: @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 ```None```s as one loops over ```seq```.  
> **username_1**: 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.  
> **username_1**: ... okay I tried predefining the list but didn't noticeably improve in my code. deque is another way in. Would make a good question.  
> **username_3**: Yes, looks like you're right. See this thread: https://stackoverflow.com/a/63600055/8861903 . I guess I must have inherited that idea from other languages (or I'm just plain wrong).  
> **username_1**: Thanks @gaspanic those timings were better than mine. I'll take another look later. It could be the data set size, this 'data set' is very small and there could be an overhead.  
> **username_3**: @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 :).  
> **username_1**: Hi @gaspanic as promised I've placed a question here: https://bioinformatics.stackexchange.com/questions/19506/improving-list-speed-glycosylation-example  
> **username_3**: Great, thx for taking the time! :)  
> **username_1**: Hi @gaspanic I've placed a 'final answer' in the thread https://bioinformatics.stackexchange.com/questions/19506/improving-list-speed-glycosylation-example . This combines an `iterator` with a `list comprehension` ... Many thanks for your input.  

---

