def mutate_seqs(seqs, nucleotides="NACTG", mutations=1):
    """Generate all sequences within `mutations` distance from a reference sequence
    
    Parameters
    ----------
    seqs : str or list of str
        Single reference sequence (a string) or a group of strings
                        
    nucleotides : list of char, optional
        Permitted nucleotide substitutions (Default: `'NACTG'`)
    
    mutations : int, optional
        Number of substitutions to make (Default: `1`)
    
    
    Returns
    -------
    set
        all sequences within `mutations` substitutions from the sequence(s)
        specified in `seqs`
    """
    if isinstance(seqs, str):
        seqs = [seqs]
    if mutations == 0:
        return set(seqs)
    else:
        seqsout = []
        for seq in seqs:
            for nuc in nucleotides:
                for i in range(len(seq)):
                    newseq = list(seq)[:]
                    newseq[i] = nuc
                    seqsout.append("".join(newseq))

        seqsout.extend(mutate_seqs(seqsout, nucleotides=nucleotides, mutations=mutations - 1))
        return set(seqsout) | set(seqs)