Python for Bioinformatics: simple Python

Showing posts with label simple Python. Show all posts

Monday, April 23, 2012

Channeling Lincoln

I found this on reddit. The post says to consider the first part of the Gettysburg address:

Four score and seven years ago our fathers brought forth upon this continent a new nation, conceived in liberty and dedicated to the proposition that all men are created equal

How long do you think it would take you to manually break the address up into 13 character lines, breaking those lines at appropriate spaces?

I'm not so concerned with the estimation process. This is a fun problem for a beginning Python coder. My solution is below, but please try it yourself first. And yes, I was pleased with my time, but I've done this before. How to handle the beginning and the end is the key. I just make them special cases.

[ UPDATE: You could also just look for the right function in the standard library! ]

> python script.py 
10 Four score
 9 and seven
13 years ago our
 7 fathers
13 brought forth
 9 upon this
11 continent a
11 new nation,
12 conceived in
11 liberty and
12 dedicated to
 3 the
11 proposition
12 that all men
11 are created
 5 equal

FH = open('data.txt','r')
data = FH.read().strip().split()
FH.close()

pL = list()
N = 13
for word in data:
    assert len(word) <= N
tL = [data.pop(0)]

for word in data:
    n = len(' '.join(tL)) 
    n += 1 + len(word)
    if n > N:
        pL.append(' '.join(tL))
        tL = [word]
    else:
        tL.append(word)

pL.append(word)

for line in pL:
    print "%2d" % len(line), line

Saturday, April 21, 2012

Testing

It seems improbable, but it is a fact that on the first exam in some courses, a few students may score lower than expected for guessing at random. We usually see some scores of less than 20 out of 100 for an exam whose questions have five possible answers ABCDE.

We can turn this into an illustration of statistical significance. How high must a score be before one can dismiss the null hypothesis (no knowledge of the subject) with a p-value of 0.05?

We have a Bernouli trial with p = 0.2.

The expected value or mean is np. There is a beautifully simple derivation of the mean and variance in the wikipedia article, the variance is np(1-p). For our example, the mean is 20 and the variance is 100 x 0.2 x 0.8, so the standard deviation is √16 = 4.

For a large class (N = 100) of randomly answering students, the distribution appraches normal, and we expect that 2.5% of such students will score at or above two standard deviations above the mean. To make our rate for a type I error less than 0.05, we should require a score of 28 or higher.

I wrote a short script to simulate this. Each student takes a test with 100 questions, and each trial consists of a class of 100 students. We count the number of students with scores at or above a threshold T. The results for 10 trials are printed, and the results for 100 trials are averaged.

A few quick points: the mean appears to fall between 20 and 21. I'm sure there's a simple explanation but it eludes me at the moment. The score that 2.5% of students regularly exceed is 28.

This is a great example of the danger of multiple hypothesis testing!

One last thing is that the code is not very smart. What I should have done was extract the statistics for each value of T for each trial, rather than repeating the trial for each value of T. But it runs fast enough, and this way I don't need a data structure to remember the results.

> python script.py 
20:  56 43 54 57 57 54 57 46 45 64 55.0
21:  39 48 42 44 39 44 48 48 50 48 43.5
22:  43 32 39 30 38 35 36 42 33 41 34.6
23:  40 29 35 31 27 19 30 20 31 25 26.0
24:  19 18 24 26 20 22 21 17 15 24 19.2
25:  11 14  8  8 13 11 12 15 17 13 12.9
26:  12  6  6  6 11  4  8  7 10  6  9.2
27:   7  4  6  5  8  7  5  7  7  7  5.9
28:   3  2  4  2  3  1  2  3  5  1  3.6
29:   4  4  2  2  2  1  3  1  0  1  2.2
30:   1  1  1  0  1  1  1  0  0  1  0.9
31:   4  0  1  0  0  0  0  0  0  0  0.6
32:   0  0  0  0  1  0  0  1  0  0  0.3
33:   0  0  0  0  0  1  0  0  0  1  0.1
34:   0  0  0  0  0  0  0  0  1  0  0.1

script.py

import random
    
def student(N=100, p=0.2):
    L = [random.random() <= p for i in range(N)]
    return sum(L)

def trial(T):
    L = [student() for i in range(100)]
    return sum([1 for n in L if n >= T])
    
def mean(L):
    return 1.0*sum(L)/len(L)

for T in range(20,35):
    L = list()
    print '%2d: ' % T,
    for i in range(100):
        L.append(trial(T))
        if i < 10:
            print '%2d' % L[-1],
    print '%4.1f' % mean(L)

Sunday, March 18, 2012

iTunes backup playlist

Since they removed the DRM from music downloads at the iTunes store, I've been buying songs from Apple using that approach. It's great for impulse buying. Also, it's nice in combination with an app for the iPhone called Shazam, which can identify a track that is playing wherever you are at the moment. I thought I was fairly modest in my purchases, but now find that I have over 500 songs after 18 months. So the question then is how to archive these purchases so they don't go **kabluie** in the night.

One strategy is to just let Apple do it, but I don't trust them.

The next idea is to manually copy all 500 songs to a backup disk. This would be easy, except that (i) I don't want the whole library, just a playlist, and (ii) iTunes uses nested folders to preserve the artist:album:songtitle information. I need to merge u/v/w with u/v/x to create a directory u/v containing both w and x. You can't do that just copying x.

The playlist info is contained in xml format, file: 'iTunes Music Library.xml' but I decided to export the playlist to disk from within iTunes ('File > Library > Export Playlist...'). Each entry in the exported text file ends with something like this as the last field (tab-separated, one entry per line):

HD:Users:Shared:iTunes Music:Bob Marley:Natty Dread:01 Lively Up Yourself.m4a

The following script finds each song on the playlist and copies it to a directory temp in the directory where the script is run. Simple. The result is 3.77 GB of .m4a files with the desired directory structure.

A detail that is (or should be) embarassing: iTunes (on OS X) uses '\r' (CR) as newline. Talk about the constraints of backward compatibility.

One last thing: please, please, please back-up and test before using this. YMMV. Caveat lector. No warranty express or implied. Don't blame me if your library vanishes.

import os, subprocess

name = 'playlist'
FH = open(name + '.txt', 'r')
data = FH.read().strip()
FH.close()

# iTunes uses '\r' (CR) as newline!
data = data.split('\r')

# data[0] is metadata (column names)
data.pop(0)

# file path is the last value
L = [item.split('\t')[-1] for item in data]

# ':' is path separator
L = [item.replace(':','/') for item in L]

for item in L:
    # remove HD name from file path
    item = item.split('/', 1)[1]
    
    # file path has spaces
    # must be quoted for shell command below
    src = '"/' + item + '"'

    artist, album, songfile = item.split('/')[-3:]
    # construct directory tree if it doesn't exist
    path = '/'.join(('temp', artist, album))
    try:
        os.stat(path)
    except OSError:
        os.makedirs(path)
    
    dst = '"' + '/'.join((path, songfile)) + '"'    
    cmd = ' '.join(('cp', src, dst))
    obj = subprocess.call(cmd,shell=True)
    if obj != 0:
        print 'e',
    else:
        print '*',
    print dst

Wednesday, August 31, 2011

Binary multiplication, again

I had a lot of fun this afternoon, but all along I knew I was "reinventing the wheel." It's easy when you know where you're going.

I implemented binary multiplication with a class that keeps its data as Python strings like '00100011'. It's probably because I have been re-reading Charles Petzold's book Code, which I think is a classic. In any event, I wrote the class "word" which allows objects to add, multiply and exponentiate themselves. The second file is test.py which exercises the class.

It's a great "simple Python" coding project. We do bit-shifting by hand. If you want to be slick, you could go further by trying to implement fast exponentiation. Here is another example by a guy who obviously knows what he's doing.

A bit of sample output:

> python test.py
x     467
y     327
x + y 00000000 00000000 00000011 00011010 = 794
check 467 + 327 = 794
x * y 00000000 00000010 01010100 10000101 = 152709
check 467 * 327 = 152709
z      3
x**z  00000110 00010010 00010010 00001011 = 101847563
check 467 ** 3 = 101847563

x     995
y     402
x + y 00000000 00000000 00000101 01110101 = 1397
check 995 + 402 = 1397
x * y 00000000 00000110 00011010 01110110 = 399990
check 995 * 402 = 399990
z      3
x**z  00111010 10110111 00001100 10111011 = 985074875
check 995 ** 3 = 985074875

x     897
y     490
x + y 00000000 00000000 00000101 01101011 = 1387
check 897 + 490 = 1387
x * y 00000000 00000110 10110100 11101010 = 439530
check 897 * 490 = 439530
z      4
x**z  10111011 11001001 10001110 00000001 = 3150548481
ran into a slight problem

word.py

class word:
    SZ = 32
    def __init__(self,n):
        if type(n) == type('a'):
            b = n.rjust(word.SZ,'0')
        elif type(n) == type(1):
            b = bin(n)[2:]
        else:
            raise ValueError, "can't do that"
        self.b = b.rjust(word.SZ,'0')
        
    def __repr__(self):
        b = self.b
        s = ' '.join([b[:8], b[8:16], b[16:24], b[24:]])
        s += ' = ' + str(eval('0b' + self.b))
        return s
        
    def __add__(self, a):
        carry = '0'
        rL = list()
        for i in range(len(self.b)-1,-1,-1):
            x = self.b[i]
            y = a.b[i]
            t = x+y
            if (t == '01' or t == '10'):
                if carry == '0':
                    rL.append('1')
                if carry == '1':
                    rL.append('0')
            elif t == '00':
                rL.append(carry)
                carry = '0'
            else:
                assert t == '11'
                if carry == '0':
                    rL.append('0')
                    carry = '1'
                else:
                    rL.append('1')
                    carry = '1'
        if carry == '1':
            raise ValueError, 'overflow'
        rL.reverse()
        return word(''.join(rL))
    
    def __mul__(self, a):
        L = list()
        for i in range(len(self.b)-1,-1,-1):
            x = self.b[i]
            if not x == '1':
                continue
            n = word.SZ - i - 1
            r = word(a.b[n:] + '0'*(n))
            L.append(r)
        res = L.pop(0)
        #print 'start:       ', res
        while L:
            next = L.pop(0)
            #print 'next:        ', next
            res = res + next
            #if L:
                #print 'intermediate:', 
            #else:
                #print 'result:      ',
            #print res
        return res

    def __pow__(self, n):
        res = self
        for i in range(n-1):
            res = res * self
        return res

test.py

import random
from word import word
        
R = range(1000)
N = 10
for i in range(N):
    x = random.choice(R)
    y = random.choice(R)    
    xw = word(x)
    print 'x    ', x
    yw = word(y)
    print 'y    ', y
    
    print 'x + y',
    r = xw + yw
    print r
    S = x+y
    assert S == int(str(r).split()[-1])
    print 'check', x, '+', y, '=', S
    
    print 'x * y',
    r = xw * yw
    print r
    P = x*y
    assert P == int(str(r).split()[-1])
    print 'check', x, '*', y, '=', P
    
    for i in range(2,N*2):
        try:
            r = xw**i
        except ValueError:
            i = i-1
            break
    
    z = i
    print 'z     ', z
    print 'x**z ', r
    E = x**z
    try:
        assert E == int(str(r).split()[-1])
    except AssertionError:
        print 'ran into a slight problem\n'
        continue
    print 'check', x,'**', z, '=', E
    print

Dissecting RSA keys in Python (4)

Still working with RSA keys. I accomplished the most important of the items left undone so far, to use the keys to encrypt a message using my own code.

In the process, I learned a new simple Python fact: the pow function (which is now a built-in), can take a modulus as a third argument.

And, it works where this doesn't seem to: (x**y) % z

(I discovered this by digging into the rsa module source).

python -m timeit -s "import rsa;  \
f = open('id_rsa');  data = f.read(); f.close(); k = rsa.PrivateKey.load_pkcs1(data)"\
 "41330915578951772302369**k.e % k.n"

10000 loops, best of 3: 79.4 usec per loop

python -m timeit -s "import rsa;  \
f = open('id_rsa');  data = f.read(); f.close(); k = rsa.PrivateKey.load_pkcs1(data)"\
 "pow(41330915578951772302369, k.e, k.n)"

10000 loops, best of 3: 70.2 usec per loop

In the example here, both work and have about the same timing. But in the code below, the first version hangs when doing decryption (with a large base). Here's the output, followed by the script.

> python encode.py 
m:   Hello, secret world!
p:   .xyz.Hello, secret world!.xyz.
a:   413309155789517723023698766343791993289928631329
c:   11434905702482726455415220687715293368262190253795 ..
i:   413309155789517723023698766343791993289928631329
r:   Hello, secret world!

encode.py

import rsa

with open('id_rsa') as f:
    data = f.read()
k = rsa.PrivateKey.load_pkcs1(data)
n = k.n
e = k.e
d = k.d

def my_atoi(s):
    L = [ord(c) for c in s]
    k = 256
    iL = L[:]
    iL.reverse()
    x = iL[0]
    for i in iL[1:]:
        x += i*k
        k *= 256
    return x

def my_itoa(i):
    rL = list()
    while i:
        rL.append(i%256)
        i = i/256
    rL.reverse()
    return ''.join([chr(n) for n in rL])

def encrypt(m):
    return m**e % n
    
def decrypt(c):
    # note:  c**d % n fails
    return pow(c,d,n)
    
if __name__ == '__main__':
    m = 'Hello, secret world!'
    pad = '.xyz.'
    p = pad + m + pad
    a = my_atoi(m)
    c = encrypt(a)
    i = decrypt(c)
    r = my_itoa(i)
    r = r.replace(pad,'')

    L = zip('mpacir',[m,p,a,c,i,r])
    N = 50
    for varname, var in L:
        s = str(var)
        print varname + ':  ', s[:N],
        if len(s) > N:  print '..'
        else:  print

Dissecting RSA keys in Python (3)

I'm working on understanding the structure of RSA keys as output by the ssh-keygen utility (previous posts here and here). At the risk of becoming redundant, I want to try to simplify what we did previously. Let's begin by rehashing old material on bytes in Python:

It's easy to get confused between a byte and its string representation, or for that matter, between an unsigned int and its string representation. Although hexadecimal is natural too, it's perhaps easiest to think of an individual byte in terms of the decimal equivalent (0..255). That's because a Python int converts easily to chr, bin, and hex.

>>> i = 35
>>> c = chr(i)
>>> c
'#'
>>> i == ord(c)
True
>>> b = bin(i)
>>> b
'0b100011'
>>> h = hex(i)
>>> h
'0x23'

Notice that hex and bin both give unpadded output:

>>> hex(1)
'0x1'
>>> bin(1)
'0b1'

Don't be fooled. h and b here are 'str' datatypes. But these string representations of bin and hex values go back to int easily (though we must specify the old base explicitly):

>>> int(b,2)
35
>>> int(h,16)
35

To interconvert hex and bin, it's easiest to go through int:

>>> bin(35)
'0b100011'
>>> hex(35)
'0x23'
>>> hex(int(bin(35),2))
'0x23'
>>> bin(int(hex(35),16))
'0b100011'

The struct module

This module performs conversions between Python values and C structs represented as Python strings. This can be used in handling binary data stored in files or from network connections, among other sources. It uses Format Strings as compact descriptions of the layout of the C structs and the intended conversion to/from Python values.

unpack works on the particular string representation of a byte in which it is a single character:

>>> i = 35
>>> c = chr(i)
>>> c
'#'
>>> unpack('B',c)
(35,)
>>> b = bin(i)
>>> unpack('B',b)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
struct.error: unpack requires a string argument of length 1
>>>

To construct a multi-byte value like a 4-byte unsigned int we can do this:

>>> h = '\x23'
>>> type(h)
<type 'str'>
>>> b_little = '\x23\x00\x00\x00'
>>> type(b_little)
<type 'str'>

With multi-byte objects (like ints) we get into endian-ness---'I' is the format for an unsigned int, and '>' means big-endian:

>>> unpack('I', b_little)
(35,)
>>> b_big = '\x00\x00\x00\x23'
>>> unpack('I', b_big)
(587202560,)
>>> unpack('>I', b_big)
(35,)

One can also use the binascii module, but this is more complicated. Let's write two bytes to a file and then read the binary data:

>>> fn = 'temp'
>>> FH = open(fn,'w')
>>> FH.write('##')
>>> FH.close()
>>> FH = open(fn,'rb')
>>> data = FH.read()
>>> FH.close()

This "data" is still a string:

>>> type(data)
<type 'str'>
>>> data
'##'
>>> len(data)
2
>>> b2a_hex(data[0])
'23'
>>> b2a_hex(data)
'2323'
>>> b2a_hex('#Az')
'23417a'

Or you can get fancy, but it seems unnecessary..

>>> import array
>>> L = array.array('B','#Az')
>>> type(L)
<type 'array.array'>
>>> L[0]
35
>>> L[1]
65
>>> L
array('B', [35, 65, 122])
>>> b2a_hex(L)
'23417a'

With all this in mind, we can redo the script from yesterday to read the base64-encoded key data in a way that is simply understandable. The output first:

> python read.py 
dlen 7
dlen 1
dlen 129

209,8,39,27 .. 169,97,145,127
0xd1,0x8,0x27,0x1b .. 0xa9,0x61,0x91,0x7f 

1467871546 .. 1722964351

We decode the data and break it into 3 parts exactly as before. Looking at the third part, we unpack it byte by byte. The result of unpacking is a list of ints. We convert that list directly into the large number n (by doing b * 256**i for each byte, moving along the list in reverse order), or we can convert the ints into hex values, assemble that list into one long hex string, and then call eval as before. It seems perfectly transparent now.

It even seems clear (in retrospect) why there is an additional null byte after the int that specifies the size of the third segment. It is probably to align the base64 encoding to begin with the first data byte.

The only remaining difficulty is that this approach fails for the private key. So that's still a mystery.

read.py

import base64
from struct import unpack

FH = open('data.txt','r')
data = FH.read()
FH.close()
b64_data = ''.join(data.strip().split())
data = base64.b64decode(b64_data)

L = list()
while data:
    dlen = unpack('>I',data[:4])[0]
    L.append(data[4:dlen+4])
    data = data[dlen+4:]
    print 'dlen', dlen
print

# let's look at n
bL = L[2]
bL = bL[1:]    # extra null first, why?
iL = [unpack('B',b)[0] for b in bL]
hL = [hex(i) for i in iL]

# look at the int and hex values
pL = [str(n) for n in iL]
print ','.join(pL[:4]), '..', ','.join(pL[-4:])
print ','.join(hL[:4]), '..',
print ','.join(hL[-4:]), '\n'

# first just do the computation ourselves
m = 256
iL.reverse()
n = iL[0]
for i in iL[1:]:
    n += i*m
    m *= 256
    
s = str(n)
print s[:10], '..', s[-10:]
    
# hex doesn't pad the output
# remove the '0x' and add the extra '0' if needed
for i in range(len(hL)):
    b = hL[i]
    if len(b) < 4:
        hL[i] = '0' + b[-1]
    else:
        hL[i] = b[-2:]  
h = '0x' + ''.join(hL)
nh = eval(h)
assert n == nh

Monday, August 29, 2011

Dissecting RSA keys in Python

I want to look at RSA keys a little bit more in this post. Let's generate a public/private pair (1024 bits) using the SSH utility keygen (no passphrase):

> ssh-keygen -b 1024 -t rsa 
Generating public/private rsa key pair.
Enter file in which to save the key (/Users/telliott_admin/.ssh/id_rsa): id_rsa
Enter passphrase (empty for no passphrase): 
Enter same passphrase again: 
Your identification has been saved in id_rsa.
Your public key has been saved in id_rsa.pub.
The key fingerprint is:
96:49:5d:40:b1:90:de:99:c5:05:43:f2:92:45:69:67 telliott_admin@c-98-236-78-154.hsd1.wv.comcast.net
The key's randomart image is:
+--[ RSA 1024]----+
|        .o=*B+.  |
|        .o B=.E  |
|       ...==.o   |
|       ..o+.     |
|        S        |
|       .         |
|                 |
|                 |
|                 |
+-----------------+

[ To be honest, at this point, I'm not actually sure this is the randomart that goes with the example in this post.. but I think it could be :) ]

Start with the "public" key. The data in the file looks like this:

id_rsa.pub

ssh-rsa AAAAB3NzaC1yc2EAAAABIwAAAIEA0QgnG+LBQosxpbkJIU0eWV9Zf7L63fy3BXp6T9RQH84yISs8SMvFV4J1NYJQtNpGopNFlvzsuFnnxaRaM8ureFIQVa8k/dRzv3r7QyIuoGyNQLDk3K8Tse/55fjuIxJFEf5jvoJs3Z/jYdPHreg8wzNCB/uMkbCccYsVt6lhkX8= my_identifier

Where my_identifier is

telliott_admin@c-myipaddress_with_dashes.hsd1.wv.comcast.net

The common representation of RSA keys (like we see in the file above) uses base64 format. I was confused about this at first, but found a great explanation here, and of course in wikipedia:

The base64, base32, and base16 encodings convert 8 bit bytes to values with 6, 5, or 4 bits of useful data per byte, allowing non-ASCII bytes to be encoded as ASCII characters for transmission over protocols that require plain ASCII, such as SMTP. The base values correspond to the length of the alphabet used in each encoding. There are also URL-safe variations of the original encodings that use slightly different results.

Consider the byte 1001-1110

We can encode this byte using hexadecimal representation as '9e'

>>> b = 0b10011110
>>> b
158
>>> hex(b)
'0x9e'

I didn't actually know you could do this, i.e. using 0b with no quotes (see SO).

The encoding of those 8 bits (one byte) as ASCII characters takes 2 bytes. The '9' and the 'e' take up that much space. An advantage of this representation of the byte is that it's easy to read, also the information can be transmitted to devices that are expecting an ASCII string, where the highest-value bit should be unset.

A disadvantage is that it doubles the size of the string that we're transmitting. As the quote indicates, the other encodings cram more bits of information into the byte to be transmitted. In base32 we could use (RFC4648) the characters 'A'..'Z' plus '2'..'7'.

In base64 we use uppercase, lowercase, digits plus '+' and '/', except that:

Using standard Base64 in URL requires encoding of '+', '/' and '=' characters into special percent-encoded hexadecimal sequences ('+' = '%2B', '/' = '%2F' and '=' = '%3D'), which makes the string unnecessarily longer.
For this reason, a modified Base64 for URL variant exists, where no padding '=' will be used, and the '+' and '/' characters of standard Base64 are respectively replaced by '-' and '_',

The wikipedia article has a great graphic showing an encoding in detail. We do conversions easily by using the Python base64 module.

>>> import base64
>>> e = base64.b64encode
>>> d = base64.b64decode
>>> e('Man')
'TWFu'
>>> d('TWFu')
'Man'
>>> s = 'Fourscore, and seven years ago'
>>> len(s)
30
>>> c = e(s)
>>> c
'Rm91cnNjb3JlLCBhbmQgc2V2ZW4geWVhcnMgYWdv'
>>> d(c)
'Fourscore, and seven years ago'

The encoding process looks at each sequence of 24 bits in the input (3 bytes) and encodes those same 24 bits spread over 4 bytes in the output. The last two characters [in the next example], the ==, are padding because the number of bits in the original string was not evenly divisible by 24 in this example.

>>> e('Ma')
'TWE='

We break this two byte (16 bit) value up into 6 + 6 + 6 = 18 bits---that makes three base64 values, and the third one is now E rather than F because the bit pattern is 000100, where the last two 00 are padding. Additional adding is added to bring us back to the byte boundary.

In Lincoln's famous phrase, there are 30 bytes, which is evenly divisible by 6, hence no padding is needed.

The second part of the problem was solved by a question I found on SO

Here is my version of the code:

import base64
import struct

data = open('id_rsa.pub').read().split(None)[1]
print data[:24]
keydata = base64.b64decode(data)
print keydata[:18]

parts = []
while keydata:
    dlen = struct.unpack('>I',keydata[:4])[0]
    data, keydata = keydata[4:dlen+4], keydata[4+dlen:]
    parts.append(data)

def f(x):
    return struct.unpack('B', x)[0]
    
e = eval('0x' + ''.join(['%02X' % f(x) for x in parts[1]]))
print 'e', e

n = eval('0x' + ''.join(['%02X' % f(x) for x in parts[2]]))
print 'n', n

I have to confess I haven't sorted this out completely yet. Nevertheless, it allows us to convert the RSA public key from ssh-keygen to an integer:

Output:

> python script.py
AAAAB3NzaC1yc2EAAAABIwAA
ssh-rsa#
e 35
n 146787154640181728129549360865206451974125178122992752499327844555082827713683060312609720092642289502088305950292885968241446393671268243539585900329449529436170858131743078225065287541399805133121852823856839448386268013284889428740476278604926908637093742329548018673369795976229837442944484597101722964351

I want to finish by showing that these values are actually correct, and tie the whole thing back to an earlier discussion of the mathematics of public key cryptography.

I found a module for Python that handles RSA cryptography.

One very nice thing about it is that it's pure Python. It might not be fast, but there's nothing to do to build it. (I used easy_install). It turns out that the rsa module can read the private RSA key that we got earlier (top of the post), but not the public key (which is why I learned everything that's come before this point):

> python
Python 2.7.1 (r271:86832, Jun 16 2011, 16:59:05) 
[GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2335.15.00)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import rsa
>>> with open('id_rsa') as f:  data = f.read()
... 
>>> data.split()[4][:50]
'MIICWgIBAAKBgQDRCCcb4sFCizGluQkhTR5ZX1l/svrd/LcFen'
>>>
>>> k = rsa.PrivateKey.load_pkcs1(data)
>>> k.n
146787154640181728129549360865206451974125178122992752499327844555082827713683060312609720092642289502088305950292885968241446393671268243539585900329449529436170858131743078225065287541399805133121852823856839448386268013284889428740476278604926908637093742329548018673369795976229837442944484597101722964351L
>>> k.e
35
>>> k.d
71296617968088267948638260991671705244575086516882194071102095926754516318074629294696149759283397758157177175856544613145845391211758861147798865874304045384903321200809536161346504709634159813307741235395231601606901597492591281205810328850837168587984349087962868939575631849929827946396907884195208774355L
>>> k.p
12214944519618901215070692893221534580187337690484280451921546594832793966376942655508370083135703324858458340368238021755353370324360760471666368429071147L
>>> k.q
12017013618393528112233953119917677248859562342149056793052943111068006176296374093622600767249370271694572460567324457033873003984957249202580919628769533L
>>> k.n == k.p * k.q
True
>>>

You can see that n equals p * q and that it matches what we eventually obtained from the public key file (id_rsa.pub) above. That's pretty cool.

Still to do: use the keys to do encryption and decryption and check that the math is right. Figure out how the hash/digest is done. And figure out how the passphrase is used to modify the private key file. We can even add a passphrase at this point to the old version:

ssh-keygen -f id_rsa -p

I used 'abcde' as a passphrase (weak!):

> cat id_rsa
-----BEGIN RSA PRIVATE KEY-----
Proc-Type: 4,ENCRYPTED
DEK-Info: AES-128-CBC,FB4D3BD371510C4B37552EF6E949CC48

E+KtbKVAwE4Dxj3vKYgfHpbgqF9Nk6O4vJ+C6wogXEHXr8vkM9SF5pz6lKNfGRbC
S93/lurvq+1+rCuldDNmEjdy36cbgrWJ1CIwB6fq5VyT3+v+AxDMqV9nQkvHjnxR
/ON2n+CjOSklnvmaBk3l8+IvCIKar2hGS/1qHqJRKJ3B1fgd5I8i9yH6hx8yl8Qm
AfVPgBJ4+91dCW2Y8odtXENgjstGtWZkYd0O+EqPuSuPHLu+J7M2KSQNe8aKOhGf
DXCEwR4aGW8EXfRNQpULbi5JruYebYlg3JILm2MaFwy+OW39q1VElvBc5wD01awW
BF7Y9x++bi9krX/G49NzkLdSVshkDMSP+eRZ4sncWvDCrkN2xlWqnNfQ06I+XCVW
FMUmyMFI87eNmjcZOgIi1S+fQtWcbl4o9Hx1z6xcMd7WxOnfIxe5zHIu8C24WONW
spcVxSsHAA57jYa8j292prtK6HDDa6UOAa32hBz1lvC8fMWPxwluTJ8eyeeK/7O6
hNI5G9gM4MsMEMau8O49fEm0RuiT/aB8IXZeEVQNbV/FdixyBBhY53DMnFQ5wau8
OJCpL51S5lj2o8eAJc61SHwX/Fc9uF+rmAFy08tthLuCF2wGdZWsE7IJ8NQzFJA+
Us2dTCjY57V/QoCC41rWh+YP+srcEK1Wp7DfvXcq98N1/+tECyCe/aDdPCIM/8QD
glKrx8Ah7Tsx6of3vCHiYG0q8Yu68NNJOiWrnGMEmP0Gq9yhLYqKrc8a3oIP0Lst
Ss+0p989df73GdM3nyqtLUwAzaCJQnDivLFT6gltEnM=
-----END RSA PRIVATE KEY-----

Wednesday, August 24, 2011

Pretty code (6)

This is absolutely the last post about html and code stuff. The whole adventure started here. It's a bit amusing to realize that if I had just understood what <br> tags were about (and <span> and <div>), this would never have happened.

I decided to write my own parser for Python code, so that I could post the results here. On the way, I learned something about "styles" and stylesheets. To organize things, I put this code for a "style" in the template for the blog at the end of <head>:

<style type='text/css'>
div.content
{
font-size:12pt;
font-family:arial;
color: #404040;
}

div.code
{ 
font-size:12pt;
font-family:courier;
margin: 20px;
border-style: solid;
padding: 5px;
border-width: 1px;
border-color:#ff00ff;
}

span.kw   { color:blue; }
span.str  { color:red; }
span.str3 { color:purple; }
span.py   { color:#00b0b0; }
span.cm   { color:green; }
</style>

The parser output looks great when run on itself. The first little bit is here:

But, alas, it was not to be. The Blogger editor is now (still) messing with my code, removing all the <br /> that I carefully put in, so that the informative whitespace is gone in the actual post (though not in the intial preview before posting). It looks hopeless at this point.

Files for the parser on Dropbox (here). The only known bugs are those visible in the comment lines above, where the strings should have been masked by being inside a comment.

Unfortunately, I do not know how to make this all work for Blogger.

What a waste! I am not happy with these guys. The only way I can see is to have my own server.

Monday, August 22, 2011

Pretty code (4)

Explanation of what this is about from last time here. I got triple-quoted strings working, and I'm really pleased with the results. More another time on other options for code-highlighting.


import sys
from keyword import kwlist
from string import punctuation as punct
from utils import load_data
import html_tags as H
pywords = ['True','False','list','dict',
           'int','float','append','extend',
           'sys','argv','pop','open','write',
           'close','string','time']

try:  
    fn = sys.argv[1]
except IndexError:  
    fn = 'example.py'
data = list(load_data(fn))

D = {'is_cm':False,
     'in_str_1':False,
     'in_str_2':False,
     'in_str_3':False }

L = list()

while data:
    c = data.pop(0)
    # comments first
    if c == '#':
        if not (D['in_str_1'] or D['in_str_2']):
            L.extend(list(H.cm_start))
            D['is_cm'] = True
    if c == "\n" and D['is_cm']:
        L.extend(list(H.cm_stop))
        D['is_cm'] = False
    L.append(c)
    
    a = '''triple-quoted-string
    with a continuation and two keywords'''
    
    # triple-quoted strings
    if c == "'" and len(data) > 1:
        if data[0] == "'" and data[1] == "'":
            if (not D['in_str_1'] and not D['in_str_2']):
                if not D['in_str_3']:
                    L.pop()
                    L.extend(list(H.str3_start))
                    L.extend(["'"] * 3)
                    data.pop(0)
                    data.pop(0)
                    D['in_str_3'] = True
                else:
                    data.pop(0)
                    data.pop(0)
                    L.extend(["'"] * 2)
                    L.extend(list(H.str3_stop))
                    D['in_str_3'] = False
                    continue
    
    # single-quoted strings
    if c == "'":
        if D['in_str_3']:
            continue        
        if (not D['in_str_1']):
            if not D['in_str_2']:
                # start a str_1
                L.pop()
                L.extend(list(H.str_start))
                L.append(c)
                D['in_str_1'] = True
            else:
                # already in str_2 or str_3
                pass
        else:
            if D['in_str_1']:
                # terminate str_1    
                L.extend(list(H.str_stop))
                D['in_str_1'] = False
            
    # double-quoted strings
    if c == '"':
        if not D['in_str_2'] and not D['in_str_3']:
            if not D['in_str_1']:
                # start a str_2
                L.pop()
                L.extend(list(H.str_start))
                L.append(c)
                D['in_str_2'] = True
            else:
                # already in str_1 or str_3
                pass
        else:
            if D['in_str_2']:
                # terminate str_2    
                L.extend(list(H.str_stop))
                D['in_str_2'] = False
s = ''.join(L)

# keywords last
pL = list()
D['is_str_3'] = False
for line in s.split('\n'):
    D['is_cm'] = False
    words = line.split()
    for w in words:
        # no kw highlighting in comments
        if w.startswith(H.cm_start):
            D['is_cm'] = True
        if w.startswith(H.str3_start):
            D['is_str_3'] = True
        if H.str3_stop in w:
            D['is_str_3'] = False
        if not D['is_cm'] and not D['is_str_3']:
            if w in kwlist:
                r = H.kw_start + w + H.kw_stop
                line = line.replace(w, r)
            for p in pywords:
                if p in w:
                    L = w.split(p)
                    if L[0] and not L[0][-1] in punct:
                        continue
                    if L[1] and not L[1][0] in punct:
                        continue
                    r = H.py_start + p + H.py_stop
                    line = line.replace(p, r)
    pL.append(line)
    
s = H.br.join(pL)
pL = [H.head, H.hr, s, H.hr, H.tail]
s = '\n'.join(pL)

fn = fn.split('.')[0] + '.html'
FH = open(fn,'w')
FH.write(s + '\n')
FH.close()

Here's a screenshot of html_tags.py:

Wednesday, August 3, 2011

e again

I came across a neat quote from Martin Gardner:

The number e is as ubiquitous as pi, turning up everywhere, and especially in probability theory. If you randomly select real numbers between 0 and 1, and continue until their sum exceeds 1, the expected number of choices is e.

-Martin Gardner
note in Calculus Made Easy, p. 153

Of course, this begs for a Python simulation:

import random

def mean(L):
    return float(sum(L))/len(L)

def oneTrial():
    counter = 1
    f = random.random()
    while f <= 1:
        counter +=1
        f += random.random()
    return counter

N = int(1E7)
L = [oneTrial() for i in range(N)]
print mean(L)

> python script.py 
2.7183671

A derivation would be nice, but that'll have to wait. I wonder about the variance.

Friday, June 24, 2011

Cramer's rule calculation

We have algebra homework that involves using Cramer's rule to solve not only 2 x 2 but also 3 x 3 systems. It seems kind of silly since this method is overkill for 2 x 2, and would never be used for 4 x 4 or larger.

(Note on the wikipedia article, start about halfway down, where it says "Explicit formulas for small systems".)

Also, and this gets closer to the point, drilling by solving 3 x 3 matrices is not really about the rule, which is pretty simple. It's about making an easy problem harder by stuffing a lot of arithmetic into it. And to me, that is symptomatic of a big difficulty with math education as I'm encountering it through my son. Computers are much better at computing sums than humans. It's just silly to drill students on arithmetic. If you want to do something complicated, why not derive Cramer's rule?

So, I decided to write a solver for 3 x 3 systems in Python. I wouldn't say it's thoroughly debugged yet, so let me know if you run into a problem. With the example shown, I did get the same answer as this online calculator.

The first code segment contains the equations explicitly entered as an array. I'm sure you know how to modify it to read input from a file.

test.py

import numpy as np
import Cramer

def test_Cramer():
    L = [2, 3, 0, 5,
         1, 1, 1, 3,
         2,-1, 3, 7]   
    A = np.array(L)
    A.shape = (3,4)
    result = Cramer.solve(A)
    if result:
        x,y,z = result
        print 'solution'
        print 'x =', x
        print 'y =', y
        print 'z =', z, '\n'
        Cramer.check(A,x,y,z)

test_Cramer()

The output looks like this:

> python test.py 
solve
[[ 2  3  0  5]
 [ 1  1  1  3]
 [ 2 -1  3  7]] 

compute 3 x 3 det of
[[ 2  3  0]
 [ 1  1  1]
 [ 2 -1  3]]
D =  5 

compute 3 x 3 det of
[[ 5  3  0]
 [ 3  1  1]
 [ 7 -1  3]]
Dx =  14 

compute 3 x 3 det of
[[2 5 0]
 [1 3 1]
 [2 7 3]]
Dy =  -1 

compute 3 x 3 det of
[[ 2  3  5]
 [ 1  1  3]
 [ 2 -1  7]]
Dz =  2 

solution
x = 2.8
y = -0.2
z = 0.4 

check
row 0 = [2 3 0 5]
2.0*2.8 + 3.0*-0.2 + 0.0*0.4  = 5.0 

row 1 = [1 1 1 3]
1.0*2.8 + 1.0*-0.2 + 1.0*0.4  = 3.0 

row 2 = [ 2 -1  3  7]
2.0*2.8 + -1.0*-0.2 + 3.0*0.4  = 7.0

Cramer.py

import numpy as np

def det2x2(A, v=False):
    if v:  print 'compute 2 x 2 det of'
    if v:  print A
    assert A.shape == (2,2)
    return A[0][0]*A[1][1] - A[0][1]*A[1][0]

def det3x3(A):
    print 'compute 3 x 3 det of'
    print A
    assert A.shape == (3,3)
    a,b,c = A[0]
    c1 = a * det2x2(A[1:3,[1,2]])
    c2 = b * det2x2(A[1:3,[0,2]])
    c3 = c * det2x2(A[1:3,[0,1]])
    return c1 - c2 + c3
         
def solve(A):
    print 'solve'
    print A, '\n'
    assert A.shape == (3,4)
    D = det3x3(A[:,:3])
    print 'D = ', D, '\n'
    if D == 0:
        print 'no solution'
        return
    Dx = det3x3(A[:,[3,1,2]])
    print 'Dx = ', Dx, '\n'
    Dy = det3x3(A[:,[0,3,2]])
    print 'Dy = ', Dy, '\n'
    Dz = det3x3(A[:,[0,1,3]])
    print 'Dz = ', Dz, '\n'
    return Dx*1.0/D, Dy*1.0/D, Dz*1.0/D
    
def check(A,x,y,z):
    print 'check'
    for i,r in enumerate(A):
        print 'row', i, '=', r
        pL = list()
        for coeff,var in zip(r[:3],(x,y,z)):
            c = str(round(coeff,2))
            v = str(round(var,2))
            pL.append(c + '*' + v)
        print ' + '.join(pL),
        print ' =', r[0]*x + r[1]*y + r[2]*z, '\n'

Wednesday, May 4, 2011

Intro to ANOVA

This is an introductory post on ANOVA (analysis of variance). We ask the question: given three (or more) groups of observations, is one or more of the group means significantly different from the others. We will compute an F-statistic, and compare that with an F-distribution (carry out an F-test). If the statistic exceeds the 95% quantile, we will reject the null hypothesis that the means are the same.

According to wikipedia:

the two-group case can be covered by a t-test (Gosset, 1908). When there are only two means to compare, the t-test and the ANOVA F-test are equivalent; the relation between ANOVA and t is given by F = t².

ANOVA is a versatile (and complex) set of methods. This is just an elementary application, where we'll use the R implementation on three simple groups of data, and then compute the result ourselves in Python to see how it works internally.

To begin with, we follow the simple example from Dalgaard. You will need R and the ISwR package (or just construct the "data frame" yourself). R code:

> library(package=ISwR)
> data(red.cell.folate)
> summary(red.cell.folate)
     folate          ventilation
 Min.   :206.0   N2O+O2,24h:8   
 1st Qu.:249.5   N2O+O2,op :9   
 Median :274.0   O2,24h    :5   
 Mean   :283.2                  
 3rd Qu.:305.5                  
 Max.   :392.0                  
> red.cell.folate
   folate ventilation
1     243  N2O+O2,24h
2     251  N2O+O2,24h
3     275  N2O+O2,24h
4     291  N2O+O2,24h
5     347  N2O+O2,24h
6     354  N2O+O2,24h
7     380  N2O+O2,24h
8     392  N2O+O2,24h
9     206   N2O+O2,op
10    210   N2O+O2,op
11    226   N2O+O2,op
12    249   N2O+O2,op
13    255   N2O+O2,op
14    273   N2O+O2,op
15    285   N2O+O2,op
16    295   N2O+O2,op
17    309   N2O+O2,op
18    241      O2,24h
19    258      O2,24h
20    270      O2,24h
21    293      O2,24h
22    328      O2,24h

We have 22 values in 3 groups.

> attach(red.cell.folate)
> anova(lm(folate~ventilation))
Analysis of Variance Table

Response: folate
            Df Sum Sq Mean Sq F value  Pr(>F)  
ventilation  2  15516  7757.9  3.7113 0.04359 *
Residuals   19  39716  2090.3                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
>

The attach gives us access to the names of the columns of values (folate) and factors (ventilation). We make a plot (shown at the top of the post):

> plot( folate ~ ventilation, data = red.cell.folate )
> stripchart(x=folate~ventilation,
  pch=16,vertical=T,add=T,col='blue')

The value 0.04359 indicates that we have P < 0.05.

We write the data to disk, and remember that the groups have (respectively) 8, 9 and 5 values.

> setwd('Desktop')
> write.table(folate,'data.txt',row.names=F,col.names=F)

We use the Python script below to compute the F-statistic:

python script.py
..
MS_W 2090.32
MS_X 7757.88
F_stat 3.71

If you look in the R output above, you'll see the same values as given here for MS_W and MS_X and the F-statistic.

To get the right F-distribution, we need to know that the degrees of freedom are:

k-1 = 2  # k = number of groups
N-k = 19  # N = total observations

Since 3.71 is just higher than the 95% quantile of this F-distribution we can reject the null hypothesis H₀.

I found a calculator online. You can see the results in the screenshot.

The underlying calculation is pretty simple. We compute sumsq, the sum of the squares of the differences from the mean for several sets of values and the relevant means.

For the within groups comparisons, using mathematical notation this is (i groups with j observations in each group):

    Σ         Σ     (x_ij - m_i)²
(over i)  (over j)

In Python, for each group we sumsq for the group compared with the group mean, and add the results for all three groups to give SSD_W.

To carry out the between groups comparisons, we first compute the grand mean m (of all of the samples). Then for each group we compute:

    Σ        Σ     (m_i - m)²
(over i)  (over j)

Since the squared value is the same within each group, this is equivalent to:

    Σ      n_i (m_i - m)²
(over i)

In the Python code this becomes:

len(g)*(mean(g)-m)**2

and sum these over all the groups to give SSD_X. This is a sum of squares of the group means.

Finally, we compute:

MS_W = SSD_W/(N-k)
MS_X = SSD_X/(k-1)
F_stat = MS_X/MS_W

As to why we do this, for now you will have to go read the article. The explanation in Dalgaard is particularly clear, indeed, the whole book is excellent.

There is a SciPy function to carry out ANOVA (stats.f_oneway), but I don't have SciPy installed right now at home, and this post is long enough already. That's for another day.

Python code:

fn = 'data.txt'
FH = open(fn,'r')
data = FH.read().strip().split()
FH.close()

data = [int(n) for n in data]
A = data[:8]
B = data[8:17]
C = data[17:]

#A = [243,251,275,291,347,354,380,392]
#B = [206,210,226,249,255,273,285,295,309]
#C = [241,258,270,293,328]

def mean(L):
    return sum(L)*1.0/len(L)

def sumsq(L):
    m = mean(L)
    print 'sumsq'
    rL = [(x-m)**2 for x in L]
    for n,p in zip(L,rL):
        print n, round(p,1)
    S = sum(rL)
    print 'total', round(S,1), '\n'
    return S

def ANOVA(G):
    # variation within groups
    SSD_W = 0
    for g in G:
        SSD_W += sumsq(g)
    
    # a bit awkward, just flattening the list of lists
    # to get the mean and N
    T = list()
    for g in G:
        T.extend(g)
    m = mean(T)
    
    # variation between groups (X for 'cross')
    SSD_X = 0
    for g in G:
        SSD_X += len(g)*(mean(g)-m)**2
        
    N = len(T)  # 22
    k = len(G)  # 3
    MS_W = SSD_W*1.0/(N-k)
    MS_X = SSD_X*1.0/(k-1)
    F_stat = MS_X/MS_W
    return MS_W, MS_X, F_stat

MS_W, MS_X, F_stat = ANOVA([A,B,C])
print 'MS_W', round(MS_W,2)
print 'MS_X', round(MS_X,2)
print 'F_stat', round(F_stat,2)

Friday, March 18, 2011

nothing in biology makes sense except in the light of evolution

I've been reading Mike Yarus's book: Life from an RNA world (Amazon). It's a very readable account of evolution from the perspective of sequences and the RNA world. Mike's a highly intelligent guy, and that and his wit inform every page. In one chapter, he gets Rob Knight and Steve Freeland to do an evolution simulation.

we let a computer write out a random string, mutate 1 in 100 characters in each generation, and select changes only if they match Dobzhansky

In my version:

in each generation we pick one position

if it already matches Dobzhansky, continue to the next generation

mutate to a random choice from the set of symbols

This isn't a very good model of evolution. It was just fun to spend a few minutes coding it.

Here is the beginning, middle and end of one run:

> python evolve.py 
    0 LVSTUknwIfGIRHFgHESZ M zthWtNQTk qgtoeMvjeJAzOidKWEZO ZwNjDyCvq
  100 LVSTUknwIfGIRHFgHESZ M zthetNQTk qgtoeMvjeJAzOidKWEZO ZwNjDyCvq
  200 LVSTUknwIfGIRHFgHESZ M zthetNQTk qgtoeMvjeJAzOidhWEZO ZwNjDyCvq
  300 LVSTiknwIfGIRHFgHESZ M zthetNQ k qgtoeMvjeJAzOidhWEZO ZwNjuyivq

 6100 notTing if biHlogy makes sense except iv the Oight of Zvolution
 6200 notTing if biHlogy makes sense except iv the Oight of Zvolution
 6300 nothing if biHlogy makes sense except iv the Oight of Zvolution
 6400 nothing if biHlogy makes sense except iv the Oight of Zvolution

13800 nothing in biHlogy makes sense except in the light of evolution
13900 nothing in biHlogy makes sense except in the light of evolution
14000 nothing in biHlogy makes sense except in the light of evolution
14100 nothing in biHlogy makes sense except in the light of evolution

14167 nothing in biology makes sense except in the light of evolution
14167 nothing in biology makes sense except in the light of evolution

import random
s = 'nothing in biology makes sense'
s += ' except in the light of evolution'
N = len(s)

symbols = 'abcdefghijklmnopqrstuvwxyz '
symbols += symbols.upper()
L = [random.choice(symbols) for i in range(N)]
print '    0 ' + ''.join(L)

i = 0
while L != list(s):
    i += 1
    v = i and not i % 100
    if v:  print str(i).rjust(5),
    j = random.choice(range(N))
    if L[j] == s[j]:  
        if v:  
            print ''.join(L)
        continue
    c = random.choice(symbols)
    if c == s[j]:  L[j] = c
    if v:  print ''.join(L)
    
print
print str(i).rjust(5), s
print str(i).rjust(5), ''.join(L)

Mutual Information (3)

We're working on a paper from Michael Laub's lab at MIT (Skerker et al 2009 PMID 18555780). Previous posts here and here.

Now it's time to analyze the data. The self comparisons (a single position in the alignment, analyzed for mutual information against itself) yield info values ranging from 4.1 to 0.01.

As you can see in the screenshot, the residue in column 16 (1-based indexing), which has a very low score, is the conserved (catalytic) histidine.

We'll filter out the self comparisons for the plots. Here is the histogram of information values that I got. It's a bit different from the paper, but not much.

In the next part, we'll try to match up the residues which look like they might interact (that have high mutual information) and see if that makes sense in terms of the protein structures.

I picked a familiar HK/RR pair to do this: EnvZ and OmpR. This screenshot of part of Fig 2 shows a section of each protein.

Searching in the alignment file (the alignments don't have titles that I recognize), I recovered a sequence that I think is probably the right one. It matches the Figure as far as I checked:

>26250004-26250005/1-457
KQLADDRTLLMAGVSHDLRTPLTRI
RLATEMMSAESINKDIEECNAIIEQ
FIDYLRTGMADLNAVLGEVIA--AE
SGYEREIETAL-YVKMHPLSIKRAV
ANMVVNAARY-GNGWIKVSSGTEAW
FQVEDDGPGIAPEQRKHLFQPFVRG
DISGTGLGLAIVQRIVDNHNGMLEL
GTSERGGLSIRAWLPNYKILVVDDD
MRLRALLERYLTEQGFQVRSVANAE
QMDRLLTRESFHLMVLDLMLPGEDG
LSICRRLRSQSPMPIIMVTAKGEEV
DRIVGLEIGADDYIPKPFNPRELLA
RIRAVLRR

I put the newlines in to help me count. The alignment is 308 residues total. From Fig 2,

EnvZ (the HK) sequence starts with:

AGVKQLADDRTLLMAGVSHDL
   KQLADDRTLLMAGVSHDL

The second line above is from the alignment. The sequence doesn't start at residue 1, however. By my calculation, residue 0 in the alignment (the K) is residue 15 in the protein (1-based index). So we'll adjust the indexes we obtain by adding 15 to compare with the actual protein sequence.

OmpR (the RR) sequence starts with

MQENYKILVVDDDMRLRALLER
   NYKILVVDDDMRLRALLER

The second line above is from the alignment, where there's an N at position 0 in this fragment. By my calculation, that N is residue 3 of OmpR in Fig 2 (1-based index), and residue 190 of the alignment, so we subtract 187 for values >= 190.

Going back to the plotting script, we filter the data for pairs in which one comes from the HK and one from the RR, and print the top 20. We do the math as outlined above. As a check, we grab the putative EnvZ/OmpR sequence from the alignment file and print the sequence starting at the position we've identified. Here are the results for the top 20 (in each pair the RR is printed first and the HK second).

[Note: to keep the code simple, I ignored the situation where a column contains '-' a gap in the alignment for some of the sequences. That's what's giving us the two strange results below at position 16 and 20.]

> python plot.py 
 14 RLRAL  42 ATEMM 0.822
 18 LLERY  42 ATEMM 0.816
 22 YLTEQ  42 ATEMM 0.769
 15 LRALL  37 TRIRL 0.730
 15 LRALL  42 ATEMM 0.694
 56 MLPGE  54 DIEEC 0.688
 14 RLRAL  54 DIEEC 0.678
 21 RYLTE  42 ATEMM 0.677
 14 RLRAL  38 RIRLA 0.668
 83 KGEEV  22 TLLMA 0.667
 22 YLTEQ  21 RTLLM 0.663
  4 YKILV  42 ATEMM 0.658
 83 KGEEV  54 DIEEC 0.657
 22 YLTEQ  18 ADDRT 0.650
 22 YLTEQ  54 DIEEC 0.646
 18 LLERY  86 --AES 0.644
 15 LRALL  54 DIEEC 0.643
 18 LLERY  38 RIRLA 0.636
 22 YLTEQ  45 MMSAE 0.633
 22 YLTEQ  86 --AES 0.631

The two residues with highest mutual information are OmpR residue 14 and EnvZ residue 42. It looks pretty good to me.

[ UPDATE: The heatmap looks pretty boring, so I'm going to skip it.
But I plotted the top 20 interactions (graphic at the top of the post), the repetition indicates we're on the right track.. ]

import sys
from utils import load_data
import matplotlib.pyplot as plt

data = load_data('results.2.txt')
data = data.strip().split('\n')
data = [e.split() for e in data]
data = [(int(t[0]), int(t[1]), float(t[2])) for t in data]

def f(t):  return float(t[2])
def part1():
    L = [t for t in data if t[0] == t[1]]
    L = sorted(L,key=f,reverse=True)
    for t in L[:10]:  print str(t[0]+1).rjust(3),round(t[2],3)
    print
    for t in L[-10:]:  print str(t[0]+1).rjust(3),round(t[2],3)
    sys.exit()
    
#part1()


L = [t[2] for t in data if t[0] != t[1]]
X = 1.0
plt.hist(L,bins=X*50)
ax = plt.axes()
ax.set_xlim(0,X)
plt.savefig('example.png')

# t[0] always > t[1]
N = 190
data = [t for t in data if t[0] >= N and t[1] < N]

aln = load_data('cell3925mmc4.fa')
aln = aln.strip().split('>')[1:]
aln = [e for e in aln if e.startswith('26250004-26250005/1-457')]
envZ_ompR = aln[0].split('\n')[1]

for t in sorted(data,key=f,reverse=True)[:20]:
    i = t[0]
    rr = i - 187
    j = t[1]
    hk = j + 15
    print str(rr).rjust(3), envZ_ompR[i:i+5],
    print str(hk).rjust(3), envZ_ompR[j:j+5],
    print '%3.3f' % t[2]

Mutual information (2)

We're working on a paper from Michael Laub's lab at MIT (Skerker et al 2009 PMID 18555780). The first post is here.

In this part we'll load the alignment (supplementary data file S4---the annotation on the page is incorrect), and crunch the numbers. I just write the results to disk.

We'll do the analysis in another post.

python info.py > results.txt

import sys
from utils import load_data
import info_helper as ih

#aln = 'AASSASSTTT\nNMWWNTTKKS\nGTSNTYRSTA\nGGGGGGGGGG'
fn = 'cell3925mmc4.fa'
data = load_data(fn)
data = data.strip().split('>')[1:]
data = [e.split('\n')[1].strip() for e in data]

def show(data):
    print 'starting:', len(data)
    for i in range(7):
        print i, 
        L = [e for e in data if e.count('-') <= i]
        print len(L)
    sys.exit()
#show(data)

def transpose(L):
    R = range(len(L[0]))
    rL = list()
    for i in R:
        rL.append(''.join([item[i] for item in L]))
    return rL

data = [e for e in data if e.count('-') <= 4]
#data = data[:100]
cols = transpose(data)
pD = ih.make_prob_dict(cols)
info = dict()

for i in range(len(cols)):
    for j in range(i+1):      
        info[(i,j)] = ih.get_info(i,j,cols,pD,v=False)
for i,j in sorted(info.keys()):
    print i,j,round(info[(i,j)],3)

Tuesday, March 15, 2011

Mutual information

I want to talk about a really nice paper from Michael Laub's lab at MIT (Skerker et al 2009 PMID 18555780). It'll give us an opportunity to exercise our matplotlib skills.

We're going to try to recreate Fig 1, which is visible in the PubMed page, or you can get the original paper from the link to Cell.

Two-component signal transduction systems are ubiquitous in bacteria (wikipedia). The canonical design consists of a membrane-bound sensor (histidine) kinase (HK) and a cytoplasmic response regulator (RR). E. coli contains about 30 such pairs. The members of each pair have substantial specificity. The HK of the ntr system has specificity for its own RR, and likewise in the phoBR system, phoB has specificity for phoR. We may speak of a HK and its cognate RR.

For our purposes the important thing is that each system comprises (in the simplest design) two protein partners with complementary surfaces. These systems (a pair of proteins) are the products of ancient gene duplication events, and have since diverged over time. Amino acids at interacting sites are constrained to co-evolve in each pair.

If this sounds too vague or too complicated, consider an even simpler example: a stem of paired RNA residues in rRNA named H15.

Here is the H15 sequence in 1D (the parentheses indicate residues involved in pairing---see the link above for details):

((((   (((((    ))))) ))))
TGCACAATGGGCGCAAGCCTGATGCA

And here is the inner stem drawn in 2D to show the base-pairing more directly:

TGGGC
GTCCG

The base-pairing of this stem is more important to rRNA function than the identities of the bases. The result is that in some bacteria the identities of the central bases have been switched:

original   co-evolved

TGGGC        TGCGC
GTCCG        GTGCG

Presumably this happened in 2 discrete steps, but I don't know of any examples where the intermediate state has been preserved. Maybe we should look for some, and it's undoubtedly been studied.

To quantify this kind of coevolution, we'll draw on a concept (and mathematical definition) called mutual information. The steps in the calculation will be:

make a multiple sequence alignment

compare column X and column Y

total number of sequences (length of each column) = c

for each residue x in column X calculate px

for each residue y in column Y calculate py

px is the probability of residue x in column X

We'll write the columns horizontally to save space.
Suppose column X and Y are:

X:  AASSASSTTT
Y:  NMWWNTTKKS

For column X we have:

pA = 0.3 (3 A out of a total of 10 residues)
pS = 0.4
pT = 0.3

For column Y:


pK = 0.2
pM = 0.1
pN = 0.2
pS = 0.1
pT = 0.2
pW = 0.2

We pre-calculate these values for each column. When we calculate the information, we'll refer to the probabilities for column Y as q rather than p, to keep them straight from the p's for column X.

Now, we consider each pair of residues, one from column X and one from column Y. This pair is made up of residues in two interacting protein surfaces or rRNA chains, that may have co-evolved.

pxy is the number of sequences with x in column X and y in column Y

divided by c, the total number of pairs:

X:  AASSASSTTT
Y:  NMWWNTTKKS


pAM = 0.1
pAN = 0.2
pST = 0.2
pSW = 0.2
pTK = 0.2
pTS = 0.1

Finally, to compute the mutual information for this pair of columns, we do this calculation for each individual pair of residues and then sum:

pAM * log (pAM / pAqM) = 0.1 * log (0.1 / (0.3 * 0.1)) = 0.0523
I would have used log₂, but Skerker et al used log₁₀, so I matched them.

Here is part of the output of the script below:

NMWWNTTKKS
AASSASSTTT
pKT 0.2 pK 0.2 qT 0.3 temp 3.33 final 0.1
pMA 0.1 pM 0.1 qA 0.3 temp 3.33 final 0.05
pNA 0.2 pN 0.2 qA 0.3 temp 3.33 final 0.1
pST 0.1 pS 0.1 qT 0.3 temp 3.33 final 0.05
pTS 0.2 pT 0.2 qS 0.4 temp 2.50 final 0.08
pWS 0.2 pW 0.2 qS 0.4 temp 2.50 final 0.08
info 0.47

temp is the result of the calculation inside the parentheses, above. Next time we'll apply this method to the data from Skerker.

from math import log
def log2(n):  return log(n)*1.0/log(2)
def log10(n):  return log(n)*1.0/log(10)

# cache character probabilities for each column
def make_prob_dict(cols):
    # input is a list of columns
    # c is the number of sequences in the alignment
    c = len(cols[0])
    pD = list()
    for col in cols:
        char_kinds = list(set(col))
        values = [col.count(k)*1.0/c for k in char_kinds]
        pD.append(dict(zip(char_kinds,values)))
    return pD

def get_info(i,j,cols,pD,v=False):
    col1, col2 = cols[i], cols[j]
    if v:  print col1 + '\n' + col2
    # as before, c is the number of sequences
    c = len(col1)
    assert c == len(col2)
    info = 0
    pairs = [col1[k] + col2[k] for k in range(c)]
    pL = sorted(list(set(pairs)))
    for p in pL: 
        pXY = pairs.count(p)*1.0/c
        pX = pD[i][p[0]]
        pY = pD[j][p[1]]
        inside = (pXY * 1.0) / (pX * pY)
        if v:  print 'p' + p, pXY,
        if v:  print 'p' + p[0], pX,
        if v:  print 'q' + p[1], pY,
        if v:  print 'temp', '%3.2f' % round(inside, 2),
        outside = pXY * log10(inside)
        if v:  print 'final', round(outside,2)
        info += outside
    if v:  print 'info', round(info,2)
    return info
  
if __name__ == '__main__':
    aln = 'AASSASSTTT\nNMWWNTTKKS\nGTSNTYRSTA\nGGGGGGGGGG'
    cols = aln.split('\n')
    pD = make_prob_dict(cols)
    info = dict()
    for i in range(len(cols)):
        for j in range(i):
            info[(i,j)] = get_info(i,j,cols,pD,v=True)

Tuesday, March 8, 2011

16S rRNA V regions, continued

As promised last time (here), this is the code to plot V regions from the sequences in the Enterobacteriales, obtained from RDP. If you compare to this post from the other day, you'll see we match pretty well.

The figure looks better with a wider window, but the limits are established more accurately with a smaller window. The threshold T needs to be adjusted based on the window size.

One detail: I did the redirect below to save the extreme values and then parsed them with the code at the end of the post.

python script.py > extreme_values.txt

script.py

import utils as ut
from info import shannon
import matplotlib.pyplot as plt

data = ut.load_data('entero.txt')
data = data.strip().split('>')[1:]
EC = [e for e in data if 'X80725' in e][0]
EC = ut.clean_fasta(EC)[1]
data = [ut.clean_fasta(e)[1] for e in data]
L = ut.make_count_list(data)

pos = 0
R = range(1,1451)
iL = list()

for i,c in enumerate(EC):
    if c == '-':  continue
    pos += 1
    if not pos in R:  continue
    cD = L[i]
    e = shannon(cD,'ACGT')
    iL.append(e)
    #print str(pos).ljust(3), c + '  ', 
    #print ''.join([str(cD[k]).ljust(5) for k in 'ACGT']),
    #print round(e,2)

aL = list()
w = 20
T = 1.8
for i in range(len(iL)):
    j, k = i - w, i + w + 1
    if j < 0:  j = 0
    if k > len(iL):  k = len(iL)
    m = ut.mean(iL[j:k])
    if m < T:  print i+1
    aL.append(m)
    
plt.plot((1,1451),(T,T),lw=2,color='r',zorder=0)
plt.scatter(R,aL)
ax = plt.axes()
ax.set_xlim(-5,1455)
ax.set_ylim(0.8,2.05)
plt.savefig('example.pdf')

analyze.py

import utils as ut
data = ut.load_data('extreme_values.txt')
L = [int(n) for n in data.strip().split('\n')]

current = L.pop(0)
print current, '-',
while L:
    next = L.pop(0)
    if next != current + 1:
        print current
        print next, '-',
    current = next
print next

16S rRNA V regions

I'm exploring ways to visualize the sequence diversity in 16S rRNA. I finally realized that to see the "V" regions clearly, I need an alignment in which the V regions can align, that is, the sequences need to be closely enough related. To get such a set, I went to RDP (Browsers), and grabbed 187 type sequences from the Order Enterobacteriales. In the next post, I'm actually going to generate the graphic above (it's a tease). The figure shows the information content of the sequence set computed for a window sliding across the length of the gene. The troughs are V regions.

But first, I thought I would post some utility code separately. Here is a script that imports the utility functions shown in the last part, below, and exercises them on the enteric sequences:

import utils as ut

data = ut.load_data('entero.txt')
data = data.strip().split('>')[1:]
seqs = [ut.clean_fasta(e)[1] for e in data]
L = ut.make_count_list(seqs,kL='ACGT')
for D in L[75:80]:
   print D

Here is what it prints:

> python rdp.py
{'A': 174, 'C': 2, 'T': 3, 'G': 1}
{'A': 51, 'C': 2, 'T': 3, 'G': 123}
{'A': 3, 'C': 172, 'T': 2, 'G': 3}
{'A': 127, 'C': 5, 'T': 1, 'G': 47}
{'A': 1, 'C': 76, 'T': 3, 'G': 99}

The first module contains a function to return the Shannon entropy for a distribution:

info.py

from math import log
def f_logf(f):  
    if f == 0:  return 0
    return f*log(f)*1.0/log(2)

def shannon(cD,kL=None):
    if not kL:
        kL = cD.keys()
    L = [cD[k] for k in kL]
    S = sum(L)
    if S == 0:
        raise ValueError('Empty list')
    L = [n*1.0/S for n in L]
    L = [f_logf(f) for f in L]
    return 2 + sum(L)

The second module contains some functions that I use all the time. I put this in my site-packages directory:

utils.py

def load_data(fn):
   FH = open(fn,'r')
   data = FH.read().strip()
   if '\r\n' in data: 
       data = s.replace('\r\n', '\n')
   FH.close()
   return data
  
def reverse_complement(seq):
   import string.maketrans
   tt = maketrans('ACGT','TGCA')
   return seq[::-1].translate(tt)

def write_data(fn,s):
   FH = open(fn,'w')
   FH.write(s)
   FH.close()
  
def clean_fasta(s):
   title, seq = s.strip().split('\n',1)
   seq = ''.join(seq.strip().split())
   return title, seq
  
def make_count_list(L,kL=None,upper=True):
   # L is a list of strings (seqs)
   # each str the same length
   # count all chars
   assert len(set([len(s) for s in L])) == 1
   iL = L[:]
   if upper:
       iL = [e.upper() for e in iL]
   if not kL:
       kL = sorted(list(set(''.join(iL))))
   rL = list()
   R = range(len(iL[0]))
   for i in R:
       temp = [e[i] for e in iL]
       counts = [temp.count(k) for k in kL]
       rL.append(dict(zip(kL,counts)))
   return rL

Saturday, March 5, 2011

The need for speed: BLAT

I started testing alternatives to BLAST for aligning Illumina reads to a genomic sequence (Haemophilus influenzae L42023.1). When first setting up the project, I should have remembered how much faster BLAT is than BLAST. But I only needed to run the search once. I just went and made myself a sandwich.

Now it's time to do a little testing. I made a sample file test.txt with 10,000 sequences. I should really do command line testing (as we discussed here), but to make things easy I wrote a script to run the shell commands and print the time elapsed.

The actual commands are:

blat Hinf.fna test.txt blat.results.txt -out=blast8

megablast -d Hinf.fna -i test.txt -e 10 -m 8  > blast.result

And the output is:

blat
Loaded 1830138 letters in 1 sequences
Searched 529941 bases in 10000 sequences
1.1
megablast
16.9

BLAT is certainly faster than BLAST. Time might be the deciding factor if we do this every day, or if we have 2e7 reads. Another important issue is accuracy. We won't address that directly, but we should at least compare the results from the two methods. I wrote a second script to compare the output. The first two lines show that BLAT found matches for 607 sequences that BLAST missed (though the converse is true for a smaller number of sequences).

We investigate sequences where the length of the alignment was different. There are lots of these (922). A common pattern is for BLAST to report the same alignment, but trim it at the end of the read (presumably due to a mismatch). That suggests we might profitably trim the ends of all the reads before doing the alignment (say, to 40 nt---it would still be plenty). At least, BLAST reported the hit. But it's messing up the quality filter.

There is a rarer pattern where BLAT reports a match that is extremely short (see the fourth and seventh items in the list).

> python compare.py 
blat  8670
blast 8099
blat  only:   607
blast only:    36
351739
blast (1818932, 1818984) 53
blat  (626892, 626840) 53
295852
blast (937875, 937843) 33
blat  (937895, 937843) 53
514095
blast (545851, 545900) 50
blat  (545851, 545903) 53
514912
blast (677240, 677292) 53
blat  (677542, 677547) 6
518007
blast (580488, 580540) 53
blat  (580493, 580540) 48
517931
blast (493879, 493929) 51
blat  (493879, 493931) 53
462351
blast (46058, 46110) 53
blat  (46152, 46157) 6

That's a real problem because we're filtering the alignments for nearly full-length. I looked at the first one, it is a repeat.

BLAST:

514912 Hinf 100.00 53 0 0 1 53 760633 760685 1e-24  105
514912 Hinf 98.11 53 1 0 1 53 677240 677292 3e-22 97.6

BLAT

514912 Hinf 100.00 53 0 0 1 53 760633 760685 7.6e-22 102.0
514912 Hinf 100.00 47 0 0 1 47 677240 677286 2.3e-18 90.0
514912 Hinf 100.00 6 0 0 48 53 677542 677547 9.9e+05 12.0

OK..
This one is my fault. It is not the repeat that's causing trouble. BLAT split a single match into two pieces. I saved the results in a dictionary, and the duplicated key caused the other result to be discarded. My parsing code needs to be smarter and look out for this situation.

Here are the alignments from NCBI's web BLAST, so you can see the mismatch:

Query  1       TAGGAAAGATAATAATCTTTATCTATCAGATAATCTCGAGCATCTTTTGTTTC  53
               |||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  760633  TAGGAAAGATAATAATCTTTATCTATCAGATAATCTCGAGCATCTTTTGTTTC  760685


Query  1       TAGGAAAGATAATAATCTTTATCTATCAGATAATCTCGAGCATCTTTTGTTTC  53
               ||||||||||||||||||||||||||||||||||||||||||||||| |||||
Sbjct  677240  TAGGAAAGATAATAATCTTTATCTATCAGATAATCTCGAGCATCTTTAGTTTC  677292

script.py

import os, time, subprocess
L = ['blat','Hinf.fna','test.txt',
      'blat.results.txt','-out=blast8']
cmd = ' '.join(L)

print L[0]
start = time.time()
p = subprocess.Popen(cmd, shell=True)
pid,ecode = os.waitpid(p.pid, 0)
print round(time.time() - start,1)

L = ['megablast','-d Hinf.fna','-i test.txt',
     '-e 10','-m 8',' > blast.results.txt']
cmd = ' '.join(L)

start = time.time()
print L[0]
p = subprocess.Popen(cmd, shell=True)
pid,ecode = os.waitpid(p.pid, 0)
print round(time.time() - start,1)

compare.py

from utils import load_data

def make_dict(fn):
    N = 51
    data = load_data(fn)
    data = data.strip().split('\n')
    D = dict()
    for item in data:
        L = item.strip().split()
        title = L[0]
        align = L[3]
        if align < N:  continue
        i,j = L[8:10]
        D[title] = (int(i),int(j))
    return D
        
blast_dict = make_dict('blast.results.txt')
blat_dict = make_dict('blat.results.txt')
print 'blat ', len(blat_dict)
print 'blast', len(blast_dict)

blat_only = [k for k in blat_dict if not k in blast_dict]
print 'blat  only:  ', len(blat_only)

blast_only = [k for k in blast_dict if not k in blat_dict]
print 'blast only:   ', len(blast_only)

L = set(blat_dict.keys()).intersection(set(blast_dict.keys()))
L = list(L)

count = 0
for e in L:
    if count > 6:  break
    if not blast_dict[e] == blat_dict[e]:
        print e
        t = blast_dict[e]
        print 'blast', t, abs(t[0]-t[1]) + 1
        t = blat_dict[e]
        print 'blat ', t, abs(t[0]-t[1]) + 1
        count += 1
#print count