<p>Here's a slightly different approach that should be more efficient. The core of the idea is to have two functions, one which changes a given nucleotide to something other than what it was, and another to perform the substitution on the input sequence. This way, you only need to build the sequence once.</p>
<pre class="lang-py prettyprint-override"><code>#!/usr/bin/env python

import random

def change_to_not_self(old, nts):
    if old not in nts:
        raise ValueError(
            f&quot;{old} is not a valid nucleotide&quot;
        )
    ## Pick the new nucleotide from the list of acceptable
    ## values minus whatever letter we have passed as 'old'.
    new = random.choice(list(set(nts) - set(old)))
    return new


def randomly_substitute(seq: str, num_mutations: int) -&gt; str:
    if num_mutations &gt; len(seq):
        raise ValueError(
            f&quot;{num_mutations} is greater than the length of the input sequence&quot;
        )
    NUCLEOTIDES = set([&quot;A&quot;, &quot;C&quot;, &quot;G&quot;, &quot;T&quot;])
    seq = list(seq)
    indices = random.sample(list(range(len(seq))), num_mutations)
    
    for index in indices:
        seq[index] = change_to_not_self(seq[index], NUCLEOTIDES)
    return ''.join(seq)
  
sequence = 'ACTGACGGTACGTACAACG'
mutated_seq = randomly_substitute(sequence, 5);
print(mutated_seq)
  
</code></pre>

CC BY-SA 4.0
2024-08-26T13:50:53.003Z
Nice! 
Although when you say `There is nothing there to check if index_choice has been used in a previous iteration of the while loop,`, that's not correct IMO, as that's what line `index_choices = [el for el in index_choices if el != index_choice]` does. Your solution is clearly an improvement though, so approving it :))
bricoletc
CC BY-SA 4.0
2024-08-26T16:14:13.260Z
@bricoletc thanks :) My point was that you were setting `index_choice = random_choice(index_choices)` and only then setting `index_choices = [el for el in index_choices if el != index_choice]`, so in the unlikely case where 3 iterations of the loop all chose the same index in `random_choice(index_choices)` you could have changed the same position 3 times. So `chosen_nucleotide` could be repeated.
terdon
CC BY-SA 4.0
2024-08-26T19:27:21.583Z
How? We start with all possible indices, say `[0, 1, 2, 3]`. Then we iterate `num_mutations` time a process where we select one index, say `0`, remove it, so `index_choices` is `[1,2,3]`, then we select another index, say `2`, and `index_choices` becomes `[1,3]`, etc. No?
bricoletc
CC BY-SA 4.0
2024-08-26T22:23:40.537Z
Ah, duh, of course! You're absolutely right, @bricoletc, my bad.
terdon
CC BY-SA 4.0
2024-07-17T15:09:35.453Z
2024-08-26T22:24:00.387Z
2024-08-26T22:24:00.387Z
terdon
<blockquote>
<p>In particular, the function maintains a list called index_choices that is of size len(seq), and keeps re-building this list after each introduced substitution. This feels not very efficient.</p>
</blockquote>
<p>Yes I agree with this. Personally, I'd simply use,</p>
<pre><code> index_choices = random.randint(1, len(seq))
</code></pre>
<p>Lists are expensive in Python and you're creating a list for every value of the sequence when <code>randint</code> is quicker and <code>random.choice</code> is if the numbers are not sequential (there's absent values). Here you just want a rectangular distribution from a straight numerical range.</p>
<p><code>mutated_seq</code> is perfectly good. A <code>numpy</code> solution is would speed this up BTW. I can see where you've got the idea of <code>random.choice</code> from, but its not needed here.</p>

CC BY-SA 4.0
2024-08-26T13:45:16.540Z
But, `index_choices` is a list of the remaining possible choices; if `num_mutations` > 1, then I can't choose any int between 1 and `len(seq)`, because that might choose the same position twice
bricoletc
CC BY-SA 4.0
2024-07-17T15:21:22.330Z
2024-07-17T15:27:09.700Z
2024-07-17T15:27:09.700Z
M__
<p>Strings in Python are immutable (can't be changed), so if you want to do something like introducing random mutations it may be better to use an intermediate mutable object like a list. Here's one way to do it via <code>sample</code> and <code>choice</code>, which will raise an error if the number of samples desired is greater than the sequence length:</p>
<pre><code>#!/usr/bin/env python
from random import sample, choice
OTHER_NUCLEOTIDES = {&quot;A&quot;: [&quot;C&quot;,&quot;G&quot;,&quot;T&quot;],
                     &quot;C&quot;: [&quot;A&quot;,&quot;G&quot;,&quot;T&quot;],
                     &quot;G&quot;: [&quot;A&quot;,&quot;C&quot;,&quot;T&quot;],
                     &quot;T&quot;: [&quot;A&quot;,&quot;C&quot;,&quot;G&quot;]}

def randomly_substitute(seq: str, num_mutations: int = 1) -&gt; str:
    seq = list(seq)
    for pos in sample(range(len(a)), num_mutations):
        seq[pos] = choice(OTHER_NUCLEOTIDES[seq[pos]])
    return ''.join(seq)
<span class="math-container">```</span>
</code></pre>

CC BY-SA 4.0
2024-07-17T22:33:46.777Z
2024-07-17T22:33:46.777Z
gringer
<p>I'm writing a function that introduces <span class="math-container">$n$</span> substitutions inside a sequence of nucleotides.
I have a working version, but I'm looking for any other ways of doing this:</p>
<pre class="lang-py prettyprint-override"><code>from random import choice as random_choice
NUCLEOTIDES = set([&quot;A&quot;, &quot;C&quot;, &quot;G&quot;, &quot;T&quot;])
def randomly_substitute(seq: str, num_mutations: int = 1) -&gt; str:
    if num_mutations &gt; len(seq):
        raise ValueError(
            f&quot;{num_mutations} is greater than the length of the input sequence&quot;
        )
    index_choices = list(range(len(seq)))
    mutated_seq = seq
    while num_mutations != 0:
        index_choice = random_choice(index_choices)
        index_choices = [el for el in index_choices if el != index_choice]
        chosen_nucleotide = seq[index_choice]
        if chosen_nucleotide not in NUCLEOTIDES:
            raise ValueError(&quot;Not a nucleotide: {chosen_nucleotide}&quot;)
        possible_mutations = [el for el in NUCLEOTIDES if el != chosen_nucleotide]
        chosen_mutation = random_choice(possible_mutations)
        mutated_seq = (
            mutated_seq[:index_choice]
            + chosen_mutation
            + mutated_seq[index_choice + 1 :]
        )
        num_mutations -= 1
    return mutated_seq
</code></pre>
<p>In particular, the function maintains a list called <code>index_choices</code> that is of size <code>len(seq)</code>, and keeps re-building this list after each introduced substitution. This feels not very efficient.</p>

<python><mutations><sequence><benchmarking>
Introducing a fixed number of random substitutions in a sequence
# Introducing a fixed number of random substitutions in a sequence

(Asked by: username_0 on 2024-07-17)

I'm writing a function that introduces $n$ substitutions inside a sequence of nucleotides.
I have a working version, but I'm looking for any other ways of doing this:
from random import choice as random_choice
NUCLEOTIDES = set(["A", "C", "G", "T"])
def randomly_substitute(seq: str, num_mutations: int = 1) -> str:
    if num_mutations > len(seq):
        raise ValueError(
            f"{num_mutations} is greater than the length of the input sequence"
        )
    index_choices = list(range(len(seq)))
    mutated_seq = seq
    while num_mutations != 0:
        index_choice = random_choice(index_choices)
        index_choices = [el for el in index_choices if el != index_choice]
        chosen_nucleotide = seq[index_choice]
        if chosen_nucleotide not in NUCLEOTIDES:
            raise ValueError("Not a nucleotide: {chosen_nucleotide}")
        possible_mutations = [el for el in NUCLEOTIDES if el != chosen_nucleotide]
        chosen_mutation = random_choice(possible_mutations)
        mutated_seq = (
            mutated_seq[:index_choice]
            + chosen_mutation
            + mutated_seq[index_choice + 1 :]
        )
        num_mutations -= 1
    return mutated_seq

In particular, the function maintains a list called index_choices that is of size len(seq), and keeps re-building this list after each introduced substitution. This feels not very efficient.

**Tags:** python, mutations, sequence, benchmarking

---

## 3 answers

### Answer by username_1 on 2024-07-17. Score: 2 | Accepted answer

Here's a slightly different approach that should be more efficient. The core of the idea is to have two functions, one which changes a given nucleotide to something other than what it was, and another to perform the substitution on the input sequence. This way, you only need to build the sequence once.
#!/usr/bin/env python

import random

def change_to_not_self(old, nts):
    if old not in nts:
        raise ValueError(
            f"{old} is not a valid nucleotide"
        )
    ## Pick the new nucleotide from the list of acceptable
    ## values minus whatever letter we have passed as 'old'.
    new = random.choice(list(set(nts) - set(old)))
    return new


def randomly_substitute(seq: str, num_mutations: int) -> str:
    if num_mutations > len(seq):
        raise ValueError(
            f"{num_mutations} is greater than the length of the input sequence"
        )
    NUCLEOTIDES = set(["A", "C", "G", "T"])
    seq = list(seq)
    indices = random.sample(list(range(len(seq))), num_mutations)
    
    for index in indices:
        seq[index] = change_to_not_self(seq[index], NUCLEOTIDES)
    return ''.join(seq)
  
sequence = 'ACTGACGGTACGTACAACG'
mutated_seq = randomly_substitute(sequence, 5);
print(mutated_seq)

> **Comments:**
> **username_0**: Nice!  Although when you say `There is nothing there to check if index_choice has been used in a previous iteration of the while loop,`, that's not correct IMO, as that's what line `index_choices = [el for el in index_choices if el != index_choice]` does. Your solution is clearly an improvement though, so approving it :))  
> **username_1**: @bricoletc thanks :) My point was that you were setting `index_choice = random_choice(index_choices)` and only then setting `index_choices = [el for el in index_choices if el != index_choice]`, so in the unlikely case where 3 iterations of the loop all chose the same index in `random_choice(index_choices)` you could have changed the same position 3 times. So `chosen_nucleotide` could be repeated.  
> **username_0**: How? We start with all possible indices, say `[0, 1, 2, 3]`. Then we iterate `num_mutations` time a process where we select one index, say `0`, remove it, so `index_choices` is `[1,2,3]`, then we select another index, say `2`, and `index_choices` becomes `[1,3]`, etc. No?  
> **username_1**: Ah, duh, of course! You're absolutely right, @bricoletc, my bad.  

---

### Answer by username_2 on 2024-07-17. Score: 1

In particular, the function maintains a list called index_choices that is of size len(seq), and keeps re-building this list after each introduced substitution. This feels not very efficient.

Yes I agree with this. Personally, I'd simply use,
 index_choices = random.randint(1, len(seq))

Lists are expensive in Python and you're creating a list for every value of the sequence when randint is quicker and random.choice is if the numbers are not sequential (there's absent values). Here you just want a rectangular distribution from a straight numerical range.
mutated_seq is perfectly good. A numpy solution is would speed this up BTW. I can see where you've got the idea of random.choice from, but its not needed here.

> **Comments:**
> **username_0**: But, `index_choices` is a list of the remaining possible choices; if `num_mutations` > 1, then I can't choose any int between 1 and `len(seq)`, because that might choose the same position twice  

---

### Answer by username_3 on 2024-07-17. Score: 0

Strings in Python are immutable (can't be changed), so if you want to do something like introducing random mutations it may be better to use an intermediate mutable object like a list. Here's one way to do it via sample and choice, which will raise an error if the number of samples desired is greater than the sequence length:
#!/usr/bin/env python
from random import sample, choice
OTHER_NUCLEOTIDES = {"A": ["C","G","T"],
                     "C": ["A","G","T"],
                     "G": ["A","C","T"],
                     "T": ["A","C","G"]}

def randomly_substitute(seq: str, num_mutations: int = 1) -> str:
    seq = list(seq)
    for pos in sample(range(len(a)), num_mutations):
        seq[pos] = choice(OTHER_NUCLEOTIDES[seq[pos]])
    return ''.join(seq)
```


---

