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.

Interleaving Two Motifs

This time we are asked to find the shortest common supersequence of two given sequences. The shortest common supersequence is the shortest sequence that contains both of the given sequences as subsequences.

Sample Dataset
ATCTGAT
TGCATA

Expected Output
ATGCATGAT

If multiple solutions exist we are told to return any of them. To solve this problem, the easiest thing to do is to pick out the longest common subsequence (lcs) of the two given sequences and then add in the nucleotides that are missing from each sequence. To find the lcs, I reused most of the code I wrote for the problem Finding a Shared Spliced Motif. The resulting code can be found on my Github as well as below:

def lcs(s,t):
    lengths = [[0 for j in range(len(t) + 1)] for i in range(len(s) + 1)]
    for i, x in enumerate(s):
        for j, y in enumerate(t):
            if x == y:
                lengths[i + 1][j + 1] = lengths[i][j] + 1
            else:
                lengths[i + 1][j + 1] = max(lengths[i + 1][j], lengths[i][j + 1])
    spliced_motif = ''
    x, y = len(s), len(t)
    while x * y != 0:
        if lengths[x][y] == lengths[x - 1][y]:
            x -= 1
        elif lengths[x][y] == lengths[x][y - 1]:
            y -= 1
        else:
            spliced_motif = s[x - 1] + spliced_motif
            x -= 1
            y -= 1
    return(spliced_motif)


def scs(s,t):
    subseq = lcs(s,t)
    superseq = ''
    i, j = 0, 0
    for nt in subseq:
        if i < len(s):
            while s[i] != nt:
                superseq += s[i]
                i += 1
            i += 1
        if j < len(t):
            while t[j] != nt:
                superseq += t[j]
                j += 1
            j += 1
        superseq += nt
    if i < len(s):
        superseq += s[i:]
    if j < len(t):
        superseq += t[j:]
    return(superseq)

s, t = [line.strip() for line in open('rosalind_scsp.txt','r')]
print(scs(s,t))

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

Maximum Matchings and RNA Secondary Structures

This problem is very similar to the problem presented in "Perfect Matchings and RNA Secondary Structures". The difference is that in this problem we don't have the same number of A as U and G as C, so we don't get any perfect matchings. Instead we are given the task to find the total possible number of maximum matchings for a given RNA string.

Sample Dataset
>Rosalind_92
AUGCUUC

Expected Output
6

As in "Perfect Matchings and RNA Secondary Structures" the programming is very simple and it all comes down to the maths. We can simplify the problem by looking at AU and GC separately, and can describe the number of maximum matches for each part as max(A,U)!/(max(A,U)-min(A,U))! and max(G,C)!/(max(G,C)-min(G,C))!, respectively. To get the total number of maximum matches, simply multiply the answers.

For the first dataset I tried, I ran the program using Python 3.5 and got the result 434412400153932848464251158517855537428692992. This was deemed incorrect by Rosalind. After seeing a tip online, I ran the program using Python 2.7, and got 434412400153932864602133769790674698240000000 instead. So I tried my program with a new dataset and ran it with Python 2.7 and this time my answer was accepted.

Here follows the code I wrote to solve this problem (works in both 2.7 and 3.5):

from Bio import SeqIO
from math import factorial
sequence = ''
with open('sampledata.fasta', 'r') as f:
    for record in SeqIO.parse(f, 'fasta'):
        sequence = str(record.seq)

A, U, G, C = 0, 0, 0, 0
for nt in sequence:
    if nt == 'A':
        A += 1
    elif nt == 'U':
        U += 1
    elif nt == 'G':
        G += 1
    elif nt == 'C':
        C += 1

AU = factorial(max(A, U)) / factorial(max(A, U) - min(A, U))
GC = factorial(max(G, C)) / factorial(max(G, C) - min(G, C))
print(int(AU * GC))

Finding a Shared Spliced Motif

In this problem we are asked to find the longest common subsequence of two DNA-sequences, s and t. As opposed to a substring, a subsequence does not need to occur contiguously in s and t. To solve this problem it might be a good idea to have a look at the Wikipedia page on the longest common subsequence problem.

Sample Dataset
>Rosalind_23
AACCTTGG
>Rosalind_64
ACACTGTGA

Expected Output
AACTGG (or any other if there are multiple subsequences of the same length)

The following is a solution to the problem using dynamic programming:

from Bio import SeqIO
sequences = []
handle = open('sampledata.fasta', 'r')
for record in SeqIO.parse(handle, 'fasta'):
    sequences.append(str(record.seq))
s = sequences[0]
t = sequences[1]

lengths = [[0 for j in range(len(t) + 1)] for i in range(len(s) + 1)]
#creates array of len(s) containing arrays of len(t) filled with 0
for i, x in enumerate(s):
    for j, y in enumerate(t):
        if x == y:
            lengths[i + 1][j + 1] = lengths[i][j] + 1
        else:
            lengths[i + 1][j + 1] = max(lengths[i + 1][j], lengths[i][j + 1])

spliced_motif = ''
x, y = len(s), len(t)
while x * y != 0:
    if lengths[x][y] == lengths[x - 1][y]:
        x -= 1
    elif lengths[x][y] == lengths[x][y - 1]:
        y -= 1
    else:
        spliced_motif = s[x - 1] + spliced_motif
        x -= 1
        y -= 1
print(spliced_motif)

Perfect Matchings and RNA Secondary Structures

This is yet another problem where the programming is really simple, but we need too figure out the maths before we can write the program.

The problem is linked to RNA folding and we are assuming that in the given RNA string, every nucleotide forms part of a base pair in the RNA molecule. We then look at the n long molecule as a graph containing n nodes where the A-nodes can form an edge (i.e base pair) with the U-nodes and the G-nodes can form an edge with the C-nodes. From this scenario we are asked to find all possible perfect matchings, where a matching is a set of edges in a graph where none of the edges include the same node. A matching is said to be perfect if for a graph of 2n nodes it contains n edges.

Given dataset:
>Rosalind_23

AGCUAGUCAU
(The sequence always contains equally many A's as U's and equally many G's as C's)

Expected Output:
12

What we need to realise is that what we are looking at can be viewed as two separate graphs, one for AU bonding and one for GC bonding. These two are complete bipartite graphs and we can describe the number of perfect matchings for each graph, if n1 = nr of A's and n2 = nr of G's, as n1! and n2!. To get the total amount of perfect matchings for the two graphs combined, we simply multiply these numbers. We thereby get:

matchings = n1! * n2!

All that is left to do now is to write a simple program to run the calculations:

from Bio import SeqIO                      

from math import factorial                 

sequence = ''                              

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

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

    sequence = str(record.seq)             

handle.close()                             

AU = 0                                     

GC = 0                                     

for nt in sequence:                        

    if nt == 'A':                          

        AU += 1                            

    elif nt == 'G':                        

        GC += 1                            

matchings = factorial(AU) * factorial(GC)  

print(matchings)                           

Longest Increasing Subsequence

This problem had me going for quite a while, and I ended up getting a lot of help from various places on the internet. But I have now managed to solve it and I will try to explain my solution here. Let's start with the description.

In this problem we are given a permutation of length n (where n is smaller than or equal to 10 000) and we are asked to find the longest increasing and decreasing subsequence of the permutation. A subsequence of a permutation is a collection of elements from the permutation in the order that they appear (note that they don't have to come directly after one another). The subsequence is said to be increasing if the elements increase (from left to right) and it is said to be decreasing if the elements decrease (again, left to right). For example, in the permutation (5, 1, 4, 2, 3) the subsequence (1, 2, 3) is increasing and the subsequence (5, 4, 2) is decreasing.

The fact that the subsequence can 'skip' numbers messed up my initial solution for this problem, but after spending some time trying to figure out how the algorithm described in this Wikipedia page works and looking through a very good explanation of it on Stack Overflow (as well as working through the algorithm with pen and paper to figure out how it operates. I know, I'm a tad bit old-school), I managed to adapt the code on Stack Overflow to also find the longest decreasing subsequence. The code below is what I ended up with. The function for finding the decreasing subsequence is talmost he same as for finding the increasing subsequence, so I have highlighted the differences. I would highly recommend having a look at the thread on Stack Overflow if you would like a good explanation of the algorithm.

data = []                                #This first bit reads the file which
with open('sampledata.txt', 'r') as f:   #contains the the length of the 
    for line in f:                       #permutation and then the permutation, 
        for nr in line.split():          #separated by spaces. The numbers are 
            data.append(int(nr))         #appended to a list as integers.
perm = data[1:]


def increasing(seq):
    P = [None] * len(seq)
    M = [None] * len(seq)

    L = 1
    M[0] = 0
    for i in range(1, len(seq)):
        lo = 0
        hi = L
        if seq[M[hi - 1]] < seq[i]:
            j = hi
        else:
            while hi - lo > 1:
                mid = (hi + lo) // 2
                if seq[M[mid - 1]] < seq[i]:
                    lo = mid
                else:
                    hi = mid

            j = lo
        P[i] = M[j - 1]
        if j == L or seq[i] < seq[M[j]]:
            M[j] = i
            L = max(L, j + 1)

    result = []
    pos = M[L - 1]
    for k in range(L):
        result.append(seq[pos])
        pos = P[pos]

    return (result[::-1])


def decreasing(seq):
    P = [None] * len(seq)
    M = [None] * len(seq)

    L = 1
    M[0] = 0
    for i in range(1, len(seq)):
        lo = 0
        hi = L
        if seq[M[hi - 1]] > seq[i]:
            j = hi
        else:
            while hi - lo > 1:
                mid = (hi + lo) // 2
                if seq[M[mid - 1]] > seq[i]:
                    lo = mid
                else:
                    hi = mid

            j = lo
        P[i] = M[j - 1]
        if j == L or seq[i] > seq[M[j]]:
            M[j] = i
            L = max(L, j + 1)

    result = []
    pos = M[L - 1]
    for k in range(L):
        result.append(seq[pos])
        pos = P[pos]

    return (result[::-1])

incr = increasing(perm)
decr = decreasing(perm)

print(*incr)
print(*decr)

Mortal Fibonacci Rabbits

In this Rosalind problem we are revisiting the Fibonacci Rabbits. This time the bunnies sadly aren't immortal and will die after m months. Apart from that they behave as in the previous problem but they only produce one pair of offspring per month in adulthood. This leaves us with two different recurrence relations to calculate the amounts of rabbits present in any given month n. The months before any rabbits start dying off, i.e. when n ≤ m can be described by the recurrence relation of the Fibonacci numbers. The months after this, i.e. when n > m, can be described by the following formula:

Fn = Fn-1 + Fn-2 - Fn-(m+1)

To calculate the number of rabbits after 96 months, when the rabbits are dying after 17 months, I wrote the following program:

n = 96                                                                         

m = 17                                                                         

bunnies = [1, 1]                                                               

months = 2                                                                     

while months < n:                                                              

    if months < m:                                                             

        bunnies.append(bunnies[-2] + bunnies[-1])                              

    elif months == m or count == m + 1:                                        

        bunnies.append(bunnies[-2] + bunnies[-1] - 1)                          

    else:                                                                      

        bunnies.append(bunnies[-2] + bunnies[-1] - bunnies[-(                  

            m + 1)])                                                           

    months += 1                                                                

print(bunnies[-1])                                                             

So, care to take a guess what the result after 96 months was? 51159459138167757395. Pairs...

Rabbits and Recurrence Relations

The fourth Rosalind problem is a bit different than the previous ones. This time we are asked to calculate how many rabbit pairs there are after n months if every grown pair produces k new pairs each month and it takes one month for a new pair to become reproductive. The problem describes the Fibonacci sequence and how it can be applied to the simplest version of the problem, the case when k=1. It's been a while since I worked with mathematical sequences, so this was a fun way to dust of things I learned during my bachelor's maths courses.

The Fibonacci sequence works on this problem when k=1 and its recurrence relation is written as follows:

Fn = Fn-1 + Fn-2

The first thing I needed to do was to figure out the  recurrence relation for any value on k. After doodling a bit and playing around with the numbers I figured out that this must be it:

Fn = Fn-1 + Fn-2 * k

Once I knew this I wrote the following program:

n = 5                      

k = 3                      

big = 1                    

small = 1                  

for months in range(1,n-1):

      bigger = big + small*k   

      small = big              

      big = bigger             

print(big)                 

I named Fn-1 and Fn-2 big and small so that I wouldn't confuse myself which was the larger number and which was the smaller of the two. Since the problem states that F1 = 1 the program is set to only iterate n-1 times to reach the answer for Fn. The first expression in the for-loop calculates Fn. The other two terms update Fn-1 and Fn-2 in preparation of the next iteration.

I really liked this problem, mainly because it felt a bit more challenging than the previous three and because I got to practice some maths skills that I haven't used in a while. It took me a bit longer than the previous ones, about an hour, but most of that time I spent working on the recurrence relation and when I got to the coding part I managed to write the program pretty quickly and with little fuss.