Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-08-26T13:50:53.003Z
Answers.Text: 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 :))
Answers.UserDisplayName: bricoletc
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-08-26T16:14:13.260Z
Answers.Text: @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.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-08-26T19:27:21.583Z
Answers.Text: 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?
Answers.UserDisplayName: bricoletc
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-08-26T22:23:40.537Z
Answers.Text: Ah, duh, of course! You're absolutely right, @bricoletc, my bad.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-07-17T15:09:35.453Z
Answers.LastActivityDate: 2024-08-26T22:24:00.387Z
Answers.LastEditDate: 2024-08-26T22:24:00.387Z
Answers.UserDisplayName: terdon
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-08-26T13:45:16.540Z
Answers.Text: 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
Answers.UserDisplayName: bricoletc
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-07-17T15:21:22.330Z
Answers.LastActivityDate: 2024-07-17T15:27:09.700Z
Answers.LastEditDate: 2024-07-17T15:27:09.700Z
Answers.UserDisplayName: M__
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-07-17T22:33:46.777Z
Answers.LastActivityDate: 2024-07-17T22:33:46.777Z
Answers.UserDisplayName: gringer
Body: <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>
  
Tags: <python><mutations><sequence><benchmarking>
Title: Introducing a fixed number of random substitutions in a sequence
ThreadText: # 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:
  
  ```py
  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.
  
  ```py
  #!/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**: @username_0 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, @username_0, 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)
  ```
  
  ```
  
  
  ---
  
  