Phylogenetics: Python Primer

From eebedia

Revision as of 14:27, 17 February 2011 by PaulLewis (Talk | contribs)
Jump to: navigation, search
Adiantum.png EEB 349: Phylogenetics
This lab represents an introduction to the Python computing language, which will be useful for upcoming homework assignments as well as for software written as Python extensions (i.e. Phycas)

Contents

What is Python?

Python logo.gif
Python is one of two programming languages that we will use this semester (the other being R). One might make the case that programs like PAUP that are associated with a unique command language also represent programming languages; however, Python and R are different in being general purpose (i.e. not written specifically for phylogenetics).

Python is one of a number of different high-level computing languages in common use. All of the software we use is written in one of these languages. PAUP is written entirely in the C language, while SplitsTree is written in Java. Knowing a little bit of computer programming can save you immense amounts of time by allowing you to automate things. You will begin realizing these savings doing homework for this class. While it is possible to do all the homework assignments by hand using a calculator, you will find that using Python will save you time and is more accurate (a mistake on a calculator early on in a calculation can be very costly in terms of time and accuracy). Python is a good language to learn first because it is relatively simple (not many words or punctuation rules to learn) and is much more forgiving than other languages. It is in a class of languages known as scripting languages because the program is interpreted as it is read - languages such as C require two additional steps (compiling and linking) before they can be run.

Installing Python


Mac logo.jpg
If you have a Mac, then Python is already installed on your computer because it is used in many parts of the Mac OS X operating system. To verify this, and to see what version of Python you have, start your Terminal program (in the Applications/Utilities folder), then type the following at the unix prompt (note: the prompt is represented by a dollar sign ($), so you should only type the "python -V" part):
$ python -V
Python 2.5.1

Win logo.png
If you have Windows, you will need to download and install Python from the http://www.python.org/download/ web site. Note that there are two current versions of Python. You should install the older one (Python 2.7.1). (We would use Python 3.1 except that some Python software that we will use later in the semester is not yet compatible with Python 3.1, so it makes sense to stick with 2.7 for now.) There are several different types of downloads: choose Python 2.7.1 Windows installer (Windows binary -- does not include source) and follow the directions once you start it up.

Once you have Python installed, you can invoke it from a console window (in Vista, use Start, All Programs, Accessories, Command Prompt) by typing python. The -V option causes Python to print out its version number and then quit:

C:\Users\Administrator>python -V
Python 2.7.1

You should download the documentation from http://docs.python.org/download.html so that it can be accessed quickly. I find the HTML form best: I unpacked the zip file and bookmarked it in my browser as follows:

file:///Users/plewis/Documents/Manuals/python-html/index.html

Python basics

This is the briefest of introductions, designed to get you just to the point where you can do the majority of your homework assignments using Python. If you get stuck trying to write a Python program, I have found the Tutorial and Global Index to be the most useful parts of the Python documentation.

Starting Python

In a console window (Windows) or in the Terminal application (Mac), start Python by typing python at the operating system's prompt. Here is what it looks like on my Mac (the part I typed is in bold):

paul-lewiss-macbook-pro:Phylogenetics plewis$ python
Python 2.5.1 (r251:54863, Apr 15 2008, 22:57:26) 
[GCC 4.0.1 (Apple Inc. build 5465)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>>

Kinds of information you can store in Python variables

Try typing in the Python code presented in the example sessions below. The >>> represents the Python prompt (you will see this when you start Python), so don't type that! The output is shown below each Python statement that generates output (don't type that either!). Finally, everything after a # character represents a comment (while you can type these in, it would probably be a waste of your time).

Integers, floats and strings

Strings are series of characters. Assign the string 'Have a nice day' to the variable s:

>>> s = 'Have a nice day'

Integers are whole positive or negative numbers. Assign the integer 9 to the variable i:

>>> i = 9

Floats are numbers with an implicit or explicit decimal point. Assign the float 9.5 to the variable f:

>>> f = 9.5

Later on you will learn how to read data from a file. When numbers are read from a file, they start out as strings. If you want these strings to be used as numbers, you will need to convert them to integers or floats. Here is how to convert a string ('5') to an integer:

>>> x = int('5')
>>> x
5

Note that putting a variable name on a line by itself causes its value to be printed to the console. Here is another attempt to convert a string to an integer. This one fails because the string '5.5' is not able to be converted to an integer:

>>> x = int('5.5')
Traceback (most recent call last):
...
ValueError: invalid literal for int() with base 10: '5.5'

Here is how to convert a string to a float:

>>> y = float('5.5')
>>> y
5.5

Lists and tuples

Lists are collections of integers, floats, strings, etc. Tuples are like lists, except that you cannot change elements of a tuple. Here is a list consisting of the integer, the float and the string we just defined:

>>> L = [i, f, s]
>>> L
[9, 9.5, 'Have a nice day']

Change the first (0th) value in the list to 5. Note that numbers begin with 0 in most programming languages, and Python is no exception.

>>> L[0] = 5
>>> L
[5, 9.5, 'Have a nice day']

Change the second value in the list to 13:

>>> L[1] = 13
>>> L
[5, 13, 'Have a nice day']

Change the last value in the list to 'Ok':

>>> L[2] = 'Ok'
>>> L
[5, 13, 'Ok']

Change the last value in the list to "G'day":

>>> L[-1] = "G'day"
>>> L
[5, 13, "G'day"]

This illustrates two things. First, the elements can be accessed by counting backward from the number of elements. Hence, L[-1] refers to the same element as L[3-1] = L[2] (i.e., the last element). Second, you can embed an apostrophe inside a string by using double quotes to surround the string, and you could embed double quotes inside a string by using single quotes to surround it.

You can find out how many elements are in the list using the len() function:

>>> len(L)
3

Here is a tuple consisting of the integer, the float and the string we just defined

>>> t = (i, f, s)
>>> t
(9, 9.5, 'Have a nice day')

Try changing the first element of the tuple (this will generate an ugly error message because tuples cannot be modified):

>>> t[0] = 5
Traceback (most recent call last):
    ...
TypeError: 'tuple' object does not support item assignment

To eliminate possible ambiguities, every tuple must have at least one comma, even if it only has just one element!

>>> t = (1.2,)

Arithmetic

Below are some examples of basic arithmetic operations in Python.

Addition

Add 1 plus 2:

>>> x = 1 + 2
>>> x
3

Multiplication

Multiply 2 and 3:

>>> x = 2*3
>>> x
6

Exponentiation

Raise 2 to the power 3:

>>> x = 2**3 
>>> x
8

Division

Divide 6 by 2:

>>> x = 6/2
>>> x
3

Order of operations

In the following expression, the multiplication is done before the division because operations with equal precedence are evaluated from left to right:

>>> x = 6*7/3
>>> x
14

If you want the division to be done first, use parentheses:

>>> x = 6*(7/3)
>>> x
12

Note that 7/3 yields 2, not 2.33333 because an integer divided by another integer yields an integer, and integers cannot have a fractional component (the value is truncated, not rounded). If this is not what you wanted, you need to use floats, which you can do by inserting a decimal point after each whole number:

>>> x = 6.*(7./3.)
>>> x
14.0

Important: This last section is extremely important! To be safe, I usually use parentheses and decimal points even if I do not need to use them.

Precedence

Some operations are evaluated before others due to their higher precedence. For example, exponentiation has higher precedence than multiplication or division, so the ** is done before the multiplcation here:

>>> x = 2*2**3
>>> x
16

But why not use parentheses to make the calculation order explicit so you can be sure what will be done first:

>>> x = 2*(2**3)
>>> x
16

Multiplication and division, in turn, have higher precedence than addition or subtraction. For the full operator precedence table, scroll to the bottom of the expressions documentation

Sum, min and max

The sum function makes it easy to sum the elements of a list:

>>> x = sum([1,2,3,4])
>>> x
10

The min and max functions make it easy to find the extreme values in a list:

>>> z = [1,2,3,4]
>>> x = max(z)
>>> x
4
>>> x = min(z)
>>> x
1

Using what you know

You already know enough Python to start using it to do your homework for this week.

Computing p-distances from numbers of differences

The first step in your homework for this week is to compute p-distances from the supplied numbers of differences. Here is one way to do this in Python:

>>> n = 3424   # store sequence length in variable n
>>> x = 293    # store number of differences between taxon1 and taxon2 in x
>>> p = x/n    # do the division
>>> p
0

Ok, you know that's not the right answer, so what went wrong? We divided one integer by another integer to produce a third integer, but instead of storing the number 0.085572429906542055 in p, it stored 0 because that is the whole number part of 0.085572429906542055. Again (this bears repeating many times), It is always a good idea to be explicit when you might get a float for an answer:

>>> p = float(x)/float(n)
>>> p
0.085572429906542055

Note that the trick of including the decimal point doesn't work on variables, only on actual numbers:

>>> p = x./n.
  File "<stdin>", line 1
    p = x./n.
          ^
SyntaxError: invalid syntax

You can, however make the variables involved into floats when you first assign values to them. The following solution uses decimal points to save some typing and will prevent you from being surprised by later calculations involving n and/or x:

>>> n = 3424.
>>> x = 293.
>>> p = x/n
>>> p
0.085572429906542055

This time p is a float because you divided one float by another.

Computing a JC distance from a p-distance

The next step is to convert your p-distance to a JC distance. Now we need some capabilities (e.g. the ability to take the logarithm of a number) that are not present by default when Python starts up. To solve this deficiency, use the import statement to bring in the needed functionality:

>>> import math
>>> jc = -0.75*math.log(1.0 - 4.0*p/3.0)
>>> jc
0.090860500068600705

In the import math statement, the math part is known as a Python module. To see a list of all modules (and documentation for them), check out the Global Module List in the Python documentation.

Loops

An important component of any programming language is the ability to do the same sort of thing many times. In the session that follows, you will create a list containing all six pairwise differences among sequences, then use a for loop to compute the JC distance for all six. Important: use a tab to indent. Python is very sensitive about indenting. If the three indented lines are not indented by exactly the same amount, Python will spit out an error message.

>>> n = 3424.     # note that you really just need the period to make it into a float
>>> x = [293.0, 277.0, 328.0, 268.0, 353.0, 353.0]
>>> jc = []       # start with an empty list
>>> for i in range(6):
...     p = x[i]/n                           # compute p-distance using the ith value in the x list
...     d = -0.75*math.log(1.0 - 4.0*p/3.0)  # compute the JC distance and store in variable d
...     jc.append(d)                         # lengthen the list jc by appending the new value d
...                                          # just hit return here to exit the loop body
>>> jc
[0.090860500068600705, 0.085604236767267153, 0.1024886399705747, 0.082663697169458011, 0.11090623674796182, 0.11090623674796182]

You may be curious about range(6). The range(n) function produces a list of integers starting with 0 and ending with n-1. Try typing the range function by itself to verify this:

>>> range(10)
[0,1,2,3,4,5,6,7,8,9]

Your for loop could therefore have been written like this and the result would be exactly the same:

for i in [0,1,2,3,4,5]:

Now, what does the line above mean? The for loop lets i take on each value in the list, in turn. Thus, the indented portion of the loop is executed 6 times, once for each value in the list.

Python script files

Now that your Python constructs are getting a little longer, it is a good time to learn about creating files for your Python programs. A file containing Python statements is called a Python script. It is a good idea to have a good text editor before starting to create scripts, so your next goal will be to download one if you do not already have one installed.

Note: Microsoft Word is not a good text editor! It is an excellent word processor, but text editors and word processors are different beasts. Text editors always save files as plain text, while Word saves the file in its proprietary file format that Python cannot read. It is possible to save Word files as plain text, but usually this is more trouble than it is worth.

Install a text editor

Fortunately, there are free (and really good) text editors for both Windows and Macs.


Mac logo.jpg If you have a Mac, download TextWrangler. If you really get into scripting, you may want to purchase the beefier version known as BBEdit, produced by the same company, but TextWrangler will suffice for this course. Once you get TextWrangler installed, start it up and create a file with the contents shown below (after the instructions for Windows users). Save the file using the name first.py in a convenient location (e.g. Documents/scripts), then navigate to that folder in your Terminal window (using the command cd $HOME/Documents/scripts)


Win logo.png If you have Windows, download Notepad++. Once you get Notepad++ installed, start it up and create a file with the contents shown below. Save the file using the name first.py in a convenient location (e.g. C:\scripts), then navigate to that folder in your command console (using the command cd C:\scripts)


Your first Python script

Here is what you should type into your new first.py file:

import math
n = 3424.
x = [293.0, 277.0, 328.0, 268.0, 353.0, 353.0]
jc = []
for i in range(6):
    p = x[i]/n
    d = -0.75*math.log(1.0 - 4.0*p/3.0)
    jc.append(d)
print jc

Note that I have added the word print on the last line. This is because our little trick of getting Python to tell us the value of a variable by placing the variable's name on a line by itself only works when you are using Python interactively. We are now switching to programming mode, where an entire script is given to the Python interpreter, and a print statement must be used in this context. The result should be the same.

Try running your script by typing the following at your operating system prompt (i.e. get out of Python if you are already in it by typing Ctrl-d if you are a Mac user or Ctrl-z if you are a Windows user):

python first.py

Fancier print statements

The program first.py spits out six JC distance values, but it is usually good to format the output so that it is clearer than this. Copy the following into a new text file named second.py and try running it:

import math
n = 3424.
x = [293.0, 277.0, 328.0, 268.0, 353.0, 353.0]
for i in range(6):
    p = x[i]/n
    d = -0.75*math.log(1.0 - 4.0*p/3.0)
    print '%12d %12.5f %12.8f' % (i, p, d)

You will note a couple of differences between first.py and second.py. First, there is no jc list anymore, we just print out values as they are computed. Second, the print statement is much more complicated now! The complexity might be daunting at first, but you will quickly become adjusted to it I think. Here is how the print statement used in second.py works. There are two parts separated by a percent symbol. The first part is a string:

'%12d %12.5f %12.8f'

while the second part is a tuple:

(i, p, d)

The string serves as a format specification. Here is a breakdown:

  • %12d says print an integer (d is code for integer here) using 12 spaces
  • %12.5f says print a float (f is code for float here) using 12 spaces total, with 5 of the 12 being devoted to the part after the decimal point
  • %12.8f says print another float using 12 spaces, this time with 8 of the 12 being after the decimal point

The tuple provides the values to insert: the integer stored in the variable i will go in the first spot, while the floats stored in p and d will go in the second and third spots. Be sure to use "d" for integers and "f" for floats in your format string, otherwise Python will complain.

Run second.py and see if Python spaced everything as you expected. Note that the second float will actually take up 13 spaces because there is a single space in the format string just before the %12.5f specification. Spaces in the format string are inserted as is into the print output.

Computing the Sum-of-Squares

The next step in your homework assignment involves computing the sum-of-squares for each of the two trees provided. The sum-of-squares formula given in lecture (for the case of power=2) is:

SS = \frac{(d_{ij} - p_{ij})^2}{d_{ij}^2}

Let's create a Python script that will spit out the numbers we need for the table labeled "SS for Tree 1" in the homework assignment. Here is the script (the parts in bold will be discussed in more detail below); save it as ss1.py:

import math
n = 3424.
rowname = ['1 vs. 2', '1 vs. 3', '1 vs. 4', '2 vs. 3', '2 vs. 4', '3 vs. 4']
path = [0.3, 0.3, 0.2, 0.2, 0.3, 0.3]
x = [293.0, 277.0, 328.0, 268.0, 353.0, 353.0]
totalSS = 0.0
for i in range(6):
    pdist = x[i]/n
    dij = -0.75*math.log(1.0 - 4.0*pdist/3.0)
    pij = path[i]
    numerator = (dij - pij)**2.0
    denominator = dij**2.0
    SS = numerator/denominator
    totalSS += SS
    print '%12s %12.5f %12.5f %12.5f' % (rowname[i], pij, dij, SS)
print 'total SS = %.5f' % totalSS

Here are the new additions:

  • A rowname variable has been added to which is assigned a list of strings. Each element of rowname will serve as the label for a row in the table that is output
  • A path variable has been added to which has been assigned a list of path lengths through Tree 1. Each element of path is the path length for the corresponding comparison named in rowname
  • A variable named totalSS was added to keep track of the sum of the individual pairwise SS values as they are calculated. Note that we have to initialize the value of totalSS to zero before the loop begins
  • I've introduced a pij variable to hold the current path length
  • The numerator and denominator variables hold the top and bottom parts of the SS calculation. The double-asterisk means "raised to the power of"
  • The SS variable holds the sum-of-squares value for pairwise comparison i
  • The totalSS += SS tells Python that we want to add the current value of SS to the running sum being stored in the variable totalSS
  • The last line prints the value of totalSS. Note that if there is only one value to be substituted in a format string, you needn't make it into a tuple (i.e. no parentheses are needed around totalSS)

Try running this script and let me know if you don't understand how it works.

Files

One of the most fundamental tasks in programming is to read information stored in an input file and write the program's output to an output file. Here are the absolute basics of doing this in Python.

Reading data from a file

Suppose you wanted to read in the six values that we've been storing in the variable x from a file. To experiment with this, we'll need a file. Create a file named ss1.txt using your text editor that contains the six numbers, one per line. It should look like this:

293.0
277.0
328.0
268.0
353.0
353.0

To read these values, add these five lines to your ss1.py file and comment out the line that currently defines x (place a pound sign # in front of the x):

lines = open('ss1.txt','r').readlines()
x = []
for line in lines:
    f = float(line)
    x.append(f)
# x = [293.0, 277.0, 328.0, 268.0, 353.0, 353.0] 

The first line above opens a file named ss1.txt for reading (this is what the 'r' means). The readlines() part causes the entire contents of the file to be read, creating a list in which each element is a separate line from the file. This list is stored in the variable lines.

The for loop visits each element in lines one at a time. This represents a different (more convenient) style of for loop that doesn't require us to know how long the lines list is (that is, this style avoids having to use range(n), which would require us to know n). Inside the for loop, the string read from one particular line in the file is stored in the variable line. The string in line is converted to a float named f, which is then appended to the growing list x. After all is said and done, the variable x will be exactly the same (a list of six float values) that it was before when we defined it directly.

Note that we had to explicitly convert each line into a float value. This is because each line of information read from the file is stored as a string, not a float. If you failed to explicitly convert each value, you would later run into trouble when you tried to do arithmetic with the values stored in x.

Writing output to a file

Writing to a file involves creating a Python variable that represents the file, then calling that variable's write function to spit out something to the file. Create a file for writing (note the 'w' for writing rather than 'r' for reading) as follows:

outfile = open('output.txt', 'w')

You can now write something to it like this:

outfile.write('Hello\n')

One big difference between writing to a file versus printing to the terminal is that the write function does not terminate the line. The end-of-line character is designated using the two-character combination \n. Thus, if you wanted two blank lines between the word "Hello" and "Goodbye" you would do this:

outfile.write('Hello\n\nGoodbye')

When you are finished writing to the file (say, as the last line of your script), be sure to close the file:

outfile.close()

A revised script

import math
n = 3424.
rowname = ['1 vs. 2', '1 vs. 3', '1 vs. 4', '2 vs. 3', '2 vs. 4', '3 vs. 4']
path = [0.3, 0.3, 0.2, 0.2, 0.3, 0.3]

# Here is where we read in the x values from a file
lines = open('ss1.txt','r').readlines()
x = []
for line in lines:
    f = float(line)
    x.append(f)

# Here we create an output file
outf = open('output.txt', 'w')

totalSS = 0.0
for i in range(6):
    pdist = x[i]/n
    dij = -0.75*math.log(1.0 - 4.0*pdist/3.0)
    pij = path[i]
    numerator = (dij - pij)**2.0
    denominator = dij**2.0
    SS = numerator/denominator
    totalSS += SS
    outf.write('%12s %12.5f %12.5f %12.5f\n' % (rowname[i], pij, dij, SS))
outf.write('total SS = %.5f\n' % totalSS)
outf.close()
Personal tools