Global Alignment with Scoring Matrix

In this problem we are asked to find the alignment score of two protein sequences using the scoring matrix BLOSUM62 and a linear gap penalty of 5. I decided to use Biopython's pairwise2 for this, as you can see in the code below and on my GitHub.

from Bio import pairwise2

from Bio import SeqIO

from Bio.SubsMat.MatrixInfo import blosum62

seqs = []

with open('rosalind_glob.fasta', 'r') as f:

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

        seqs.append(record.seq)

s = seqs[0]

t = seqs[1]

alignments = pairwise2.align.globalds(s, t, blosum62, -5, -5)

print(pairwise2.format_alignment(*alignments[0]))

Note that the gap open penalty and the gap extension penalty are both set to -5, as we are told in the problem description to use a linear gap penalty of 5. 

Newick Format with Edge Weights

If you have already solved "Distances in Trees", then this problem will probably not take that long to complete. When I solved distances in trees, I used the Biopython module Phylo. The problem I had then was that the given trees did not contain any branch lengths, so I had to add this in the program. However in this problem the branch lengths are already included in the trees, so the new program is in fact even shorter than the original. You can find my new version of the code below and on my GitHub.

import sys
from Bio import Phylo
import io

f = open('rosalind_nkew.txt', 'r')
pairs = [i.split('\n') for i in f.read().strip().split('\n\n')]

for i, line in pairs:
    x, y = line.split()
    tree = Phylo.read(io.StringIO(i), 'newick')
    sys.stdout.write('%s' % round(tree.distance(x,y)) + ' ')
sys.stdout.write('\n')

Creating a Character Table (Finished!)

The approach I wrote about in my last post turned out to be a lot harder to implement than I initially thought, and I actually never made it past the step of splitting the tree into the possible sub-trees. After bashing my head against the problem a bit too long, I decided to try a different approach, using Biopython. The Phylo module of Biopython has a method for parsing trees in Newick format, and by using this I was also able to extract internal and external nodes of the tree in a much easier way than in my initial approach. The final version of my code can be found here on my GitHub, along with the problem description and some sample data files.

Distances in Trees

In this problem we are looking att the Newick format and how to find the distance between two nodes in a phylogenetic tree. We are given a file containing trees in Newick format and two nodes for each tree, and are asked to find the distance between those nodes.

Sample Dataset
(cat)dog;
dog cat

(dog,cat);
dog cat

Expected Output
1 2

I remember working with the Newick format before, in one of the bioinformatics courses I took, so when I started working on this problem I recalled that there were functions for the format available in Biopython. So I had a look at the documentation, and sure enough, there is a function called distance that would be suitable. However, as always, there was a slight problem. The trees given by Rosalind did not contain any branch lengths, which is what the distance function uses to calculate the distance between two nodes. To enable using this function I therefor had to assign the branches a length of 1 (done on rows 18-21). The following code (also available on Github here) yielded a result accepted by Rosalind:

import sys
from Bio import Phylo
import io

#open file and parse data
f = open('rosalind_nwck.txt','r')
pairs = [i.split('\n') for i in f.read().strip().split('\n\n')]

#for each pair:
#-parse data further with biopython
#-add branch length 1 to all branches
#-use bioputhons Phylo distance funktion to get distances
#-print result on requested format

for i, line in pairs:
    x,y = line.split()
    tree = Phylo.read(io.StringIO(i),'newick')
    clades = tree.find_clades()
    for clade in clades:
        clade.branch_length = 1
    sys.stdout.write('%s' % tree.distance(x,y) + ' ')
sys.stdout.write('\n')

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!

Calculating Protein Mass

In this problem we are asked to calculate the molecular weight of a peptide. The peptide is assumed to come from the middle of a protein so we are to assume that it consists entirely of amino acid residues (meaning we don't have to account for the extra weight of the "water molecule" present when we include the two edges of the protein).

I wrote two programs for this problem, the first one uses Biopython and the second doesn't. The following code is the one using Biopython:

from Bio.Seq import Seq                  

from Bio.Alphabet import generic_protein 

from Bio.SeqUtils import molecular_weight

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

    for line in f:                       

        prot_seq = line.strip('\n')      

print('%0.3f' % (molecular_weight(       

    Seq(prot_seq, generic_protein),      

    monoisotopic=True) - 18.01056))      

This program uses the Biopython function molecular_weight from SeqUtils. The function sets monoisotopic=False by default, but because the problem specified that we should use the monoisotopic weitghts we need to set it to monoisotopic=True. Also, the function includes the extra weight of one water molecule, so we need to remove that manually (hence the - 18.01056). Come to think of it there is an option called circular that states that the sequence has no ends. Maybe that would work as well?

Anyhow, I thought I'd try to make this program without the use of Biopython as well, and the following is what I came up with:

weights = {'A': 71.03711,             

           'C': 103.00919,            

           'D': 115.02694,            

           'E': 129.04259,            

           'F': 147.06841,            

           'G': 57.02146,             

           'H': 137.05891,            

           'I': 113.08406,            

           'K': 128.09496,            

           'L': 113.08406,            

           'M': 131.04049,            

           'N': 114.04293,            

           'P': 97.05276,             

           'Q': 128.05858,            

           'R': 156.10111,            

           'S': 87.03203,             

           'T': 101.04768,            

           'V': 99.06841,             

           'W': 186.07931,            

           'Y': 163.06333}            

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

    for line in f:                    

        prot_seq = line.strip('\n')   

weight = 0                            

for aa in prot_seq:                   

    weight += weights[aa]             

print('%0.3f' % weight)               

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)                                              

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.

Translating RNA into Protein

In this task we are asked to translate a RNA-sequence into its corresponding protein sequence. To write this program 'from scratch' would probably take me a lot of time, mainly because there are 64 different codons available when translating RNA into proteins. Instead, there is an excellent function in Biopython for doing this, so I decided to save myself some time and effort and just use this instead. The following is the final code that i came up with:

from Bio.Seq import Seq              

from Bio.Alphabet import generic_rna 

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

    data = f.read()                  

rnaseq = Seq(data, generic_rna)      

protein = rnaseq.translate()         

print(protein)                       

I think this approach was a lot faster than if I had set to write this program from scratch. At this point I feel like I have a good grasp of the basics of Python, and I will continue working more with Biopython when applicable. 

Computing GC Content

This problem felt like a good opportunity to recap my skills in using Biopython. It asks you to calculate the GC-content of some sequences in a FASTA file and print the highest CG-content along with the corresponding sequence ID. To do this, the program first needs to parse the FASTA file, then it needs to calculate the GC-content of the strings, and finally it needs to print the largest GC-content, coupled with the correct ID.

I initially had some trouble getting Biopython to work with Python 3.4, but it turned out that I had installed it to Python 2.7, and not 3.4. When this problem was solved, I set to writing my program.

As I mentioned, the first thing I needed my program to do was to parse the FASTA file. To do this, Biopython has a function called SeqIO. I read up on it here. I then wrote this piece of code which reads the file, calculates the GC-content of the sequences in the file and prints it together with the corresponding ID.

from Bio import SeqIO                      

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

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

    count = 0                              

    totalcount = 0                         

    print(record.id)                       

    for nt in record.seq:                  

        totalcount = totalcount + 1        

        if nt == 'G' or nt == 'C':         

            count = count + 1              

    percent = count/totalcount*100         

    print(percent)                         

handle.close()                             

This worked fine, but what I needed to do now was to somehow save the highest GC-content, coupled with its ID, and print only that at the end of the program. I thought this was going to be tricky, but it turned out to be quite simple. In the following code, which is the final version of the program, I simply added an if-statement to overwrite the variables 'GC' and 'ID' every time a higher GC-content is found.

from Bio import SeqIO                      

GC = 0                                     

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

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

    count = 0                              

    totalcount = 0                         

    for nt in record.seq:                  

        totalcount = totalcount + 1        

        if nt == 'G' or nt == 'C':         

            count = count + 1              

    percent = count/totalcount*100         

    if percent > GC:                       

        GC = percent                       

        ID = record.id                     

print(ID)                                  

print(GC)                                  

handle.close()                             

Of course, I could have just used the built-in function GC in Biopython instead, but where is the fun in that?