Beginning Python for Bioinformatics

We are currently following Chapter 4 of Beginning Perl for Bioinformatics, which is the first chapter of the book that actually has code snippets and real programming. The last exercises in this chapter deal with the ability to read files and operate with information extracted from these files, to create arrays and scalar list in perl.

We are going to check how to read files in python. The book tells you how to read protein sequences. Here we are going to read DNA and protein sequences from files and change them.

Let say you have a file with a DNA sequence in some directory in your hard disk. The file cannot be a FASTA type (we will see later how to handle FASTA files), just pure sequence, something like this:

GTGACTTTGTTCAACGGCCGCGGTATCCTAACCGTGCGAAGGTAGCGTAATCACTTGTTC
TTTAAATAAGGACTAGTATGAATGGCATCACGAGGGCTTTACTGTCTCCTTTTTCTAATC
AGTGAAACTAATCTCCCGTGAAGAAGCGGGAATTAACTTATAAGACGAGAAGACCCTATG
GAGCTTTAAACCAAATAACATTTGCTATTTTACAACATTCAGATATCTAATCTTTATAGC
ACTATGATTACAAGTTTTAGGTTGGGGTGACCGCGGAGTAAAAATTAACCTCCACATTGA
AGGAATTTCTAAGCAAAAAGCTACAACTTTAAGCATCAACAAATTGACACTTATTGACCC
AATATTTTGATCAACGAACCATTACCCTAGGGATAACAGCGCAATCCATTATGAGAGCTA
TTATCGACAAGTGGGCTTACGACCTCGATGTTGGATCAGGG

You can download the file here. This is a partial sequence of a mitochondrial gene from a South American frog species called Hylodes ornatus. For our purposes you can save the file in the same directory you are going to run the script from, or if you are using the Python interpreter start it in the directory that contains the file.

The file name is not important but we will use AY162388.seq from now on. The first thing we have to do is to open the file for reading. We define a string variable that will contain the file name.

dnafile = "AY162388.seq"

In order to open the file, we can use the command open, that receives two strings: the first is the file name (it can be the whole location too) to be opened and the mode to be used, which is what you want to do with the file. The mode can be one or more letters that tell the interpreter what to do. For now we are going to use the r mode , which tells Python to read the file, and only do that. So our code is

file = open(dnafile, 'r')

file is a file object that contains the directives to read our DNA sequence file. Now, we have actually read the contents of the file but they are stored in a file object and we did not accessed it yet. We can achieve that by using a myriad of commands. We will start with the commonest one: read the file line by line.

file is our file object. We have just opened it, but Python already knows that any file contains lines (remember that this is a regular ASCII file). We are going to use a loop to read each line of the file, one by one

[sourcecode language='python'>for line in file:
print line[/sourcecode]

In Python, the for loop/statement iterates over items in a sequence of items, usually a string or a list (we will see Python’s list soon), instead of iterating over a progression of numbers. Basically, our for above will iterate over each line in the file until EOF (end-of-file) is reached. Our simple script to read a DNA sequence from a file and output to the screen is

#! /usr/bin/env python
dnafile = "AY162388.seq"
file = open(dnafile, 'r')
for line in file:
    print line

You can download the above script here. To run it have the AY162388.seq in the same directory.

We have seen how to transcribe DNA using regular expression, even though the regex we used cannot be considered a real one. Now we are going to simplify our small script even more and take advantage of some string capabilities of Python. Instead of using two lines, we are going to use only one. And also we won’t need to import anything.

Let’s start again with the same DNA sequence

myDNA = 'ACGTTGCAACGTTGCAACGTTGCA'

This time we are going to use replace. This is one of the Python’s methods to manipulate strings. Basically, we are asking the interpreter to replace a certain string by another. The method returns a new copy of your string:

myRNA = myDNA.replace('T', 'U')

This tells Python: myRNA will receive a copy of myDNA where all Ts were changed by Us. the “dot” after myDNA means that the method replace will get that variable as input on that variable.

So our code from above would like this

#! /usr/bin/env python

myDNA = 'ACGTTGCAACGTTGCAACGTTGCA'
myRNA = myDNA.replace('T', 'U')
print myRNA

Simple and efficient. Next we will use the same approach on generating the reverse complement of a DNA sequence, with no regex pattern.

As mentioned above, regex in Python are provided by the re module, which provides an interface for the regular expression engine. First thing we have to do is to tell the interpreter what to do and what expression to use.

Let’s start with a DNA sequence.

myDNA = 'ACGTTGCAACGTTGCAACGTTGCA'

How to transcribe it to RNA? Transcription creates a single-strand RNA molecule from the double-strand DNA; basically the final result is a similar sequence, with all T’s changed to U’s. So our regular expression has to find all T nucleotides in the above sequence and then replace them.

Regular expressions in Python need to be compiled into a RegexObject, that contains all possible regular expression operations. In our case we need to search and replace, what can be done by using the sub() method. According to the Python’s Regular Expression HOWTO sub() returns the string obtained by replacing the leftmost non-overlapping occurrences of the RE in string by the replacement replacement. If the pattern isn’t found, string is returned unchanged.

Let’s put everything above in real code. First we need to compile a new RegexObject that will search for all thymines in our sequence. It can be achieved by using this:

regexp = re.compile('T')

Simple as that. This line of code tells the Python interpreter that our “regular expression” is every T in our string. Now, we have to make replace those Ts with Us. In order to do that we just tell the interpreter:

myRNA = regexp.sub('U', myDNA)

Let’s look at the last two lines of code. On the first line we created a new RegexObject, regexp (that could have any name, as any variable) and compiled it, making our regular expression to be every T in our string. On the second line, we assigned our soon to be created RNA sequence to a new string (remember that strings in Python are immutable) and used the command sub to replace in the Ts by Us present in our original DNA string. Putting all together our transcription code will be

#! /usr/bin/env python

import re
myDNA = 'ACGTTGCAACGTTGCAACGTTGCA'
regexp = re.compile('T')
myRNA = regexp.sub('U', myDNA)
print myRNA

You can download the resulting script here.

Just an apart from the bioinformatics aspect of programming: Python’s print statement.

As in most computer languages Python allows an easy way to write to the standard output. Python’s print only accepts output of strings, and if the variable sent to it is not a string it is first converted and then output.

The print always put a linebreak ('\n' or "\n") at the end of the expression to be output, except when the print statement ends with a comma. For example:

print "This is a"
print "test"

will print

This is a

test

while

print "This is a",
print "test",

will print

This is a test

Of course Python’s print statement allows any programming escape characte, such as '\n' and '\t'.

Concatenating strings on output

To concatenate two strings on output there are two possible ways in Python. You can either separate the strings with a comma, like we did here

print myDNA, myDNA2

or you can use the “+” sign in roder to obtain almost the same result. This is similar to what was used here

myDNA3 = myDNA + myDNA2

but instead we would use the print command as

print myDNA3 + myDNA

In the latter case, both strings will not be separated by a space and will be merged. To get the same result you would have to concatenate an extra space between the strings like

print myDNA3 + " " + myDNA

Tisdall’s book on Perl introduces next the ability to transcripts DNA sequences into RNA. In order to do that we need to check a different aspect of programming: regular expressions (or regex). Regular expressions is a pattern/string expression that works matching/describing/filtering other strings.

Let’s say you want to examine or extract all vowels contained in one phrase, one page, one word. Another example would be to remove all html tags from a downloaded webpage. As HTML tags are encapsulated between < and > signs we can create a regex that will search for any characters in between the signs and remove (parse) them from our page. We will deal very briefly with regex, and if you are interested in learning more about it you can search for countless references on the internet (such as this one).

In order to use regular expression in Python we need to check another concept present in the language: importing modules. Python core functionality provides most of the usual operations and also comes with a built-in library of functions that "access to operations that are not part of the core of the language but are nevertheless built in". One of this operations is the ability to interpret regular expression that in Python is located in the re module. Apart from the language core, built in modules, Python can be further extended by using third-party modules imported into the language. Anyone can create a module and distribute to every Python user and programmer.

So, in order to have regex capabilities we literally have to import the regex module. We do that by entering the line:

import re

Python's code style guide suggests that import statements should be on separate lines

import sys
import re
 ...

and to be put on the top of the file (usually below the line that tells your OS to use Python's interpreter on the script).

So the first two lines of our new script would be

#! /usr/bin/env python
import re

Now that we have the ability to use regex, we need to create one expression that will transcribe our DNA sequences into RNA.

Another important task for many biologists is to merge/concatenate different strings of DNA in one unique sequence. We can modify the previous script to concatenate two distinct DNA sequences in one.

We start using code_01 structure, adding some extra elements (line 3):

#! /usr/bin/env python
myDNA = "ACGTACGTACGTACGTACGTACGT"
myDNA2 = "TCGATCGATCGATCGATCGA"
print myDNA, myDNA2

So far, we added a new string containing an extra DNA sequence and we print both sequences. In Python print statement automatically adds a new line at the end of the string to be printed, unless you add a comma (,) to the end. The comma is also needed if you are going to print more than one string in order to separate them (try removing the comma from the code above).

Now, how do we merge myDNA and myDNA2? Easy in Python: just sum them with a plus signal:

myDNA3 = myDNA + myDNA2
print myDNA3

Notice that in Python strings are immutable, meaning they cannot be changed. This immutability confers some advantages to the code where strings (in Python strings are not variables) cannot be modified anywhere in the program and also allowing some performance gain in the interpreter. So, in order to have our sequences merged we create a third sequence that receives both strings. Finally our code will be (some captions were added):

#! /usr/bin/env python
myDNA = "ACGTACGTACGTACGTACGTACGT"
myDNA2 = "TCGATCGATCGATCGATCGA"
print "First and Second sequences"
print myDNA, myDNA2
myDNA3 = myDNA + myDNA2
print "Concatenated sequence"
print myDNA3

Easy, eh? Of course these two simple scripts do no scratch the surface of Python programming, but they are a start.

The above code can be downloaded from the repository.

As pointed in Beginning Perl for Bioinformatics, a large percentage of bioinformatics procedures deals with strings, especially DNA and amino acids sequence data. As is largely known DNA is composed of four different nucleotides: A, C, T and G and proteins can contain up to 20 amino acids. Each one of these elements have one letter of the alphabet assigned to them. In the DNA case some letters represent one or more nucleotides that can be present at some sequence position (click here for more ).

So, as the amino acid is the basic building block of proteins (AKA polypeptides), strings containing sequence is our most basic block, from where all the bioinformatics magic will work on.

Usually in Perl a string is represented by the dollar sign in front of the variable name, like this $sequence. Python is dynamically typed, meaning variable types are assigned/discovered by the interpreter at run time. This means that the value after the equal sign will tell the interpreter what variable type you are declaring. So in Python if you want to store a DNA sequence you can just enter

mydna="ACGTACGTACGTACGTACGTACGT"

a quick note: Python can be used with the interpreter command line or by previously saved scripts. I will try to use the latter in the code examples.

OK, we are ready to create our first Bioinformatics Python Hello World script. Let’s get the sequence above and print it on the screen. The first line will tell the operating system what to use and where to find the Python interpreter

#! /usr/bin/env python

Next we will create the variable myDNA and assign the corresponding sequence

myDNA = "ACGTACGTACGTACGTACGTACGT"

And finally, we will print the contents of the variable to the screen:

print myDNA

As mentioned above, Python mandates that you have your code indented, but in our final script this is not needed:

#! /usr/bin/env python
myDNA = "ACGTACGTACGTACGTACGTACGT"
print myDNA

The first line tells your operating system that this is a Python script and to use the interpreter located in that directory; line two declares a variable called myDNA and assigns the sequence string to it and the last line simply output this variable to the screen. That simple!

To run this (extremely simple) script you can copy and paste the code above to your favourite text editor save the file with a .py extension (recommended but not necessary). To run the script, as long as you have Python installed, just open a shell and type on the command line:

> python code_01.py

To download the above script check the repository.