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