Introduction to Set Operations

In this problem we are introduced to set operations. After finishing a course in algebra and geometry almost six years ago, I thought that I'd never have to work with set theory again. Apparently I was wrong. Fortunately, at least in this problem, it turned out to be much easier than I remembered.

For a set A and B, that are both subsets of U (where U is the set {1, 2, ..., n} for a given n), we are asked to find the following six sets: A∪B, A∩B, A−B, B−A, the set complement of A and the set complement of B (with respect to U).

Sample Dataset
10
{1, 2, 3, 4, 5}
{2, 8, 5, 10}

Expected Output
{1, 2, 3, 4, 5, 8, 10}
{2, 5}
{1, 3, 4}
{8, 10}
{8, 9, 10, 6, 7}
{1, 3, 4, 6, 7, 9}

As it turns out, Python has built in set operators, which made this problem a lot easier to solve. The following code is what I came up with. You can also find it on my GitHub.

#read file and format data into sets of integers
data = []
with open('rosalind_seto.txt','r') as f:
    for line in f:
        data.append(line.strip('\n'))
U = set(range(1, int(data[0])+1))
A = set([int(x) for x in data[1].strip('{}').split(',')])
B = set([int(x) for x in data[2].strip('{}').split(',')])

#perform set operations
union = A.union(B)
intersection = set(A).intersection(B)
diff_AB = A-B
diff_BA = B-A
A_comp = U-A
B_comp = U-B

#print answers
print(union)
print(intersection)
print(diff_AB)
print(diff_BA)
print(A_comp)
print(B_comp)

Learning Git

Learning Git is something that has been on my to-do-list for quite a while now, and I have finally gotten around to it. I have been following this tutorial, which includes some basic commands and also some more advanced stuff towards the end. I also created an account on GitHub where I will start adding my solutions to the Rosalind problems.

To try it out, I made a repository and pushed my solution to "Introduction to Alternative Splicing" into it. And immediately ran into a problem. I got the following error message:

! [rejected]        master -> master (fetch first)
error: failed to push some refs to 'https://github.com/ssjunnebo/Rosalind.git'
hint: Updates were rejected because the remote contains work that you do
hint: not have locally. This is usually caused by another repository pushing
hint: to the same ref. You may want to first integrate the remote changes
hint: (e.g., 'git pull ...') before pushing again.
hint: See the 'Note about fast-forwards' in 'git push --help' for details

A quick bit of Google-ing told me that the problem was that when I created the repository on GitHub, I added a README, and since this wasn't available locally on my computer I first needed to merge the GitHub repository with the folder on my computer. I did this using the following command:

> git pull origin master

 I was then able to push my first piece of code onto GitHub!

Expected Number of Restriction Sites

This problem is very similar to "Introduction to Random Strings". As in that problem, we are given a string, s, and an array, A, containing some GC contents. However, in this case we are also given an integer, n, representing the length of a second string, t. t is the random string formed with each given GC content and we are asked to find the probability of finding s as a substring of each t. What we need to realise to solve this problem is that the number of opportunities for finding s in t is equal to n-len(s)+1. To get the overall probability of finding s in t, we simply need to add all the individual probabilities of randomly forming s with the specified GC contents.

Sample Dataset
10
AG
0.25 0.5 0.75

Expected Output
0.422 0.563 0.422

The following is the code I wrote to solve the problem. When I wrote it I became aware that there is a difference in how Python3 and Python2 handles rounding of 0.5. In Python2 it is rounded up, but in Python3 it is rounded down. On an earlier problem I had to run my program using Python2 because the answer I got from Python3 wasn't accepted by Rosalind. In this case however, the sample data set I got yielded the same answer regardless of which version I used (suggesting it no cases of having to round 0.5 occurred in this data set). (Note: if the following code is to be run with Python2 the formating of the output needs to be rewritten for it to work).

data = []
with open('rosalind_eval.txt', 'r') as f:
    for line in f:
        data.append(line.strip('\n'))
n = int(data[0])
s = data[1]
A = [float(x) for x in data[2].split()]

AT, GC = 0, 0
for nt in s:
    if nt == 'A' or nt == 'T':
        AT += 1
    elif nt == 'G' or nt == 'C':
        GC += 1

B = [None]*len(A)
for i, j in enumerate(A):
    P = (((1 - j)/2)**AT)*((j/2)**GC)*(n - len(s)+1)
    B[i] = '%0.3f' % P
print(*B, sep=' ')

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

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)      

Matching Random Motifs

This problem is very similar to “Introduction to Random Strings”. In this case however, we are given a sting s, and are asked to calculate the probability of randomly forming s when generating N nr of stings of equal length of s and a GC content of x.

Sample Dataset
90000 0.6
ATAGCCGA

Sample Output
0.689

To calculate the probability of randomly getting a string that equals s when constructing it with GC content x, we can use the same equation as in introduction to random strings. To go from this probability (lets call it s_prob) to the probability of getting at least one string that is equal to s we use the following equation:

P(at least 1 match of s) = 1 − P(no matches out of N strings) = 1 − [1 - s_prob]^N

The following program makes the above calculations and outputs the answer with three significant figures:

N = 90000
x = 0.6
s = 'ATAGCCGA'
AT = 0
GC = 0
for nt in s:
    if nt == 'A' or nt == 'T':
        AT += 1
    elif nt == 'G' or nt == 'C':
        GC += 1

s_prob = (((1 - x) / 2)**AT) * (((x) / 2)**GC)
prob = 1 - (1 - s_prob)**N
print('%0.3f' % prob)