<jupyter_start><jupyter_text>Package installsIf you are using jupyter lab online, all packages will be available. If you are running this on your local computer, you may need to install some packages. Run the cell below if using jupyter lab locally.<jupyter_code>!pip install numpy
!pip install scipy
!pip install pandas
!pip install scikit-learn
!pip install seaborn<jupyter_output>Requirement already satisfied: numpy in /home/hp/anaconda3/lib/python3.7/site-packages (1.19.1)
Requirement already satisfied: scipy in /home/hp/anaconda3/lib/python3.7/site-packages (1.5.2)
Requirement already satisfied: numpy>=1.14.5 in /home/hp/anaconda3/lib/python3.7/site-packages (from scipy) (1.19.1)
Requirement already satisfied: pandas in /home/hp/anaconda3/lib/python3.7/site-packages (1.1.1)
Requirement already satisfied: pytz>=2017.2 in /home/hp/anaconda3/lib/python3.7/site-packages (from pandas) (2020.1)
Requirement already satisfied: python-dateutil>=2.7.3 in /home/hp/anaconda3/lib/python3.7/site-packages (from pandas) (2.8.1)
Requirement already satisfied: numpy>=1.15.4 in /home/hp/anaconda3/lib/python3.7/site-packages (from pandas) (1.19.1)
Requirement already satisfied: six>=1.5 in /home/hp/anaconda3/lib/python3.7/site-packages (from python-dateutil>=2.7.3->pandas) (1.15.0)
Requirement already satisfied: scikit-learn in /home/hp/anaconda3/lib/python3.7/site-packages (0.2[...]<jupyter_text>Importing dataTo begin, we need to understand the data.The ribosome genes are available in a .fasta file called 'ribosome_genes.fasta'. You can have a look if you like.These genes will be imported as classes (RibosomeGene). Each RibosomeGene object has a name, accession, sequence and length.You can access these properties using '.' syntax. (see below). Try to think of each gene as something physical, rather than code.In real life, each gene has a .length, its organism has a .name, and it has a .sequence. We can write code in this way too. we will import these into **ribosome_genes**, which is a list of our genes.<jupyter_code>import warnings
warnings.filterwarnings('ignore')

from utilities import import_16s_sequences
ribosome_genes = import_16s_sequences()
print('{:<20}{:<30}{:<15}{:<10}'.format('gene.accession', 'gene.name', 'gene.length', 'gene.sequence'))
for gene in ribosome_genes:
    print('{:<20}{:<30}{:<15}{:<10}'.format(gene.accession, gene.name[:27], gene.length, gene.sequence[:8] + '...'))<jupyter_output>gene.accession      gene.name                     gene.length    gene.sequence
EC_A                Escherichia coli DP101        100            CCUUCGGG...
EC_B                Escherichia coli BEBOP        100            CCUUCGGG...
EC_C                Escherichia coli SAMURAI_CH   100            CCUUCGGG...
EC_D                Escherichia coli SAO          100            CCUUCGGG...
EC_E                Escherichia coli code_GEASS   100            CCUUCGGG...
EC_F                Escherichia coli DEATHNOTE    100            CCUUCGGG...
EC_G                Escherichia coli STEINS_GAT   100            CCUUCGGG...
EF_A                Enterococcus faecalis ATCC    100            CCUUUUUG...
EF_B                Enterococcus faecalis EnGen   100            CCUUUUUG...
EF_C                Enterococcus faecalis ATCC    100            CCUUUUUG...
EF_D                Enterococcus faecalis GAH94   100            CCUUUUUG...
HS_A                Homo sapiens (barack obamas   100            CCCGCCGG.[...]<jupyter_text>SECTION 1: PAIRWISE DISTANCESTo be able to compare organisms via their sequences, we need a way to measure their difference as a distance.  **K-mer distance**The kmer distance between two sequences is defined here as the total number of k-mers that are unique to either sequence.eg: If seq1 has 3 unique kmers not found in seq2 (copy number difference also matters), and seq2 has 2 unique kmers, the kmer distance is 5.<jupyter_code>def create_kmer_dictionary(seq, k):
    kmer_dict = {}
    num = len(seq) - k + 1
    for i in range(num):
        kmer = seq[i:i+k]
        if kmer not in kmer_dict:
            kmer_dict[kmer] = 0
        kmer_dict[kmer] += 1
    return kmer_dict
create_kmer_dictionary('CCUUCGGG', 2)
def calculate_total_unique_kmers(kmers1, kmers2):
    unique_kmers = 0
    c1 = c2 = c3 = c4 = c5 = 0
    for k in kmers1:
        if k not in kmers2:
            c1 += kmers1[k]
        elif k in kmers1 and k in kmers2:
            c2 = c2+ abs(kmers1[k] - kmers2[k])
            
    for k2 in kmers2:
        if k2 not in kmers1:
            c3 += kmers2[k2]
    unique_kmers = c1+c2+c3
            
    return unique_kmers
kmers1 = {'CCUUCGGG':1}
kmers2 = {'CCUUUUUG':2}
calculate_total_unique_kmers(kmers1, kmers2)
def kmer_distance(seq1, seq2, k): 
    kmers1 = create_kmer_dictionary(seq1,k)
    kmers2 = create_kmer_dictionary(seq2,k)
    distance = calculate_total_unique_kmers(kmers1, kmers2)
    return distance<jupyter_output><empty_output><jupyter_text>Let's check our function. We can use the first two entries in the 'ribosome_genes' list.If implemented correctly, the following should return 24<jupyter_code>distance = kmer_distance(ribosome_genes[1].sequence, ribosome_genes[3].sequence, 8)
print(distance)<jupyter_output>24<jupyter_text>**smith-waterman alignment**Another way to compare the similarity of two sequences is through alignment.The alignment score of two sequences will be high when they are similar, and low when they are distinct.  Keep in mind the matrix must be 1 element larger than the sequence lengths. Consider whether indel scores for the first row and column need to be filled in.<jupyter_code>import numpy as np
def init_scoregrid(seq1, seq2, indel_score=-4):
    rs = len(seq1) +1
    cs = len(seq2) +1
    scoregrid = np.zeros((rs, cs), np.int)
    return scoregrid<jupyter_output><empty_output><jupyter_text>Let's do a sanity check that the grid has been initialised properly. The following should print the initialised scoregrid<jupyter_code>print(init_scoregrid('hello', 'kittycat'))<jupyter_output>[[0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]]<jupyter_text>Write a function that calculates the initialised scoregrid. It accepts two sequences, a scoregrid and match/mismatch and indel scores.<jupyter_code>import itertools
def calculate_scoregrid(seq1, seq2, scoregrid, match_score=1, mismatch_score=-4, indel_score=-4):
    for i, j in itertools.product(range(1, scoregrid.shape[0]), range(1, scoregrid.shape[1])):
        match = scoregrid[i - 1, j - 1] + (match_score if seq1[i - 1] == seq2[j - 1] else + mismatch_score)
        delete = scoregrid[i - 1, j] + indel_score
        insert = scoregrid[i, j - 1] + indel_score
        scoregrid[i, j] = max(match, delete, insert, 0)

    return scoregrid<jupyter_output><empty_output><jupyter_text>Let's do another sanity check. The following should print a calculated scoregrid, with the these numbers in the bottom right corner: 2  0 0  3<jupyter_code>scoregrid = init_scoregrid('hello', 'helllo')
print(calculate_scoregrid('hello', 'helllo', scoregrid))
def report_alignment_score(scoregrid):
    # given a completed scoregrid, return the smith-waterman alignment score.
    sw_alignment_score = scoregrid.max()
    return sw_alignment_score<jupyter_output><empty_output><jupyter_text>Final sanity check. Should return 4.<jupyter_code>scoregrid = init_scoregrid('hello', 'helllo')
calculated_scoregrid = calculate_scoregrid('hello', 'helllo', scoregrid)
print(report_alignment_score(calculated_scoregrid))<jupyter_output>4<jupyter_text>Ok! now we're ready to put it all together. Fill in the function below with the three functions you wrote to calculate the alignment score of two sequences<jupyter_code>def smith_waterman(seq1, seq2):
    matrix = init_scoregrid(seq1, seq2, indel_score=-4)
    element_scores = calculate_scoregrid(seq1, seq2, scoregrid, match_score=1, mismatch_score=-4, indel_score=-4)
    alignment_score = report_alignment_score(scoregrid)
    return alignment_score<jupyter_output><empty_output><jupyter_text>The following should print 4<jupyter_code>print(smith_waterman('hello', 'helllo'))<jupyter_output>4<jupyter_text>**pairwise distances** We have now written two functions which can calculate the distance of two sequences. We can calculate the k-mer distance, and the smith-waterman alignment score. lets use these two methods to calculate the pairwise distance of our genes.<jupyter_code>import numpy as np
def init_distance_matrix(genes):
    values=[]
    for gene in genes:
        s=gene.accession
        values.append(s)
        values.append(0)
    
    distance_matrix = {}
    for gene in ribosome_genes:
        key = gene.accession        
        distance_matrix[key]={values[i]: values[i + 1] for i in range(0, len(values), 2)} 
    return distance_matrix<jupyter_output><empty_output><jupyter_text>Let's print the distance matrix to make sure it worked.<jupyter_code>from utilities import print_distance_matrix
distance_matrix = init_distance_matrix(ribosome_genes)
print_distance_matrix(distance_matrix)<jupyter_output>EC_A  EC_B  EC_C  EC_D  EC_E  EC_F  EC_G  EF_A  EF_B  EF_C  EF_D  HS_A  MI_A  SD_A  SF_A  SS_A  SS_B
EC_A     0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
EC_B     0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
EC_C     0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
EC_D     0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
EC_E     0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
EC_F     0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
EC_G     0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
EF_A     0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
EF_B     0.0   0.0 [...]<jupyter_text>Time to fill in the matrix with distances. Write a function which calculates the pairwise distance of genes using kmer distance.you will need to call the 'kmer_distance' function you have written above.<jupyter_code>def calculate_kmer_distance_matrix(genes, matrix, k):
    for gene1 in genes:
        key1=gene1.accession
        for gene2 in genes:
            key2=gene2.accession
            matrix[key1][key2]=kmer_distance(gene1.sequence,gene2.sequence,k)

    return matrix<jupyter_output><empty_output><jupyter_text>Let's do the same as above, but this time use the 'smith_waterman' alignment distance function you wrote.<jupyter_code>def calculate_sw_alignment_distance_matrix(genes, matrix): 
    for gene1 in genes:
        key1=gene1.accession
        for gene2 in genes:
            key2=gene2.accession
            matrix[key1][key2]=smith_waterman(gene1.sequence,gene2.sequence)
    return matrix<jupyter_output><empty_output><jupyter_text>Let's test them out. The two cells below will use your calculate_kmer_distance_matrix, and calculate_sw_alignment_distance_matrix functions to add distances to the matrix. **NOTE:** the smith-waterman distance calculations can take time. Give it a minute.<jupyter_code>distance_matrix = init_distance_matrix(ribosome_genes)
kmer_distance_matrix = calculate_kmer_distance_matrix(ribosome_genes, distance_matrix, 8)
print('\nkmer distance matrix')
print_distance_matrix(kmer_distance_matrix)
distance_matrix = init_distance_matrix(ribosome_genes)
sw_alignment_distance_matrix = calculate_sw_alignment_distance_matrix(ribosome_genes, distance_matrix)
print('\nsmith waterman alignment score distance matrix')
print_distance_matrix(sw_alignment_distance_matrix)<jupyter_output>smith waterman alignment score distance matrix
        EC_A  EC_B  EC_C  EC_D  EC_E  EC_F  EC_G  EF_A  EF_B  EF_C  EF_D  HS_A  MI_A  SD_A  SF_A  SS_A  SS_B
EC_A     5.0   5.0   5.0   5.0   5.0   5.0   5.0   4.0   4.0   4.0   4.0   2.0   5.0   5.0   5.0   5.0   4.0
EC_B     5.0   5.0   5.0   5.0   5.0   5.0   5.0   4.0   4.0   4.0   4.0   2.0   5.0   5.0   5.0   5.0   4.0
EC_C     5.0   5.0   5.0   5.0   5.0   5.0   5.0   4.0   4.0   4.0   4.0   2.0   5.0   5.0   5.0   5.0   4.0
EC_D     5.0   5.0   5.0   5.0   5.0   5.0   5.0   4.0   4.0   4.0   4.0   2.0   5.0   5.0   5.0   5.0   4.0
EC_E     5.0   5.0   5.0   5.0   5.0   5.0   5.0   4.0   4.0   4.0   4.0   2.0   5.0   5.0   5.0   5.0   4.0
EC_F     5.0   5.0   5.0   5.0   5.0   5.0   5.0   4.0   4.0   4.0   4.0   2.0   5.0   5.0   5.0   5.0   4.0
EC_G     5.0   5.0   5.0   5.0   5.0   5.0   5.0   4.0   4.0   4.0   4.0   2.0   5.0   5.0   5.0   5.0   4.0
EF_A     4.0   4.0   4.0   4.0   4.0   4.0   4.0   5.0   5.0   5.0   5.0   2.0  [...]<jupyter_text>Let's visualise those in a better manner for human eyes. The cell below will plot heatmaps instead of raw numbers.<jupyter_code>from utilities import heatmap
heatmap(kmer_distance_matrix, sw_alignment_distance_matrix)<jupyter_output><empty_output><jupyter_text>SECTION 2: CLUSTERING From the heatmaps, it seems like there are a few clusters in the data. First, lets convert the pairwise distances to 2D coordinates. This is possible using Multidimensional scaling (MDS).After we have transformed the distance matrix to 2D coordinates, we can plot it to see if any clusters are evident.<jupyter_code>from utilities import mds_scatterplot, distance_matrix_to_coordinates_MDS
kmer_distances_xy = distance_matrix_to_coordinates_MDS(kmer_distance_matrix)
sw_distances_xy = distance_matrix_to_coordinates_MDS(sw_alignment_distance_matrix)
mds_scatterplot(kmer_distances_xy)
mds_scatterplot(sw_distances_xy)<jupyter_output><empty_output><jupyter_text>Seems like there is some clustering happening. Let's use some clustering algorithms to define the clusters. in this manner, we can have an objective way to talk about the patterns in the data. Let's implement the k-means algorithm.<jupyter_code>from utilities import initialise_centroids, average_point, assign_points, plot_kmeans, points_equal, euclidean_distance
def calculate_mean_centroids(data, assignments, k):
    centroids = []
    for cluster in range(k):
        points = [point for point, assignment in zip(data, assignments) if assignment == cluster]
        centroids.append(average_point(points))
    return centroids<jupyter_output><empty_output><jupyter_text>Place calculate_mean_centroids() in the kmeans function below to complete kmeans<jupyter_code>def kmeans(data, k):
    centroids=initialise_centroids(data,k)    
    cluster_assignments=assign_points(centroids,data)
    centroids=calculate_mean_centroids(data,cluster_assignments,k)
    return centroids, cluster_assignments<jupyter_output><empty_output><jupyter_text>You can check your implementation using the cell below:<jupyter_code>centroids, cluster_assignments = kmeans(kmer_distances_xy, 3)
plot_kmeans(kmer_distances_xy, centroids, cluster_assignments, 3)<jupyter_output><empty_output><jupyter_text>Let's also implement k-medoids while we're at it. The only difference between k-means and k-medoids is the calculate_mean_centroids() step, which will instead be calculate_median_centroids()the median can be taken here as the point in the cluster which has smallest cumulative distance to the other points in the clusterYou can use the provided euclidean_distance() function to calculate distances between pointswrite a function which calculates new centroid locations (using the median)<jupyter_code>def calculate_median_centroids(data, assignments, k):
    centroids = []
    for cluster in range(k):
        points = [point for point, assignment in zip(data, assignments) if assignment == cluster]
        centroids.append(tuple(np.median(np.array(points), axis=0)))
    return centroids<jupyter_output><empty_output><jupyter_text>Place calculate_median_centroids() in the kmedoids function below to complete kmedoids<jupyter_code>def kmedoids(data, k):
    centroids=initialise_centroids(data,k)    
    cluster_assignments=assign_points(centroids,data)
    centroids=calculate_median_centroids(data,cluster_assignments,k)
    return centroids, cluster_assignments<jupyter_output><empty_output><jupyter_text>Here is another check cell, for kmedoids this time:<jupyter_code>centroids,cluster_assignments = kmedoids(kmer_distances_xy, 3)
plot_kmeans(kmer_distances_xy, centroids, cluster_assignments, 3)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08955133631636232, 'fraction_numerical': 0.052847986455435966, 'mean_word_length': 3.300312012480499, 'pattern_counts': {'":': 0, '<': 97, '<?xml version=': 0, '>': 96, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '13281', 'license': 'MIT', 'n_tokens_mistral': 6150, 'n_tokens_neox': 5801, 'n_words': 1680, 'path': 'Machine Learning/kmeans from scratch/Kmeans.ipynb', 'repo_name': 'iqbalikku/Hacktoberfest2020-Expert', 'size': '16538'}