# Exercises - "Introduction to Python" :snake: ## Python - II ## Assignment - Similarity of Sequences Write in an editor the program, which calculates the distance between two sequences.  python3 seq1 = "ACGT" seq2 = "AGGT"  A simple program (without function and modules) is sufficient. :::spoiler Solution  python3 # define sequences seq1 = "ACGT" seq2 = "AGGT" # initiate counter distance_score = 0 # for each letter in the sequences for i in range(len(seq1)): # if they don't match, add a distance point if (seq1[i] != seq2[i]): distance_score += 1 # print result to the terminal print("Distance between A and B: ", distance_score)   Distance between A and B: 1  ::: </br> 1. **Calculate the distance between the following sequences and print out the result**. Since the following sequences are already aligned, we can calculate the distance between them. Change your program so that it can read two aligned sequences from the command line. Test your program with the following sequences.  a) ACGT and A-GT b) AC-GT and AGT-- c) AC-CGT and AGT--- d) ACCGT and TGCCA e) GATT-ACA and TACCATAC f) --GA--TT--AC-A and TA--CC--AT--CA  :::spoiler Solution  python3 import sys # read sequences from command line arguments seq1 = sys.argv[1] seq2 = sys.argv[2] # initiate counter distance_score = 0 # for each letter in the sequences for i in range(len(seq1)): # if they don't match, add a distance point if (seq1[i] != seq2[i]): distance_score += 1 # print result to terminal print("Distance between seq1 and seq2: ", distance_score)   Distance between seq1 and seq2: 1 Distance between seq1 and seq2: 4 Distance between seq1 and seq2: 5 Distance between seq1 and seq2: 4 Distance between seq1 and seq2: 7 Distance between seq1 and seq2: 13  ::: </br> 2. **Extend the program that the aligned sequences are printed out additionally to their distance**. :::spoiler Solution  python3 print("Sequence seq1: ", seq1) print("Sequence seq2: ", seq2) print("Distance between seq1 and seq2: ", distance_score)   Sequence seq1: ACGT Sequence seq2: A-GT Distance between seq1 and seq2: 1 Sequence seq1: AC-GT Sequence seq2: AGT-- Distance between seq1 and seq2: 4  ::: </br> 3. Extend the program that the distance between two sequences is only calculated when both sequences have the same length. Test your program with the input sequences:  a) ACGT and AGT b) ACCGT and TGCCA  :::spoiler Solution  python3 import sys # read sequences from command line arguments seq1 = sys.argv[1] seq2 = sys.argv[2] # check if the sequence lengths match if (len(seq1) == len(seq2)): # initiate the counter distance_score = 0 # for each letter in the sequences for i in range(len(seq1)): if (seq1[i] != seq2[i]): distance_score += 1 # print the result print("Sequence seq1: ", seq1) print("Sequence seq2: ", seq2) print("Distance between seq1 and seq2: ", distance_score) # tell the user the sequence length don't match else: print("Sequences seq1 and seq2 are of different length.")   Sequences seq1 and seq2 are of different length. Sequence seq1: ACCGT Sequence seq2: TGCCA Distance between seq1 and seq2: 4  ::: </br> 4. **Extend the program that the second sequence is inverted and assigned to a third sequence**. Please, read the first and second sequence from the command line. Calculate the distances between the first and the second and between the first and the third sequence. Compare the distance between the first and the second and the first and the third sequence and print the alignment with the smaller distance. If the distances are equal, then print the alignment of the first and second sequence. Test your program with the following sequences:  a) ACGT and A-GT b) AC-GT and AGT-- c) ACCGT and TGCCA d) GATT-ACA and TACCATAC  :::spoiler Solution  python3 import sys # get sequences from the command line arguments seq1 = sys.argv[1] seq2 = sys.argv[2] ### reverse the seq2 string and save as seq2_rev # initiate variable seq2_rev = "" # for each letter in seq2 for i in range(len(seq2)): # add the next character to the reversed seq2 string seq2_rev += seq2[len(seq2) - i - 1] # only run the calculation if the sequences have the same length if (len(seq1) == len(seq2)): # initiate counters dist_1_2 = 0 dist_1_2rev = 0 # for each letter in the sequences for i in range(len(seq1)): # if they don't match, add a distance point if (seq1[i] != seq2[i]): dist_1_2 += 1 # if they don't match, add a distance point if (seq1[i] != seq2_rev[i]): dist_1_2rev += 1 # if the distance seq1seq2 is less or eq to distance seq1seq2_rev if (dist_1_2 <= dist_1_2rev): # print the seq1seq2 sequences and distance score print("Sequence seq1: ", seq1) print("Sequence seq2: ", seq2) print("Distance between seq1 and seq2: ", dist_1_2) # else, if the seq1seq2_rev distance is less than seq1seq2 else: # print the seq1seq2_rev sequences and distance score print("Sequence seq1: ", seq1) print("Sequence seq2_rev: ", seq2_rev) print("Distance between seq1 and seq2_rev: ", dist_1_2rev) # tell the user the lengths differ else: print("Sequences seq1 and seq2 are of different length.")   Sequence seq1: ACGT Sequence seq2: A-GT Distance between seq1 and seq2: 1 Sequence seq1: AC-GT Sequence seq2: AGT-- Distance between seq1 and seq2: 4 # and so on  ::: </br> ## :red_circle: Bonus Exercises ### 1. Functions Open an editor and save your new program. In this program we will create a few functions. 1.1 Define the two functions similarity and distance: \begin{equation*} similarity(a,b) =\left\{ \begin{array}{rl} 1,& \mathrm{if} \quad a = b\\ 0.5,& \mathrm{if} \quad a \neq b, \mbox{a and b are both purines or pyrimidines}\\ 0,& \mathrm{if} \quad a \neq b, \mbox{a and b are not the same kind} \end{array}\right. \end{equation*} \begin{equation*} distance(a,b) =\left\{ \begin{array}{rl} 0,& \mathrm{if} \quad a = b\\ 0.5,& \mathrm{if} \quad a \neq b, \mbox{a and b are both purines or pyrimidines}\\ 1,& \mathrm{if} \quad a \neq b, \mbox{a and b are not the same kind} \end{array}\right. \end{equation*} **Note**: Purines are A and G, pyrimidines are C and T. :::spoiler Solution  python3 # define which bases are purines and pyrimidines pur = ["A", "G"] pyr = ["C", "T"] # define the similarity function for two single bases def similarity(base1, base2): # if they match, return 1 if (base1 == base2): return 1 # else,if they dont match but are of the same kind elif (((base1 in pur) and (base2 in pur)) or ((base1 in pyr) and (base2 in pyr))) return 0.5 # if they neither matches or are of the same kind, return 0 else: return 0 # define the distance function for two single bases def distance(base1, base2): # if they match, return 0 if (base1 == base2): return 0 # else,if they dont match but are of the same kind elif (((base1 in pur) and (base2 in pur)) or ((base1 in pyr) and (base2 in pyr))) return 0.5 # if they neither matches or are of the same kind, return 1 else: return 1  ::: </br> 1.2 **Write two functions sequence_similarity and sequence_distance, which calculates the similarity and distance of two whole sequences**. :::spoiler Solution  python3 # define the similarity function for whole sequences def sequence_similarity (seq1, seq2): # initiate counter similarity_score = 0.0 # go through all bases in seq1 for i in range(len(seq1)): # calculate their similarity and add to the score similarity_score = similarity_score + similarity(seq1[i], seq2[i]) # return the final score return similarity_score # define the distance function for whole sequences def sequence_distance(seq1, seq2): # initiate counter distance_score = 0.0 # go through all bases in seq1 for i in range(len(seq1)): # calculate the distance and add to the score distance_score = distance_score + distance(seq1[i], seq2[i]) # return the final score return distance_score  ::: </br> 1.3 **Calculate the similarity and distance for the following sequences**. Read these sequences from the command line and print out their similarity and distance.  a) ACGT and TGCA b) ATAG and ACAC c) ACGC and ATTT d) AGTT and ACTT e) TCGC and AGAG  :::spoiler Solution  python3 import sys # read the sequences from command line arguments seq1 = sys.argv[1] seq2 = sys.argv[2] # print the similarity and distance print("Similarity: ", sequence_similarity(seq1, seq2)) print("Distance: ", sequence_distance(seq1, seq2))   Similarity: 0.0 Distance: 4.0 Similarity: 2.5 Distance: 1.5 # and so on  ::: </br> ### 2. Modules In this exercise we will write three different programs. 2.1 **Write a new Python file (module) called sequence_tools.py which contain both the two functions similarity and distance as defined previously**. :::spoiler Solution  python3 ######################### ### sequence_tools.py ### ######################### # define which bases are purines and pyrimidines pur = ["A", "G"] pyr = ["C", "T"] # define the similarity function for two single bases def similarity(base1, base2): # if they match, return 1 if (base1 == base2): return 1 # else,if they dont match but are of the same kind elif (((base1 in pur) and (base2 in pur)) or ((base1 in pyr) and (base2 in pyr))) return 0.5 # if they neither matches or are of the same kind, return 0 else: return 0 # define the distance function for two single bases def distance(base1, base2): # if they match, return 0 if (base1 == base2): return 0 # else,if they dont match but are of the same kind elif (((base1 in pur) and (base2 in pur)) or ((base1 in pyr) and (base2 in pyr))) return 0.5 # if they neither matches or are of the same kind, return 1 else: return 1 # define the similarity function for whole sequences def sequence_similarity (seq1, seq2): # initiate counter similarity_score = 0.0 # go through all bases in seq1 for i in range(len(seq1)): # calculate their similarity and add to the score similarity_score = similarity_score + similarity(seq1[i], seq2[i]) # return the final score return similarity_score # define the distance function for whole sequences def sequence_distance(seq1, seq2): # initiate counter distance_score = 0.0 # go through all bases in seq1 for i in range(len(seq1)): # calculate the distance and add to the score distance_score = distance_score + distance(seq1[i], seq2[i]) # return the final score return distance_score  ::: </br> 2.2 **Write another Python file that calculates for each combination of two sequences stored in list seq_list the similarity and distance using the module defined previously**. l = ["ATCCGGT", "GCGTTAC", "CTACTGC", "TTGCAGT", "AGTCACC"] :::spoiler Solution  python3 from sequence_tools import * # define sequences seq_list = ["ATCCGGT", "GCGTTAC", "CTACTGC", "TTGCAGT", "AGTCACC"] # loop over each sequence in seq_list for i in range(len(seq_list)): # loop over the remaining sequences in seq_list for j in range(i+1, len(seq_list)): # calculate the similarity and distance similarity_score = sequence_similarity(seq_list[i], seq_list[j]) distance_score = sequence_distance(seq_list[i], seq_list[j]) # print the result for this comparison print(seq_list[i], seq_list[j], " Similarity: ", similarity_score, " Distance: ", distance_score)   ATCCGGT GCGTTAC Similarity: 2.5 Distance: 4.5 ATCCGGT CTACTGC Similarity: 3.5 Distance: 3.5 ATCCGGT TTGCAGT Similarity: 4.5 Distance: 2.5 ATCCGGT AGTCACC Similarity: 3.5 Distance: 3.5 GCGTTAC CTACTGC Similarity: 4.0 Distance: 3.0 GCGTTAC TTGCAGT Similarity: 3.0 Distance: 4.0 GCGTTAC AGTCACC Similarity: 2.0 Distance: 5.0 CTACTGC TTGCAGT Similarity: 4.5 Distance: 2.5 CTACTGC AGTCACC Similarity: 2.0 Distance: 5.0 TTGCAGT AGTCACC Similarity: 2.5 Distance: 4.5  ::: </br> 2.3 **Extend your program. Determine the combination of sequences with the highest similarity of all sequences stored in list l. Write these two sequences and the alignment into a new file, called similar_sequences.txt.** For example for two given sequences: "ATC" and "ACC" The alignment would be:  ATC | | ACC  And this alignment should be written to a new output file. **Hint**: A line-break in Python can be made by adding ’\n’ to the end of the line. :::spoiler Solution  python3 from sequence_tools import * # define sequences seq_list = ["ATCCGGT", "GCGTTAC", "CTACTGC", "TTGCAGT", "AGTCACC"] # define variables similarity_highscore = 0 best_seq1 = "" best_seq2 = "" # loop over each sequence in seq_list for i in range(len(seq_list)): # compare the sequence to all remaining sequences in seq_list for j in range(i+1, len(seq_list)): # calculate the similarity similarity_score = sequence_similarity(seq_list[i], seq_list[j]) # check if it's a new similarity highscore if (similarity_score > similarity_highscore): # if it is, save this as the new highscore similarity_highscore = similarity_score best_seq1 = seq_list[i] best_seq2 = seq_list[j] # create an empty string to add the alignment to alignment_matches = "" # go through each letter the best aligned pair for i in range(len(best_seq1)): # find places where they match if (best_seq1[i] == best_seq2[i]): alignment_matches = alignment_matches + "|" # and places they don't match else: alignment_matches = alignment_matches + " " # write the sequences and the match symbols to file outfile = open("similar_sequences.txt", "w") outfile.write(best_seq1 + "\n") outfile.write(alignment_matches + "\n") outfile.write(best_seq2 + "\n")   ATCCGGT | | || TTGCAGT  ###### tags: UPPMAX Intro course