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:

Error Correction in Reads

In this problem we are looking at the errors that occur during genome sequencing. We are given a set of reads of equal length, some of which contain an error (one nt is exchanged). The reads are either of the following:

1. The read was correctly sequenced and appears in the dataset at least twice (possibly as a reverse complement)
2. The read is incorrect, it appears in the dataset exactly once, and its Hamming distance is 1 with respect to exactly one correct read in the dataset (or its reverse complement)

From this dataset we are expected to find the incorrect reads and correct them.

Sample Dataset

>Rosalind_52

TCATC

>Rosalind_44

TTCAT

>Rosalind_68

TCATC

>Rosalind_28

TGAAA

>Rosalind_95

GAGGA

>Rosalind_66

TTTCA

>Rosalind_33

ATCAA

>Rosalind_21

TTGAT

>Rosalind_18

TTTCC

Expected Output

TTCAT->TTGAT

GAGGA->GATGA

TTTCC->TTTCA

I really enjoyed working on this problem and I managed to solve it pretty quickly. A problem I had though was that I would get the correct result for the sample dataset, but when I tried with a 'real' dataset I got way too many reads in my answer, something like 4 times as many as the input reads. I looked at my output and could conclude that a lot of the reads occurred multiple times. I figured something was wrong in the part of my program where I build the list of correct reads, and sure enough, I managed to find the error and fix it (highlighted below). The following code is what I finally ended up with:

from Bio import SeqIO                                                 

from Bio.Seq import Seq                                               

from Bio.Alphabet import generic_dna                                  

reads = []                                                            

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

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

    reads.append(str(record.seq))                                     

handle.close()                                                        

right = []                                                            

wrong = []                                                            

for i, j in enumerate(reads):                                         

    read = Seq(j, generic_dna)                                        

    rev_read = read.reverse_complement()                              

    for k in range(i + 1, len(reads)):                                

        if read == reads[k] or rev_read == reads[k]:                  

            if read not in right and rev_read not in right:           

                right.append(str(read))                               

                right.append(str(rev_read))                           

for l in reads:                                                       

    if l not in right:                                                

        wrong.append(l)                                               

for incorrect in wrong:                                               

    for correct in right:                                             

        hamming = 0                                                   

        for nt1, nt2 in zip(incorrect, correct):                      

            if nt1 != nt2:                                            

                hamming += 1                                          

                if hamming > 2:                                       

                    break                                             

        if hamming == 1:                                              

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

                print(incorrect, '->', correct, sep='', file=textfile)

Completing a Tree

In this problem we are asked to find out the minimum number of edges we need to add to a graph of n nodes described by a given adjacency list in order to produce a tree. I quite quickly managed to overthink the problem and got caught up trying to build the graphs described by the adjacency list, until I finally realised there is a much simpler solution.

Sample dataset:

10

1 2

2 8

4 10

5 9

6 10

7 9

Expected output:

3

The key to this problem is to realise that a tree with n nodes contains n-1 edges. If we then look at the adjacency list, every element of this list represents a node. Thus, to obtain the number of nodes needed to produce a tree, we simply need to make the calculation n-1-elements in the list. The following Python code does just that:

data = []                                          

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

    for line in f:                                 

        split_data = [int(x) for x in line.split()]

        data.append(split_data)                    

n = data[0][0]                                     

edges = data[1:]                                   

print(n - len(edges) - 1)                          

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

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

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

Introduction to Random Strings

We know that genomes are not just a random collection of A,T,C and G, and that across a genome there are loads of different motifs, some of which are similar in many different species. But since the human genome is huge, we need to account for the possibility that the subsequence we are looking at was created at random. In this problem we look at a simplified way to calculate the probability that a certain subsequence occurs randomly. We are given a substring of at most 80 bp and an array containing the GC-content of up to 20 random strings, and we are asked to find log10 of the probability that the substring will match the random strings exactly, and return the result in an array of equal length of the given array.

Sample dataset:

ACGATACAA

0.129 0.287 0.423 0.476 0.641 0.742 0.783

Expected output:

-5.737 -5.217 -5.263 -5.360 -5.958 -6.628 -7.009

The probability, P, of the subsequence occuring in a sequence of a GC content x, can be simplified and written as follows, where AC is the total nr of A and C in the subsequence and GC is the total number of G and C:

To solve this problem, we need to write a program that extracts the GC contents of the array and counts AC and GC of the subsequence. Then it's just a simple matter of iterating over the list of GC contents to calculate the probabilities using the above equation. Here is my final version:

import math

AT = 0

GC = 0

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

    for line in f:

        if line[0] != 'A' and line[0] != 'T' and line[0] != 'G' and line[

                0] != 'C':

            numbers = line.split()

            GC_contents = [float(x) for x in numbers]

        for i in line:

            if i == 'A' or i == 'T':

                AT += 1

            elif i == 'G' or i == 'C':

                GC += 1

probabilities = []

for j in range(len(GC_contents)):

    prob = math.log10((((1 - GC_contents[j]) / 2)**AT) * (GC_contents[j] / 2)

                      **GC)

    probabilities.append('%0.3f' % prob)

print(*probabilities, sep=' ')