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

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

Ordering Strings of Varying Length Lexicographically

In this problem we are given a string representing an alphabet. We are also given an integer n ≤ 4. We are to return a list of all strings of length at most n that can be formed from the alphabet string. Repeats should be included and the strings should be sorted lexicographically according to the given alphabet (and this time they really should, as opposed to in “Enumerating k-mers Lexicographically”).

Sample Dataset
D N A
3

Expected Output
D
DD
DDD
DDN
DDA
DN
DND
DNN
DNA
DA
DAD
DAN
DAA
N
ND
NDD
NDN
NDA
NN
NND
NNN
NNA
NA
NAD
NAN
NAA
A
AD
ADD
ADN
ADA
AN
AND
ANN
ANA
AA
AAD
AAN
AAA

To begin with this problem is fairly similar to “Enumerating k-mers Lexicographically”, so I started by having a look at what I wrote there. As in that problem, I decided to use itertools, but this time, since the answer should include repeats, I used product instead. I also put the creation of the permutations into a for-loop in order to get all the different lengths, as you can see below. This generated a nested list, which I then flattened into a single list using the itertools function chain(). Then all that was left to do was to sort and print the permutations to a file. In order to sort them, I made use of the function sorted() and used a lambda function as key. Here is an explanation of the lambda function in Python. The final code can be seen below:

import itertools
n = 3
alphabet = 'DNA'
perm = []
for i in range(1, n + 1):
    perm.append(list(map(''.join, (itertools.product(alphabet, repeat=i)))))
permutations = list(itertools.chain(*perm))
srt_perm = sorted(permutations,
                  key=lambda word: [alphabet.index(c) for c in word])
with open('answer.txt', 'a') as f:
    for j in srt_perm:
        f.write('%s\n' % j)

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)

Speeding Up Motif Finding

In this problem we are looking at how to speed up motif finding using the Knuth-Morris-Pratt algorithm. A good explanation of the algorithm can be found here. We are given a DNA-string and are expected to return the failure array for the sting.

Sample Dataset
>Rosalind_87
CAGCATGGTATCACAGCAGAG

Expected Output
0 0 0 1 2 0 0 0 0 0 0 1 2 1 2 3 4 5 3 0 0

The following code is what I ended up with:

from Bio import SeqIO
record = SeqIO.read('sampledata.fasta', 'fasta')
sequence = list(record.seq)

F_array = [0] * len(sequence)
k = 0
for i in range(2, len(sequence) + 1):
    while k > 0 and sequence[k] != sequence[i - 1]:
        k = F_array[k - 1]
    if sequence[k] == sequence[i - 1]:
        k += 1
    F_array[i - 1] = k

with open('array.txt', 'w') as answer:
    answer.write(' '.join(map(str, F_array)))

k-Mer Composition

In this problem we are given a DNA-string for which we are to return a matrix consisting of the frequencies of all possible 4-mers in the string, ordered alphabetically.

Sample Dataset
>Rosalind_6431
CTTCGAAAGTTTGGGCCGAGTCTTACAGTCGGTCTTGAAGCAAAGTAACGAACTCCACGGCCCTGACTACCGAACCAGTTGTGAGTACTCAACTGGGTGAGAGTGCAGTCCCTATTGAGTTTCCGAGACTCACCGGGATTTTCGATCCAGCCTCAGTCCAGTCTTGTGGCCAACTCACCAAATGACGTTGGAATATCCCTGTCTAGCTCACGCAGTACTTAGTAAGAGGTCGCTGCAGCGGGGCAAGGAGATCGGAAAATGTGCTCTATATGCGACTAAAGCTCCTAACTTACACGTAGACTTGCCCGTGTTAAAAACTCGGCTCACATGCTGTCTGCGGCTGGCTGTATACAGTATCTA
CCTAATACCCTTCAGTTCGCCGCACAAAAGCTGGGAGTTACCGCGGAAATCACAG

Expected Output
4 1 4 3 0 1 1 5 1 3 1 2 2 1 2 0 1 1 3 1 2 1 3 1 1 1 1 2 2 5 1 3 0 2 2 1 1 1 1 3 1 0 0 1 5 5 1 5 0 2 0 2 1 2 1 1 1 2 0 1 0 0 1 1 3 2 1 0 3 2 3 0 0 2 0 8 0 0 1 0 2 1 3 0 0 0 1 4 3 2 1 1 3 1 2 1 3 1 2 1 2 1 1 1 2 3 2 1 1 0 1 1 3 2 1 2 6 2 1 1 1 2 3 3 3 2 3 0 3 2 1 1 0 0 1 4 3 0 1 5 0 2 0 1 2 1 3 0 1 2 2 1 1 0 3 0 0 4 5 0 3 0 2 1 1 3 0 3 2 2 1 1 0 2 1 0 2 2 1 2 0 2 2 5 2 2 1 1 2 1 2 2 2 2 1 1 3 4 0 2 1 1 0 1 2 2 1 1 1 5 2 0 3 2 1 1 2 2 3 0 3 0 1 3 1 2 3 0 2 1 2 2 1 2 3 0 1 2 3 1 1 3 1 0 1 1 3 0 2 1 2 2 0 2 1 1

Regardless of the given sequence, the possible 4-mers are always the same and we can generate them with the following code:

import itertools
nt = 'ACGT'        #Use this order of nt to get correct order later without sorting
permutations = itertools.product(nt, repeat=4)

kmers = []
for i, j in enumerate(list(permutations)):
    kmer = ''
    for item in j:
        kmer += str(item)
    kmers.append(kmer)

This gives us a list of all the possible 4-mers in alphabetical order (note that the funktion sorts the permutations in the order the letters are listed in. If nt is not entered alphabetically you will need to sort kmers).

Now we can extrakt the sequence from the FASTA file and use regex to find all the occurrences of the k-mers. Remember to use ?= in the pattern to include overlapping k-mers:

import re
from Bio import SeqIO

record = SeqIO.read('sampledata.fasta', 'fasta')
sequence = record.seq

A = []
for k in kmers:
    occurence = 0
    pattern = re.compile(r'(?=(' + k + '))')
    for l in re.findall(pattern, str(sequence)):
        occurence += 1
    A.append(occurence)
print(*A, sep=' ')

Finding a Spliced Motif

This time we are once again asked to find the position of a given subsequence for a given sequence. However, this time we should take into account that the subsequence can be spliced, i.e. it can be split up in the sequence and there can be other nucleotides between the parts. There could be multiple ways that the subsequence can be found in the sequence, but we only have to return one of them in the form of the positions each letter of the subsequence has in the sequence.

Sample dataset:
>Rosalind_14
ACGTACGTGACG
>Rosalind_18

GTA

Expected output:
3 8 10
(or any of the other possible combinations)

My first thought was to look at my solution for Finding a Motif in DNA, but in that problem I used Biopython to find the motifs and I wasn't able to find a way to adapt it to finding spliced motifs. Another thought was to use regular expressions. However, I quite quickly managed to come up with this rather simple solution instead:

from Bio import SeqIO                      

sequences = []                             

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

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

    sequences.append(str(record.seq))      

handle.close()                             

s = sequences[0]                           

t = sequences[1]                           

pos = 0                                    

positions = []                             

for i in range(len(t)):                    

    for j in range(pos, len(s)):           

        pos += 1                           

        if len(positions) < len(t):        

            if t[i] == s[j]:               

                positions.append(pos)      

                break                      

print(*positions, sep=' ')                 

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)                           

Enumerating k-mers Lexicographically

In this problem we  are given a set of characters and a positive integer n and we are asked to return all the strings of length n that can be formed from the characters (including duplicates). The answer should be sorted lexicographically according to the order that the characters are presented in the sample data (i.e. not alphabetically). Once again the built-in Python library itertools comes in handy. This time I used the function products because the problem asks for duplicates as well (that is, each letter may occur up to n times in a string an not only once). The following is the code I came up with:

import itertools                                  

n = 4                                             

s = ['H', 'U', 'P', 'M']                          

perm = itertools.product(s, repeat=n)             

for i, j in enumerate(list(perm)):                

    permutation = ''                              

    for item in j:                                

        permutation += str(item)                  

    with open('answer.txt', 'a') as text_file:    

        print(permutation.strip(), file=text_file)

Now, this code works just fine, and when I look at the results they match what is described in the problem description, but I can't seem to get my answers to pass. I have tried seven different datasets now, either pasting the answer into the browser or uploading the answer file. I have made sure there are no additional blank spaces or new lines and even tried concatenating the answer into one long string, but nothing I do seems to work. I have noticed that several people have had the same problem rather recently, which has lead me to believe that the problem is probably not from my side, but rather from Rosalind's side. Either way, I think I will just let this problem rest for a while and press on with the other problems instead.

EDIT:

The problem has finally been solved! Apparently, the answer SHOULD be sorted alphabetically, even though it says not to in the note below the problem statement. Also, according to someone in the questions section, the permutations should not be separated with new lines but instead with blank spaces. This is rather annoying because this is exactly what my first attempt at the program did, but then I read the note and changed it to the code above... Anyhow, I can now move on with my life. Below is the new code that produces results that are accepted by Rosalind, although that isn't what they asked for...

import itertools                                  

n = 4                                             

s = ['H', 'U', 'P', 'M']                          

perm = itertools.product(s, repeat=n)

answer = []                                       

for i, j in enumerate(list(perm)):                

    permutation = ''                              

    for item in j:                                

        permutation += str(item)                  

    answer.append(permutation)                    

sorted_answer = sorted(answer)                    

with open('answer.txt', 'a') as text_file:        

    print(*sorted_answer, sep=' ', file=text_file)

RNA Splicing

This time we are looking at splicing. When mRNA has been transcribed from a DNA strand it often contains introns that have to be spliced away before the mRNA is translated into a protein. In this task we are given a FASTA file containing the sequences of one DNA strand and several introns. The goal is to remove the introns and print the resulting protein.

So, I wrote the following program that does exactly that:

from Bio import SeqIO                          

from Bio.Seq import Seq                        

from Bio.Alphabet import generic_dna           

sequences = []                                 

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

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

    sequence = ''                              

    for nt in record.seq:                      

        sequence += nt                         

    sequences.append(sequence)                 

handle.close()                                 

long_seq = sequences[0]                        

introns = sequences[1:]                        

for i in range(len(introns)):                  

    long_seq = long_seq.replace(introns[i], '')

long_seq = Seq(long_seq)                       

print(long_seq.translate(to_stop=True))        

Locating Restriction Sites

In this problem we are asked to find the reverse palindromes of a given DNA sequence. A piece of DNA is said to be a reverse palindrome when it is equal to it's reverse complement.

This problem was quite similar to the previous problem "Finding a shared motif" so I adapted that code slightly and ended up with this:

from Bio import SeqIO                           

record = SeqIO.read('sampledata.fasta', 'fasta')

frw_seq = str(record.seq)                       

rev_seq = str(record.seq.complement())          

for i in range(len(frw_seq)):                   

    for j in range(i, len(frw_seq)):            

        m = frw_seq[i:j + 1]                    

        rev_m = rev_seq[i:j + 1]                

        if len(m) >= 4 and len(m) <= 12:        

            if m == rev_m[::-1]:                

                print(i + 1, len(m))            

The program works through the forward string and its complementary string and picks out all pieces that are  longer or equal to 4 but shorter or equal to 12. Then it compares each forward piece with its reverse complement and if they are equal, prints the position and the size of the palindrome. There are probably much more efficient ways to do this, but it works and I managed to write the program a lot quicker than I thought I would!

Finding a Shared Motif

When I first looked at this problem I thought that solving it would be pretty straightforward and that I could somehow reuse some of the code I wrote for Finding a Motif in DNA, but I was wrong. In this problem we are asked to find the longest common substring of a set of DNA-sequences in a FASTA-file. The problem sounds similar to the one presented in Finding a Motif in DNA, but this time we don't know anything about the motif, so using the same approach as last time won't work. Time to rethink.

The first thing however, was as usual to parse the FASTA file and get the data into a suitable format. I chose to put the sequences as strings in a list. The following piece of code does that:

from Bio import SeqIO                      

sequences = []                             

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

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

    sequence = []                          

    seq = ''                               

    for nt in record.seq:                  

        seq += nt                          

    sequences.append(seq)                  

handle.close()                             

What I then wanted to do was to compare all the possible motifs in the shortest sequence to the remaining sequences. To do this I first sorted the list containing the sequences and picked out the shortest one. I could then iterate over all the possible motifs in that sequence and test if they were also present in all of the other sequences. The longest of the motifs that are present in all of the sequences is then saved and printed. This was the easiest solution to the problem that I could think of, but I'll admit it's probably not the neatest or the most efficient solution.

srt_seq = sorted(sequences, key=len)     

short_seq = srt_seq[0]                   

comp_seq = srt_seq[1:]                   

motif = ''                               

for i in range(len(short_seq)):          

    for j in range(i, len(short_seq)):   

        m = short_seq[i:j + 1]           

        found = False                    

        for sequ in comp_seq:            

            if m in sequ:                

                found = True             

            else:                        

                found = False            

                break                    

        if found and len(m) > len(motif):

            motif = m                    

print(motif)                             

Consensus and Profile

For this problem we are asked to find a sequence that is as similar to the given set of sequences as possible (the consensus sequence). To complete the task the problem states that we should make a profile matrix containing the number of times each nucleotide (A, T, G or C) occurs at each position of the sequences. We should then be able to build the consensus sequence from this.

I started by breaking down the problem into four different parts. First off, I wanted the program to parse the FASTA-file and add the sequences to a nested list. This way I could then easily iterate over the sequences in the list in order to count the occurrence of each nucleotide in each position. The following code opens the FASTA-file and adds each sequence to a separate list inside the list sequences. 

from Bio import SeqIO                      

sequences = []                             

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

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

     sequence = []                         

     for nt in record.seq:                 

          sequence.extend(nt)              

     sequences.append(sequence)            

handle.close()                             

For the second part of the program I needed to count the occurrence of each nucleotide at each position and summarize the results in a profile matrix. I must admit that I had some problems with this part of the program. All those problems were finally explained by the way I had initialized the profile matrix:

profile = [[0]*len(sequences)]*4

What I didn't realize when I wrote this is that the *4 operation makes dependent copies, so while this gave me the look I wanted for my profile matrix, each time I added something to a specific position of one of the nested lists, the same was added to the corresponding position of the other three. When I finally realized this I changed the profile matrix into a numpy matrix and all problems were resolved. The following code iterates over each nucleotide in each of the sequences and counts the number of times they occur at a specific position. The results are added continually to the profile matrix.

import numpy                                                  

profile = numpy.zeros((4, len(sequences[0])), dtype=numpy.int)

for i,line in enumerate(sequences):                           

     for j, nt in enumerate(line):                            

          if nt == 'A':                                       

               profile[0][j] += 1                             

          elif nt == 'C':                                     

               profile[1][j] += 1                             

          elif nt == 'G':                                     

               profile[2][j] += 1                             

          elif nt == 'T':                                     

               profile[3][j] += 1                             

For the third part the program needed to evaluate the profile matrix and create the consensus sequence (any one of them would suffice in the case that there were more than one). The following code does just that. If two nucleotides occur equally many times in a position than only the first one in the order given in the program will be added to the consensus sequence.

consensus = ''                                                  

for A,C,G,T in zip(profile[0],profile[1],profile[2],profile[3]):

     if A >= C and A >= G and A >= T:                           

          consensus += 'A'                                      

     elif C >= A and C >= G and C >= T:                         

          consensus += 'C'                                      

     elif G >= A and G >= C and G >= T:                         

          consensus += 'G'                                      

     elif T >= A and T >= C and T >= G:                         

          consensus += 'T'                                      

In the fourth and final part of the program I just needed it to print the consensus sequence and the profile matrix in the format specified in the problem description. The final part looks like this:

print(consensus)                                   

print('A: ' + ' '.join(str(e) for e in profile[0]))

print('C: ' + ' '.join(str(e) for e in profile[1]))

print('G: ' + ' '.join(str(e) for e in profile[2]))

print('T: ' + ' '.join(str(e) for e in profile[3]))

I used the .join command to get rid of the brackets of the matrix that were otherwise printed.

Finding a Motif in DNA

This time we are asked to find the positions of a given motif in a given DNA-sequence. Having just read about Seq Objects in Biopython, I had noticed a function called find, which can be used to find the position of motifs in sequences. However, this function only seems to output the position of the first occurrence of the motif. Instead I decided to go with the Biopython function motifs (have a look at chapter 14.6). Here is my code:

from Bio import motifs                                      

from Bio.Seq import Seq                                     

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

instances =[Seq(data[1])]                                   

m = motifs.create(instances)                                

sequence = Seq(data[0])                                     

positions = ''                                              

for pos, seq in m.instances.search(sequence):               

    positions += str(pos+1) + ' '                           

print(positions)                                            

The function returns all the positions of the motifs, but it seems to have a different definition of the position than Rosalind. In order to receive the correct result I had to add 1 to all the positions, as you can see in line 9. I also added a blank space in this step to achieve the correct formatting.