Introduction to Alternative Splicing

In this problem we continue looking at subsets. This time we imagine them as exons and want to know how many different combinations we can form from subsets of the set. We are given two integers n and m, which are the length of the set ([1,2,3,4,5,6] in the sample case) and the minimum length of the subsets to be formed. The order of the exons cannot be altered, so for n=3 and m=2 we would get the four subsets [1,2], [1,3], [2,3] and [1,2,3]. The answer should be given modulo 1000000.

Sample Dataset
6 3

Expected Output
42

The problem can quite easily be solved by counting all the different combinations of each length from m to  n, and then summarising them to form the final answer. The number of combinations can be calculated as n!/(k!(n-k)!), where k is the length of the subset. In python3 code this can look like:

from math import factorial
n, m = 6, 3
subsets = 0
for k in range(m, n + 1):
    subsets += (factorial(n) // (factorial(k) * factorial(n - k)))
print(subsets % 1000000)

At first I used / instead of //, which worked fine for the sample dataset, but when I downloaded another dataset the numbers were bigger and therefore they got too big for float division and I got the error message "OverflowError: integer division result too large for a float". A quick Google search took me to this thread on Stack Overflow that solved my problem.

Counting Subsets

In this problem we are introduced to set theory (causing flashbacks from the linear algebra course...). For a set {1,2,3,...,n}, we are asked to return the number of subsets of the set. What we need to remember is that the empty set {} and the whole set {1,2,3,...,n} are also subsets (aptly explained here). So, to find the total number of subsets, all we need to do is calculate 2^n. Since this number quickly gets very large, the problem states that we should return the answer modulo 1 000 000.

Sample Dataset
3

Expected Output
8

The following program is probably the shortest one yet:

n = 3                      

print(2**n % 1000000)      

Or, if you like, a fancier version where the user can input n in the terminal:

n = int(input('Enter n: '))

print(2**n % 1000000)      

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

Counting Phylogenetic Ancestors

This was probably the easiest problem so far. In fact, you don't even need to do any programming because the calculation is so simple. The only thing you need to realise is that for any unrooted binary tree with n leaves, the number of internal nodes is equal to n-2.

An easy way to convince yourself of this is to doodle the trees on some paper. Then you can see that a tree with 3 leaves has 1 internal node, 4 leaves have 2 internal nodes, 5 leaves have 3 internal nodes, 6 has 4, 7 has 5, 8 has 6 and so on:

Enumerating Oriented Gene Orderings

This is yet another problem with permutations. As in the previous problems of this type we get a positive integer n and are supposed to find some permutations. This time we are looking at signed permutations, that is, for all the positive integers in the range 1 to n we should also take into account their corresponding negative value. As in Enumerating Gene Orders, we should print the total number of permutations followed by each of the permutations.

Sample dataset:

2

Expected output:

8

-1 -2

-1 2

1 -2

1 2

-2 -1

-2 1

2 -1

2 1

Note that we should not count permutations of the same integer (i. e. -1 1, 2 -1 and so on should not be included).

For this problem we can once again utilise the itertools function permutations. However, this time we need to pair it to the itertools function product. The following is my code:

import itertools                                                           

n = 3                                                                      

permutation = []                                                           

nr = 0                                                                     

for i in itertools.permutations(list(range(1, n + 1))):                    

    for j in itertools.product([-1, 1], repeat=len(list(range(1, n + 1)))):

        perm = [a * sign for a, sign in zip(i, j)]                         

        permutation.append(perm)                                           

        nr += 1                                                            

print(nr)                                                                  

for i in range(len(permutation)):                                          

    print(*permutation[i], sep=' ')                                        

Partial Permutations

We are back with the permutations once again. This time we are set the task to find the number of partial permutations that can be formed from a set of numbers. We are given two integers, n and k, where the range of 1 to n is the set of numbers and k is the length of the partial mutations and should return the number of partial mutations modulo 1 000 000.

Sample data:
21 7

Expected output:
51200

I wrote two programs that solve this problem. The first uses the math module permutations to find all possible permutations and count them. You can have a look at the code below. However, this program is inefficient which is why I chose to write a second, faster program.

Slow approach:

import itertools
n = 21
k = 7
count = 0
for p in itertools.permutations(range(1, n + 1), k):
    count += 1
print(count % 1000000)

Faster approach:

n = 21
k = 7

partial_perm = 1

for i in range(k):
    partial_perm *= (n - i)
print(partial_perm % 1000000)

Both of the programs work and complete the problem within the time limits of Rosalind, but the second program is much faster than the first one. 

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 Gene Orders

In this problem we are given the task to find all the possible permutations of a list from 1 to n, where n ≤ 7. We are given n and are supposed to return the total number of permutations and a list of all the different permutations. To clarify, if we are told that n = 3, then the total number of permutations is 3! which is equal to 6 and we should print the following:

6      

1 2 3

1 3 2

2 1 3

2 3 1

3 1 2

3 2 1

It was fairly easy to make this program since Python3 has a built in library for factorials and another built in library for making permutations. The following is the final code:

import math                                         

import itertools                                    

n = 3                                               

print(math.factorial(n))                            

perm = itertools.permutations(list(range(1, n + 1)))

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

    permutation = ''                                

    for item in j:                                  

        permutation += str(item) + ' '              

    print(permutation)                              

Open Reading Frames

This time we are asked to find all the open reading frames of a DNA-sequence and to translate these into the different proteins they encode. Only unique proteins should be printed and we need to remember that there are in total six different reading frames ( 3 for each DNA-strand).

My first thought was to use Biopython for this, because surely there is a built in function for this sort of problem. This turned out not to be the case, but in the Biopython documentation there is a description on how to identify open reading frames. I tried this, and managed to print some protein sequences with it, but they were completely wrong compared to the ones shown in the Rosalind example. I tried to fix the code, but finally I gave up and decided to write my own code instead, using regular expressions.

The first thing the program needs to do after reading the given FASTA-file is to save the forward strand and then also pick out the reverse complement of it. Then we need to define the pattern that should be searched for using regular expression. In this case, we want to find the parts of the sequence that start with a start codon and ends with a stop codon. For regex this can be written as:

(?=(ATG(?:...)*?)(?=TAG|TGA|TAA))

The program should then look for this pattern in both the forward strand and the reverse complement. If an ORF is found it should be translated into the corresponding protein sequence and if that one is unique it should be printed. The following code is what I ended up with. Again, not the prettiest program ever written, but I managed to write it within an hour.

import re                                                 

from Bio import SeqIO                                     

from Bio.Seq import Seq                                   

from Bio.Alphabet import generic_dna                      

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

pattern = re.compile(r'(?=(ATG(?:...)*?)(?=TAG|TGA|TAA))')

frw_seq = record.seq                                      

rev_seq = frw_seq.reverse_complement()                    

sequences = []                                            

for m in re.findall(pattern, str(frw_seq)):               

    dna_seq = Seq(m, generic_dna)                         

    prot_seq = dna_seq.translate()                        

    if prot_seq not in sequences:                         

        sequences.append(prot_seq)                        

for n in re.findall(pattern, str(rev_seq)):               

    rev_dna_seq = Seq(n, generic_dna)                     

    rev_prot_seq = rev_dna_seq.translate()                

    if rev_prot_seq not in sequences:                     

        sequences.append(rev_prot_seq)                    

for i, s in enumerate(sequences):                         

    print(s)                                              

Inferring mRNA from Protein

This problem was fairly straightforward, so I thought I'd try my hand at making my code a bit more readable by dividing it into functions. I have done this many times before during courses, but it has been a while since then.

For this problem we are asked to find all possible mRNA-sequences that a given protein sequence may have come from. Because an amino acid can be encoded by several different codons, the number of sequences that the protein may have come from rapidly gets very large as the protein sequence gets longer, so we are asked for the answer modulo 1000000.

For every possible sequence (s) that encodes the given protein, there are three possible stop codons. This means that we get three variants (s-stop1, s-stop2, s-stop3) out of each possible sequence. So, the total amount of possible sequences becomes three times the number of codons that each amino acid in the protein could be encoded by (lets call it the frequency of the amino acid):

3*Freq(aa1)*Freq(aa2)... and so on

The program I wrote looks like this:

codons = {                                                     

    'UUU': 'F',     'CUU': 'L',     'AUU': 'I',     'GUU': 'V',

    'UUC': 'F',     'CUC': 'L',     'AUC': 'I',     'GUC': 'V',

    'UUA': 'L',     'CUA': 'L',     'AUA': 'I',     'GUA': 'V',

    'UUG': 'L',     'CUG': 'L',     'AUG': 'M',     'GUG': 'V',

    'UCU': 'S',     'CCU': 'P',     'ACU': 'T',     'GCU': 'A',

    'UCC': 'S',     'CCC': 'P',     'ACC': 'T',     'GCC': 'A',

    'UCA': 'S',     'CCA': 'P',     'ACA': 'T',     'GCA': 'A',

    'UCG': 'S',     'CCG': 'P',     'ACG': 'T',     'GCG': 'A',

    'UAU': 'Y',     'CAU': 'H',     'AAU': 'N',     'GAU': 'D',

    'UAC': 'Y',     'CAC': 'H',     'AAC': 'N',     'GAC': 'D',

    'UAA': 'Stop',  'CAA': 'Q',     'AAA': 'K',     'GAA': 'E',

    'UAG': 'Stop',  'CAG': 'Q',     'AAG': 'K',     'GAG': 'E',

    'UGU': 'C',     'CGU': 'R',     'AGU': 'S',     'GGU': 'G',

    'UGC': 'C',     'CGC': 'R',     'AGC': 'S',     'GGC': 'G',

    'UGA': 'Stop',  'CGA': 'R',     'AGA': 'R',     'GGA': 'G',

    'UGG': 'W',     'CGG': 'R',     'AGG': 'R',     'GGG': 'G' 

}                                                              

def codon_frequence():                                         

    frequence = {}                                             

    for k, v in codons.items():                                

        if v not in frequence:                                 

            frequence[v] = 0                                   

        frequence[v] += 1                                      

    return (frequence)                                         

def possible_sequences(sequence):                              

    f = codon_frequence()                                      

    n = f['Stop']                                              

    for aa in sequence:                                        

        n *= f[aa]                                             

    return (n % 1000000)                                       

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

    prot_seq = f.read().strip()                                

print(possible_sequences(prot_seq))                            

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.