Global Alignment with Scoring Matrix

In this problem we are asked to find the alignment score of two protein sequences using the scoring matrix BLOSUM62 and a linear gap penalty of 5. I decided to use Biopython's pairwise2 for this, as you can see in the code below and on my GitHub.

from Bio import pairwise2

from Bio import SeqIO

from Bio.SubsMat.MatrixInfo import blosum62

seqs = []

with open('rosalind_glob.fasta', 'r') as f:

    for record in SeqIO.parse(f, 'fasta'):

        seqs.append(record.seq)

s = seqs[0]

t = seqs[1]

alignments = pairwise2.align.globalds(s, t, blosum62, -5, -5)

print(pairwise2.format_alignment(*alignments[0]))

Note that the gap open penalty and the gap extension penalty are both set to -5, as we are told in the problem description to use a linear gap penalty of 5. 

Edit Distance Alignment

This problem builds on a previous one called Edit Distance, for which you can find a solution here. In that problem we were asked to find the minimum number of changes needed to transform a string s into a string t. In this problem, we are asked once again to find this number (the edit distance), but this time we should also print the two augmented strings representing an optimal alignment of s and t.

To do this I went back to my solution of the previous problem and added some code to pick out one of the possible optimal alignments from the matrix formed by the existing program. The code is a bit messy, but for those who are interested it it is available on my GitHub.

Edit Distance

In this problem we are looking for the Levenshtein distance between two protein sequences. The distance is calculated as the number of single character changes (insertion, deletion or substitution) needed to turn one sequence into the other. One way to do this is described in the Wikipedia page and involves building a matrix of distances between the prefixes of the two strings. The last value to be calculated (in the bottom right position of the matrix) then equals the distance between the two strings.

Sample Dataset
>Rosalind_39
PLEASANTLY
>Rosalind_11
MEANLY

Expected Output
5

The following python3 code initiates a matrix M with len(s) rows and len(t) columns. It then fills the matrix based on if and what type of change is needed between the two strings. The last value is then printed as the answer.

from Bio import SeqIO
seqs = []
with open('sampledata.fasta', 'r') as f:
    for record in SeqIO.parse(f, 'fasta'):
        seqs.append(str(record.seq))
s = seqs[0]
t = seqs[1]
M = [[0 for j in range(len(t) + 1)] for i in range(len(s) + 1)]
for i in range(1, len(s) + 1):
    M[i][0] = i
for i in range(1, len(t) + 1):
    M[0][i] = i

for i in range(1, len(s) + 1):
    for j in range(1, len(t) + 1):
        if s[i - 1] == t[j - 1]:
            M[i][j] = M[i - 1][j - 1]
        else:
            M[i][j] = min(M[i - 1][j] + 1, M[i][j - 1] + 1,
                          M[i - 1][j - 1] + 1)
print(M[len(s)][len(t)])

Creating a Distance Matrix

This problem was fairly similar to "Counting Point Mutations" and "Error Correction in Reads". We are given a set of DNA strings and are asked to return the distance matrix of the strings. The distance between two given strings can be calculated by dividing the Hamming distance (i.e. the number of nucleotides that differ between the stings) with the length of the sequences (given that all sequences are of equal length). These values should then be printed on matrix form with 5 significant figures.

Sample Dataset

>Rosalind_9499

TTTCCATTTA

>Rosalind_0942

GATTCATTTC

>Rosalind_6568

TTTCCATTTT

>Rosalind_1833

GTTCCATTTA

Expected Output

0.00000 0.40000 0.10000 0.10000

0.40000 0.00000 0.40000 0.30000

0.10000 0.40000 0.00000 0.20000

0.10000 0.30000 0.20000 0.00000

To write this program I reused parts of the code from "Counting Point Mutations" and "Error Correction in Reads" and modified it to suit this problem. The following is the final code, which took me only 20 minutes to write. Hurray!

from Bio import SeqIO

reads = []

with open('sampledata.fasta', 'r') as f:

    for record in SeqIO.parse(f, 'fasta'):

        reads.append(str(record.seq))

read_len = len(reads[0])

for curr_read in reads:

    distance = []

    for comp_read in reads:

        hamming = 0

        for nt1, nt2 in zip(curr_read, comp_read):

            if nt1 != nt2:

                hamming += 1

        distance.append(str.format('{0:.5f}', hamming / read_len))

    print(*distance, sep=' ')

Transitions and Transversions

This problem was a really quick one. It took me less than 20 minutes to solve! Hurray!

We are asked to compare two sequences of equal length and classify the mutations as either transitions (substituting a purine to another purine or a pyrimidine to another pyrimidine) or transversions (substituting a purine to a pyrimidine or vice versa). We should then return the transition/transversion ratio for the sequences.

Sample dataset:
>Rosalind_0209
GCAACGCACAACGAAAACCCTTAGGGACTGGATTATTTCGTGATCGTTGTAGTTATTGGA
AGTACGGGCATCAACCCAGTT
>Rosalind_2200
TTATCTGACAAAGAAAGCCGTCAACGGCTGGATAATTTCGCGATCGTGCTGGTTACTGGC
GGTACGAGTGTTCCTTTGGGT

Expected output:
1.21428571429

The problem is very similar to Counting Point Mutations in which we calculated the Hamming distance. I used my code from that problem as a starting point and this is the altered code:

from Bio import SeqIO                        

sequences = []                               

handle = open('sampledata.fasta', 'r')       

for record in SeqIO.parse(handle, 'fasta'):  

    sequences.append(str(record.seq))        

handle.close()                               

s1 = sequences[0]                            

s2 = sequences[1]                            

transition = 0                               

transversion = 0                             

AG = ['A', 'G']                              

CT = ['C', 'T']                              

for nt1, nt2 in zip(s1, s2):                 

    if nt1 != nt2:                           

        if nt1 in AG and nt2 in AG:          

            transition += 1                  

        elif nt1 in CT and nt2 in CT:        

            transition += 1                  

        else:                                

            transversion += 1                

print('%0.11f' % (transition / transversion))

Counting Point Mutations

In this problem we are asked to calculate the Hamming distance of two sequences. For this task, my program would need to read the file containing the two sequences and then compare each position of the sequences to see if the nucleotides were the same. Each time there is a difference the program should add 1 to the Hamming distance.

The first thing i wrote was a piece of code to load the two sequences into a list. Because the sequences are separated by a new line in the sample file, I used the command line.strip to remove this from the resulting list. To be able to compare each position of the two sequences I used the command zip(), which iterates over each element of the strings. Then it was just a simple matter of comparing the nucleotides and adding 1 to the variable called "hamming" each time they differed. The following is the final code that I came up with:

seq = [line.strip('\n') for line in open('sampledata.txt')]

hamming = 0                                                

for nt1, nt2 in zip(seq[0],seq[1]):                        

    if nt1 != nt2:                                         

        hamming += 1                                       

print(hamming)                                             

Note that I haven't added a check to make sure that the length of the two sequences are the same. For real applications this would definitely be a good idea, but with the data set given from Rosalind I didn't feel it was necessary this time.

A simpler way to write this program would be to use the already existing package Distance which includes a function for calculating the Hamming distance.