Exploring ctypes (3)

Continuing with ctypes---here is an example of scoring a genome for potential sites for a DNA-binding protein according to the methods we talked about here.

We have a single function in a C module that will be compiled into a shared dynamic library. Then we load that library from Python, using ctypes, and call the function. The heavy lifting is done by the function score. The arguments are:

• the DNA sequence as a string (a str on the Python side)
• an array of floats which holds the scoring table
(in order acgt for each position in the site)
• the length of the site (n = 2 in the example)
• an array to hold the result for each position in the sequence

other variables:

• loop counters / index variables i, j, k
• the current character, c, so we can print it
• the current score for this char, d, also so we can print it
• the current score for this position in the sequence, r

The logic is straightforward and elementary. (Although perhaps there are faster approaches that people might suggest?)

mylib.c

#include <stdio.h> #include <string.h> void score(const char* dna, double p[], int n, double result[]) { char c; int i,j,k,N; double d,r; N = strlen(dna)-n+1; for (i=0; i<N; i++){ r = 0; for (j=0; j<n; j++){ c = dna[i+j]; switch (c) { case 'a': { k=0; break; } case 'c': { k=1; break; } case 'g': { k=2; break; } case 't': { k=3; break; } } printf("%c ", c); d = p[j*4 + k]; printf("%3.3f ",d); r += d; } printf("%3.3f\n", r); result[i] = r; } }

We build the shared library using clang (or gcc, a quick test showed no speed difference):

clang -g -Wall -c mylib.c clang -dynamiclib -current_version 1.0 mylib.o -o mylib.dylib

The Python script is a little more unusual. The array initialization was just adapted from the docs (here), I don't yet understand how it works. (It seems like a pointer would be better).

Here is the second of those:

NFloats = ctypes.c_double * N result = [0.0] * N result = NFloats(*result)

And I was surprised that we don't have to worry about the types for arguments to our score function---that's why those lines are commented out at the top.

For this toy example, the DNA sequence is 14 nt, and the scoring table is given in the list L. Here is the script:

script.py

import ctypes pre = '/Users/telliott/Desktop' mylib = ctypes.CDLL(pre + '/mylib.dylib', ctypes.RTLD_GLOBAL) #mylib.score.argtypes = [ctypes.c_char_p] #ctypes.c_float_p, #ctypes.c_int] dna = 'actgtcgactcgag' L = [ 0.567, -1.603, -0.2245, 0.3605, -0.1175, -0.4655, -0.5326, 0.6898 ] EightFloats = ctypes.c_double * len(L) ff = EightFloats(*L) n = int(len(L)/4) N = len(dna) - n + 1 NFloats = ctypes.c_double * N result = [0.0] * N result = NFloats(*result) mylib.score(dna, ff, n, result) for f in result[:6]: print round(f,3)

Here's the output. Each position is evaluated on a separate line. The last six values printed are the results received back in Python.

> python script.py a 0.567 c -0.466 0.101 c -1.603 t 0.690 -0.913 t 0.360 g -0.533 -0.172 g -0.225 t 0.690 0.465 t 0.360 c -0.466 -0.105 c -1.603 g -0.533 -2.136 g -0.225 a -0.117 -0.342 a 0.567 c -0.466 0.101 c -1.603 t 0.690 -0.913 t 0.360 c -0.466 -0.105 c -1.603 g -0.533 -2.136 g -0.225 a -0.117 -0.342 a 0.567 g -0.533 0.034 0.101 -0.913 -0.172 0.465 -0.105 -2.136

Not bad for 50 lines of code and a little more than an hour of work! If you want to see it go fast, comment out the printf statements in the C code, and change:

dna = 'actgtcgactcgag'*400000

Now the DNA length is 5.6E6, and it finishes in about 5 seconds. It's impressive.

DNA binding sites 5

Here is an example that shows how difficult the problem of searching for sites is, at least by the method we've used. I include two scripts in the zipped files: search.EC.py and plot.py. The site we want to search for is hard-coded in the file as crp.

As its name implies, the first script searches the E. coli sequence. You'll need to get the sequence first (and save it in the right place---see the script) before this will run. We slide a window over the length of the sequence, shifting one base at a time, and score what's showing in the window under the site's scoring scheme. Scores were multiplied by 10 (to convert to ints) and then saved to disk (22 MB or so).

In the second phase, the search.EC.py script also randomizes the sequence and re-runs the same search. The plot.py script filters the results for those above a threshhold and plots a histogram of the results (see the above graphic). The red bars are the results obtained with the authentic sequence and the yellow bars are with the randomized sequence.

The point is that although the extreme high values are clearly higher in the real sequence, for a scoring range like 110 - 120 (that is 11.0 - 12.0 in the original scheme), the ratio of the likelihoods for the two models (real E. coli v. random) is not even greater than 2. So, upon observing a site with a value of 11.5, say, its significance isn't clear.

One idea that would improve the significance is to note that the genome is subject to selection. For the Crp system to work properly, there has likely been selection against randomly placed sites, so the random model is not really the appropriate one to test against.

Again, zipped files here.

DNA binding sites 4

I grabbed the data for crp and purR sites in E. coli from George Church's server (here). This is the first part of the crp data set:

>aldB -18->4 attcgtgatagctgtcgtaaag >ansB 103->125 ttttgttacctgcctctaactt >araB1 109->131 aagtgtgacgccgtgcaaataa >araB2 147->169 tgccgtgattatagacactttt

The graphics above are representations of the information analysis for crp and purR produced by site_score.py. Notice the different patterns. The important bases for crp are two short pentamers separated by one turn of the helix. Not so for purR, which suggests a different mode of binding.

According to (Schumacher 1994 PMID 7973627):

The DNA-binding domain contains a helix-turn-helix motif that makes base-specific contacts in the major groove of the DNA. Base contacts are also made by residues of symmetry-related alpha helices, the "hinge" helices, which bind deeply in the minor groove. Critical to hinge helix-minor groove binding is the intercalation of the side chains of Leu54 and its symmetry-related mate, Leu54', into the central CpG-base pair step. These residues thereby act as "leucine levers" to pry open the minor groove and kink the purF operator by 45 degrees.

It's that minor groove interaction that is giving the strong signal in the "middle" of the site.

Here is what we calculate for site scores for crp. Notice that the lac site is a relatively poor one:

$ python site_utils.py crp tnaL 20.2 nupG2 18.7 lac 17.6 cdd 17.3 deoP2 16.7 malT 16.5 .. cya 12.8 .. crp 12.0 .. lac 8.8 .. avg for 100000 random seqs: -15.64 13.16 12.18 11.87 11.51 10.78 10.66 10.58 10.56 10.2 9.87 158 sites in random seq for cutoff = 5

DNA binding sites 3

Continuing with binding site analysis, first two posts here and here.

Tom Schneider also invented Sequence Logos, which display the information for binding sites in an intuitive, graphical way (Schneider 1990 PMID 2172928). In that paper, (following Shannon) they define an uncertainty measure for each position in an alignment:



where H(l) equals minus the sum over the four nucleotides of the frequency of each base b at that position times the log₂(freq). Then, the information is:



where e(n) is a small sample correction factor. Thus, uncertainty plus information is constant, and approximately equal to 2 (bits). It's no coincidence that lacking any information about which nucleotide is present at some position in a sequence, you need to ask me two yes-no questions to obtain the identity. For example: is it a purine? Yes. Then is is adenine? Yes. Two questions, two bits.

The script site_score.py does this calculation for the fis sites example, and we plot our home-grown version of the logo as the graphic below.


You can compare that output to what is in the paper:



The colors are switched for the central position because the values for A and T are exactly equal, and we sorted to plot T on top, while Schneider did the reverse.

There is also a site on the web for making logos. To use that, we need to strip the names out of the sequence file.

FH = open('fis.sites.txt','r') data = FH.read() FH.close() for line in data.strip().split('\n')[1:]: print line.split()[1]

The only significant difference is at the middle position (11). We didn't use the reversed sequences, so we see mainly A at that position. This is an artifact of the web site's approach.



Zipped project files on Dropbox (here).

DNA binding sites 2

The first few sequences in fis.sites.txt are:

#format Schneider 1 tttgccgattatttacgcaaa 2 agtgactaaaatttacactca 3 gtggtgcgataattactcata 4 attgcatttaaaatgagcgtg 5 attggtcaaagtttggccttt

The first few lines of output from site_utils.py:

ttatgtacaaatagtaagaaatgtctgaga.. [45, 10, 26, 39] 0.567 -1.603 -0.2245 0.3605 ttggtatatatactatacacctatatttga.. [28, 22, 21, 49] -0.1175 -0.4655 -0.5326 0.6898

For Schneider's approach each site is present in the alignment in both forward and reverse complement orientation. In the output we are looking at each column of the alignment. The first base 't' of the top line of output is from sequence #1, the second base 't' is from the reverse of sequence #1, the third base 'a' is from sequence #2, and so on. The counts are given next, and sum to N = 120 = 60 * 2. The scores are calculated as given previously:

2 + log2(freq) - 0.018

The scores for the first 11 sequences are:

1 12.2 2 11.8 3 9.0 4 6.5 5 12.2 6 8.5 7 8.4 8 4.6 9 12.0 10 5.3 11 10.4

If you look at the graphic from last time (or the paper) you'll see that we match. So I think we're doing things correctly.

The last thing we do when running this basic site_utils.py file as __main__ is to look at some random sequences. For starters, we calculate

avg for 100000 random seqs: -10.85

which is roughly -0.5 for each position in the 21 nucleotide alignment. According to my understanding the purpose of the correction term was to make this be zero, and I'm not sure why it isn't.

The top 10 sites from the random sequences had scores:

11.72 10.73 10.73 10.52 10.48 10.28 10.08 10.08 10.07 9.86

There were a total of 9 sequences with scores > 10, while only 14 of the 60 authentic sites scored that high. So depending on where our cutoff is for an authentic site we'll get a lot of false positives. For this example, with a cutoff of 5, I found 267 of them.

DNA binding sites 1

I haven't written much about classification of DNA binding sites or motifs on the blog. (I did a series on Gibbs sampling to find new motifs starting here).

But a couple of years ago, before the blog started, I did some work on this and there is a modest writeup on my pages at mac.com (here). At that time, I put a lot of effort into providing a Quartz GUI for the basic functionality, but I haven't tested it recently and I'd be pretty surprised if it still works. Still, that discussion (and see here) is reasonable. Let's see if we can improve upon it.

The reason the subject came up is in considering Cython (here). One use case that's definitely appealing is to predict new members of the family in a bacterial genome (given a set of known binding sites), where there are (just about) as many sites to be tested as base pairs. That can take a substantial time using Python (a few minutes).

To start working on this problem, we need some data. A repository of alignments for various transcription factors of E. coli is available from George Church's lab (here), but I'm going to use Tom Schneider's method (website) for constructing the scoring system, so I'll use the example from his paper about fis (Hengen 1997 PMID 9396807). Here is a graphic from that paper showing a few of the 60 sites they analyzed and the bit scores for each site. We're going try to to recreate this. I've got a file with the 60 sequences; you can get them and the code we'll be using from Dropbox (here).



As usual for these problems, we consider each position in the alignment to be independent, and add up the scores for all the columns to obtain a total score. The score for a single column of the alignment depends on the the count of nucleotides at that position. Suppose we have N = 20 sequences and this distribution:

A:5 C:5 G:5 T:5

We calculate the frequency as 5/20 = 0.25 for each base, and then the score is 2 + log₂(freq) - correction. The correction is for the small sample size of known sites (see here to begin) and it's about 0.018 in this particular case, but we'll neglect it for this brief discussion. The method essentially computes a log odds score for the competing hypotheses of a real binding site versus random sequence.

The score for each base in this alignment, which seems to be random, is calculated as 2 log₂(0.25) = 0. On the other hand, suppose we have:

A:18 C:1 G:1 T:0

Then the scores are:

A = 2 + log2(0.9) = 2 - 0.152 = 1.848 C = G = 2 + log2(0.05) = 2 - 4.322 = -2.322

We assign a pseudocount of 1 for T, even though no site had T, so the score for T is the same as for C and G. In this way we end up with an m x n matrix of floats, where m = 4 (ACGT) and n = the length of the alignment.

To score a candidate site, we observe each nucleotide in turn, and retrieve the corresponding score for that nucleotide and position in the alignment. For example, in the non-random case if the site contains A at that position we add 1.848 to the score, while if it has C we subtract 2.322 from the score. Here we see the great strength of this approach: we don't just reward sequences for being close to the "consensus", we penalize them for having a nucleotide that is rarely observed at the corresponding position in authentic sites.

That's the simple, basic idea. As I say, the files will be on Dropbox (here). Next time we'll show some results.