The Perl code had a bug, but I found it!

The bells-and-whistles version of the Perl model of USS evolution still had a bug, which became apparent once I fiddled the fragment scoring system to strongly favour good matches, and turned off mutation of the genome sequence (so only the fragments mutated). The bug manifested itself in the program cycles stopping, at fairly random points in the run (never stopping twice at the same cycle number or genome score, as far as I could tell).

After a LOT of careful detective work on my part, entirely unencumbered by knowledge of any Perl debugging tools, I found that a 'while' counter was being incremented at the wrong place (inside an 'if' instruction that was inside its 'while' loop, instead of just inside its 'while' loop). I still don't understand why this would cause the runs to stick at random points, but maybe the undergrad can explain it to me tomorrow.

(Confession added later: Solving the problem was not just the result of my careful detective work. The final discovery was helped by luck. I had added an 'else' statement to print a report that the next step had happened, but had incorrectly inserted one too many } brackets. In solving this I accidentally removed a different bracket than the one I had incorrectly inserted, which moved the while counter outside of the 'if' loop and, I discovered, eliminated the stopping problem.)

As usual, the Perl problem was the line feeds

Somehow, in being emailed to me, the USSv4a.pl and settings.txt files both acquired nasty Mac carriage returns instead of nice well-behaved Unix line feeds. This kind of problem has arisen often enough in the past that I knew to suspect it, but I had to figure out how Komodo deals with line feed issues before I could confirm that this was the problem and correct it.

Now the program runs fine, but it prints some reporting lines probably created by the undergrad while he was debugging the new bells and whistles....

...Three hours later.... I've found, understood and removed the unwanted reporting lines, and found and fixed a big mistake in how the program decided whether it was time to print out a report. And found and fixed several minor problems....

Not as simple a task as it should have been

Before he left last night the undergrad sent me what should be a fully functional version of our USS model, with the desired new features all working (they're described here).  But I just tried to run it on my home computer and absolutely nothing happens.  I enter "perl USSv4a.pl" at the Terminal prompt, and I just get another terminal prompt.

The problem is likely with the computer, not the program.  This is the first program I've tried to run on it since I replaced the hard drive a few months ago - as everything on the old drive was lost I was effectively starting from scratch.  I don't think I need to install perl - I've never had to do that before.  The USSv4a.pl file is in the Rosie directory, as is the settings.txt file it calls.  I could create a 'Hello World' perl file to test that perl is working but my perl book is at my office; maybe I can find one online ('Hello World' file, not perl book).

Yes I can (it's just one line), but now I discover that TextEdit will only save files as rtf, html, word and wordxml!  Wtf?  Guess I'd better install Komodo Edit.  OK, I remember that finding a Komodo Edit download site took a bit of looking (Komodo IDE was easy to find but costs a bundle), but Google found a site for downloading Komodo Edit.  Unfortunately what this has given me has a '.msi' suffix, and my Mac doesn't know what to do with it.  On looking more carefully at the source site, I see that it says something about bittorrent, which I don't know how to use.  Back to Google - OK, a site offering a Mac version, which is coming in (very slowly) as '.dmg'.

OK, got Komodo Edit, in Applications folder, opening for first time involves 'pre-loading' various bells and whistles, Mac gets VERY HOT (Fan Control maxed out).

OK, I created the Hello World file, and it ran fine.  Hmm, maybe there IS something wrong with the USSv4a.pl program.  Well, now I have Komodo Edit so I can take a look at it.  It looks fine, so I created a 'test' version with Hello World print lines at various places.  None of them gets printed.  I see that the undergrad has inserted a 'die' reporting step if the program can't open its settings file, so I think that's not the problem.  Instead Terminal acts as if the USSv4a.pl file and the new USSv4atest.pl file don't exist.

Back to the original USSv4a.pl file; now Komodo is claiming that every line fails to end with the expected "LF" end-of-line marker, but when I have it show the markers they all look fine.  The USSv4atest.pl file markers look exactly the same, and Komodo has no problem with them.

Time to quit working on this for now.  Maybe I'll send an email to the undergrad, or maybe I'll just wait till I get to my office (maybe not till tomorrow) and try to run it there.

Time to start preparing some talks

I have three talks to give in the next month and a half, so I need to start preparing them now.

First a 20- or 25-minute one at the annual workshop of the new CIfAR Program in Integrated Microbial Diversity, held somewhere not far from here, sometime close to the end of May. The guy in the next office invited me but he's out of town so I can't recheck the details. This talk will describe what we know about how natural selection has acted on the genes that cause bacterial genetic exchange. I think I can probably do this with slides I already have prepared.

Next, a 20-minute talk at a conference titled "Sex and Recombination: In Theory and In Practice", at the University of Iowa in mid-June. This talk will begin by introducing everything that the above talk will take 20 minutes to cover, and will then go on to explain how we are using computer simulations to understand how uptake sequences can accumulate in the genomes of competent bacteria. I hope to discuss results from our Perl model; the undergraduate is adding necessary embellishments (an oxymoron). (He thought he would have them all in place by late today but he left without passing the improved model on to me, so I suspect they're not quite debugged yet.) I don't yet have any model-specific slides for the talk, nor even a good idea of what the results will be.

And a few days after that, a 15-minute talk at the big Evolution meeting at the University of Minnesota. This will be on my work with the bioinformatician, on how the accumulation of uptake sequences in bacterial genomes has affected the ability of their genes to code for well-adapted proteins. Almost all the work is done here, and we have nice figures prepared for our manuscript. Unfortunately my collaborator has just redone some of the analysis and sent me new figures I don't completely understand. And her email response time has gotten very slow (though I guess this is only fair payback for my long silence while I was swamped with teaching).

My very old 'analytical' USS evolution model

I found my old USS evolution model in the files with my 1999 NIH grant proposal. It's not a computer simulation of evolution but analysis of equations describing an equilibrium. It uses only very simple algebra, so calling it an 'analytical' model is probably giving it more credit than it deserves. The introductory text gives an excellent description of the background, so here it is:

This model starts with the following assumptions:

1. H. influenzae cells have a preexisting DNA uptake system that preferentially transports fragments containing a USS0 (the 9bp core: 5'AAGTGCGGT). Fragments with imperfect (USS1) sites are not favoured.
2. Fragments of H. influenzae DNA are frequently brought into cells by this system.
3. Once inside the cell these fragments recombine with and replace homologous regions of the chromosome.
4. The DNA in the cells' environment comes from lysed cells ('donors') having the same distribution of USS0 and USS1 sites as the cells that are taking up DNA (the 'recipients').
5. Random mutation acts on both USS0 and USS1 sites to destroy and create new USS.

How these assumptions cause USS0 to accumulate:
Because the DNA uptake mechanism is biased in favour of fragments containing USS0 sites, any donor DNA fragment containing a new mutation that has created a USS0 will be taken up more efficiently than the wildtype version of the fragment. Similarly, donor fragments containing new mutations that have eliminated a USS0 will be taken up less efficiently than their wildtype counterparts. Consequently the DNA fragments within cells will be enriched for USSs relative to the external DNA, and recombination between these fragments and the resident chromosome will increase the number of genomic USS0s more frequently than it will decrease it. The bias will similarly affect the fate of new mutations arising within the recipient cell; mutations removing a USS0 will often be replaced by donor fragments carrying the wildtype USS0, whereas mutations creating new USS0s in the recipient will less frequently be replaced. This biased gene conversion can thus both compensate for mutational loss and amplify mutational gain of USS0s, and will cause them to be maintained at a higher equilibrium number than would be expected for a random sequence of the same base composition.

How the model works:
This model addresses processes acting within a single genome. The model uses the observed equilibrium numbers of perfect and imperfect USS sites to derive an equation relating transformation frequency to the bias of the DNA uptake system. This equation tells us the values that the transformation frequency and uptake bias parameters would have to take in order to be responsible for the maintenance of perfect and imperfect USS at the high equilibrium ratio that we observe.

The model then defines its terms ('variables'? 'parameters'?). It assumes that the frequencies of USS0 and USS1 observed in the real H. influenzae genome represent an equilibrium between forces that create USS0s and forces that convert them into USS1s. It then derives an equation relating the bias of the DNA uptake machinery (enrichment of USS0 over USS1) to the transformation frequency (the probability that a USS site will be replaced a fragment brought in by the uptake machinery in each generation). Conveniently, the mutation rate drops out of the equation.

Again quoting from what I wrote ten years ago:

What the final equation means: This equation tells us how frequent transformation must be (T), given a specified bias B of the DNA uptake system in favor of the USS, in order to fully account for the observed numbers of USS0 and USS1 sites in the Rd genome. The range of values is given below.












This is a pleasingly tidy result. It makes good biological sense, and the values are not unreasonable. Estimates of the actual bias vary, no doubt partly because they have been determined using USSs with different flanking sequences, but are usually between 10 and 100. We have no good estimate of actual transformation frequencies for H. influenzae in its natural environment, but if cells grow slowly, and frequently take up DNA as food then an average of 1% transformation per generation seems plausible, and even 10% not impossible.

I'm not very sure that this work is sound. The mutation rate dropping out is a bit suspicious, and I'm not sure how to interpret the transformation frequency when there are no other terms depending on generations (such as mutation rate).

Gene Transfer Agent evolution

My GTA colleague suggests that GTA may persist because of species-level selection. 'Species' is a tricky concept because these are bacteria, but we can simplify this to consider selection of lineages.

The basic idea is like that proposed to explain the surprisingly high frequency of bacterial lineages with defective mismatch-repair genes. Like most mutations, most GTA-mediated recombinational events will probably be deleterious. But some will be beneficial. Each time a beneficial recombination event occurs the lineage of cells descending from it will all contain GTA as well as the new combination. Provided the short-term and long-term costs of GTA don't cause the new lineage to go extinct or lose its GTA genes before the GTA-mediated beneficial change, lineages with GTA could take over.

Evolution-of-sex theorists have produced theoretical evaluations of such 'hitchhiking' processes, treating the allele in question as a 'recombination-modifier' (in this case the 'allele' would be the block of GTA genes). I haven't looked at the literature recently, but I think the general conclusion is that hitchhiking is common but weak; it is very unlikely to push an otherwise harmful allele to fixation. But these analyses weren't done for bacteria but for sexually reproducing eukaryotes, with the modifier controlling the frequency of meiotic crossovers between two other loci. I don't know how they would apply to bacteria.

Nevertheless, we know the ability of relatively rare beneficial events to preserve the GTA gene cluster must depend on how frequent the beneficial events are, how beneficial they are, how often GTA genes themselves undergo mutations that block their function, and how much (if any) harm the GTA genes cause. For example, if beneficial events are very rare, functional GTA genes may be lost by random mutation in the interval between beneficial events. Subsequent selection might cause the lineage to then go extinct, but it wouldn't bring the GTA genes back.

The important question is, how could we tell if this sort of thing is responsible for the persistence of the GTA genes? Short answer: I don't know.

USS don't accumulate because the bias is much too weak

In our very-preliminary version of the USS evolution model, we've been using a very simple scheme to score the similarity of DNA sequences to the USS motif. We just count the number of matches to the 10bp core USS sequence. Right now I'm keeping everything simple by running the model with DNA fragments that are only 13 bp long (and a genome that's 200-1000 bp long).

So a fragment with a perfect 10 bp match to the motif is only twice as likely to recombine back into the genome as a fragment with only 5 bp matching. We know from our earlier model, and from a calculation I did years ago, that the bias favouring USS needs to be much stronger than this if it is to overcome the randomizing effects of mutation. For example, (ignoring the effects of base composition) a sequence that matches the motif at 9 positions has 27 different ways to mutate into a fragment that matches at only 8 positions, and only one way to mutate into a fragment that matches at 10 positions. To overcome this disproportion, the bias favouring 9 over 8 (or is it 10 over 9?) has to be proportionally strong (i.e. 27-fold).

Now I need to find that old calculation - I think it might be with my 1999 NIH proposal.

Why don't USS accumulate in our model?

The undergraduate is working at home, making improvements to the use-ability of our Perl model of USS evolution, while I'm here doing some test runs with the slightly unwieldy version we have now.

This is the version that incorporates the first major improvement. Instead of a single "best-score" fragment recombining with the genome in each cycle, each of the scored fragments can recombine with the genome, with a probability proportional to its score. In principle this should allow the recombination to introduce new USS-like sequences faster than they are lost from the genome by mutation. But in practice the USS-score of the genome drifts up and down but doesn't consistently increase.

"When in doubt, run more controls."

So I've made a modified version of this program that has the same feature as our PositiveControl.pl program (it's name is PositiveControlv2.pl). Instead of mutating random fragments and recombining them back into the genome, it replaces their sequences with perfect USS cores. Of course the fragments need to be the same length as the core.

OK, this positive control tells me that the program is doing what it should. 10 fragments are recombining into the genome each cycle, which is what should happen to the 10 fragments the program considers each generation. Now I'll test PositiveControlv3, which puts in an imperfect USS instead of a perfect one. If the program is doing what I want, on average 8 of the 10 fragments will recombine each generation, and the genome score will not get as high as with perfect USSs.

OK, that worked too. So I think the problem may not be with the implementation of the code, but with the design of this version of the model. Back to the drawing board... (literally, back to the big whiteboard in the hall outside my office).

Komodo rocks!

Today we started using Komodo Edit as our editor for Perl.  It's much better than Mi:  better colouring (the numbers are all in red), more reliable indentation, and it catches syntax errors on the fly.  And, very cute, when I typed "if (" , it automatically added the second bracket ")" on the far side of my cursor, to make sure I remembered that I needed to close the (if condition) brackets.

Our most immediate goal is to transfer some of the features of our old/abandoned program into this new one.  This includes some basic/sensible features: reading the parameter settings for each run from a separate file, rather than changing the code of the program, and having a mechanism to neatly end the program and save the interim work when it's interrupted with a control-C.

There are also a couple of quite clever features specific to the issues our models raises.  One is when to print interim reports on the model genome's status.  Depending on the parameters being tested, runs can take anywhere from a few hundred to a hundred thousand or more cycles, with rapid change only at the beginning.  Rather than printing reports at a fixed interval, we have the reports printed whenever the number of cycles increases by a specified percentage.  This gives very frequent updates at the beginning of a run, and reports at increasingly long intervals as the number of cycles gets large.  

The other clever feature is how the model decides that the evolving features of the genome have reached an equilibrium.  This is checked by comparing the recent state of the genome (mean of some attribute in the most recent print interval) with the mean of the same attribute in the previous interval and over the whole run.  By using the increasing print interval as the unit of measure, we get sensitivity appropriate to the rate of change. 

No time for Perl today

But I did spend some time with a post-doc working on the Discussion of her manuscript. I had forgotten that last time we worked on it we tore the existing Discussion into shreds and came up with a new organization. So this morning I was discouraged to see that we had improved our Discussion out of existence, but our new organization is so much better that we soon had at least half of the text in place.

On the Perl side of things, a very helpful commenter (Neil) pointed out that finding the missing/extraneous curly bracket would have been easy if we were using an editor that highlights and validates syntax. We're using an editor called 'mi'. It lets us specify that our text is Perl, and uses colours to distinguish between different kinds of text (comments are red, text to be printed is grey, functions are green, operators are sort of purplish, 'while's and 'if's and 'my's are blue), but it doesn't sort out hierarchical stuff like the levels of brackets, and the levels of indentation keep going to hell (possibly my own fault). I'd love to hear about a better Perl editor for Macs, if any reader knows of one.

Who knew that 'while' loops can't be nested...

Substantial progress on the Perl model of USS evolution.

First the undergrad and I added the code that tallies up the scores of every sliding-window position in the genome. It didn't take us long to get it running (a few stray semicolons, etc.) However after now Keith's comment on yesterday's post I think there are more efficient ways to do what we've done.

Then we created code that does the recombination a completely different way, so each of the fragments being tested has its own chance to recombine (probability of recombination depends on its score). The code wasn't tricky, but getting it running took ages of tracking down stray curly brackets and discovering that I can't nest two 'while' loops (I changed one of them to an internal 'if' test).

I still need to add in the feature I came up with yesterday, that writes out the genome sequence and the tally of scores only at specified intervals, but now I have to dash off to a 'visioning education' meeting imposed on us by the administration.

Positive control progress on the USS model

The simplest version of our new Perl model of USS evolution has progressed to the state where it runs correctly. This afternoon I've been doing lots of runs, both with a 'positive control' version that replaces a random genome position with a single perfect USS core in every cycle, and with a test version that mutates random fragments and scores them for goodness of match to the USS motif, and then recombined the best-matched one back into the genome. Tomorrow the undergrad and I are going to create a modified version, to try a different way of having the fragments' scores determine whether they recombine with the genome.

With the positive-control version I've been examining the effect of changing the genomic mutation rate. If the mutation rate is zero, the only limit to USS accumulation is the way insertion of new USS disrupts existing USSs. (This happens only because each 10bp fragment is changed to a perfect USS before recombination, and so bears no relation to the original sequence at that position.) Not surprisingly, more USSs are present at equilibrium when the mutation rate is zero, and fewer when the mutation rate is 0.01 or 0.05 changes per position per cycle. The rate of increase in genome score is largely independent of the mutation rate. Because only a single USS is inserted per cycle, the number of cycles to equilibrium depends on the length of the genome.

Wait - good idea! I think we need to add some code to give us the frequency of each sliding-window score at the end of the run. This would let us make a histogram of how many USS-like sequences the genome has at the beginning of the run, and how many it ahs accumulated at the end. Basically, as the sliding-window is scoring match to the motif at each position, it should record the score in a tally (number of sites scoring 0, number scoring 1, number scoring 2, ...... number scoring 10). I could write some inefficient code to do this (barring about a thousand syntax errors - I really should go back and reread the first few chapters of Beginning Perl for Bioinformatics), but this sounds like something the undergrad might have learned an efficient way to do.

Did I learn anything else from the positive control runs? If the genome is very short the program runs very fast but the scores are noisy (no surprise there). I learned that I have no practical insight into how a sequence's USS 'score' reflects the quality of its matches to the motif - that's why we need the tally. I played around with the 'threshold' we use as a cutoff for insignificant matches to the USS consensus, but I think we can't really understand what this accomplishes until we have the score tally.

I also did some runs with the test version (not the positive control). The results of these mostly served to reinforced the importance of the genomic mutations. Under the present recombination system, USS can't accumulate in the genome because they mutate away faster than they're improved by recombination. I tried turning off the mutation of the genome, so that mutation only happens to fragments that are about to be scored for possible uptake. Even with this 'cheating', the genome's USS score crept up slowly and then plateaued at what looks (without the tally) to be a genome with only weak USS sites.

Model systems in evolutionary biology

I spend the weekend, with the rest of my lab, at the Evo-WIBO meeting of evolutionary biologists.  Over breakfast we got into a discussion of the role (non-role?) of 'model systems'.  I did my usual rant about how evolutionary biologists don't even understand the concept of model systems, but I'm wondering whether the problem is partly just the nature of evolution research.

Molecular biology, biochemistry, cell biology and physiology have made dramatic advances, largely because many of them work on the same organisms, so that the findings of one study can be directly used as the groundwork for more studies.  Evolutionary biologists (and ecologists) almost always work on different organisms, and although they publish lots of nice papers these rarely can be applied to studies by other research groups.  

This is partly tradition - a mark of academic independence in evolutionary biology seems to be choosing your own research system (organism+field site+questions of interest), but it might also partly arise from the nature of the field.  The process of evolution is intrinsically tied more to variation than to shared properties - natural selection acts on differences, not similarities.  So maybe choosing to work on different systems just looks like the sensible thing to do.

But it's consequences are unfortunate, because although every nice bit of research claims to have big-picture implications, the lack of transferability means we haven't really gotten anywhere.

Outline of the perl program

(in response to good advice in the comments)

Below is just a list of the main sections of the program, in its present 'test' incarnation.

MAIN PROGRAM:

1.  Get parameter settings from a file (except it doesn't, the settings are hard-coded in this version).

2. Create a random-sequence 'genome' of the specified length and base composition.

3.  Simulate a set number of cycles (presently 100), each consisting of genome mutation, fragment creation, mutagenesis and scoring, and recombination.    

     3A.  Mutate the genome by randomly changing bases with a specified probability.  (This step should be later in the cycle, not here.)

     3B.  Select a specified number of segments of the same lengths, from random positions in the genome.  This will represent fragments in the external gene pool.

     3C.  Record each fragment's sequence and 5' end position.

     3D.  Mutate each fragment's sequence.

     3E.  Score each fragment's sequence for goodness of match to the uptake sequence motif, using a sliding window.  I think the sliding window scores are not being correctly cumulated.

     3F.  Choose the fragment with the highest score.  Put its sequence at the corresponding genome position, replacing the original genome sequence of this fragment.  (This is a simulated form of recombination by gene conversion.)  Any mutations of this fragment that occurred in step 3D will thus become changes in the genome sequence.

     3G.  Score the genome for how well its sequence matches the USS motif.

4.  At the end of all the specified cycles, stop and report.

-----------------------------------------------------

SUBROUTINES:

I.  Creating the original random genome sequence:  This is pretty simple; it just picks bases randomly, with probabilities specified by the base composition.  

II. Mutating the genome or fragment sequence:  This is more complex, partly because the mutations need to maintain the base composition (see subroutine III), but mainly because it does it a relatively non-obvious but more efficient way.  It first decides how many mutations to make, by dividing the genome length by the mutation rate and taking the integer value.  (Oops, this will only work if the genome or fragment is big enough to get more than one mutation per cycle. The 'test' version has only a 100nt genome and a specified mutation rate of 0.001, so it has a real mutation rate of zero.)   The subroutine then randomly chooses positions for this number of mutations, and makes the mutations at these positions.

III. Doing the calculations for the mutagenesis probabilities:  This creates arrays holding the mutation probabilities for each base (A or G or C or T) to mutate to each other base.

I've got to stop this and work on my course's final exam for a while.

Still reading code

I'm working my way through the undergraduate's Perl code, annotating it with detailed comments to explain to myself what I think is going on.

I'm only about 30% through the code, but I think I've found a number of problems, one of them big.  Tomorrow morning I'll tie down the undergraduate and go over everything with him.  Or maybe tomorrow early afternoon - tomorrow morning I'm supposed to help one of the post-docs polish her talk for this weekend's Evo-WIBO meeting.

Reading code

Our Perl-programming undergrad has just sent me a copy of the latest version of his program simulation the evolution of uptake sequences by molecular drive.  So far I've gotten to about line 100 and found several trivial typos and one source of confusion (to me).  I had thought that the order of steps in each cycle was as follows:

1. Choose random fragments of a specified length from genome and mutate them (as if they came from different daughter cells).
2. Score each fragment for its match to the USS consensus.
3. Mutate the original genome according to the same procedure used on the fragments.
4. Replace the corresponding segment of the genome with the fragment that has the best USS score.

But the standard version of this code seems to instead do the following:

1. Mutate the whole genome.
2. Chose random fragments and mutate them (again).
3. Score each fragment for its match to the USS consensus.
4. Replace the corresponding segment of the genome with the fragment that has the best USS score.

So the fragments are getting mutated twice.

In actuality, this 'test' version of the code has a couple of steps commented out, and short-circuits the fragment-generation and mutation steps by simply specifying the sequence of every fragment (as a perfect USS).  I think this makes it a lot easier to confirm the the code that does the scoring is working as intended.  Tomorrow I'll sit down with the undergrad and go over it.

Why "Expelled" is a bad movie

The guy in the next office just got off the phone with a colleague, arguing whether atheists should rise up and express their views. Triggered by the pending visit of Richard Dawkins to UBC. This horrible movie is one of the reasons I think we should:

Expelled

Classes are over!

Let's see if I can get back to posting every day.  First let's see if I can remember what's been on the back burner.

The USS-evolution computer simulation is up and running, and the programming assistant will be able to spend more time working with me and the post-doc to get it doing what we want.  It's close, so I'm hoping for lots of advances and discoveries.  I've signed up to give a short talk about this at an upcoming meeting on sex and recombination (in Iowa, just before the big evolution meeting in Minnesota at the end of June).

The "How do USS constrain genome function" project with the out-of-town bioinformaticist is ready for its final polishing.  A couple of weeks ago she sent me an email which (I think) contains the final data, and it really shouldn't take long now to have the manuscript ready for submission.  I've signed up to give a short talk on this work at the evolution meeting.

The post-doc who's been analysing the variation in competence in a diverse set of H. influenzae strains now has a manuscript that doesn't need a lot more work.  She's going to give a short talk on this work at next week's Evo-WIBO meeting (evolutionary biologists in Washington, Idaho,British Columbia and Oregon).

Both of the molecular biology post-docs have data that has yet to be put into manuscripts (at least I have yet to see the manuscripts).  One has been analysing how CRP binds to recognition sites, and the other has been doing microarrays to find out how CRP and Sxy regulate genes in E. coli.

On the teaching front, I still have to:  fix up the final exam so its a valid assessment tool for our homework-research project as well as for students' understanding; grade 17 term papers that are evaluating intelligent design as a scientific alternative to natural selection; help the graders grade the other ~75 project reports; help our homework grader finish grading the last two homeworks, and prepare detailed keys for these; analyse and post the marks for the 'clicker' questions the students have been answering in classes; administer and help grade final exams for about 360 students; and get all the grades analysed and submitted.

And on the homework-research project front, I (and the wonderful teaching fellow I'm working with) still have to finish a proposal for a small grant to hire assistants to assess the quality of writing in papers and exams by this year's and last year's students; find someone with sufficient Excel skills to transform our clicker-collected survey data into something we can work with; read the literature (I'm hoping the teaching fellow will point me to the appropriate papers); analyse the data; and write the paper.

Insights from a visitor

We are blessed this week by a visit from a potential collaborator - a computer scientist who's done work on uptake sequences.  When we described our computer-simulation model of uptake sequence evolution, he quickly discovered a serious problem arising from the way we have set up the uptake bias function to work.

In our model, a pool of DNA fragments is created each cycle, and the fragment with the best uptake sequence score gets to replace the homologous segment of the genome whose evolution is being simulated.  The uptake sequence score is determined by using a sliding window to compare the fragment's sequence to the designated uptake sequence.  For a long fragment this score is expected to reflect the combined effect of multiple uptake sequences.  But at present the model is using short fragments, so the score is likely to just be that of a single uptake sequence. 

This means that, once the genome contains a 'perfect' uptake sequence at one position, a fragment homologous to that position is expected to always out-score any other fragments in the pool, and thus always be the one that replaces the resident sequence.  Thus one good match prevents the gradual evolution of other not-quite-as-good matches at other positions.

There's a different way to do the competition that doesn't have this problem.  Rather than having multiple fragments 'compete' for the best score, with the winner taking the only opportunity to recombine with the genome, we can have each fragment in the pool independently challenge the odds of being taken up.  Uptake of one fragment would not affect the chance of any other fragment being taken up in the same cycle.  

This will require some rewriting of the code. But it shouldn't be a big deal, and luckily our undergraduate programmer will have more time to work for us once exams are over at the end of the month.

Can I remember how our USS-evolution model works?

Today we're having the first lab meeting in weeks. (When my turn came around a few weeks ago I just kept canceling them, but now I'm starting to see the light at the end of the teaching tunnel.)  We're going to discuss an issue that's arisen in the USS-modeling work being done by an undergrad research assistant, but first I promised to introduce this project.  What can I remember (or rediscover by reading my old blog posts about it)?

The big goal is to simulate how uptake sequences accumulate in genomes of competent bacteria, under the combination of mutation pressure (a randomizing force) and biased uptake preferring fragments containing these sequences.  The model follows a single genome-sized sequence through repeated cycles in which 

1. Random segments of the genome are treated as if they were fragments in an external DNA pool released by descendants of the 'index genome'
2. These fragments are scored for quality of their match to the ideal uptake sequence.  The best fragment is chosen for the uptake step
   
3. In the conceptual meantime, the index genome itself undergoes random mutation, becoming the descendant index genome.
   
4. The chosen fragment's sequence replaces the homologous sequence in the descendant  index genome.
5. This recombinant sequence becomes the new index sequence and the cycle starts again at step 1.

There have been lots of issues to resolve along the way (how the mutation steps maintain the base composition of the sequence, how the uptake sequences are scored), but we finally have a program that runs.  It seems to be working correctly, but the undergraduate who's done most of the work tells me that it isn't causing any uptake sequences to accumulate.  He's quite a sophisticated undergraduate - he has a Biochemistry degree under his belt and is nearly finished a second degree, in computer science - and he's done a lot of statistical analysis to look for the expected accumulation.

I suspect that the problem is that so far the model is using inappropriate parameters (mutation rates too low or too high? uptake bias settings too weak or too fussy? numbers of cycles too short?).  Today's goal is to figure out what these might be.