Estimating equilibrium with noisy data

I'm trying to find a better way to identify an 'equilibrium' genome score for our simulations of uptake sequence evolution. Ideally I would just let the runs go until the score was stable, but two factors complicate this. First, the scores are very noisy due to genetic drift, so there is no stable core. In principle this can be smoothed out by taking the mean score over a sufficiently long interval, but in practice the runs take a long time.

I have been working on complicated averaging strategies to deal with this, but for many runs I thought I could just stop trying to be economical of computer time. So I queue'd up a few very long runs, but this didn't really solve the problem, because the mean scores just drift up and down.

Finally I tried just plotting the 'grand mean' scores of the up-from-random-sequence-genome and down-from-seeded-sequence-genome runs ('up' and 'down' runs) on the same graph, and taking the average of the final values. The 'grand mean' score is the geometric mean of the genome scores at every cycle from the start of the run. For up runs it's a underestimate of the recent mean, and for the down runs it's an overestimate, but in both cases its much smoother than the recent mean. I'll paste a graph in later, but these give what look like pretty good stable estimates of equilibrium values. So I think I'll just state in the methods that runs were continued until the up and down grand mean scores differed by no more than two-fold, and then the final scores were averaged to give an equilibrium score.

Simulating uptake of degenerate USSs

In addition to using our uptake simulation to generate surrogate data for uptake of chromosomal fragments, I think I can also use it to generate surrogate data for the uptake experiments that start with degenerate pools based on perfect uptake sequences. In these experiments we'll generate a pool of synthetic DNA fragments (~70 bp in the latest plan, I think, although our CIHR proposal says 200 bp) with a 30 bp internal USS that is 12% degenerate at each position (maybe it's degenerate throughout its length. We plan to reisolate fragments from this pool that cells have taken up, sequence these, and analyze the patterns of variation.

I can give the uptake simulation a 'genome' consisting of just this fragment sequence, and instruct it to take up fragments of this full length (so every fragment will start out the same). I can then set the fragment-specific mutation rate to 12%, so that each fragment will be independently mutagenized to this level before being scored for possible uptake. I can test the effects of different uptake biases by using different matrices for the scoring, and different ways of using the matrices (additive, multiplicative). (I'll also ask my mathematical colleagues if there is any function that gives effects intermediate between additive and multiplicative.) If our proposed fragment is degenerate only in its USS, I'll use a USS-length fragment; if it's degenerate throughout I'll use the full length sequence.

Then I'll give the pool of simulated-uptake sequences (and a control pool of degenerate but unselected sequences) to the post-doc for analysis. He's already done a very nice preliminary analysis of the simulated chromosomal fragments.

I've spend the last 48 hours watching simulations of uptake sequence evolution run (very slowly) and fiddling with code and queue-submission files and settings, but now it all seems to be falling into place.

First the former post-doc said that she thought the int was indeed an error, and that it probably arose when the coding-expert undergrad copied the code from another probabilistic step to the recombination step. She also made the excellent suggestion that we use the additive matrix in the US-variation paper.

Now that the erroneous int has been eliminated from the simulation code, runs with additive scoring do indeed work. So do runs with the less-extreme tenfold bias, and runs using the Gibbs position weight matrices derived from the H. influenzae and N. meningitidis genomes. None of this give as strong accumulation as the 100-fold bias matrix I have been using, but that's OK.

Testing the additive scoring led to several other improvements. I increased the resolution of the recombination decision by introducing a lower-cutoff for the random number range, as I mentioned in the previous post. I also found that when using the additive scoring we no longer needed to have the bias gradually decrease during each cycle, but could just leave it fixed. And I rediscovered a problem with scoring the genome additively. [Aside: scoring is done in two situations - first to decide how well the different positions in an individual fragment correspond to the matrix, so a decision about recombination can be made, and second to track the accumulation of preferred sequences in the whole genome, to follow the progress of the simulation.] When individual fragments are scored, only the highest position-score is used, and the scores of all other positions are discarded. But when the genome is scored we add up the scores of all the positions in the genome. With multiplicative scoring the low scores are many orders of magnitude lower than the good scores, so including them in the sum has little effect. But with additive scoring even really lousy matches to the matrix have significant scores, and the sum of these swamps out the effect of the rarer good scores. The solution to this problem is to only use additive scoring for the individual fragments, and to always score the whole genomes by using the matrix multiplicatively.

I also spent some time relearning how to run the simulation on the computer cluster in Nova Scotia that a colleague kindly gave me access to, and modifying the program to run properly in this different situation (just trivial fiddly stuff). I kept queueing runs on this system, thinking they were running OK, and then looking at the error file )"just in case" and finding that they had aborted at the start because forgot to specify a file, or specified it wrong, or specified it in the wrong order. But now I have about 15 runs going overnight!

My new plan for this part of the US-variation manuscript (it's most of the Results) is to start with the additively scored matrix, both because that's how position-weight matrices are commonly used in other situations and because this means we don't have to explain about reducing the bias during the cycles. I'll start with the 100-fold matrix, and show how we detect equilibrium, and that genome length doesn't affect the outcome (score per kb). Then I'll explain about mutation rates (genomic and fragment-specific). The preliminary runs I did today indicated that the mutation rate of the fragments doesn't affect the equilibrium score, which I continue to find surprising. I had suspected that this was somehow due to using the bias-reduction feature with the multiplicative scoring, but now I see it with additive scoring and fixed bias. But I have a more detailed set of runs going overnight which should make the situation clearer. Then I'll consider how the amount of recombination (amount of genome recombined in each cycle) affects the outcome, and the effect of fragment length on number and spacing of uptake sequences

Then I'll consider the effects of different matrices, first the ten-fold matrix I tried today, and then the multiplicatively scored 100-fold and ten-fold matrices, and then the Gibbs matrices. This last step links the model back to the Gibbs analyses that the Results starts with.

Two more ideas about scoring sequences with a matrix

1. While the post-doc and I were fussing with the Perl simulation, figuring out the best settings for scoring sequences with the Gibbs position-weight matrix, we also tested the effect of calculating the score in each window of the sequence by adding the scores of the individual positions in the window, rather than by multiplying them. In fact the original versions of the simulation did add the positions cores, but we changed that because it didn't seem to give enough resolution (not enough difference in score between what our experiments had told us were good and bad uptake sequences).

Below I think I'll be able to explain how to get better resolution, but when the post-doc and I were discussing the additive and multiplicative versions we realized that of course neither additive nor multiplicative has any chance of being realistic. That's because the proteins that interact with the sequence don't use arithmetic to decide whether to act on a particular sequence. Instead each base at each position will affect the probability of uptake in its own way. Thus the additive and multiplicative ways of applying the matrix should probably be seen as the extreme outliers of how the real interactions might work.

So when we present this analysis in our proposal to NIH, I think we should include the additive case as well as the multiplicative case, with the no-selection sequences as control. The additive case will show only modest selection for fragments containing good matches to the matrix, and the multiplicative version wills how very strong selection. And if we're lucky enough to have some real uptake data by then, we can compare that to the two kinds of simulated data; it should fall somewhere in between. (Ideally I'd put a graphic here showing what I expect these to look like, but I'll leave that for the post-doc to do once he's analyzed the simulated data.)

2. Here's a way to improve the resolution of the additive selection: In the previous post I said that the simple selection method started by choosing a random number from the range between zero and the maximum possible score. But the range doesn't need to start at zero, and for the additive matrix it certainly shouldn't. Instead maybe it should start at the minimum possible score, or at the average score, or at some higher-than average score.

Why would we have used "int(rand(" rather than "rand("?

On Tuesday night I successfully modified our Perl simulation model so it printed out lists of sequences of simulated fragments that it had found acceptable, based on their scores using a position-weight matrix. But I couldn't get it to choose the best fragments scored with the Gibbs-derived matrix; instead it just gave me a list of randomly chosen fragment sequences. Scoring with a fake matrix worked OK.

So tonight I was reading through the code, trying to remind myself of how the model uses each fragment's score and Perl's random number generator to decide whether or not to accept the fragment for recombination (or, in this case, for printing). The logic is simple, but its implementation surprised me.

Here's the logic: First choose a random number from the range between zero and the maximum possible score. Then compare it to the fragment's score. If the score is lower, discard the fragment. If the score is higher, accept the fragment. So a fragment with the best possible score will always be accepted, and a fragment whose score is 50% of the best will get accepted half the time.

But the code the simulation actually used added a step. When it got each random number, it converted it to its integer value (instead of rand($max) it did int(rand($max). This didn't create a problem with our fake matrix, because the preferred bases at each of its 10 positions all had value 10, giving the total matrix a maximum score of 10^10, and taking the integer values just stripped off irrelevant decimal points. But the Gibbs matrix has true probabilities in its elements, so its maximum score is only 1.9 x 10^-7. This means that all the random numbers will be ≤ 1.9 x 10^-7, so their integer values will all be zero. And this means that all the scores will be larger than the random number, so all the fragments will be accepted.

I don't know why the int function was included in the original program, but we tried deleting it and the program still ran fine. And now the program works as it should with the Gibbs matrix, preferentially accepting fragments with higher scores. So I used it to collect the positions of 200 preferred fragments from each of the forward and reverse-complement strands of the first 50 kb of the H. influenzae genome, and an equal number of control fragments (with the program set to accept every fragment it scored.

As a bonus, solving this problem may also have solved the problem I was having with the real version of the simulation - the one we are using to model the evolution of uptake sequences. In the past I couldn't see any evolution of uptake sequences when I used a Gibbs-derived matrix, only with the fake matrix. At the time I though this must be because the Gibbs matrix didn't incorporate strong enough selectivity, and that was worrisome for my understanding of how uptake sequences evolve. But now I suspect that the problem may just have been that int function!

But I still would feel better if I knew why the student who did the coding for us included the int step. Did we discuss it and decide that this was the right thing to do, or did he just do it on his own initiative? I've checked the old versions of the program and it's there from the beginning. Maybe it's some sort of Perl good-practice thing.

Preliminary data goals for the NIH proposal

I took a turn at lab meeting this afternoon, and used the time to brainstorm about what experiments and analyses we should try to get done before the NIH deadline in Early February. The work for the first NIH Aim will also be preliminary data for a resubmission of the CIHR proposal I just submitted, should that become necessary. I had hoped to also go over other work for the CIHR resubmission but didn't get to it before my time was up, so I'm going to do that next week (remembering to also update the lab meeting schedule).

Goals for Aim 1: Here we have two goals. The big one is to reisolate enough clean DNA from the periplasm that we can either send it for Illumina sequencing (if we have the opportunity to do this - the post-doc has a contact) or characterize it with conventional sequencing, and then do some analysis.

The second goal is to generate and analyze some surrogate sequence data, just to show that we know how to analyze it. By first doing this surrogate analysis, we'll be all set to quickly analyze whatever read sequence data we are able to get before the NIH deadline. Made bold by my recent success in adapting our Perl simulation of uptake sequence evolution to score DNA sequences, I've claimed that I can easily modify the simulation code so it generates sets of DNA sequences from a genome, a random-fragment set and a high-uptake set. Ideally these sets would contain fragments with a range of lengths, but for now I'll just generate several sets of each type (e.g. a 0.5 kb random set, a 1 kb random set set, a 2 kb random set) and mix them together.

Goals for Aim 2: The big issue here is whether we can successfully reisolate input DNA from the cytoplasm. This DNA is expected to be single-stranded, and to be only 5-10% of the total DNA in the cytoplasm (the rest is double-stranded chromosomal DNA). Traditional methods of separating ssDNA from dsDNA rely on their differing physical properties; both hydroxyapatite and BND-cellulose can be used. An alternative I just thought of is to use random hexamers to prime DNA synthesis on the single-stranded fraction of total cytoplasmic DNA, including a biotinylated nucleotide in the dNTP mix. The new DNA can then be separated from the rest using streptavidin agarose or beads, and the template strands eluted by denaturing the DNA.

One complication is that the chromosomes will contain short single-stranded segments, but these are likely to be a small contribution relative to the input DNA. Another complication is that the DNA in the cytoplasm is transient. We'll use a rec-1 mutant so the DNA doesn't recombine with the chromosome, but degradation by cellular nucleases is certainly a concern.

Goals for Aim 3: Here we plan to characterize the results of chromosomal recombination, by deeply sequencing a very large pool of DNA from cells transformed with DNA of another strain. Because there is no DNA reisolation step, I think the best preliminary data we could have is the sequence of one transformed strain (lab strain transformed with DNA from a diverged strain). We can easily get this information from one Illumina lane, and we may be able to get one lane done well before the NIH deadline.

The plan is for the post-doc to first transform a cloned novobiocin-resistance (novR) allele into the most transformable of the three sequenced diverged strains. This can be done very quickly because the former post-doc left a tube the novR DNA, and the transformation can be done by overnight co-culture of cells and DNA. The next step is to grow up and prepare DNA from one of the novR transformants, and to use that DNA to transform the standard strain KW20. Then we (i.e. the post-doc) just need to prep DNA from one of the transformants and pass it on to the Illumina specialists. They will shear the DNA to the appropriate length for Illumina sequencing, ligate the Illumina adapters onto the ends, and sequence it for us. We can use this sequence to specify the expected density of recombinant segments in the pooled chromosomes we propose to sequence, and the length distribution of recombination tracts.

Juggling data into the US-variation manuscript

I have two uptake sequence analyses I want to publish, and I'm trying to decide whether they would fit into the US-variation manuscript. One is the analysis of covariation between different USS positions and between different DUS positions (see here for example), and the other is BLAST analysis of the amounts of within-species sequence variation seen at different USS positions (see here). Both are certainly about variation in uptake sequences, and in previous blog posts I've described them as parts of this manuscript, but since then I've gotten bogged down in thinking that they don't really belong.
But let's assume that they do belong in this manuscript, and consider the next question of where would they fit best and how they would be connected with the other parts. They have to go after the first major section, which describes the analyses of the H. influenzae and N. meningitidis genomes with the Gibbs Motif Sampler, because they use datasets this generated. Should they go before the other major section, which describes the Perl model of uptake sequence evolution, or after it? These two major sections fit well together, because the model makes predictions that can be tested with the Gibbs datasets.
What if the model went first, using the generic uptake sequence? It could then be followed by the Gibbs analyses, applied both to the model's output and to the real genomes. And then by more versions of the model, using the position-weight matrices generated from the real genomes. We could then do the other analyses of the genomic data....
Here's an attempt at an outline using this order:
Intro:
Emphasize the goal of understanding the cause of uptake sequences. Summarize what we know of their properties and the evidence that they arise by point mutation and spread by transformation more efficiently than other sequences because they are preferentially taken up (cartoon figure contrasting drive and beneficial-variation models). We have developed a model of this evolutionary process, which we describe and test below.
Results:

1. The model (how it works, with a cartoon figure).
2. We run the simulations to equilibrium (evaluated as score or US count), using a 10 bp generic position-weight matrix.
3. Characterization of model: Effects of (i) genome length (each cycle takes longer); (ii) mutation rate (getting to equilibrium takes more cycles); and (iii) recombination rate (determines equilibrium score, below saturation).
4. Properties of equilibrium sequences (i): Proportions of perfect and mismatched occurrences. Use Gibbs to find them? Compare with direct counts of perfects, one-offs and two-offs? Explain why Gibbs is, in principle, better?
   
5. Properties of equilibrium sequences (ii): Spacings between occurrences (found by Gibbs. These are more even than random, as are the spacings between real uptake sequences. The spacing depends on the length of the recombining fragments.
   
6. Use Gibb to reanalyze the Perl output sequences.
7. Use Gibbs to reanalyze the genome sequences.
8. Repeat the key Perl runs using the position-weight matrices from the Gibbs analyses of real genomes.
   
9. Compare what we find.
10. Do other reanalyses and new analyses of the Gibbs output of the genome analyses. (Variation in subsets of genome sequences (no interesting results but worthwhile anyway) This would include the covariation analysis and the BLAST analysis of within-species sequence variation.

That seems workable, and as sensible as the previous order. Now I need to decide whether it's sufficiently better than the old order to be worth the trouble of rearranging the whole Results section.
A couple of hours later: I don't think it is. For now I'm going to treat the manuscript as two distinct parts. First, the Gibbs analysis of the genome and sub-genome sequences, and analyses using the occurrences these identify, including the covariation and within-species variation analyses. Second, the Perl simulation of uptake sequence evolution, and comparison of its results with the genome analyses.

More about Gibbs scores

Yesterday the post-doc and I were again trying to understand the scoring system used by the Gibbs Motif Sampler. It seems peculiar to us because, in a search where most occurrences of the motif get scores close or equal to 1.00, it produces a lot of sites that score zero, as illustrated in the histograms above. In general, sites that are scored higher also better match the position-weight matrix that describes the motif, but the scaling seems weird. All of these sites, even the ones scoring 0.00, are pretty good matches to the motif (see below, copied from earlier post), and we don't understand why there are a lot of sites with 0.00 scores but few or (usually) none with scores between 0.01 and o.50.

The Gibbs documentation isn't very helpful, describing the score as "the probability of the motif belonging to the alignment". And the legend provided with the scores at the end of the run only says "Column 6 : Probability of Element". I decided to email my helpful Gibbs expert how the score is calculated, but while searching for his email address I found that I'd already asked him in early 2007 (!) Here's what he said:

Well, it’s a little confusing. For the near optimal solution (the one labeled NEAR OPTIMAL RESULTS), the probability given is the probability of the particular site according to the motif model for that motif type. It is calculated by computing the posterior probability of each of the sites according to the motif probability model. In your first example, its calculated as P(A in position1)*P(C in position 2)…. The probability is normalized so that the sum of the probability of the particular site belonging to each possible motif type and background sums to 1. The first site below is pretty definitely a good fit for the model, the second is not so great and might be better described as a background site. Both of the statements you list below are true in the sense that the probability correlates with how well the site fits the model.

Why are there low probability sites in the model? When the MAP is calculated, the program tries to maximize the MAP. Sometimes, there is a marginal improvement to the MAP by adding these low probability sites. The frequency solution will usually not include these sites, so in some ways it gives a better picture of the motif model.

First about the 'normalization': I don't understand what he means, but I think I don't need to worry about it. I've drawn a sketch of what I think might be going on - the output first calculates a score for every position in the genome, and then puts only those above some criterion into its alignment. But before it reports those aligned sequences it recalculates their scores, maybe using some normalization procedure that stretches out the middle, cramming most of the scores that would have been in the lower half down to zero.

Now about the 'background': In one sense his statement that the low-scoring sites in the alignment should be treated as background seems to be wrong, because even the sites with scores of 0.00 have significant matches to the matrix (see figure below). But maybe sites with this strength of match to the matrix are expected to occur by chance in a genome sequence of this length and base composition, so we can't assume that they result from selection for the process/effect of interest. In that sense they might indeed be considered background. I would think that the decision about background status would also have to take into account how often such sites are expected by chance and how often they're seen. For example, if we found 100 occurrences where theory predicts 10, we would conclude that most are due to non-background forces (though we wouldn't know which ones).

Silane-treated coverslips?

I'm hoping to get started learning how to work with the optical-tweezers apparatus soon. One of the things I want to test is adhesion of H. influenzae cells to silane-coated coverslips, as this will eliminate the need to restrain the cells by suction from a micropipette. The RA has already found that cells don't stick to uncoated coverslips, nor to ones coated with BSA. I had been hoping that I could just order them from Sigma or somewhere, but now it looks like preparing the coverslips might be the first of many very fussy procedures.

A google search didn't find anyone selling them, but it found lots of protocols for making them. They all involve precleaning the coverslips by boiling in strong acid or soaking in strong base or scrubbing with bathtub cleanser (Bon Ami was used), and then dipping or soaking them in solutions containing 3-aminopropyltriethoxy-silane/epoxy-silae/TMA silane/trichloro-octenyl-silane/methacryloxypropyltrimethoxysilane, and then air-drying or baking or ??? and storing in a dessicator or under unspecified conditions. I have no idea what kind of silane I should use or how to prepare the coverslips, so I've emailed the author of the Bacillus subtilis paper that did tweezers analysis, to ask what she did. I've also emailed my physicist collaborator (whose apparatus I'll be using) to ask if she has experience with this.

I think I'll also email around our institute, to find out if anyone here is routinely using silane-treated cover slips. Maybe they can (1) give me one to test, and (2) show me how to make them.

To do list

Last night's post began to list the things I need to get done in the next ~6 months. I'll try to elaborate on it now, adding other goals.

Finish and submit the US-variation manuscript: This is the only manuscript we have in the works right now (except one on teaching that's fallen by the wayside). Hmmm, that in itself is probably a cause for concern.

Write the NIH proposal: Now that NIH has cut its proposal page limit to 12 pages, including figures but not references (?), perhaps we can just revise (and revise and revise) the 10-page proposal the post-doc submitted for his fellowship application. Yesterday I discovered that HERRO, the UBC group that supports research in the health sciences, is offering a half-day workshop on writing NIH grants at the end of the month. They don't seem to have advertised this to the UBC community, just announced it on their web page, but a colleague and I tracked it down and signed up.

Get preliminary data for the NIH proposal: The proposal has 3 (or maybe 4) parts, each consisting of sequencing specific very diverse DNA samples we've purified from transformed cells. For the first and third of these parts we can do a small-scale experiment to show that the components all work, if necessary doing the preliminary sequencing by conventional methods rather than Illumina. The second part is the most difficult, but we should be able to get preliminary data there too.

Get more preliminary data for the CIHR proposal: This will be critical if this proposal is rejected. The first part of this proposal is the same as the first part of the NIH proposal, so that new data will do double duty. But we also need more data for the second part, making and analyzing mutants and testing candidate DNA-binding proteins, and for the third part, optical tweezers measurements and studies of USS conformation and flexibility. The optical tweezers part is the most risky but that makes it the one most in need of evidence that it can be done. It's also my own project and I want to get started having fun with it soon.

Grant proposal's done!

I clicked "Submit" on Sunday night. It was wonderful writing it with the research associate and post-doc - we made a good team. If the proposal doesn't get funded it will be only because we don't have much preliminary data yet.

But I didn't get to relax because I headed off to Chicago on Monday morning to give a talk at Northwestern School of Medicine. I learned lots of interesting things in my talks with various faculty and students. One was something I should have considered for the proposal. It discussed the problem of fitting DNA through the narrow secretin pore in the outer membrane, especially if the DNA has to share the pore with a type IV pilus. The same pore is used to pass proteins out of the cell, and the same size problem arises - some of the proteins are probably too big to fit through the pore. So maybe we're wrong in attributing to the pore the size it appears to be when purified - maybe it's stretchy. Another particularly interesting thing was some unexpected properties of the DNA uptake specificity of Neisseria gonorrhoeae. I won't spell them out because they're not published yet, but we'll want to do some experiments to see if the same things happen in H. influenzae.

Now it's time for some careful thinking about what to try to get done in the next 6 months. I plan to submit a proposal to NIH in early February, and we need to obtain preliminary data. So we need to decide what data will be the most compelling. I also want to be prepared with more preliminary data in case CIHR rejects our DNA uptake proposal - that would need to be resubmitted by March 1. Some of the data we'd be getting for NIH will also be good for CIHR, but others will be different. And I need to get the US-variation manuscript finished. As I recall, I need to run some more simulations and then the former post-doc and I need to finish the writing.

A better analysis of the data

The previous post had a rather silly way of handling the matrix-score data. Now I've just graphed the actual scores on a log scale.

One line of Perl code changed = new data!

This afternoon the post-doc and I were discussing the uptake-sequence analysis we're proposing for the CIHR grant. We wanted a way to compare the actual uptake of genome fragments with that predicted by their scores with the matrix I'd derived from the genome sequence using the Gibbs motif sampler. I realized that we might be able to use my existing Perl program (the one that evolves uptake sequences in simulated genomes) to score these.

So, after I burned out on rewriting paragraphs in the proposal, I opened up the giant Perl program US.pl. In the 965 lines of code I found just the spot I wanted, and told the program to print the score it had just calculated. I put the first 10 kb of the H. influenzae genome sequence into a text file, and changed the program's settings file so it would score this sequence but wouldn't try to evolve it. And I created a scoring matrix using the output from one of my old Gibbs runs on the H. influenzae genome.

It ran! I pasted the output into a Word file and converted the '-' signs of the exponents of the scores (range -32 to -8) to a tab mark, and reopened the file in Excel. I then subtracted these values from 32, so now the lowest pseudoscore value was zero and the highest was 24. Then I plotted this 'score' against the row numbers

And here's the result (black bars). To check that I hadn't screwed up, I looked up the position of the highest score in the Excel file (score = 24, at position 5319) and checked the sequence I'd started with. Sure enough, at position 5319 I found AAAGTGCGGTCAATTTTTCTGGTATTTTTT, a near-perfect match to the USS consensus.

What this graph doesn't show very well is the large number of positions with very low scores. So here's another graph (blue bars) to the same scale, but showing only the first 300 positions. Now we see that most positions have scores less than 10.

Versions of these figures will go in the Appendix of the proposal, as preliminary data.

H. influenzae and influenza

Haemophilus influenzae got its name because (1) it likes to grow on medium that contains blood (it needs heme and NAD, which blood has lots of) and (2) it was identified as commonly present in lings and sputum of patients ill with influenza in the 1918 global epidemic. It was initially suspected of being the causative agent of influenza, before the virus was discovered. If you're interested in the history of medicine, I strongly recommend the book The Great Influenza, by John M. Barry.

Anyway, the new H1N1 influenza epidemic probably means we'll see a rise in H. influenzae pneumonia. I'm mentioning this in my grant application, in the paragraph on the medical relevance of H. influenzae.

A terminology problem

I've only just realized how strongly the terms we use have been affecting our thinking about uptake sequences. I'll try to write about the problem here, but the nature of terminology problems makes it hard to write about them...

First some history of uptake sequences: It's been known since 1974 (Scocca, Poland and Zoon) that competent H. influenzae cells take up H. influenzae DNA fragments in preference to fragments from unrelated species, and since 1979 that this is a sequence preference (i.e. not a modification preference). In 1980 Danner et al. showed that H. influenzae cells preferentially take up fragments containing the 11bp sequence AAGTGCGGTCA, and that ethylation of positions in and close to this sequence interferes with uptake. They correctly concluded that the species preference is seen because H. influenzae DNA is greatly enriched for this sequence (or for similarly-acting sequences). [They presciently also proposed that the uptake bias might be at least partly responsible for this enrichment, a hypothesis we've since formalized in our molecular drive model.] Later Goodgal showed that the 9 bp sequence AAGTGCGGT was probably sufficient to give strong uptake. In 1995 Smith et al. analyzed the newly available H. influenzae genome sequence for the 9 bp sequence and found 1465 (later 1471) perfect matches. Alignment of these also showed a strong consensus for two 6 bp AT-rich segments on one side of the 9 bp sequence. And in the last five or ten years we showed that these sequences are the consensus of a diverse motif that accumulates in the genome by point mutation and homologous recombination, and attempted to measure the preference of the uptake machinery for different variants of this consensus.

The terminology problem: So we need terminology for several different phenomena, all of which have at least sometimes been called uptake sequences or uptake signal sequences (USSs). There's the short sequence initially identified as the uptake sequence; the 9 bp version of this is commonly referred to as the uptake signal sequence or USS and is often the only sequence considered in analyses of uptake. Then there's the consensus derived from analysis of the genome sequence, with the 9 bp sequence called the core USS and the flanking AT-rich segments called segments 2 and 3. Then there's the position-weight matrix description of the genomic enrichment, which we refer to only in convoluted ways. And finally there's the true bias of the DNA uptake machinery, which we have only very poor estimates of and no term for.

Until recently nobody knew enough about DNA uptake and the genome to clearly distinguish between these, and this confusion has led to a lot of muddled thinking. We've been among the main perpetrators of this confusion, using 'USS' interchangeably for genomic sequences and uptake biases. But now we're proposing to obtain a very detailed description of the 'true' uptake bias, and to compare it to the genomic motif, so I need to carefully choose the terms I use, and to explicitly draw the reader's attention to them.

I could call the genomic motif the genomic USS motif, and its consensus the genomic USS consensus. We think it really has evolved as a motif, just like those associated with the preferences of known DNA-binding proteins, but we haven't yet identified the proteins and we aren't defining it by the proteins' preferences (which we don't know) but by its abundance in the genome. Although this consensus is useful in some experiments and simulations, it isn't much use when considering what's actually in the genome. If the genomic USS motif does approximate the real bias of the uptake machinery, then counting the number of occurrences of the core consensus (1471) is quite misleading about the actual uptake of different parts of the genome. (If it doesn't, counting is even more misleading.)

Should I call the real uptake bias just that (real uptake bias)? Or the real uptake motif? I could leave off the real, and just use uptake motif. Like the genomic USS motif, it would be described by a position-weight matrix. What would I call the consensus of this, which I am proposing to use in additional experiments in place of the genomic USS consensus? Should I avoid using the USS abbreviation when referring to the matrix and consensus of the real uptake bias? This would help the reader keep things straight, but would mean getting rid of the customary description of this bias as the USS.

I could keep it simple, using genomic motif and genomic consensus, and uptake motif and uptake consensus. Or I could insert USS into each of these terms (genomic USS motif etc.). I think I'll try the latter for now - it's a bit wordy but has the virtue of being unambiguous, and if I later come up with something better I can easily make the needed changes.

Later: I overlooked one more terminology problem. How do I refer to individual sequences in the genome that might (or might not) promote uptake, but that haven't ever been tested? We used to distinguish between fragments that contain 'USS' and those that don't, but the motif/molecular drive perspective predicts a smooth continuum between sequences, with those that strongly promote uptake at one end and those that resist uptake at the other. It gives no justification for distinguishing between 'real' USS and sequences that look a bit like USS. Should I call them all USS-like sequences? Candidate USS sequences? Fragments with different degrees of correspondence to the genomic USS motif?

Proposals plan

Last week I went down to Bellingham to collect pearls of wisdom from a friend and colleague who has lots of experience with NIH. I came away with lots of good ideas for strengthening our NIH proposal, but the most important was the decision to delay submission to the next deadline, in February. This will not only let us generate more preliminary data but will let me get the US-variation manuscript submitted. This submission is important because I now realize that its contents will play an important role in the proposal.

So for the next two weeks I'll only have to think about the CIHR proposal. It's due Sept. 15, but I'll need to click the submit button Sept. 13 (Sunday), because I'll be traveling on the 14th and 15th. (The paperwork (electronic and real) needs to be in a few days before that.)

Although the reviewers thought that the previous version of the proposal was well written, I now see that the overall structure was weak. In particular, the Introduction (a short section, only 3 paragraphs) introduced competence but failed to introduce the problem/s I want to address, the hypothesis/hypotheses I was testing, or the approach I will take to test the hypothesis/hypotheses. That was because I hadn't done a good enough job of spelling these out even to myself. So here's a new attempt:

The phenomenon: Bacteria actively transport DNA molecules into their cells, across their outer and inner membranes.

The specific problems this proposal addresses: How does H. influenzae transport DNA across its outer membrane: (a) How is this uptake initiated? (b) How are long fragments pulled into the periplasm? (c) Why is uptake so strongly biased towards specific sequences?

The hypotheses: (1) Initiation and progression of DNA uptake are mechanistically distinct processes. (2) Uptake sequences are favoured because they facilitate the DNA deformation needed to initiate uptake. (3) Progression of uptake uses an ATP-driven ratchet driven by assembly and disassembly of type four pseudopili (this hypothesis may well be wrong, but that's OK).

The Specific Aims: (1) We will fully characterize the uptake specificity, using both synthetic sequences and the full H. influenzae genome. (2) We will identify the proteins that contact DNA during uptake and those responsible for the sequence specificity. (3) We will test whether the uptake sequence motif polarizes uptake and whether it promotes DNA kinking. (4) We will characterize the forces acting on DNA during binding and uptake.

I'm going to paste this summary on the first page of my draft proposal, so I incorporate it into the Introduction and use the Background to establish the problems and build to the hypotheses.

Learning from the experts, part 2

The other expert I consulted this morning works on type four pili, the filaments we think are responsible for pulling DNA across the outer membrane. Here are some things I learned from her.

We shouldn't put too much faith in the present measurements of, for example, secretin structure (secretin is the protein pore the filaments pass through). It may be wider or more flexible, which would allow a pilus plus DNA to pass through.

She agreed that, if H. influenzae uses short pseudopili to bring in DNA, the pilus has to function as a ratchet that (i) binds DNA at the cell surface, (ii) retracts and pulls the DNA across only the thickness of the periplasm, (iii) releases the DNA (perhaps because the pseudopilus has become too short to hold on to it), (iv) elongates, and (v) binds the same DNA at a more distant position before retracting and pulling again. She thought this quite unlikely.

We need to revisit the question of whether competent cells might make long pili, rather than just short pseudopilus stubs that don't protrude beyond the outer membrane. We can look for pilin subunits in the culture medium of competent cells, but we will need a good anti-pilin antibody. We may be able to get this from a colleague, or she can help us design a peptide for a commercial antibody. We also may need to add BSA to the mix to prevent the intrinsically stickly pili/pilin from being lost by binding to the tubes and tips. And to add EDTA to inhibit any secreted metalloproteases that would degrade the pili (Vibrio has these but I don't know if H. influenzae does).

If cells make long pili, the DNA might wrap around them, eliminating the need for the improbable pseudopilus ratchet mechanism.

She agrees with a reviewer of our earlier proposal that secretin is quite stable, so that we may be able to purify it to test whether it binds DNA. But this would be a part of the proposal that I don't think would be very strong, so we might instead eliminate it in favour of the search for long pili.

She thinks that H. influenzae (and B. subtilis?) might have an as-yet-unrecognized protein that does the same job as PilT, using ATP hydrolysis to power pilus retraction. She mainly works on the toxin-coregulated pili of Vibrio cholerae - these don't themselves retract, and they have no PilT, but they are responsible for pulling filamentous phage into the cell so must be able to exert a similar force.

She thinks that the pilC gene encodes a crucial inner membrane protein that makes the connection between PilB/PilT (on the cytoplasmic side of the inner membrane) and the pilin subunits at the base of the pilus. This is likely to be the weak link in the strength of the pilus, as the pilin subunits are tightly braided together.

Learning from the experts, part I

This morning I went to the university across town to consult with a couple of experts about components of my CIHR proposal. My first meeting was with the physicist whose optical tweezers apparatus expertise I plan to rely on. Here are some issues we discussed:

To get evidence for USS-dependent binding, we should collect data for a histogram of rupture forces (i.e. how strong a pull is needed to break the bead free of the cell), for DNAs with and without USS. If binding is distinct from uptake, this distribution should be bimodal, probably with a much higher peak for the low 'binding' rupture force. She said this would be quite a lot of work (weeks, at least).

Can we stick H. influenzae cells onto a glass coverslip, (as the B. subtilis people did) rather than restraining them with a micropipette? This would make the manipulations much easier. She recommended trying BSA-coated coverslips as well as silane-coated ones. The coating will also help reduce non-specific sticking of DNA to the glass.

To test whether uptake is due to 'diffusive motor' or a 'power stroke' (I'll need to read up on these concepts) we should create a force-velocity curve. If we find that, below a critical stalling point, the velocity of DNA uptake is independent of the opposing force pulling the bead back into the center of the laser trap, we would conclude that uptake was due to machinery exerting a power stroke, whereas if we observed a gradual slowing of uptake with increasing opposing force, it would be due to a diffusive motor. I think Tfp use a power stroke, and a process driven by flow of ions (maybe the PMF (= proton motive force?) would be diffusive. The B. subtilis rate is largely independent of the opposing force on the bead, although it appears to somehow be generated by the PMF.

If uptake is indeed mediated by a short pseudopilus that acts as a ratchet on the DNA, we might be able to see the ratcheting in a velocity vs time (or force vs time?) graph. This would depend on the resolution of the velocity (or force?) measurements. If the pseudopilus is only, say 20 nm long (a guess of the thickness of the periplasm), we would need to detect fluctuations on this scale.

In principle optical traps can exert forces greater than 100 pN on trapped beads. However DNA deforms at 65 pN, allowing it to stretch. Not only will such stretching complicate length/velocity measurements but the deformation may perturb interactions with the uptake proteins, causing them to lose their grip on the DNA. So we should keep our measurements below about 60 pN.

She recommended making a flow chart to lay out the tests we will do and how the results will affect our hypotheses and guide our next steps.

Six sentences in Science!

Looking to Bacteria for Clues

IN HIS NEWS FOCUS STORY “ON THE ORIGIN of sexual reproduction” (5 June, p. 1254), C. Zimmer highlights the importance of the phylogenetic perspective championed by John Logsdon, but by considering only eukaryotes he overlooks an important bacterial clue to the evolution-of-sex puzzle.

Until recently, bacteria were thought to be sexual; they have well-characterized processes that cause recombination of chromosomal alleles, and these parasexual processes were assumed to have evolved for recombination in the same way as meiotic sex in eukaryotes. However, a more critical analysis of the genes responsible for the parasexual processes suggests that they did not evolve for sex after all. Instead, the chromosomal recombination they cause appears to arise as unselected effects of related processes, the evolutionary functions of which are well established (1).

The fact that bacteria lack genes evolved for recombination indicates that meiotic sex must have evolved in eukaryotes to solve a problem that bacteria don’t have. Bacteria apparently get whatever recombination they need by accident—why do eukaryotes need so much more?

ROSEMARY J. REDFIELD
Department of Zoology, University of British Columbia, Vancouver, BC V6T 3Z4, Canada. E-mail: redfield@interchange.ubc.ca

We need a hypothesis.

We have lots of hypotheses, Lord knows, but I've just read the NIH advice for preparing a strong proposal, and the first thing they emphasize is the importance of starting with a clear hypothesis that the proposed work will test. Of course at one level I know this, but at another level I've been ignoring it. The series of experiments I would like to propose fit together nicely, building a coherent picture of the limits to recombination at both molecular and population levels, but it's not obvious to me that they test a specific hypothesis.

I can come up with a hypothesis for the first parts (the ones covered by the current post-doc's NIH fellowship proposal). "Every stage of transformational recombination is limited by sequence biases." This hypothesis might well be false - we only know of biases at two steps: DNA uptake and mismatch repair. But it provides a clear explicit framework for finding out where and what the biases are.

But I would like to extend the research to include the biases due to the differences in different isolates. Well, that's a way of saying it that sort-of fits the hypothesis but doesn't really capture what I want to propose. Instead I should first describe what I would like to propose, and then try to find a way to describe it that can be part of a unified hypothesis.

I would like to propose to find out the molecular reasons why different isolates of H. influenzae differ so much in their ability to transform. This should be fairly straightforward, given the molecular and genetic/genomic tools we have. We can also try to estimate how long these strains have had their specific transformation phenotypes - did a very recent ancestor get a mutation (or recombination!) that inactivated a competence/transformation gene, or is it descended from a long lineage of non-competent ancestors.

But I would also like to find out how different are the recombination histories of strains that can and can't transform - looking not just at the individual sequences they have acquired, replaced, or lost, but at the pattern underlying these events. Do strains that don't transform at all in the lab have a history of greatly reduced recombination, relative to strains that transform very well in the lab?

I don't know how do-able this is. If we get genome sequences of enough strains, can we build up a detailed picture of their recombination histories? The previous post-doc used a program called Mauve to identify recombination tracts...

I also don't know to what extent this analysis will be facilitated by the information we'll get from the first parts described above. We'll first build up a complete picture of the relationship between (i) a pair of donor and recipient genomes and (ii) the recombinant genomes that result from DNA uptake. (If resources permit we can do more than one pair of genomes.) Then we try to use this knowledge to infer the recombinational histories of other genomes. Is there a hypothesis in there? Maybe the hypothesis is that this will be possible. We could hypothesize that "Understanding the limits to transformational recombination in lab cultures will let us identify the recombination histories of natural isolates."

But as science this is what I consider a 'pseudo-hypothesis'. It's not a hypothesis about the nature of reality, but a hypothesis about the nature of our abilities. How about "The limits to transformational recombination in lab cultures explain the recombination histories of natural isolates."

Hmm, if that's our hypothesis then maybe this proposal does belong in the Evolution of Infections Disease program after all, and not in the Prokaryotic Cell Growth, Differentiation and Adaptation program. I think I'll handle this question by putting together a decent summary (once we have a single hypothesis) and asking both Program Directors where they think it belongs.