<jupyter_start><jupyter_text>DNA mutation simulation via Jukes-Cantor method Below are the main Jukes-Cantor simulation functions. Sequences are represented by a list of values {1, 2, 3, 4} corresponding to bases {A, C, G, T} respectively. Start by generating ancestor sequence A, filled with randomly selected bases with each base being equally likely. Descendent sequences B and C are mutated following the Jules-Cantor model as explained in the assignment brief.<jupyter_code>import numpy

def length_mutation_simulation(sequence_length, t, mu):
    """
    Generates an ancestor DNA sequence of given length and then mutatates two
    other DNA sequences descended from that ancestor using Jukes-Cantor model.
    Sequences are represented as a list of values {1, 2, 3, 4}
    """
    A = generate_new_sequence(sequence_length)
    B = A.copy()
    mutate_sequence(B, t, mu)
    C = A.copy()
    mutate_sequence(C, t, mu)
    return A, B, C
   

def ancestor_mutation_simulation(A, t, mu):
    """
    Mutatates two DNA sequences descended from given ancestor sequence, A, 
    using Jukes-Cantor model. Sequences are represented as a list of values {1, 2, 3, 4}
    """
    B = A.copy()
    mutate_sequence(B, t, mu)
    C = A.copy()
    mutate_sequence(C, t, mu)
    return B, C
    
def generate_new_sequence(sequence_length):
    """
    Randomly generates a DNA sequence represented as a list of values {1, 2, 3, 4}
    of a given length. Probability of a given site being any base is equally likely
    for all bases.
    """
    bases = [1, 2, 3, 4]
    cumulative_probabilities = numpy.cumsum([0.25, 0.25, 0.25, 0.25])
    sequence = []
    for _ in range(sequence_length):
        u = numpy.random.random()
        sequence.append(bases[numpy.argmax(cumulative_probabilities > u)])
    return sequence

def randexp(lambd):
    """
    Generates a random exponential with parameter lambd
    """
    return - numpy.log(numpy.random.random()) / lambd

def mutate_sequence(sequence, T, mu):
    """
    Takes a sequence represented by a list of values {1, 2, 3, 4} and randomly mutates 
    each site in the sequence to another base with rate (3/4)*mu within time T.
    """
    for site in range(len(sequence)):
        t = randexp(3/4*mu)
        while t < T:
            sequence[site] = mutate_base(sequence[site])
            t += randexp(3/4*mu)
            
def mutate_base(base):
    """
    Generates a mutation of a base represented by one of the values {1, 2, 3, 4}
    to one of the other bases with each other base being equally likely.
    """
    cumulative_probabilities = numpy.cumsum([1/3, 1/3, 1/3])
    u = numpy.random.random()
    if base == 1:
        bases = [2, 3, 4]
        base = bases[numpy.argmax(cumulative_probabilities > u)]
    elif base == 2:
        bases = [1, 3, 4]
        base = bases[numpy.argmax(cumulative_probabilities > u)]
    elif base == 3:
        bases = [1, 2, 4]
        base = bases[numpy.argmax(cumulative_probabilities > u)]
    elif base == 4:
        bases = [1, 2, 3]
        base = bases[numpy.argmax(cumulative_probabilities > u)]
    return base

def calculate_differences(object1, object2):
    """
    Calculates difference between two iterable objects of same length
    """
    count = 0
    for i in range(len(object1)):
        if object1[i] != object2[i]:
            count += 1
    return count

def convert_sequence_to_string(sequence):
    """
    Converts a DNA sequence represented as a list of values {1, 2, 3, 4}
    to a string of bases {A, C, G T}
    """
    string = ""
    for value in sequence:
        if value == 1:
            string += "A"
        elif value == 2:
            string += "C"
        elif value == 3:
            string += "G"
        elif value == 4:
            string += "T"
    return string<jupyter_output><empty_output><jupyter_text>Example simulation with sequence length 50, time 10 and mutation parameter 0.01.<jupyter_code>A1, B1, C1 = length_mutation_simulation(50, 10, 0.01)
diff_A_and_B = calculate_differences(A1, B1)
diff_A_and_C = calculate_differences(A1, C1)
diff_B_and_C = calculate_differences(B1, C1)
print("A: {0}\nB: {1}\nC: {2}".format(convert_sequence_to_string(A1), convert_sequence_to_string(B1), \
                                      convert_sequence_to_string(C1)))
print("Differences between B and ancestor A: {0}".format(diff_A_and_B))
print("Differences between C and ancestor A: {0}".format(diff_A_and_C))
print("Differences between siblings B and C: {0}".format(diff_B_and_C))<jupyter_output>A: AATTGCAGTCGCATTCACAACTGTCTCCTACAAACCGTCATGTGATCACT
B: AATTGCAGTCGTATTCACAACTGTCTCCTACAAACAGTCATGTGATCACT
C: AATTGCAGTCTCATTAACAACTGTCTCCTACAAACTGCCATGTAATCACT
Differences between B and ancestor A: 2
Differences between C and ancestor A: 5
Differences between siblings B and C: 6<jupyter_text>Below is an expanded version of the simulation with insertions and deletions modeled, representing gap characters by 5.<jupyter_code>def length_mutation_simulation_id(sequence_length, t, mu):
    """
    Generates an ancestor DNA sequence of given length and then mutatates two
    other DNA sequences descended from that ancestor using Jukes-Cantor model
    Also simulates insertions and deletions.
    Sequences are represented as a list of values {1, 2, 3, 4, 5}
    """
    A = generate_new_sequence(sequence_length)
    B = A.copy()
    mutate_sequence(B, t, mu)
    C = A.copy()
    mutate_sequence(C, t, mu)
    B, C = simulate_insertions(B, C, len(A), t, mu)
    B, C = simulate_deletions(B, C, len(A), t, mu)
    return A, B, C

def ancestor_mutation_simulation_id(A, t, mu):
    """
    Mutatates two DNA sequences descended from given ancestor sequence, A, 
    using Jukes-Cantor model. 
    Also simulates insertions and deletions.
    Sequences are represented as a list of values {1, 2, 3, 4, 5}
    """
    B = A.copy()
    mutate_sequence(B, t, mu)
    C = A.copy()
    mutate_sequence(C, t, mu)
    B, C = simulate_insertions(B, C, len(A), t, mu)
    B, C = simulate_deletions(B, C, len(A), t, mu)
    return B, C

def poisson_process(lambd, T):
    """
    Generates a Poisson random variable with parameter lambda over duration T
    """
    sum = 0
    t = randexp(lambd)
    while t < T:
        sum += t
        t += randexp(lambd)
    return sum

def convert_sequence_to_string(sequence):
    """
    Converts a DNA sequence represented as a list of values {1, 2, 3, 4, 5}
    to a string of bases {A, C, G, T} and gap character -
    """
    string = ""
    for value in sequence:
        if value == 1:
            string += "A"
        elif value == 2:
            string += "C"
        elif value == 3:
            string += "G"
        elif value == 4:
            string += "T"
        elif value == 5:
            string += "-"
    return string

def simulate_insertions(B, C, L, t, mu):
    """
    Simulates insertions on DNA sequences B and C. Sequences are represented as a list of values {1, 2, 3, 4}
    Gap characters, 5, are inserted into other sequence to ensure homologous sites are aligned
    """
    hI = L * t * mu / 10
    B_insertions = poisson_process(hI, 1)
    i = 0
    while B_insertions > i:
        cumulative_probabilities = numpy.cumsum([1/len(B)]*len(B))
        v = numpy.random.random()
        index = numpy.argmax(cumulative_probabilities > v)
        if index + 1 != len(B):
            if B[index] != 5:
                B = B[:index] + generate_new_sequence(3) + B[index-1:]
                C = C[:index] + [5, 5, 5] + C[index-1:]
                i += 1
        else:
            B += generate_new_sequence(3)
            C += [5, 5, 5]
            i += 1
    C_insertions = poisson_process(hI, 1)
    i = 0
    while C_insertions > i:
        cumulative_probabilities = numpy.cumsum([1/len(C)]*len(C))
        v = numpy.random.random()
        index = numpy.argmax(cumulative_probabilities > v)
        if index + 1 != len(C):
            if C[index] != 5:
                C = C[:index + 1] + generate_new_sequence(3) + C[index + 2:]
                B = B[:index + 1] + [5, 5, 5] + B[index + 2:]
                i += 1
        else:
            C += generate_new_sequence(3)
            B += [5, 5, 5]
            i += 1
    return B, C

def simulate_deletions(B, C, L, t, mu):
    """
    Simulates deletions on sequences B and C. Sequences are represented as a list of values {1, 2, 3, 4}
    Gap characters, 5, are inserted into other sequence to ensure homologous sites are aligned
    """
    hD = L * t * mu / 10
    B_deletions = poisson_process(hD, 1)
    i = 0
    while B_deletions > i:
        cumulative_probabilities = numpy.cumsum([1/len(B)]*len(B))
        v = numpy.random.random()
        deletions = 3
        index = numpy.argmax(cumulative_probabilities > v)
        if B[index] != 5:
            while deletions > 0:
                if index == len(B):
                    deletions = 0
                elif B[index] != 5:
                    B[index] = 5
                    deletions -= 1
                index += 1
            i += 1
    C_deletions = poisson_process(hD, 1)
    i = 0
    while C_deletions > i:
        cumulative_probabilities = numpy.cumsum([1/len(C)]*len(C))
        v = numpy.random.random()
        deletions = 3
        index = numpy.argmax(cumulative_probabilities > v)
        if C[index] != 5:
            while deletions > 0:
                if index == len(C):
                    deletions = 0
                elif C[index] != 5:
                    C[index] = 5
                    deletions -= 1
                index += 1
            i += 1
    return B, C<jupyter_output><empty_output><jupyter_text>Example simulation with sequence length 50, time 10 and mutation parameter 0.01 using insertions and deletions.<jupyter_code>A2, B2, C2 = length_mutation_simulation_id(50, 20, 0.01)
print("B: {0}\nC: {1}".format(convert_sequence_to_string(B2), convert_sequence_to_string(C2)))<jupyter_output>B: TCCGCTAGAAGCTGACACCGCTTCGGTAATTGTTTTCATGATAGGTGAGTCGGT
C: TCCACTAGTAGCTGACCCAGAGTCGGGAATTGTATTCATCGT---TGAG---GT<jupyter_text>Below is a global alignment function for sequences. Remove gaps function is provided so sequences simulated with insertions and deletions can be aligned via gloab alignment system in addition to viewing true alignment from simulation.<jupyter_code>def globally_align(A, B, score, gap_penalty):
    """
    Globally aligns sequences A and B.
    Sequences are represented as a list of values {1, 2, 3, 4}
    """
    F = []
    for i in range(len(A)+1):
        F.append([])
        for j in range(len(B)+1):
            F[i].append(0)
    for i in range(len(A)+1):
        F[i][0] = i * gap_penalty
    for j in range(len(B)+1):
        F[0][j] = j * gap_penalty
    for i in range(1, len(A)+1):
        for j in range(1, len(B)+1):
            match = F[i-1][j-1] + score[A[i-1]-1][B[j-1]-1]
            delete = F[i-1][j] + gap_penalty
            insert = F[i][j-1] + gap_penalty
            F[i][j] = max(match, delete, insert)
    A_aligned = []
    B_aligned = []
    i = len(A)
    j = len(B)
    while (i > 0 or j > 0):
        if (i > 0 and j > 0 and F[i][j] == F[i-1][j-1] + score[A[i-1]-1][B[j-1]-1]):
            A_aligned = [A[i-1]] + A_aligned
            B_aligned = [B[j-1]] + B_aligned
            i -= 1
            j -= 1
        elif (i > 0 and F[i][j] == F[i-1][j] + gap_penalty):
            A_aligned = [A[i-1]] + A_aligned
            B_aligned = [5] + B_aligned
            i -= 1
        else:
            A_aligned = [5] + A_aligned
            B_aligned = [B[j-1]] + B_aligned
            j -= 1
    return A_aligned, B_aligned

def remove_gaps(sequence):
    """
    Converts a DNA sequence represented as a list of values {1, 2, 3, 4, 5}
    representing bases {A, C, G, T} and gap character - to an equivalent sequence 
    representation without gap characters.
    """
    new_sequence = []
    for base in sequence:
        if base != 5:
            new_sequence.append(base)
    return new_sequence<jupyter_output><empty_output><jupyter_text>Example simulation of sequences using insertions and deletions followed by global alignment of the first 35 charcters of those sequences.<jupyter_code>A3, B3, C3 = length_mutation_simulation_id(50, 20, 0.01)
print("B: {0}\nC: {1}".format(convert_sequence_to_string(B3), convert_sequence_to_string(C3)))

score_matrix = [[2, -2, -2, -2], [-2, 2, -2, -2], [-2, -2, 2, -2], [-2, -2, -2, 2]]
B3 = remove_gaps(B3)[:35]
C3 = remove_gaps(C3)[:35]
B3_aligned, C3_aligned = globally_align(B3, C3, score_matrix, -3)
print("Gap Penalty -3\nB aligned: {0}\nC aligned: {1}".format(convert_sequence_to_string(B3_aligned), \
                                                              convert_sequence_to_string(C3_aligned)))<jupyter_output>B: TCTGTC---TGCTTTGTCCCTAAAACGAAACCATTCCTTGGTCGCTGATAAC---C
C: TC---CGTTTGCTATGTCCTCAACACGAAAGCATTCCTTCGTCGTTGATAACTGTA
Gap Penalty -3
B aligned: TCTGTCTGCTTTGTCCCTAAAACGAAACCATTCCT
C aligned: TCCGTTTGCTATGTCCTCAACACGAAAGCATTCCT
{'fraction_non_alphanumeric': 0.11300738007380073, 'fraction_numerical': 0.025907134071340713, 'mean_word_length': 2.507414397411701, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 42, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '276394', 'license': 'Apache-2.0', 'n_tokens_mistral': 4655, 'n_tokens_neox': 4318, 'n_words': 1324, 'path': 'DNA Mutation Simulation.ipynb', 'repo_name': 'MoanaStirling/DNA-Mutation-Simulation', 'size': '13008'}