# Introducing a fixed number of random substitutions in a sequence

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.

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.

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. Plus, your version didn't protect against changing the same position more than once:

    while num_mutations != 0:
index_choice = random_choice(index_choices)


There is nothing there to check if index_choice has been used in a previous iteration of the while loop, so your version could return a sequence with just a single change in the admittedly unlikely case where every index_choice turned out to be the same index. In my version, we ensure the exact number of substitutions.

#!/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)



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)
$$$$
`