Whence Gibbs?

I successfully worked out how to command the Gibbs Motif Sampler to analyze the new genome sequences. I've only done it for two of them, because a better option has appeared.

A new version of the Gibbs motif software is available. It gives the option of using a 'centroid' sampling method that combines the best sites found in different runs (runs initiated with different random-number seeds), rather than simply taking all the sites identified in the run that had the best score. This has the big advantage of eliminating most of the weakly-matched 'false positive' sites.

It took me a few days to work out how to get it running on the computer cluster (the helpful administrator reset some permissions for me). The new release includes a version that runs in the Mac terminal, and I now have that working too. But it didn't take long to discover that it runs about 100-fold (no, I'm not exaggerating) slower than the usual (non-centroid) version. This means that a good run analyzing a whole genome would take several weeks (or more?); getting rid of the false positives isn't worth that big an investment.

But the very helpful Gibbs expert has again offered to help - he says the centroid version shouldn't be slower at all. So I've sent him the test file I've been using (2% of the genome) plus examples of the output I get. He's going to see if he can find the problem and fix it.

Lots of Gibbs search progress

Yesterday I worked out a way to nudge the Gibbs motif sampler into finding the Neisseria meningitidis DUS (their term for their uptake signal sequence). Even though the DUS is present in Neisserial genomes even more frequently than the H. influenzae USS is in its genome, the sampler couldn't find it without prompting. This may be because it's much shorter than the USS (only 12 contiguous bp vs 22 bp spread over 29 positions), or for some other reason I don't understand.

I didn't want to give the sampler a prior file specifying the pattern to look for, so instead I added two lines of fake sequence with a very high frequency of the DUS to the start of the genome file. This 'seed' was enough to get the sampler started on the right motif. Once it's started it has no trouble finding the DUS, and I can later delete the seeded DUSs from the list it generates.

This morning I obtained the A. pleuropneumoniae genome sequence my collaborators have been working with, split it into pieces, and generated reverse complements of both it and the N. meningitidis genome, and combined each genome's forward and reverse-complement sequences into single 'F+RC' files for searching. I did this because I need to have the sampler search both strands, and (I think) I have better control if I tell it to search just the sequence I've given it. The A. pleuropneumoniae USS is very similar to but even longer than the H. influenzae USS, so I did test runs with the 'prior' masking file I'd used for H. influenzae to make sure everything worked.

I did this and all my other tests using only 10% of the genome and only one orientation, because I wanted them to run very fast and because the guys who manage the computer cluster want all long runs to be entered through their 'Fair Share' queueing system. And now I've successfully queue'd requests for full-genome searches. I don't expect to get the results until tonight or tomorrow.

I also emailed my collaborators to let them know I'm finally back working on this project. The PI is on vacation, but the bioinformatician has been taking advantage of his absence to work full time on it! She's going to send me her new data and rewrite in a few days, so I'm not going to do any work on the manuscript until then. I could go ahead and do Gibbs analysis of all the genomes we might want to consider, but I think I should wait to see how the three main foci of our work (H. influenzae, A. pleuropneumoniae and N. meningitidis) fit into the manuscript.

Do USS constrain protein-coding?

About eight years ago I started collaborating with people with bioinformatics skills on an analysis of how USSs constrain (or don't constrain) the ability of the genome to code for proteins. The project is still unfinished, but it's made a lot of progress. I think I should make one last push to get it done before I dive into doing experiments.

Originally I was working with a biophysicist in Taiwan and his excellent graduate student. The student did a lot of nice analysis and was coauthor on one paper with us, but we never got the second part of the work finished (or published). He's since moved on to other things, and the analysis is being redone (new data is available, and we now realize the flaws in the original analysis) by a bioinformatician (bioinformaticist?) working at the National Research Council labs in Ottawa. I'm the main idea person, and the main manuscript-writer, and she's the person who can write Perl scripts and deal efficiently with databases.

Last summer, when we finished our first joint manuscript on USS evolution (pdf here), I did the sensible thing of writing a rough draft of this second manuscript before we'd done most of the analysis. I even put in mock-up figures of the results I expected (based on the previous work by the grad student). I'm told that the best scientists always have a pretty good idea what the paper will say before they do (get their students and post-docs to do) the experiments, but I'm rarely that far ahead.

Since then the bioinformatician has done quite a bit of the work, but I've mostly let my contribution slide while I did more urgent things. Today I read through my rough draft of the manuscript (a very nice aid to my lousy memory) and realized that I'm far from clear about what analysis has been done and what still needs to be done, by her and by me. I'm pretty sure that I just need to read back through our emails and associated attachments to get this clear.

One new bit of analysis will be Gibbs motif searches on the various genomes. I've already done this to death for H. influenzae but now want to do it with other genomes. Unfortunately I've forgotten such important basic information as how to connect to the computer cluster where I run the searches (username? password?), and how to format the search instructions. Not a big deal - I just wish I had been more organized in storing this useful information.

20 minutes later: I was unduly pessimistic. I started up Fugu (the program I'd used to interface with the computer cluster) and found that it had not only remembered the name and location of the cluster server, it knew my password. And I found a folder called "run pbs scripts" in the Gibbs folder on my computer, and this turns out to contain the instruction files I used when queueing my searches on the cluster server! (I had even forgotten that they needed to be queued.) Now I just need to get the genome sequences, and get them into the right format (Fasta, in big fragments?), and reread the pages explaining what my previous instructions meant. Then I can set up the new searches and put the files onto the server. I'll still need to log on to the server using Mac's Terminal interface, to put the searches into the queue, but I know those instructions are around somewhere....

New BLAST book!

Thanks to a suggestion from a reader, I ordered the O'Reilly Press book on BLAST. It just arrived and looks to be exactly what I and the rest of my lab need.

We all use BLAST all the time, but we've never really had any understanding of how our search query sequence became the search results. We sort-of knew that this was asking for trouble, but haven't taken the time to learn more. Probably this was partly because doing a BLAST search is so fast and easy that you want to use the results right away, not 'waste time' reading the manual.

The new book has a Glossary! (No more using Google to find hints of what the terms mean.) It has a detailed index! Chapter 2 has a section on Evolution, which opens with the wonderful statement that "BLAST works because evolution is happening."

Yesterday I used my newly gained ability to do local BLAST searches to set up a search for one of the post-docs. We blundered around a bit because I couldn't remember what the different letters controlling the parameter settings did. Now I have the book, all the information I need is at hand.

The only problem is that the book was published in 2002, and some details have changed. Right now I only notice that the BLAST web interface has changed a lot. The available version of BLAST has also changed, from 2.2.6 (new when the book was written) to 2.2.16. I suspect I'll need to read the book before I'll have the background to let me understand the changes.

Should we write a proposal to NIH?

Yesterday one of the post-docs and I discussed whether we should submit a proposal to NIH. Yes, we did just get one proposal funded, to work on the regulation of competence genes in H. influenzae and E. coli. But her project is completely different, and some parts of it are going to be expensive.

She's studying when, how and why different lineages of H. influenzae lose (or maybe gain) the ability to take up DNA and recombine it into the chromosome. We already knew that this occurs in various bacteria, and she's now completing a thorough analysis of the variation in DNA uptake and recombination ability in a broad selection of H. influenzae strains.

Some of these strains were chosen because their genomes have been or are in the process of being completely sequenced (one of the benefits of working on a sometime pathogen), and her next goal is to analyze these sequences for differences that could explain their different phenotypes. This bioinformatics work won't cost much except her time; we've already bought a nice fast computer for it. And her time doesn't cost the lab anything, because she's supported by a lovely post-doctoral fellowship from NIH.

But the next steps will be expensive. She wants to use the bioinformatics information to design investigations into the genetic differences of strains that haven't been sequenced. Her original plan was to develop a microarray chip containing all the genes and alleles that the bioinformatics and other work suggested might be involved. This still seems like a good approach, but the field is changing so fast that better ways to survey genomes are becoming available faster than we can keep up. One thing they have in common is that they'll all cost a bundle.

Subsequent work will also be pricey. We'll probably want to follow up the H. influenzae findings with investigation into related bacteria. This will be beyond the scope of the present post-doc, so we'll need new post-docs or grad students or technicians, as well as money for the tools and techniques.

The other reason to write a proposal to NIH is that proposal-writing is the best framework I know of for clear thinking about research plans.

Gibbs motif search progress continues unabated

The sxy manuscript has been on hold, partly because one of the two grad students involved in it was in the far north. But he's back, and the manuscript is close enough to being finished that I'm hopeful it will be done soon ('soon' being an elastic term here). So I need to switch my attention to it and away from the motif searches for the USS-defined manuscript that have been consuming my brain power lately.

But before I stop I'm seeing how much I can get finished. I ran and analyzed the leading-strand and lagging-strand searches - their motifs are indeed identical to the composite one I posted.

The more-stringent and less-stringent searches gave the results I expected (fewer and more sites with the motif, with stronger and poorer mean scores, respectively). I used the run that gave the most sites and the worst mean scores to do a correlation analysis. (Having more sites that are imperfectly matched to the consensus increases the power of this analysis to detect weak interactions.)

The goal of the correlation analysis is to find out whether the bases at different positions of the USS interact. For example, the most common bases at position 17 are A and G, and the most common bases at position 21 are T and A. If we find that the individual USSs that have A at position 17 usually have T at position 21, and those that have G at position 17 usually have A at position 21, we would conclude that the bases at these positions interact during DNA uptake. Said another way, we'd conclude that USSs with a G at position 17 function better if they have a A at position 21.

Results: MatrixPlot found only weak correlations between only a few adjacent positions in two clusters. A colleague has kindly used software he wrote to also analyze a preliminary data set for us; I'm going to ask him if he can test the big set. Before doing this, I realized that I only have half the data, as I only did the low stringency searches on the forward strand. So I've queue'd up more searches, with the same and even lower stringency, on both forward and reverse-complement strands.

And, finally, some of the gene searches are working, thanks to fine-tuning advice from the helpful Gibbs expert. These runs are searching the sequences of only the parts of the genome that code for proteins, to see if the direction of coding affects the motif. I had to split the gene set into four parts, and two of these managed to find their motifs. So I've queue'd up more replicates, using more seeds, and also runs looking for the reverse-strand motif.

And last night I read over the Introduction and improved parts of it, though it still needs more work.

Not as random as it should be

On Saturday I described a key step in the motif-search computer program thusly ("thusly"? Pomposity alert!):

The program uses a random-number 'seed' to start searching the genome sequence for sequences whose consensus motif fits this pattern. ... Once the score is stable for a 'plateau' number of cycles the search ends and the program stores the final set of sequences that gave this pattern and goes on to try another random number seed.

Well. the seeds are supposed to be random numbers. The program uses them to get approximately random starting patterns, from some of which it converge on very similar final motif patterns. In any one run only a few of any 100 seeds lead to any approximation of the USS pattern, with the rest stuck muddling around in poorly-matched pseudo-patterns. But of course starting with identical seeds will give identical outcomes

This morning I noticed that two replicate runs that had been put into the computer queue within a few seconds of each other had given exactly the same results for all 100 seeds they claimed to have used. The program doesn't report each of the seeds used, but it does give you the initial seed it started the whole run with. Surprise surprise, the identical runs had started with identical 10-digit seeds (1162169861). And the next run queue's had the number 1162169821, differing only in its second-last digit. And the next two runs also had identical seeds (1162169911) Not only are these numbers much too similar to be considered 'random' the identity of all the results implies that each of the other 99 seeds is completely determined by the first seed.

I had assumed that any respectable computer program would draw whatever random numbers it needed from the random-number generator that's built in to most programming language compilers. (That's what our Perl programs do.) These actually generate pseudo-random numbers, and a great deal of analysis has gone into designing and evaluating the algorithms they use (they're also used in computer-security applications).

But I was wrong. I just went back to the original program notes, and found that one of the many settings users can specify is a 'random number generator seed'. This feature isn't discussed anywhere in the notes, and none of the examples provided in the notes use it. But maybe it prevents exactly the problem I'm having.

I've sent an email to my helpful advisor, and in the meantime I've queue'd up substitute runs, spaced far enough apart in time that they should have distinct seeds.

Deciding how much (computer) data is enough

I need to decide how many replicates of the motif sampler analysis I should do for our 'defining the USS' paper. And at what level of stringency to do each replicate. For real experiments, such decisions are made by factoring in how much work each experiment will be, and how expensive the needed resources will be. But the computer cluster we're using is free to us, and starting a run takes very little of my time. (It's deciding what to do that takes the time.)

I'm starting with the first analysis: searching the whole genome sequence for sites fitting the USS pattern whose spacing is specified by my 'fragmentation mask' prior (a 10bp motif followed by two 6bp motifs separated by 1 and 6 bp. The prior file uses this line:
++++++++++x++++++-xxxx-++++++
where ‘+’s specify positions whose consensus is to be included in the motif, and ‘x’ and ‘-’ positions that are not to be included and that are optional, respectively.

The program uses a random-number 'seed' to start searching the genome sequence for sequences whose consensus motif fits this pattern. Initially it finds only sequences whose consensus motifs very weakly match this pattern, but it goes through many cycles trying to improve the motif its found to get a better match.

It scores the quality of the matches with a MAP number. The first cycles always give MAP scores much less than zero (-4000 is typical); scores improve in later cycles cycles. Once the score is stable for a 'plateau' number of cycles the search ends and the program stores the final set of sequences that gave this pattern and goes on to try another random number seed.

With genome-sized sequences, most of the time the score never gets much better than -4000 because the search hasn't found the real USS motif. But sometimes it gets lucky and finds the motif, and then the scores on successive cycles rapidly increase to values around +4000. Once the program has completed trying the specified number of seeds, it checks the results from all the seeds to find the one that gave the best MAP score. It then polishes this set of sequences to see if it can get an even better score, and then reports this as the final set of sequences and consensus.

So I control stringency by specifying two things. First, I specify how long the 'plateau period should be (how long the program keeps trying to see if it can find a better motif). Second, I specify how how many different random number seeds the program tries before it picks the best results for polishing.

The default is a plateau of 20 cycles and 10 seeds. But I've been advised to use a plateau of 100 or 200, and lots more seeds. One complication with running analyses on the computer cluster is that you have to specify in advance how much time your run will need (the cluster uses this info to set your place in the queue). I don't have a good idea how long my runs will take on this cluster, so I submitted a test run using plateau = 200 and seeds = 100. I set the time limit to 12hr because a much shorter run had taken only 2 minutes, but it's been going for 7.5 hours now.... If it runs out of time before it's done I suspect that I won't get any interim results, just a blank file. So a few hours later I queue'd up an identical run with a 48hr time limit. I hope that's long enough.

On Wednesday one of the post-docs and I are meeting with the local experts for this computer cluster. We sent them an email confessing that we don't really know what we're doing, and I suspect they feel they'd better straighten us out before we seriously mess up their system.

Motif analysis update, part 2: coding direction

The second of my three questions about the USS motif was whether, for the large subset of USSs that are in genes, the orientation of USSs with respect to the proteins they help code for affects their motif consensus. So my plan was to assemble all the coding sequences of the genome, all oriented in the direction their proteins are specified (not in the direction their DNA is replicated), and to then compare the motifs of the USSs in the two possible USS orientations.

Assembling the sequences seemed straightforward (download in one file from TIGR, remove unwanted characters). But the motif-search program couldn't find the USSs I knew were there (see last week's post). I spent a week or more trying more tests and variations, to try to figure out what was going wrong, because I didn't feel that I understood the problem well enough to clearly explain it in an email to the helpful expert. Was the number of sequences over the limit? Were the 'N's I'd had to insert causing problems? Were the sequences too short? Was the problem dominant or recessive to a well-behaved sequence?

Yesterday the same problem appeared in some new sequence files, and then was corrected (see previous post). I wasn't entirely sure what I'd done that made the difference, but this did give me confidence that the problem with my gene sequence files was in the formatting, and my prime suspect was the hated carriage returns. These are a nightmare for Unix beginners like me - they're often invisible, they come in several incompatible flavours (Mac vs PC vs Unix), and Unix/Linux is very fussy about them. I can't remember exactly what I did, but I think it involved global search-and-destroy missions against carriage returns in both Word and Unix, then global restoration of the important returns in Word, then a passage through the text editor Mi to convert any Mac-style returns to Unix ones. And presto, the problem seems to be solved!

So while I've been sleeping the program has been busy searching the gene sequence file for USS motifs, and later this morning I hope to be able to compare the forward- and reverse-direction motifs. We know that protein coding constraints do affect the reading frame that USSs are found in - for each USS orientation there's a preferred reading frame that USSs are best tolerated in. So it's reasonable to suspect that the USS consensuses might also differ between the orientations. If they do, we'll understand a bit more about how natural selection acts on USSs.

Later: No, I was overly optimistic. The program is able to find a short version of the USS motif (10bp) but it can't find anything when asked to search for the full-length motif (22bp). I suspect it needs to be given a stronger prior expectation than just the spacing I'm giving it. Maybe I'll try suggesting the consensus sequence.