Introducing a fixed number of random substitutions in a sequence

Question

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.

Answer 1

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)
  

Answer 2

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.

Answer 3

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