Exercises - "Introduction to Python"

Python - II

Assignment - Similarity of Sequences

Write in an editor the program, which calculates the distance between two sequences.

seq1 = "ACGT"
seq2 = "AGGT"

A simple program (without function and modules) is sufficient.

Solution

# 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

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

Solution

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

2. Extend the program that the aligned sequences are printed out additionally to their distance.

Solution

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

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

Solution

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

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

Solution

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

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:

$$similarity(a,b)=\{\begin{array}{rl}1,& \mathrm{if}a=b\\ 0.5,& \mathrm{if}a\ne b,\text{a and b are both purines or pyrimidines}\\ 0,& \mathrm{if}a\ne b,\text{a and b are not the same kind}\end{array}$$

$$distance(a,b)=\{\begin{array}{rl}0,& \mathrm{if}a=b\\ 0.5,& \mathrm{if}a\ne b,\text{a and b are both purines or pyrimidines}\\ 1,& \mathrm{if}a\ne b,\text{a and b are not the same kind}\end{array}$$

Note: Purines are A and G, pyrimidines are C and T.

Solution

# 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

1.2 Write two functions sequence_similarity and sequence_distance, which calculates the similarity and distance of two whole sequences.

Solution

# 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

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

Solution

import sys

### Paste here the code for the functions you wrote in 1.2 and 1.3 ###

# 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

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.

Solution

#########################
### 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

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"]

Solution

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

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.

Solution

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