No DNA Control

Phase II

2011-03-16T00:15:00.000-07:00

The Human Microbiome meeting was awesome...

Here's a 16S phylogeny of bacteria that live in the human airway:

here's some other places on the body

Anyways, cool stuff.

So my ears are itching, because yesterday the local genome centre was doing quality control on our libraries of Haemophilus influenzae recombinants (and controls)… If all goes well, we'll get scads of data in a couple weeks!

But I also need to get the work of our "Phase II" plan for Genome BC off the ground quickly, so we'll have samples to submit within a couple months of getting this incoming "Phase I" data.

The goal is to make maps of transformation frequency across the H.flu genome, when our standard lab strain Rd is transformed by DNA from clinical isolates with varying distributions of genetic variation. We hope to be able to investigate the sequence factors controlling the chance of homologous recombination, not only the effects of divergence, but also the role of chromosomal position, for example proximity to transformation hotspots.

There are two parts to making these transformation maps, the experiment itself and the sequence analysis. For the latter, we have worked out a strategy, but I'll just give a brief synopsis here: The challenge is basically to be able to make frequency estimates that are nearly as low as the sequencing error rate. We plan to handle this by using the very high co-transformation rate of adjacent SNPs and the high density of SNPs between H.flu isolates. Thus we will use information within reads and between paired reads to distinguish true transformation from spurious errors.

The basic experiment is really straightforward: (1) Incubate donor DNA with competent cells, allowing for uptake and recombination. (2) Purify the recipient chromosomes. (3) Sequence these to extremely high genomic coverage and measure the donor-specific allele frequency at each polymorphic site.

While that sounds simple, there will be a bunch of challenges ahead in getting this done properly (and on-time). We've got a pretty got a pretty good plan for handling the issues, so I think we'll be able to make appropriate samples for sequencing:

We need as highly competent of cultures as possible. Our sensitivity will depend on sequencing chromosomal DNA fragments from as few non-competent cells as possible.
We need very high complexity libraries, since we will be sequencing to several thousand fold genomic coverage of each library.
We need replication.
We need more than one clinical isolate to donate its DNA (so that we can actually make this into a real experiment).
We need controls.
We need to show that we don't have to worry about donor DNA contamination (i.e. donor sequences that weren't recombined into the chromosome).

So the first thing to do is make the competent-est competent cultures I can. Our sensitivity will depend on sequencing DNA fragments from as few non-competent cells as possible.

This week, I'm going to test MIV transformation of wild-type Rd and a couple of our hypercompetent mutants (murE and sxy-1) with/without cyclic AMP to see which gives the highest transformation frequencies (and least congression). (I'll start just with our regular old MAP7 DNA). We reason that this will not alter the transformation frequency map in an important way, because the hypercompetent mutants and addition of cAMP affect the regulation of competence, rather than its consequences.

I am tempted to try transforming H.flu circa 1961 and make the crazy blood brain heart media and do an anoxic transformation. Their transformation frequencies were so high in the sixties! It might be worth doing, though the difference could perhaps be something other than the culture conditions (perhaps our strain has changed?)… It would also require, ahem, some fairly odious preparation...

Anyways, I'll try and keep the blog updated on this experiment's progress (and later the sequence analysis) as things move along…

Uptake specificity as a tool for microbiome research

2011-03-01T19:20:00.000-08:00

image credit

Sorry for neglecting you, blog.

So Rosie and I are headed to the Human Microbiome Meeting here in Vancouver in just over a week. This means there's two posters to make! AAccK! We'll be a little out-of-place at the meeting, since our work is mostly in the lab, but we're hoping to learn more about the normal environment of H. influenzae and maybe even find people with access to mucus interested in bacterial genetic exchange.

So our posters won't be directly about the microbiome, but instead one will be about our genomics of natural transformation experiment, and the other will be about the uptake specificity measurements using deep sequencing. As I was thinking about the future of the uptake specificity work, it occurred to me that we could use H. influenzae as a way to filter Haemophilus and other Pasterellaceaen chromosomal fragments from a lung mucus DNA sample, for example. This would at least be a good way to screen out the gobs and gobs of human DNA likely there (at least most of it) and maybe focus a microbiome project down to a narrowish group... save a lot of money that way...

Sequence submission, ho!

2010-12-11T12:49:00.000-08:00

link

Hi blog! Long time, no update…

On Thursday, I submitted 91 genomic DNA samples to our local genome centre (center for you American reader). They’ll be sheared up, converted into indexed (barcoded) libraries, and sequenced on an Illumina GA2 over 4 “lanes”. Whew! That should have been easier.

I won’t get the data until probably around February, so with that out-of-the-way for now, I’d best get my act together. But before that, here’s what this first round of sequencing for our Genome BC grant is about:

We want to get a basic genetic picture of what happens when we do our standard transformation experiments; we lack simple rules of genetic transmission by natural transformation. How many DNA fragments do competent cells add to their chromosomes by homologous recombination? How long are the segments? How does genetic divergence affect transformation?

link

Haemophilus influenzae KW20 cells can be made naturally competent by starvation of log-phase cells in M-IV media. A fraction of the cells in the culture becomes competent, able to take up added DNA fragments into their cytosol (in H. influenzae, uptake is biased to fragments containing USS, abundant motifs in H. influenzae chromosomes). Taken up DNA fragments can recombine with the cells’ chromosomes, if sufficient sequence identity between the fragment and part of the chromsome allows it.

We added genomic DNA from an antibiotic-resistant derivative of a divergent clinical isolate NP (86-028NP NovR NalR) to competent KW20 cultures in three separate experiments. In each case, we isolated antibiotic-resistant and –sensitive clones, grew these up in overnight cultures, saved frozen stocks, and extracted genomic DNA from the rest of the cultures for sequencing.

In these experiment, we are not getting a particularly “natural” view of natural transformation, since the donor DNA we’re using is purified chromosomal DNA from a single divergent isolate. Since we know effectively nothing about the actual DNA in the natural environment of naturally competent H. influenzae—its source, size-distribution, or associated muck—this seems like a good place to start. We know the genome sequence of the donor and can use the genetic differences between donor and recipient to identify transforming DNA segments.

Based on our preliminary work, we think that these transformed clones will contain long(ish) segments from NP, replacing segments of the KW20 chromosome. We’re using multiplexed Illumina GA2 paired-end sequencing to identify recombinant segments in a large set of clones from these three experiments. We won’t get perfect genome sequences out the other end; instead we expect a median read depth of ~60 per clone, enough to probably capture the vast majority of diagnostic polymorphisms in each clone. Then we’ll look at the distribution of recombinant segments across the independent transformants. That’ll be the Science part.

Here’s what the 91 samples were:

Selected set: 72
Because we estimate that only about 10% of the cells in each culture was actually competent, we selected for either novobiocin or nalidixic acid resistant clones (3x12 of each) to ensure that the clone was derived from a competent cell. To allow for selection for transformed cells, the NP donor DNA was purified from a clone made NovR and NalR by PCR-mediated transformation. This means that for a given selected clone, a recombinant segment is always predicted that spans the appropriate locus. At the two selected loci, we will be able to ask about “LD decay” away from the selected positions. We also expect a large number Independent recombinant segments in these clones; these will provide basic information about the distribution, size, and breakpoints of recombination tracts in natural transformants.

Unselected set: 8 pools of 2
While we estimate percent-competence by evaluating the “congression” (unexpectedly high co-transformation) of “unlinked” markers (not on the same DNA fragment), this calculation relies on an assumption about competent cultures, namely that cells come in only two flavors: either non-competent, unable to take up DNA, or competent, able to take up several long fragments of DNA. It could also be that “congression” arises from a more quantitative distribution of states, where cells are more or less competent; i.e. able to take up and recombine a variable number of fragments.
For the most part, the use of “congression” to evaluate the overall competence of a culture is probably valid, but in terms of generalizing from the above large set of selected transformants, we’d like to know how well this assumption holds true. Because we expect that most unselected clones will look just like the recipient KW20 chromosome, we pooled 16 unselected clones into pairs and include these pools as 8 of our submitted samples. Since we predict that only 10% of cells in the cultures was competent, our null hypothesis is that 1-3 clones will look like selected clones (i.e. a handful of independent recombination tracts) and the rest will look just like KW20. We may, however, find that many more clones show evidence of transformation, but containing only one or only very short recombination tracts.

Selected late-log set: 8
Similarly, we also included 8 NovR clones selected from a late-log transformation. Cultures of cells in late-log are considerably less competent than MIV cultures (~1-2 orders of magnitude), and Rosie showed me some old data in her notebook suggesting that this could roughly be accounted for by percent-competence. The suggestion is that the low transformation frequency of late-log cultures vs MIV cultures is only because fewer cells are competent, rather than late-log competent cells being less transformable. Sequencing a handful of transformants from a late-log culture will help sort this out, in a similar fashion as sequencing unselected clones from MIV cultures.

Controls: 3
We included both the donor and recipient genomes as samples; while we’ve already obtained ridiculous coverage of these samples, including them here acts as an internal “coverage control”. We also threw in some MAP7 DNA, which is a KW20-derivative with a bunch of antibiotic-resistances. While it’s going to look almost like KW20, we use it on a practically daily basis, so it’d be nice to see its sequence.

So that’s the set of 91 clones we’re having sequenced with our first round of Genome BC cash. I’m keeping my fingers crossed that the data is good and abundant and reveals new things…

It's more delicious because it's more nutritious

2010-08-21T00:37:00.000-07:00

I think I got somewhere with my degenerate sequencing data... possibly even discovered something unexpected. The first thing I had to do was more or less start over my analysis, and temporarily forget everything I’ve read about “information”, “bits”, “logos”, etc. and just go back to thinking about the experiment from the get-go...

When it comes to DNA, H. influenzae are picky eaters. They seem to find some DNA fragments tastier than others. When I came to the lab, Rosie provided me with a DNA sequence: AAAGTGCGGTTAATTTTTACAGTATTTTTGG. She told me that DNA fragments containing this sequence (the “consensus USS”) are especially delicious to cells, and she wanted to know why.

Well, before I can aspire to helping with that, it'd be helpful to know which parts (or combinations of parts) of this special sequence actually make it delicious. Are all positions equally important? Do different positions make independent contributions to deliciousness, or do they work together to make an especially tasty morsel? Perhaps such insight into H.flu’s palate will get me closer to the bigger project of understanding why H.flu is so picky.

I did the following experiment last month to this end:

(1) I provided a massive smorgasbord of DNA to a whole bunch of cells... hundreds of culinary choices per cell to a billion cells. This tremendous buffet also had ridiculous variety. Every DNA snack was based on Rosie’s extra-tasty sequence, but there were differences in every bite (due to incorporating 24% degeneracy at each base during DNA synthesis). Indeed, there were hundreds of billions of different taste sensations in the input DNA pool.

(2) I harvested the undigested meals from the cells (facilitated by a mutation that kept the food from getting digested). The cells probably had terrible heartburn at that point, but each on average managed to stuff down dozens of DNA molecules.

(3) I sent two DNA samples to my friend to get sequenced: a bit of the original DNA buffet, and the DNA that I purified from the cells’ guts. By comparing the relative abundance of different DNA molecules, we can learn which positions are important to making the sequence tasty and, hopefully, how the different positions work together to impart deliciousness.

(4) Now I have two giant lists of >10 million sequences each from my friend. So far, I’m working with the first 100,000 of each dataset so I can figure out how to process the data. For the below analysis, this is more than sufficient to make accurate frequency measures.

First, I asked how the uptake sample (the H.flu guts) compares to the input sample, treating each position in the sequence independently. To keep things extra simple, I only considered whether a given position on a particular sequence read is matched or mismatched from the original sequence.

So, for each position, I counted how many of the 100,000 reads in each set were mismatched. To calculate the relative enrichment of the mismatched sequence in the uptake sample compared to the input sample, I simply took the ratio of uptake:input, then took the log of that ratio. This normalization makes it so that positive numbers indicate enrichment between uptake and input, while negative numbers indicate depletion. Here’s what I get:
Lovely. It looks like some parts of the original sequence are more important than other parts. For example, sequences that didn’t have a C at position 7 were depleted more than 64-fold in the uptake sample compared with the input sample. This must be an especially important flavor. On the other hand, some of the other positions look like they contribute less to tastiness, and a few positions don’t seem to contribute at all.

Great! We could now go and compare this to some other information we have, but for now, I’m going to ignore all of that and move on to a pairwise mismatch analysis. This will tell us how pairwise combinations of mismatches affect tastiness.

First, here’s that first barplot again, with the bars re-colored by the amount of depletion. Mismatches at red positions make the DNA taste bad, while mismatches at blue positions don’t.
That’s the legend for this next image, which shows the log2(uptake/input) for each pair of positions. To make the dataset, I counted how many of the 100,000 reads in each set had mismatches at a particular pair of positions. As with the above, I did not condition on the total number of mismatches in the read, just that at least those two were. On the diagonal is shown the consensus sequence and values corresponding to the above barplot. The plot is symmetric about this diagonal.
The big cross of red is like what we saw in the 1D barplot above... Those same important positions are still important. The cross tells us that, for example, not having a C at position 7 means you taste terrible, independent of where else you have a mismatch.

BUT, Hot damn! I think there’s a discovery here beyond that delectable flavor packet of GCG at positions 6-8 (and possibly a few flanking bases).

Note the 3x3 grid of regularly-spaced "mismatch interactions. Note the spacing between these 3 interacting parts of the sequence are just about a helical turn of DNA apart (10-12 base pairs).

I think this picture has some implications for our molecular models...

brain whirring...

(Notes: After getting the number-crunching steps nice and tidy, I spent an inordinate amount of time making these plots. I started using heatmap(), graduated to heatmap.2(), but ultimately took it to the next level and went fully customized with the image() function. It also took a ridiculously long time to figure out how to have fine-control over the color-scaling. This was necessary to be able to make the legend. I also made a more normal-looking legend (below), which was quite tricky to do right... If I wasn't so distracted by helical turns at the moment, I might explain how I made these in more detail.)

Direct uptake specificity measurements

2010-08-17T00:44:00.001-07:00

So I got a bit distracted from our plans with the Genome BC dough, since I got some more sequencing data! Woo! I'll return to the design of our 92 clones worth of sequencing later. For now, I am trying to work through some more grant writing, incorporating our findings from the new sequence data.

Below, I'll (try to) quickly summarize what I've found with my (very) preliminary data analysis and then try to describe the thing I'm having a hard time writing clearly about...

First off: It worked! My scheming panned out quite nicely. I got ~15 million single-end sequence reads from each of two DNA pools. One was from the input 24% degenerate USS-containing fragment. The other was from the periplasm-purified DNA fragments. The experiment and its outcome were described here and here.

Filtering the data:

The first thing I had to do was somehow parse the raw data file. For my first pass, I used a quick and dirty, but probably quite stringent, quality filter:

I excluded any read where the four first and six last bases were not exactly what they should be. (These sequenced positions were supposed to be non-degenerate, so should have zero mismatches from the expected input sequence.) This excludes reads with known "bad" bases and forces the alignment to be correct. The unfiltered data contained a bunch of reads that clearly had insertions or deletions within the USS, probably due to the limitations of the original oligo synthesis.
I excluded all reads that contained any Ns. This makes the analysis a bit simpler, especially when calculating the number of mismatches in the degenerate USS from consensus.

This left me with ~10 million reads in each dataset. That's a lot of loss, and I certainly will try and return to those sequence reads, but for now, leaving them out is a decent filter that will simplify the analysis (mainly by forcing a "true" alignment).

Details:

I piped a few greps together to pull all the sequence reads matching my criteria (the first four bases ATGC, the last six bases GGTCGA, and no Ns). From a UNIX terminal:

grep ^ATGC sequence.txt | grep GGTCGA$ | grep -v N > newfile.txt

That's it! (Someday, I should probably do something with all the quality scores, but for now, I'm ignoring them. This operation also lost all the read IDs.)

I soon found that when analyzing these new files that my poor computer-foo is making very slow-running scripts, so I took the first 100,000 reads from each filtered dataset onwards with a simple:

head -n 100000 newfile.txt

Making a Position Weight matrix:

I next calculated the frequency of each base at each position for the two datasets (generating two position weight matrices, or PWMs). Luckily, I have previously worked out how to use the input PWM to background-correct the periplasm-purified PWM here and here.

Details: On to some R...

I read in the data with scan. I've usually used thing like read.table or read.csv, but scan looks like the way to go with a big giant vector:

sequences <- scan("newfile.txt", what="character" )

Next, I calculated the per position base frequencies:

basesread <- 42
newmatrix <- matrix(, basesread, 4)
for (i in 1:basesread) { newmatrix[i, ] <-
table(substr(sequences, i, i))
}

Hmmm... Would've been nice to avoid that loop...

From the two matrices, I got a background-corrected periplasmic uptake specificity motif that looks like this:
This excludes the fixed bases and one other anomalous base at the beginning of the read.
(Sorry I can't use Weblogo or other things that make nicer logos, as none accomodate position-specific base composition.)

Hey, look at that C at position 7! It just sticks way out! That must be a particularly important residue for uptake. A previous postdoc and graduate student both had measured uptake efficiency for fragments with point mutations in a consensus USS. This data offers a great independent validation of this basic result. Here's what that data looks like:
It looks like the two approaches are giving comparable results. Fabulous!

The next step:

So what else can we learn from this much larger deep sequencing dataset? I'm sure all sorts of things. For now, the only other thing I've done is calculate the number of mismatches between each of the 100,000 filtered sequence reads and the consensus USS the degenerate construct was based on.

Here's what the distribution of mismatch number looks like for the two datasets:

Awesome! Shifted, just as previously simulated!

Details:

At first I tried a bunch of extremely laborious and slow ways of doing this. Then I found a package called "cba" with a function called "sdists", which can do edit distances with vectors of character strings. (The sdists step is how I murdered the computers when I tried using the full dataset.)

install.packages("cba")
library(cba) consensus <- "ATGCCAAAGTGCGGTTAATTTTTACAGTATTTTTGGGTTCGA" distances <- as.vector(sdists(consensus, sequences, weight=c(1,0,1) ) )

I guess I could just do it in chunks for the whole dataset to save on RAM...

Anyways, this vector of distances is going to serve as an important index to the sequence reads for the next set of analyses (which will look at co-variation between positions and how specific mismatches affect the motif).

Surprise:

Okay, now to try and describe the thing I'm having trouble saying in a succinct way. These information content measurements, as in the above motif, are in bits. This is "information content", but can also be thought of as "surprise". This is why the USS motif derived from the genome has so much higher information content than the new experimentally determined uptake motif. (Here's the genomic USS:)
Because most of the bases at a given position in the degenerate USS pool were already the preferred consensus base, little "surprise" was elicited when the preferred base was enriched in the periplasm-purified pool. As measure in bits. That's why the scale is so different between the two motifs.

Anyways, the issue arises when thinking about some other possible experiments. For example, if uptake specificity were partially modified (as we think we can accomplish using engineered strains), the positions with a new specificity would become extremely informative relative to the other when using the normal degenerate fragments. But if a new construct were designed based on the altered specificity, these positions would no longer have such elevated information content relative to the other bases...

This isn't actually a bad thing--in fact we can exploit it--it's just something hard to describe clearly...

How to spend some dough?

2010-08-11T15:14:00.000-07:00

We have some money to spend (kindly provided by Genome BC) for identifying recombination tracts in individual transformants, where donor genomic DNA from the clinical isolate NP is used to transform Rd competent cells (NP-into-Rd transformations). This will help us answer some very basic questions about natural transformation: What are the numbers and sizes of recombination tracts in individual transformants? Are different parts of the genome equally transformable? Do mismatch repair and other degradative mechanisms limit transformation?

We already have a few clones worth of data (which is now mostly analyzed), but we need more clone sequences to do any worthwhile statistics. Actually, I've made nicer pictures and done some re-analysis for the poster I took to the Evolution meeting. I should post some of that sometime...

We’d originally had a cost-saving plan of sequencing overlapping pools of clones to obtain data from 64 clones. We even came up with a slicker plan (since sequence yields keep increasing), in which we'd get 128 clones worth of data, and 64 of these would be "de-convoluted":
But now we’ve discovered that our local genome center provides inexpensive indexing of DNA samples (i.e. barcoding). This means we can obtain 92 clones worth of data for only ~$15K. Wow!
So now, I need to collect the clones for sequencing and extract their DNA. The question is what clones? I will be sequencing 4 sets of 23 clones (plus a control), so might think of the data collection in that manner. (This number of indexed clones per pool will be sufficient to obtain ~50-100 fold median coverage per clone, given current estimated yields.) Three are two things to consider: (a) technical concerns, and (b) using some clones to look at different kinds of transformations.

This post will cover the basic technical concerns, and I’ll use the next post to discuss some alternative things we might do with some of this sequencing capacity.

Replication:
To deal with potential criticism later, I should collect the NP-into-Rd transformants from three separate transformation experiments using three separate competent cell preps. This will allow me to make statements about how reproducible the distribution of recombination tracts between transformations is. Because each replicate will only have a couple dozen clones, this will likely not be sufficient for detecting subtle differences between the independent replicates, but will be sufficient in determining overall consistency of different transformations

Selection for transformants:
Since only a fraction of cells in a competent culture appear to be competent (as measured by co-transformation of unlinked marker loci), I will again select for either the NovR or NalR alleles of gyrB and gyA (respectively) that I’ve already added to the donor NP strain. This selection step ensures that every sequenced clone has at least one recombination tract, thereby selecting against clones derived from non-competent cells.

With indexing, it shouldn’t matter how I organize the clones submitted for sequencing, but this data can also be analyzed without paying attention to the index (for which I have my reasons), so it would make sense to organize a plate like this:

We’d originally decided 64 additional clones (atop the 4 we’ve got) would be sufficient for a basic picture. Now we can do 92 for much less money. Should we then do 92 NP->Rd selected transformants instead of 64, or should we use some of this new sequencing space to investigate new things?

Doing more clones would give us greater statistical power, but the more we do the more get will get diminishing returns. Since we estimate that ~60-70 clones will be sufficient for a nice first look, maybe we don’t gain much by doing 92 of the “same” thing and would be served better by using some of the extra space for other endeavours.

For example, if we reserved one set (23 clones) for one or more other DNA samples, we’d still get 69 clones worth of recombination tracts. What could we do with this extra space (keeping in mind that these should be transformations that can be done immediately)?

The most important factor to consider with each of these is whether the sample size (23 clones) would be sufficient to obtain useful data, despite having given up some statistical power for the primary set of 69 clones. (In my mind this means what we get should publishable as is alone or with the primary dataset, but could also mean several small pilot experiments for “preliminary data”.)

There are a whole lot of things we might do with this space… (stay tuned)

The anticipation is palpable...

2010-08-10T13:32:00.001-07:00

We’re resubmitting a grant soon, in part of which we propose to measure the specificity of DNA uptake in H. influenzae for the “USS motif”, a ~28 bp sequence motif that is sufficient to mobilize efficient DNA uptake. One of the things we want to know is just how the motif lends itself to efficient uptake by competent cells. In previous work, the lab has shown that point mutations at positions in the motif with very high consensus sometimes affect uptake efficiency, but other times do not.

Just before leaving for vacation, I managed to submit exciting samples to my friend for sequencing. These DNA samples are being sequenced AS WE SPEAK!!!!!!! RIGHT NOW!!!

The two samples were: (1) a complex mix of 200 bp fragments containing a degenerate USS with a 24% chance of a mismatch at each base of a consensus USS 32mer; and (2) an enriched fraction of these fragments taken up by competent cells into the periplasm. My previous post on my pilot experiment describes the experiment in more detail with additional links to other posts on the degenerate USS.

In brief, I am trying to compare a complex input pool of DNA fragments with what is actually taken up by cells to make a new uptake motif, based on the uptake process itself, rather than inferring it from genome sequence analysis.

This is the uptake saturation experiment I did with the consensus USS (USS-C) and the degenerate USS (USS-24D):
As I've seen before, USS-24D is not taken up as well as USS-C, since there are many suboptimal sequences that are inefficiently taken up. At very high DNA concentrations, similar amounts of USS-C and USS-24D are taken up. This is presumably because there are enough optimal fragments to saturate the uptake machinery at very high USS-24D concentrations.

I purified the periplasmic DNA from three of these USS-24D uptake samples: SUB (10 ng/ml), MID (70 ng/ml), and SAT (508 ng/ml). Because the yields were quite low, I used PCR to produce more material from these periplasmic DNA preps. It took a little effort to optimize this PCR. Perhaps I will get back to the weird artifacts I discovered at some later point, but suffice it to say, I ensured that I had a clean PCR and did few rounds of PCR, so that I wouldn't alter the complexity of the library too much.

I then took this periplasmic DNA (which should be enriched for optimal sequences) and re-incubated it with fresh competent cells. If the experiment worked (so the DNA samples are worth sequencing), then the periplasmic USS-24D prep should be taken up better than the original input. Indeed this was the case:
Note, this experiment did not take the saturation curve as far out as the original, due to limits on the number of samples I could process in a single experiment and limits on available material before leaving for vacation.

Unfortunately, I didn't see an appreciable difference between the three different periplasmic preps (SUB, MID, and SAT), which I'd hoped to. I had originally thought that periplasmic DNA recovered for unsaturated uptake experiments would include more suboptimal sequences, while saturated experiments would yield only optimal sequences. I hoped that this contrast would allow me to investigate competition between different sequences. Oh well.

So I decided to send only two samples for sequencing: the input and MID. These are the ones BEING SEQUENCED RIGHT NOW!!! WOO!!!

It's still possible the different purifications would behave differently at the higher end of the saturation curve (i.e. perhaps SAT would saturate sooner than SUB). So it's probably worth running some more saturation curves to see if its worth sequencing the other samples. It'd be really nice to have an actual experimental condition changing to get the most accurate depiction of uptake specificity.

SO, what do I need to do next?

(1) Prepare for the coming data deluge: I think I know how to start the data analysis once I've processed the raw read data into a comprehensible ~2e7 rows X 32 columns (though probably in a computationally slow fashion... meh). What I'm less prepared for is the initial processing of the data. I designed the experiment to sequence 4 non-degenerate bases upstream and 6 non-degenerate bases downstream of the USS, so I will probably do a simple crude quality filter, where I demand the first 4 and last 6 bases are exactly correct. This will tend to eliminate poor reads and force the alignment of the 32 degenerate bases in the middle to already be correct. This filter will exclude erroneous oligo synthesis or fragment construction that introduced indels into the USS, which will simplify things initially. At this point, I will also need to apply a base quality filter to ensure I ignore base calls that have low confidence. Sequencing error is a problem for this analysis, so its important I use stringent filters. Even if I only ended up with 10% of the raw reads, I'd have an enormous amount of data for sequence motif analysis.

(2) Prepare for abject failure. It's possible I screwed up the design of the constructs or that there is some unanticipated challenge sequencing through those 32 bases or using this approach. I'll know in a few days, but need to think about what the next step would be, if indeed I misplaced some bases in my reverse engineering scheme.

(3) Prepare for a seeming failure that isn't really. I may have screwed nothing up but get back data from Illumina's pipeline that says something's terribly wrong. I don't fully understand the details of this, but the base-calling step in Illumina's pipeline (which is reading the raw image files from each sequencing step) may be screwball, because of the extremely skewed base composition my constructs will have during each cycle. E.g. the first base for every single cluster should be "A", whereas 76% of clusters should have "A" at the fifth base. Apparently, this can create some weird artifacts in the data processing step, which I need to be prepared for (and not despair when I find out "it didn't work". Mostly this will involve working with my friend and his sequencing facility to re-run the base-calling with an alternate pipeline.

(4) Work on updating the grant for re-submission. There are several paragraphs of our grant application that need to be modified to show our progress. To some extent, this will involve the soon-to-come data, but I can begin by identifying the parts that will need changes an including some different figures.

(5) Work out the molecular biology. I've done a bunch of uptake experiments in a bunch of settings, and I should now have enough to put together some kind of little story, even without the sequence data. With a mind towards a paper, I need to work out just what I need to do. Do I need to do another sequencing experiment with a different construct (more or less degenerate) or under different conditions? Or is a single periplasmic enrichment enough? If the latter, i will certainly want to show a bunch of other experimental data, but which experiments? And which need to be repeated. My first step in this direction is to go through my notebook and figure out what I have...

(6) Work out a large-scale periplasmic prep. I circumvented this for the degenerate USS experiment by doing PCR. SInce I know my exact construct, I could use PCR to amplify it and didn't need to have a good yield, nor particularly pure of a prep (since genomic DNA won't amplify using the construct's primers). However, if I want to look at uptake across a chromosome, I will need to fractionate periplasmic DNA away from chromosomal DNA both to a high level of purity and with a high yield. I've accomplished each of these individually, but so far have not managed a large-scale periplasmic prep that leaves me with enough pure DNA to make sequencing libraries reliably. I refuse to use whole-genome amplification for this experiment.

Hi blog!

2010-08-10T13:18:00.000-07:00

So I’ve returned to work after a 3.5 week vacation to Bali! Woo! Time to catch up on some long non-blogging. But before that...

Here’s where my crazy Uncle Jimmy took me and the missus for a few days… Gili Gede (meaning “Little Island Big”) off the coast of Lombok (which appears to be experiencing a gold-rush, so go now before its all built up!).
The ferry between Bali and Lombok took us across the legendary "Wallace Line", though my total ignorance of systematics means I couldn't really tell the difference in flora and fauna between the islands. (Though we did see a HUGE monitor lizard near a brick-making operation.)

I stayed out at the end of this pier on the lower left at the “Secret Island Resort”. (Note the vessel that took us snorkeling to the immediate South... the "Scorpio".
And here's an actual shot of the infamous "Rocky Docky" with the lovely Heather in the foreground:
The reef was right over the edge of the pier, so snorkeling to see all the corals and fishies was optional. You could just look over the edge! Excellent!

Anyways, I've been getting my head back into things and am writing a few blog posts about what’s going on in my world-of-science. So expect a deluge of catch-up posts during the next day or two…

Manuscript plans?

2010-06-22T09:17:00.000-07:00

So I’m a bit over one year into my postdoc. What have I got to show for myself? Well, plenty of work, but not any papers, or even written manuscripts, so that’s a bit of a problem. Can I turn my first set of genome sequencing data into a manuscript?

Seems likely. I collected >5 Gigabases of Illumina sequence data from several Haemophilus influenzae chromosomes, and this could be used as the basis of a manuscript. I obtained data from a donor strain (86-028NP NovR NalR) and a recipient strain (Rd, RR722) as control data (in order to evaluate the ability of the sequencing and read alignment to correctly identify polymorphisms). I also obtained data from two individual transformants and a pool of four transformants to identify donor alleles in transformed recipient chromosomes. I even found some things out.

Does this a paper make? One outstanding issue is that, in spite of being a lot of data, which has required a fair amount of work to get a handle on, there is not a tremendous amount of biologically relevant data. Yes, I obtained extremely accurate and comprehensive data for the four transformants sequenced. But it was still only four transformants. There are some biologically meaningful results; they just aren’t terribly novel or statistically robust. The bigger biologically meaningful results will have to wait until we can collect more data.

So to turn this into something publishable, the approach and method need to be important enough (and made explicit enough) to be of value to others. So far, I have not done anything in my analysis that is truly novel, but I have managed to produce the bare-bones of a “pipeline” for measuring allele frequencies from pools, and identifying recombination tracts in transformants. The data we got was also extremely high coverage, so we were able to see the limits of the technology fairly well: i.e. depth-of-coverage variation, errors, and issues with read alignment.

Though everything I’ve done so far uses “off-the-shelf” bioinformatics tools, there are so many people trying to do similar things, it might be useful to write a paper that is sort of an “application” of the technology and tools I’ve been using. It took me months to piece everything together, so maybe I could save someone else some time by having everything in one place. But with each passing day, the value of such a paper is probably diminishing, so I’d best get started!

There are still a few analyses I’d like to do that would give the paper a little more spice:

Structural variant analysis: This is something that will involve our collaborators at UVA, who are experts. We can see these pretty well (at least the larger ones), but something systematic has yet to be done.
Reciprocal read mapping: I’ve mapped all the data to both the donor and recipient genomes, but I have not really fully leveraged this fact. The read alignment artifacts that arise mapping data from one strain onto the other could be handled much better, if I was able to assign individual reads to either of the two reference genomes, based on the mapping quality. I’d really like to do this. It’d be novel, mainly because most people sequencing are doing SNP discovery. I already have all my SNPs discovered, so doing an extra good job at calling SNP frequencies using reciprocal alignment would be at least something new. This will take a bit of work, however, and I’ll need to figure out the best computational way to do it. Aside from doing uber-detailed error analysis for a technical paper, I think this is really the best chance to make a novel contribution bioinformatically.

Here is a rough outline of the manuscript, as I’m viewing it so far:

1. Introduction

Many bacteria become naturally competent. Natural transformation is important in evolution.
Previous sequencing studies have focused on only a handful of defined constructs. (Bacillus, Helicobacter, Actinobacillus).
For any organism, the total extent of recombined fragments in individual transformants has never been directly evaluated, and the factors dictating the chance of transformation are only poorly understood as a result.
Haemophilus influenzae is a model system for natural transformation. The mechanism is well-defined. Transformation is efficient in the lab strain.
The extensive natural genetic variation between H. influenzae strains provides tens of thousands of markers to identify recombination tracts in individual transformants. Not only single-nucleotide, but structural variation in the form of indels and other rearrangments. The “supragenome” hypothesis.
We investigated the use of massively parallel sequencing (or “next generation sequencing”, NGS) to characterize natural transformation at a whole-genome scale.
Our results show the Illumina platform to be an excellent method to obtain nearly exhaustive information on recombination tracts in individual transformants. Our approach uses the alignment of sequence reads to both donor and recipient reference sequences. We obtained donor and recipient genome sequence as controls for evaluating sequencing error, depth of coverage, and polymorphism identification. We also obtained the sequence data from two individual transformants and a pool of four transformants.
Individual recombination tracts are longer than previously appreciated and can bring hundreds of polymorphisms from donor to recipient chromosomes (both single-nucleotide, insertion, deletion, and insertional deletion). However, recombination tracts often appear interrupted by or terminated at sites of structural variation between the two genomes. This shows that such variation are barriers to strand exchange and/or are preferred mismatch repair substrates.

2. Materials and Methods

Strains
DNA
Transformations
Library preparation
Illumina sequencing and initial data processing pipeline
Reference genome alignment by MUMmer, MAUVE
Reciprocal read alignment by BWA
SAMtools pileup
Galaxy pileup parser
Variant frequency analysis
Assignment of reads (unimplemented)
Donor segment calling
Analysis of structural variation by HYDRA

3. Results

Genetic transformation of competent cells: Marker to marker variation. Dependence on sequence identity. Congression and linkage
Illumina sequencing and read alignment: Table of sequencing results and fraction of mapped reads. Variation in depth-of-coverage. Sources of sequencing error and read mapping artifacts. Reciprocal read alignment? Varying alignment stringency?
Comparison of donor and recipient strains. Identification of SNPs and structural variants between the donor and recipient strains. Comparison to whole-genome alignment methods.
Identification of donor alleles in transformed recipient chromosomes. Accounting for SV alleles. Identifying novel alleles.
Identification of allele frequencies in a pool of four transformants.
Identification of donor segments and putative recombination tracts
Enrichment of SVs at donor segment breakpoints

4. Discussion

A first look at transformation… still few transformants
Excellent method. Limitations are circumvented by very high coverage, knowledge of both donor and recipient genome sequences, and the use of reciprocal read alignment (unimplemented)
Big recombination tracts. Evidence of mismatch repair. SVs as blocks to recombination tract progression.
Speculations: Hotspots? Role of uptake specificity? Supragenome transfer?
Future: Aside from collecting more transformants, making a transformation frequency map to investigate the “cis-acting” factors controlling the efficiency of transformation. Long-term utility in understanding the population genetics of human pathogens.

Still a rough outline, but something to start with....

Degenerate Uptake: Pilot Study

2010-06-11T03:48:00.000-07:00

Something to blog about! (Wow; it's been over a month... sorry to my three loyal blog readers.)

I’ve gotten around to doing a pilot-scale experiment on the specificity of H. influenzae DNA uptake for the “uptake signal sequence” (USS). The USS is a ~29 base pair motif highly abundant in the H. influenzae genome, and sites that match the consensus USS are known to be preferred substrates for DNA uptake by competent cells. The presence of many USS in the chromosome is presumed to be why H. influenzae competent cells prefer H. influenzae DNA over DNA from other organisms.

However, little is known about how the structure of USS contributes to uptake of USS-containing fragments: Limited analyses of mutations of a DNA fragment containing a consensus USS suggests that some but not all informative positions in the USS motif are important to uptake, indicating that other forces (perhaps later steps in transformation) contribute to the structure of the USS motif.

To carefully dissect uptake specificity for the USS motif, we have devised an enrichment experiment:
(1) A complex pool of DNA fragments containing a degenerate USS library is incubated with competent cells.
(2) The fragments preferentially taken up by cells are purified from the periplasm.
(3) DNA sequencing is used to compare the input and periplasm-purified pools of sequences.

Details and Pilot-scale Results:

I’ve previously discussed the design of the input DNA pools. The control 200 bp construct is designed to already contain the sequences needed for Illumina single-end sequencing, along with a 32 bp consensus USS site near the middle of the fragment. The test construct is the same, except the USS is degenerate, having a 24% chance of a non-consensus base at each position. Thus in the degenerate-USS pool, the average site has ~7-8 mismatches from the consensus sequence.

The expectation is that, while the consensus-USS construct (USS-C) will be taken up by cells well, the degenerate-USS construct (USS-D) will be taken up more poorly, since it contains many suboptimal sequences (i.e. it is less uniformly delicious). Indeed this is the case, with USS-C being taken up about 10 times better than USS-D at sub-saturating DNA concentrations (see below). The notion is that comparing the USS-D input to that taken up by cells will provide a precise measurement of uptake specificity for the USS (i.e. which sequences are tastiest). We think this will tell us a lot about the mechanism of uptake.

It occurred to me a couple weeks ago that before moving on to the data collection (i.e. the DNA sequencing), I should first make sure that the USS-D fragments recovered from the periplasmic purification are taken up better than the original USS-D input (i.e. the competent cells selected more delicious sequences). This would provide the clearest indication that the experiment worked and the material is worth sequencing. It is!

I compared the uptake of USS-C and USS-D before and after periplasmic purification of taken up DNA from rec-2 competent cells across a range of DNA concentrations. Here are the results:
A and B show the % DNA uptake for USS-C and USS-D, respectively, for different amounts of added DNA (to 200 ml competent cultures). C and D show the same data: C is a dose-response curve, and D is a double-reciprocal plot (since I used 2 ng of hot label, along with an additional amount of cold label for these experiments).

Input USS-C and periplasm-purified USS-C were quite similar, while periplasm-purified USS-D was taken up substantially better than input USS-D.

Notably, at low (sub-saturating) concentrations of DNA, periplasm-purified USS-D is taken up less well than USS-C, while at high (saturating) concentrations similar amounts of DNA are taken up. Also of note is that the input USS-D does not saturate until higher concentrations than the other three samples.

This is all good news. I left out a fair number of details, but this pilot-scale experiments is extremely encouraging. Next week, I plan to repeat the experiment, but this time on an appropriate scale for recovering samples for sequencing. I will also investigate how periplasm-purified USS-D samples behave when recovered from uptake experiments with varying amounts of DNA. I expect that at sub-saturating concentration, the cells will be less “picky”, such that periplasm-purified USS-D will be taken up less well than that purified from saturating concentration. This would provide a useful experimental condition, as in the sequence analysis we would be able to investigate the role of competition in shaping USS specificity.

I think this might end up working swimmingly... Onward!

Mining some old array data

2010-04-13T12:35:00.000-07:00

So in an effort to re-examine some of the lab’s old array data, I made a fairly simple R script to plot the change in expression of competence genes, putative purR-regulated genes, and genes involved in utilizing secondary sugars. We no longer have our expensive license for fancy-pants software, but all I needed to do was some arithematic to columns, and then find the rows of interest, so it’s R-tastic!

I looked at a time course dataset, in which expression was monitored over the course of growth in sBHI and after transfer to MIV. I also looked at a single one-off array comparing purR- to purR+ strains growing in late-log +cAMP.

Here’s the results for the time course. All values are normalized to the first time point. Blue are sBHI timepoints, and red are MIV timepoints. MIV cultures were split from the sBHI cultures at t=0 minutes.

It’s pretty clear that the competence genes are strongly induced in MIV, but are also induced in late-log phase, as expected. Putative PurR-regulated genes are strongly and quickly induced in MIV, indicating that purine pools are quickly depleted, and the purine biosynthetic pathway is activated quite quickly (much faster than the competence genes, it appears). The “non-PTS” genes (several genes induced by CRP when cAMP levels are high) appear to be briefly weakly induced in MIV, as well as being weakly induced in late-log.

Here’s the same sets of genes plotted as the ratio of expression in purR- vs purR+ cultures (late-log, induced with cAMP). Here, I plot the ratios from both array elements for each gene (open and closed circles) and colored them just so they’d be easy to see. Also note, I normalized everything to the median ratio to account for dye effects (under the assumption that the median gene is not PurR regulated). Again, strong induction of the putative purine-regulated genes, a weak repression of the competence genes (presumably due to purine repression), and not much happening with the non-PTS sugars.

Conclusion: Nothing we didn’t already suspect, but it’s good to see that things are behaved as expected. One point of note is that the hypothesized regulation of rec2 by PurR isn’t something that jumps out of this, but if purine repression acts upstream of rec2, we wouldn’t be able to see the effects of deleting PurR here anyways…

SNP densities

2010-04-02T13:25:00.001-07:00

So I’ve been writing yet another grant, which has been distracting me from blogging (this isn't supposed to be a monthly blog, but this will hopefully be the last grant application for a while).

But I’ve also been doing several analyses lately. Here’s one. I took the sequences of an ~300 kb restriction fragment from three H. influenzae isolates (Rd, 86-028NP, and PittGG). They’re all similarly divergent from each other (~2.5%), and I wondered how well the level of divergence of Rd vs NP and Rd vs GG correlated along the chromosome...

So I aligned the sequences in Mauve, took its SNP calling output, and did a couple simple sliding window analyses inside R (using the zoo package for rolling means). Here’s what divergence looked like averaged over 5 kb windows (click to enlarge):
The divergence between Rd and the two other isolates are quite well correlated (r2= 0.8, using linear modeling). But since NP and GG are similarly divergent, I made two other plots.

First, here’s a comparison of the density of SNPs that are shared by NP and GG and those that are unique to either NP or GG:
The correlation is a lot worse (r2=0.4).

And if I further break the “unshared” line into NP and GG-specific SNPs (i.e. positions are different between Rd and NP but not GG, and vice versa).
The correlation is worse still (r2=0.2)

Similar results applied to smaller windows, but the plots looked a lot messier. Note that it’s not exactly totally straightforward to measure SNP density... What does one do at indels?? I just ignored them, so the results above are rough. Part of the reason I focused on only a co-linear segment of chromosome was to minimize this problem, but there are still several indels between each of the three strains.

Indels aside, what’s this mean? One of the goals of my transformation frequency mapping is to be able to distinguish the effects of sequence divergence on transformation from the effects of other local chromosomal properties (base composition, sequence motifs, etc.). Since NP and GG have correlated SNP densities relative to Rd, transformation frequencies across the Rd chromosome are expected to also be correlated. Discrepencies in transformation frequency by NP and GG donors could indicate that SNPs specific to the isolates are somehow modulating transformation independent of divergence per se.

Distinguishing chromosome “position effects” from sequence divergence will probably require a third donor DNA. Deciding what this would be requires some thought. All of the sequence H. influenzae are similarly divergent from Rd (and for the most part each other), and phylogeny poorly distinguishes separate clades (i.e. they kind of give a star phylogeny).

So I should use either a strain much more closely related to Rd or one more distantly related (perhaps another species). Using a closely related strain has the advantage that transformation frequencies are expected to be higher and divergence will play less of a role, making the focus more on divergence-independent factors, but I would also have far fewer markers.

Based on MLST comparisons, several strains are sisters of Rd (RM7033, RM7429, RM7271). These assignments are made in several phylogenetic and put the three at ~0.5% divergent from Rd. So I would expect that RM7033 (for example) would have ~6000 SNPs from Rd (far more than our Rd or the other sequenced Rd), ample to have markers across the chromosome...

Chromosome "Position Effect"

2010-03-16T10:25:00.000-07:00

I got some data for the transformation frequency at five different markers. They varied. This isn’t anything ground-breaking; reports of a "position effect" for transformation go back decades and a couple of recent studies in other organisms bear it out. The underlying cause of variation in transformation rate at different positions likely stem from two sources: the physical structure of the chromosome and the sequence composition of the recombination substrates. The former case is reasonably well-worked out for analogous processes in eukaryotes: For example, in yeast, heterochromatic regions are recalcitrant to recombination, but these sites become recombinogenic in mutants with defective heterochromatin assembly. Sequence composition has also been shown to affect the efficiency of recombinational strand exchange in several different contexts, both genetic and biochemical. This latter type of variation is not traditionally considered a "position effect", but is difficult to distinguish from the former.

Anyways, I wanted preliminary data showing that I can, in fact, detect differences in transformation at different genomic positions, since a big part of my proposed work will involve measuring to very high resolution this position effect...

I used MAP7 donor DNA to transform three independent Rd competent cell preps. MAP7 is highly similar to Rd, except that it carries several point mutations that confer antibiotic resistance. There are likely other unselected differences between Rd and MAP7, but few. Thus, differences between marker transformation rates are likely to predominantly reflect chromosome position effects, rather than sequence divergence between donor and recipient.

This latter point isn’t strictly true: in order to see transformation, a genetic change has to be made, and the selected MAP7 point mutations are genetic differences. But because in our preliminary sequencing data, we saw long stretches of donor-specific DNA with dozens to hundreds of SNPs, I don’t think these single-nucleotide differences are contributing too hugely to the observed variation in transformation rate.

Here’s the data for the five markers individually. Vertical bars indicate the mean transformation frequency per viable cell to the indicated antibiotic resistance allele. The inset circle shows a rough map of the location of the MAP7 markers. (Sorry about the lack of an origin.)
Indeed, I see a ~5-fold range of transformation frequencies, from ~1/500 to ~1/100. Since this is only an arbitrary sampling of five sites, the range of variation across the chromosome could be much higher.

As previously discussed, these values underestimate the transformation frequency per competent cell. Competent cultures typically have both competent and non-competent cells, and the “fraction competence” is typically measured by looking at co-transformation frequencies. These are often higher than expected, even for unlinked markers, a phenomenon termed “congression” and interpreted as a binary distinction between competent and non-competent cells in the culture.

The technical value of this is that I can elevate the observed transformation frequency at one locus by selecting for transformation at another unlinked locus, since this eliminates all non-competent cells from the culture, providing potentially higher sensitivity on our proposed sequencing experiments. It also dampens differences in culture-to-culture variation caused by big differences in fraction competence (not shown).

However, there is at least one old report using Bacillus that suggests congression does not simply reflect a binary distinction between competent and non-competent cells . If cells only came in those two flavors, we would predict that any pair of unlinked markers would show the same level of congression, yet they report that different pairs of markers had different congression frequencies (aside from due to linkage). They go on to suggest an interesting model for their observations, but my concern is more technical:

Does selecting for transformation at different loci affect the tranformation rate at a second unlinked locus?

If the answer is yes, then selection for transformants at a locus would be a poor way to elevate the transformation rate at other unlinked loci, since it would be biased in an unknown way. I also measured co-transformation of Nal resistance and each of the other four. Nal is “unlinked” from all the others (i.e. DNA fragments from standard DNA preps will always be too short to contain the NalR allele with another antibiotic resistance allele), so I can measure “congression” four times.

Here is the data. So, for example, the first bar was calculated as: f(kanR nalR) / f(kanR). This normalizes each bar to the nalR rate (i.e. “the frequency of nalR among kanR transformants”).
The first thing to note is that the scale bar has changed relative to the transformation/cfu. For each of the 3 cultures, there was ~3.5 fold increase in the observed transformation rate, which would be expected if ~1/3 of cells in each culture were competent.

The second thing to note is that selecting for any of the four markers had no effect on the NalR transformation frequency. So the answer to the above question is no. Phew! The Bacillus result was cool, but I’m glad it isn’t the case here. A binary competent/non-competent model is perfectly reasonable in our system (though this does not exclude the possibility of variation among competent cells). With this in hand, I can now plot the co-transformation data with respect to NalR. If selecting for NalR only eliminates non-competent cells but does not change the underlying transformation frequencies per competent cell at the other unlinked markers, then life is good.

Here’s the data. So for example, the first bar was calculated as f(kanR nalR) / f(nalR). This normalizes each bar to its own rate (i.e. “the frequency of kanR among nalR transformants”). For the nalR/competent cells, I used the average of all 12 points in the previous plot.
This data closely resembles that of the first figure, except all the values are ~3.5 -fold higher.

Woo! Next I should probably repeat congression data for linked markers, and repeat experiments with more divergent donor DNA.

Repression of competence induction by purines

2010-03-14T13:01:00.000-07:00

My illustrious colleagues have been re-examining some of the lab’s old data regarding the repression of competence by purines. The work has been slowly ongoing for years, and it may be close to being a complete story. I want to try and express what I think their model is and what seem to be its predictions, so they can tell me whether my understanding is straight…

First, a schematic depiction of how I interpret what we already know the induction of the competence regulon:
What we knew: Transferring cells growing in rich medium to competence medium induces 15 operons driven from the novel CRP-S promoter, and then cells become naturally transformable. Cells in competence medium have elevated cyclic AMP levels, directing the CRP protein to induce expression of genes with canonical CRP-N promoters, including the sxy gene. Sxy protein alters the binding specificity of CRP to also bind at CRP-S promoters, thereby inducing competence gene expression.

But Sxy levels are also regulated at translation, in addition to at transcription. The wild-type sxy mRNA transcript contains a stem-loop structure that inhibits its translation. Mutations that disrupt the stem-loop structure in the 5’-UTR are hypercompetent (e.g. the sxy-1 mutation). In wild-type cells, unknown factor(s) disrupt the stem-loop to induce the translation of sxy transcript in competence medium.

Now a schematic depiction of how I interpret the model for purine repression of competence:
The observation: Addition of purines to competence medium represses competence. Purine biosynthesis is repressed by the PurR protein when cellular pools of purines are high. Deletion of the purR gene reduces competence (presumably indirectly, by increasing cellular pools of purine), but mutations disrupting the sxy 5’UTR’s stem-loop suppress the purR mutant defect (Rosie’s last post).

A hypothesis: The sxy transcript stem-loop is stabilized in the presence of purines (either directly or indirectly), blocking the production of Sxy protein and thus the activation of the competence regulon. When purine pools are depleted, the stem-loop is disrupted. This predicts that addition of purines and purR mutations will inhibit sxy translation more than sxy transcription.

A corollary hypothesis: Purines block DNA translocation by PurR-dependent repression of the rec-2 gene, whose promoter contains a putative PurR binding site. A potential test of this hypothesis would be to treat sxy-1 competent cultures with purines. We would predict that if PurR directly represses rec-2, DNA translocation would be inhibited (but DNA uptake would not). Obviously, checking rec-2 transcription relative to other competence genes would make sense here as well, but the functional test would be most compelling.

Is that the basic notion? I know there’s a bunch of other experiments that have been done that I need to find out about…

What have I got?

2010-03-11T09:30:00.000-08:00

Okay, grant planning part 2... Below is a dense description of the preliminary data I have/will have for writing this next grant...

PRELIMINARY DATA

Transformation frequency depends on chromosome position:
DNA from a multiply marked derivative of Rd (MAP7) was briefly incubated with competent Rd cultures. The resulting transformation frequency at each of four loci was evaluated by selecting for cells that acquired the corresponding MAP7-specific antibiotic resistance allele. MAP7 DNA transformed each Rd locus at a different frequency. (Repeat experiment in progress… stay tuned but looks good.)

Sequence divergence decreases transformation frequency:
DNA from an antibiotic-resistant derivative of NP (1350NN) transformed Rd competent cells less efficiently than did DNA from MAP7, and vice versa. NP differs from Rd by ~2.4% per alignable base position (and an additional 10% of each genome is absent from the other, contained in indel polymorphisms) while the transformation frequencies at two loci were affected ~2 to 4-fold. (Data in hand.)

Co-transformation frequencies are non-random due to congression and linkage:
(Wish I didn’t have to describe this, but it’s too fundamental. Data mostly in hand.)

Transformants acquire hundreds of donor-specific alleles:
Several large DNA fragments recombined into the chromosomes of four individual Rd competent cells, as revealed by genome sequencing. Each of the four transformants was selected for resistance to one of two antibiotics encoded in the 1350NN strain (two NalR and two NovR), and the corresponding donor-specific allele was present in each of the four. In all, 24 donor segments (contiguous stretches of donor-specific alleles) were found across the 4 transformants, with an average of 1.4% of each recipient chromosome replaced with donor DNA (~25 kb and ~600 SNPs each). Mismatch repair is likely responsible for the disruption of contiguous stretches of donor-specific alleles in the transformants; assuming that for closely adjoined segments this was true, a total of 10 (instead of 24) independent transformation events occurred across the four transformants (6 of which were unselected; notably two of these were overlapping in independent transformants). (Data in hand; re-analysis in progress.)

DNA uptake signal sequences (USS) are densely distributed in the two genomes:
Both the Rd and NP chromosomes contain USSs nearly every kilobase and most are syntenic. (Cursory data only. Need a better analysis.)

Sequence preferences in DNA uptake can be captured by periplasmic DNA purification:
DNA fragments containing uptake signal sequences are efficiently taken up into cells, and taken up fragments can be cleanly purified away from both free DNA and chromosomal DNA. The use of rec-2 and rec-1 mutations will facilitate separating sequence biases at different stages of natural transformation. (Data in hand, except rec-1.)

LIST OF FIGURES:

Four/five marker transformation rates
Rd vs NP transformation rates
SNP spacing histogram with embedded SV table
Genome sequencing figure (pool data)
USS analysis
Molecular biology figure (uptake data)

Another day, another attempt to get a dollar

2010-03-09T09:28:00.000-08:00

Sigh... another grant due soon; this time, it's my last attempt to get an NIH postdoctoral fellowship. My last reviews mainly took issue with my proposal, which they found to be overly ambitious and somewhat unfocused. So below, is my first attempt at a summary/specific aims page, followed by a couple of preliminary data collection things I'd like to do before it's due (on April 8)...

The introduction: Naturally competent bacteria take up intact DNA from their surroundings and can incorporate it into their chromosomes by homologous recombination. Akin to sexual recombination in eukaryotes, this natural transformation pathway moves alleles and genes between otherwise clonal lineages; and human bacterial pathogens have used this pathway to share antibiotic resistance genes, antigenic determinants, and virulence factors. To better elucidate the mechanism of transformation and to inform population/epidemiological studies, the proposed work will use the opportunistic Gram-negative bacterium Haemophilus influenzae to disentangle the sequence biases intrinsic to the DNA uptake and DNA recombination mechanisms by combining classical microbiology with modern DNA sequencing.

The specific aims:

Define the genetic consequences of natural competence to H. influenzae. Transformation frequencies vary for different sequences and at different chromosomal locations, and this could strongly influence the rate of sequence evolution and adaptation along the genome. I will transform competent cultures of the standard lab strain with the genomic DNA of a clinical isolate and use deep sequencing to measure transformation across the lab strain’s chromosome for all the ~40,000 sites differing in the clinical isolate. This will provide an unparalleled dataset for investigating the sequence factors that promote and limit genetic exchange between bacterial cells.
Measure the contribution of DNA uptake specificity to natural transformation. In several human pathogens, including H. influenzae, the uptake machinery prefers DNA fragments containing short “uptake sequences”, and abundant sequence motifs in many bacterial chromosomes suggest that biased DNA uptake has had a profound influence on genome evolution. I will purify the intact DNA molecules taken up into the periplasm and cytosol of competent cultures and use deep sequencing to measure the sequence biases of the uptake machinery. In combination with (a), this will disentangle the contributions of DNA uptake from those of DNA recombination during natural transformation.

The platitudes: The proposed work will link molecular studies of transformation to the growing genome sequence data being collected from many isolates of many bacterial species. By establishing my approach with completely sequenced chromosomes and using a well-defined experimental system, later studies could include a greater diversity of sequences or mimic more and more natural conditions. As a directly applicable outcome, the work will also produce the beginnings of a new type of genetic resource for mapping traits that differ between natural bacterial isolates (as in eukaryotic quantitative genetics) by generating fully genotyped recombinants. In the future, such studies will give empirical underpinnings to population genomic studies of bacterial genetic exchange, as well as provide new testable hypotheses for investigating the molecular mechanism of transformation.

Preliminary data I would like: (besides what I’ve got)

Properly replicated transformation frequencies for several markers. (I did this before, but it hasn’t been properly replicated.)
Follow a few molecules through uptake and recombination?
Population genetic inferences of “recombination” in H. influenzae? (I did this before, but it sucked.)

Multiplexing sans barcodes

2010-03-05T16:01:00.001-08:00

Previously, I’d said I wanted to go over how we might obtain many recombinant genotypes by deep sequencing pools of recombinants, since our tiny genome is TOO EASILY SEQUENCED using modern methods, making the sequencing of individual clones inefficient. The challenge is then in assigning donor DNA segments in the pools to individual clones. To a first approximation, this isn’t really necessary, since one of our main motivations for sequencing recombinants is simply to determine whether the locations and endpoints of donor DNA segments are biased: i.e. whether there are recombination hotspots or whether certain types of donor-recipient differences are recalcitrant to recombination.

However, our preliminary data showed that donor segments were often clustered in individual recombinants, probably due to mismatch repair disrupting larger donor fragments during transformation. We were only able to pin this down, because we individually sequenced 2 of the 4 transformants that we’d pooled. To illustrate, here’s a zoom of the region containing one of our selected sites at gyrB. 2 of 4 clones carry the causal allele (the red dot). But the pool data indicates several additional segments:

Are they in different clones? The same clone? How do we disentangle, without sequencing individuals (as was done here; shown as sets of colored bars at the top)?
Several methods for handling pooled data exist. The one typically referred to is “barcoding” where samples are processed individually and have unique sequence codes added during library construction, so that individual sequence reads can be assigned to individual clones. This is powerful method, but extremely expensive and labor-intensive. It surely has useful contexts, but for our purposes, we don’t really need to assign every read to every clone… only donor segments.

An alternate approach, outlined below, would simply ensure that any given clone appears in two different otherwise non-overlapping pools. In its simplest form this would simply be to pool by rows and also by columns (other more involved ways are here and here). I recently did a transformation experiment, where afterwards I grew up independent transformants in 64 wells of a 96-well culture plate.

They were arrayed in a checkerboard grid… 8X8 clones (yellow = NalR, and blue=NovR). If I prep DNA from all these clones, I could then produce Row Pools 1-8 and Column Pools A-H and each would have four clones of each resistant type. One issue would be distinguishing which endpoints belong together when segments are overlapping; another issue would be deciding which segments belong in the same clone.

If a donor segment appeared in clone 3C, for example, and it had unique endpoints (i.e. that donor segment is present only in clone 3C), then we would see those unique endpoints solely in pool 3 and pool C.

So we would have no difficulty assigning the segment to clone 3C.

On the other hand, if the segment was NOT unique, but present in, say clones 3C and 7E, we’d be unable to assign the segment to a particular clone due to "ghost" signals, but would instead know that there were two identical segments, but either in 3C and 7E, or in 3E and 7C.

(We’d be able to do this, since we’d still know the frequency of the segment in the different pools.)

So this is a good plan. We could first sequence by rows, giving us 64 more clones worth of data. And as long as there aren’t a whole bunch of identical endpoints for independent donor segments, we could then sequence pooled columns to assign segments to clones. If there were tons of identical endpoints, this would be such a shocking result, we’d need to re-think our next step anyways…

grantgrantgrant

2010-02-05T12:16:00.000-08:00

Whew! So another grant out-of-the way for now; another one almost done; and my own postdoc re-application in the works… A brief respite… Maybe it’s time to blog.

Last time (two months ago), I started showing some pictures from IGV showing our raw sequence data aligned to the Rd genome. Here, I’ll do yet another summary of our preliminary experiment, as we pitched it in our grant application, which will lead nicely into what I’d like to do next time… talk about alternatives to multiplexing DNA samples by barcoding…

So we’re studying transformational recombination in H. influenzae, where cells take up DNA from the media and incorporate it into their chromosomes. We think we have a decent model of the mechanism from studies in H. influenzae and other organisms:

(My figures here might have been a little degraded on their journey into the blog)
But until our little sequencing experiment, we could only infer the extent of transformational recombination of a chromosome based on the transformation and co-transformation frequencies of phenotypic markers. We’ve learned a lot just from obtaining four recombinant genotypes. Here’s what the experiment looked like in overview:
DNA from one isolate (NP NovR NalR) was incubated with competent cells of another (Rd), and transformants were selected. Two were NovR and two were NalR. We got sequence data from our collaborator for all four of these in a pool, two of them individually, and each parent (Rd and NP) individually. Here’s the figure we used to illustrate what our data looked like:
Hmm… probably that’s a low-resolution picture, but working from the bottom of the figure:
The lower panel shows the frequency of NP-specific SNP alleles across the Rd chromosome for the pool of four chromosomes. Blue dots at 25% indicate that 1 of 4 recombinants contained the donor-specific allele, while blue dots at 50% indicate that 2 of 4 recombinants did. The two red dots indicate the two selected markers (NovR and NalR), which as expected are at 50%.

In the upper panel, a zoomed view around the NovR-containing region is shown. The blue dots clearly define the donor DNA segments, but since there are overlapping donor segments, their appropriate assignment to different recombinants is unclear:
But because we also sequenced one of the NovR recombinants, the assignment of all the segments is made apparent. The green bars at the top of the figure show the donor DNA segments in Recombinant A, and so the donor segment spanning NovR in Recombinant B is unambiguously inferred.

Notably, there are several clustered donor segments in Recombinant A. This suggests that processes like mismatch repair may be disrupting larger original DNA fragments during recombination. For example in the upper panel of Figure 3 above, the area shown by the small purple circle appears to be a mismatch repair event around an insertional deletion difference between Rd and NP. Here is what that region looks like in IGV:
This IGV picture is showing our sequencing reads against the NP genome (the donor). The top track shows our Rd reads mapped to NP; the middle track show NP reads mapped to NP, and the bottom shows Recombinant A reads onto NP. I looked at the whole-genome alignment in this interval and found that the structural variation here is due to an insertional deletion: the alignment breaks and NP has 128 bp that doesn’t align with 52 bp of Rd.

Here is how I interpreted this event in the context of the larger NP donor segment:
Okay! Cool!

It’s going to take a while to fully parse this data, but more important is how we should go about collecting more. We certainly think we can increase our pool size, but as it is now, we can’t obtain “linkage” information from the pool. The obvious solution, barcoding individual DNA samples, presents monetary, technical, and computational problems. However, there may be another way…

Using IGV to look at Illumina data

2010-01-07T13:51:00.000-08:00

Out with the oughts in with the tens... Out with Maq, in with BWA!

Recap: I got a bunch of Illumina GA2 data from 5 DNA samples: Our recipient chromosome, Rd; our donor chromosome, NP; two transformed chromosomes; and a pool of four transformed chromosomes. I’d previously done a bunch of work with this data using the Maq alignment algorithm and SNP caller and hacked away at the SNPs in R.

Since then, I re-mapped our sequence reads using the BWA alignment algorithm, which was pretty much as easily installed and used (from the command line) as was Maq, though there are far more settings that might be manipulated. BWA has two big advantages over Maq:

1) It performs gapped alignment, so reads containing short indels can still be mapped to a reference sequence. This partially helped to overcome variation in read depth (coverage) due to mapping artifacts.
2) It outputs data in the SAM format, which is a newly-minted standard format for reference mappings of deep sequencing data. Thus I can use several downstream analysis tools that have been developed to work with SAM files and their binary equivalent BAM. Namely, the SAMtools package and the Broad Institute’s Integrated Genome Viewer (IGV) work with this format.

It still has one similar disadvantage as Maq, in that when a read maps to multiple locations in the genome, only one arbitrary mapping is kept, while all alternate mappings are discarded. This can create some odd artifacts at multicopy sequences in the genome. Largely this isn’t such a problem for me, but does cause some challenges in interpreting the correct base at every position.

To make a long story short, as expected, this re-analysis ended up yielding the same recombinant donor segments in our transformed chromosomes at a gross level, but I was also able to examine the genomes we sequenced in much greater detail by manually scanning the raw data in IGV.

To make the short story long...

I had an odd problem running IGV initially: To conserve memory, IGV only loads in a portion of the genome at a time. For large genomes with low coverage data, it works great with little RAM, so launching the application directly from the website is no problem. But because our dataset has such tremendous coverage of the genome, IGV was continuously crashing when I tried to run the 2Gb web version, and the lab computer only has 6Gb, so I couldn’t run the 10Gb web version.

So instead, I had to download the package and run it myself from the command-line. First I modified the relevant shell script for my computer called igv_mac-intel.sh by changing the switch -Xmx750m to -Xmx4800m using nano. This changed the memory used by IGV from 750 Mb to 4800 Mb. Then, after making the script executable (using chmod +x), I could run the script from the IGV directory by typing:
./igv_mac-intel.sh.

I did confront one other issue, which was that my FastA header was not identical to the name of the chromosome used in the SAM format file, so before loading the genome, I had to modify the FastA header from Genbank’s version to what I had used in naming the chromosome in the SAM file.

From there it was simple to load in a genome and then load in BAM files to look at the BWA alignments. (While this worked like a charm, I completely failed to add annotation tracks; IGV can apparently load GFF3 and BED formatted files, but mine wouldn’t load properly for whatever reason).

IGV’s visualization of read alignments really helped me understand the nature of this huge data set. I have now scanned through all the datasets against the Rd sequence reference and am part way through scanning through them all using the NP sequence reference. It’s been a rather cumbersome and slow chore, which I’ve spent the last three days doing, but it is also rewarding to look at the raw data, rather than something that I just have the computer spit out. I’ll get to that next.

I will use the rest of this post to go over some basic things about looking at this kind of data and will later show some interesting bits of our transformants.

Rd versus Rd

Here is what mapping our Rd sequence reads against the Rd reference genome looks like for a tiny segment of the genome. I picked a location that seems to have an “error-prone” region. There’s a whole slew of locations that look like this.
The top panel shows the portion of the Rd chromosome in view.

The next panel shows coverage, or “read depth” per position (so for this window, coverage ranges from ~800-1000 sequence reads mapped per position).

The three colored columns in the read depth panel are indicating potential sequence variants, based on a user-defined threshold percentage (I am using the early access version of IGV to have this control. I set it at 10%.) Just below the coverage plot is a running string of colored boxes, which indicate the bases in the Rd reference sequence.

So at the central base (position 904,862), the reference has a C (blue), whereas in our dataset, this position was covered 770 times, where 78% of reads were C, while ~21% were T (along with a few other stragglers).

In the lower panel is a pileup of the individual sequence reads. Grey indicates a base matching the reference, whereas colors indicate mismatches. As Ilumina sequencing is rather error-prone, as scattering of color is seen all over the place, but for most positions, the vast majority of reads match the reference base. The orientation of the read is also obvious in this view.

However, this is only showing a handful of the reads across this region. I can also compress the reads by right-clicking and selecting the “Collapse Track” option:
Now, the high degree of “errors” for several bases becomes evident. Other features also become apparent, namely the red and black labeled reads:

Red reads have an “orphaned” pair; i.e. the other read from the same molecule wasn’t mapped. This could be because the mate pair was a low-quality read, or could be because the mate pair was in sequence absent from the Rd genome.

Black reads have a paired read outside the user-defined maximum DNA fragment size (which I set at 300 bp). That is, the mate pair maps more than 300 bases away. As with most of the sequencing errors outside of the error-prone region, the red and black lines are fairly rare and evenly scattered. Manually checking where the black mate pairs were showed that most were just outside my threshold, so probably indicate no indel or rearrangement.

So what’s going on with the “error-prone” region? Because there’s not a lot of strangely mapping paired reads, I doubt this is a BWA artifact, but more likely is a sequencing artifact, possibly due to a bit of DNA that is “slippery” to the polymerase.

The interesting possibility that these are due to “clonal variation” (especially since the three variants shown in the coverage column are all ~25%) is unlikely in this instance, since there seem to be error-prone regions immediately adjacent to these three that just didn’t make my threshold.

NP versus Rd

The next image adds the NP sequencing data to the mix, still using Rd as the reference. Again, the error-prone bit shows up, further arguing that this isn’t due to clonal variation, since this is an independent culture and strain. However, on the left, several bona fide NP-specific SNPs were quite unambiguously identified.

The next image illustrates a couple things:
First, the region I previously showed had rather even coverage, whereas many, if not most, regions show much more variation in coverage. Here, coverage varies from ~200 to ~1800. In most instances, all DNA samples showed similar coverage for the same region, unless there is a structural variant.

Second, there is clearly a structural variant of some kind here between Rd and NP. The black discordant reads mostly map ~375 bp away, suggesting an ~75 bp insertion in Rd relative to NP. The red orphaned reads (mostly flanking the black reads) suggest that there is also an insertion in NP relative to Rd... so an insertional deletion then.

Okay, that’s it for now. Later, I’ll put up some more interesting stuff...

Vacation Post: Uncoiled!

2009-12-28T11:05:00.001-08:00

Number 6...

...in the "Top 10 New Species of 2009"

Number 10 is pretty mind-boggling...

Update on problems with analysis

2009-12-11T17:58:00.000-08:00

Ugh. Too much to blog about. The pace of computer work is totally different than with lab work (though in the end they seem equally labor-intensive), and I've made so many figures and gotten so many numbers in the past week, I barely know where to start...

Well, I guess I'll just blather about a couple problems I've been working through, but leave out all the charts and graphs and figures for the time being:

500-fold sequence coverage still "misses" parts of the genome:

We got a ton of data. For our controls, it was surely massive overkill. Nevertheless, "read depth" (how many read mappings cover a particular position) still varies by a substantial amount. There are numerous sources of this variation (%GC being one that is quite apparent), but I am most worried about variation in "read depth" due to the alignment algorithm I'm using to map reads to the genome.

As I try to root out artifacts in my discovery of putative recombination-associated mutations, I confront the fact that "read depth" is on average reduced when recombinant donor segments are mapped back to the recipient genome, so the novel mutations I found in these strains are on average supported by far fewer sequence reads than the average base... Most of them still look pretty solid to me (though there are several obvious artifacts), but I don't have a good rationale for whether or not to trust them.

I'm trying several things to work this out, namely by examinining the reciprocal mappings (using the donor chromosome as the reference).

So far, my analysis has a huge risk of false negatives:

Several problems here.

(a) Part of this problem and the last one is that I am using an alignment package that does not account for gaps (Maq). This means even a single nucleotide indel reduces "read depth" dramatically on either side (out to ~42 bases, or read length). See above.

(b) Another issue I'm facing with several of Maq's downstream outputs is that "read depth" is capped at 255. Presumably, they were conserving memory and only assigned a byte to this number. But what I haven't quite figured out is whether the SNP output (for example) is ignoring any possible SNPs where coverage exceeded 255. My cursory look at the more raw output (the "pileup") suggests this might well be the case. This could mean that I'm missing a lot, since the mean "read depth" per genome position in our datasets is ~500.

(c) Finally, I've been ignoring all Maq's self-SNP and "heterozygous" SNP calls in my downstream analysis using R. I presume that SNPs called in my mapping of the recipient genome to the complete recipient sequence are simply mutations between our wild-type Rd strain and the sequenced one. (As an aside, several hundred of the SNPs called by Maq were actually giving the correct base for an ambiguous base in the "complete" genome. I'd like to find a way to somehow revise the archived Rd sequence to get rid of all the ambiguous bases.) And I don't have a solid plan on how to deal with the "heterozygous" calls. Because the Maq assembly program can only have greater than or equal to two haplotypes, positions with mixed base signals are called heterozygotes. These is actually pretty cool and could reflect cool stuff like clonal variation, but largely these are probably due to multiply mapping reads and/or persistent sequencing errors.

Solutions: The solutions to these problems will initially largely be a matter of doing everything all again with different alignment software. My plan is to use the BWA aligner and SAMtools. BWA allows for gaps (so the "read depth" issue should be partially solved), and SAMtools not only keeps everything in a new agreed-upon standard format, but has several other tools, including what looks to be a better SNP caller (at least it has more modifiable settings). I would also like to try to do some de novo assembly, perhaps with Velvet, since we have such absurd coverage and a simple enough genome.

In the meantime, my R-fu has been improving, though I am convinced that I am missing some really basic principles that would make my code run a lot faster.

Mutagenic recombination?

2009-12-02T13:44:00.000-08:00

Okay, this is pretty cool. I will probably discover that it’s just some artifact of the mapping, but digging into the transformant data some more reveals what appears to be a high number of mutations within recombined segments (alleles that have neither donor nor recipient identity)...

For this analysis, I used much more stringent criteria for calling SNPs, so that low quality SNP calls would not contaminate the result. This is particularly important here, since we might expect to get lower quality SNP calls in the recombinant segments, due to the relatively high divergence between the donor and recipient genomes and the limitations of current mapping algorithms.

For TfA, there were 802 unambiguous donor alleles and 19 high-quality novel alleles, while for TfB, there were 902 unambiguous donor alleles and 21 high-quality novel alleles.

The two plots below indicate the presence of unambiguous donor alleles in blue bars going to 1 (which defined the recombinant segments), and the presence of unambiguous mutant alleles in red bars going to -1. (Click to enlarge)

That looks pretty striking! Mutations are clearly clustered into the recombined segments!

A few of the "novel alleles" in the two genomes are shared. 4 of 6 are in the first overlapping donor segment, and the other two are outside the donor segments. It is still early to be too confident in this result, but still! It is very suggestive.

I’ve never really taken the supposed causal connection between recombination and mutation too seriously, since the evidence mostly seems correlative to me, but if this result holds up, I think it will be a uniquely clear-cut example of mutations induced by recombination.

Before being confident in the result, I need to map the data back to the donor genome, and cross-check the result. If that works, some simple PCR and traditional sequencing should readily confirm or refute this deep sequencing result.

Recombinant Genomes: First Pass

2009-11-30T21:08:00.001-08:00

We’ve now got deep sequencing data of donor, recipient, and transformant genomes. And it is, indeed, “deep”, at least in quantity. Here’s the post where I illustrate the details of how the sequenced DNA was obtained. There were five “lanes” of sequencing obtained, of which I will talk about the first four here:

Lane 1: Rd (the recipient genome)

Lane 2: 1350NN (the donor genome 86-028NP plus NovR and NalR alleles of gyrB and gyrA, respectively)

Lane 3: Transformant A (the genome of a NovR transformed clone)

Lane 4: Transformant B (the genome of a NalR transformed clone)

So I’ve done a first-pass analysis of the sequencing data using Maq as the mapping algorithm. I left everything at default settings, and so far have only analyzed the datasets with respect to the recipient genome, Rd.

The expectation is that when I map Lane 1 to Rd, there will be few differences detected by Maq; when Lane 2 is mapped to Rd, the bulk of SNPs between Rd and 86-028NP will be detected; and when Lanes 3 and 4 are mapped to Rd, the donor DNA segments transformed into the recipient background will be identified. For lanes 3 and 4, since we know the locus that we selected for donor information, we have controls for whether or not an appropriate donor segment was detected.

Before continuing, I must specify that this is a very preliminary analysis: I do not consider the quality of the SNP calls, beyond whatever Maq does to make the calls (and culling ambiguous SNPs, as described below). I have not mapped anything with respect to the donor genome. I have not considered any polymorphism between the strains other than simple single-nucleotide substitutions (since Maq only does ungapped alignment). I also missed regions of very high SNP density, since the Maq default will not map any read with >2 mismatches from the reference. Finally, I have only cursorily examined depth of coverage across each genome (it is ~500-fold on average, ranging from ~100 to ~1000 over each dataset).

However, even with these caveats, the approximate sizes and locations of transformed DNA segments were pretty clear...

Here are the number of “SNPs” called by Maq between each dataset and Rd:

Lane 1: 933 (Rd)
Lane 2: 30,284 (1350NN)
Lane 3: 1,870 (TfA-NovR)
Lane 4: 1,881 (TfB-NalR)

Rd versus Rd

The first obvious issue is that when the Rd DNA was mapped against the Rd genome sequence, Maq called 933 “SNPs”. What is all this supposed “variation”?

108 had “N” as the reference base
432 had an ambiguous base as the query

So 540 / 933 “SNPs” are easily explained artifacts-- either ambiguous positions in the complete Rd genome sequence reference, or ambiguous base calls by Maq from the Illumina GA2 dataset.

The remaining “SNPs” may also be persistent sequencing/mapping artifacts, or they may be true genetic differences between the “Rd” we sequenced (our lab’s wild type) and the original DNA sample sequenced back in the 1990s.

To simplify matters I culled any position called a “SNP” between Rd and Rd, as well as all other ambiguous positions, from the remaining datasets.

Rd versus 1350NN

Before turning to the tranformants, I used Lane 2 to make a list of detectable SNP positions between the donor DNA and the recipient chromosome. Of the 30,284 “SNPs” detected by Maq, 29,002 were neither in the Rd SNP set nor had an ambiguous base in either the reference or query.

Note that I am not using SNPs identified by comparison of the two complete genome sequences, but rather those that were unambiguously determined by this sequencing experiment. Rd remains the only reference genome I have used.

I used this set of SNPs as the set of “Donor-specific alleles” to map transforming DNA segments.

The transformants

To identify donor-specific alleles in the transformants, I took the intersection of the “Donor-specific alleles” and the unambiguous SNPs identified by Maq for the two transformants, yielding the following number of SNPs in each transformant:

TfA: 890
TfB: 975

This suggests that about 3.0-3.5% of each transformant genome consists of donor DNA. This value is consistent with what we might expect, based on the co-transformation frequency of the two donor markers into the recipient genome when I did the original transformation that produced the sequenced clones.

Here are plots of the TfA and TfB genomes (genotype 0 = recipient allele, genotype 1 = donor allele):
(Note that images can be enlarged by clicking)

This makes me happy. I was kind of hoping our thinking was wrong and that there’d be all kinds of kookiness going on, but in many ways, having our expectations met is a vastly better situation, since it means that the designs of our other planned experiments are probably sound.

Note that the right-most "donor segment" represents only a single "donor-specific allele" that is identical in both recombinants, as well as being surrounded by Rd vs. Rd SNPs. It is highly likely that this singleton is an artifact. All other donor segments are supported by many SNPs, including the overlapping segments ~200,000 bp. This latter shared segment may suggest a hotspot of recombination.

Control loci

The following two plots zoom in on the segments containing the control selected loci. The first has red lines bounding the gyrB gene in TfA, while the second has blue lines bounding the gyrA genome in TfB. Also in each plot, the masked positions (those that were left out due to ambiguity or presence in the Rd vs. Rd comparison) are show (for the first, in orange; for the second, in grey):
Within each marker gene, there are a few recipient-specific SNPs. These have to do with the fact that the PCR fragment I used to add the NovR and NalR alleles of gyrB and gyrA contained recipient SNPs and some of these ended up in the donor genome.

Okay! That’s almost it for now. There’s a looong way to go, but happily I suspect that there is indeed biology to learn even from these “meagre” couple of genome sequences.

My next task will be to account better for where I am blind. I used stringent criteria to determine allele identity here. I am quite confident in what was found, but I’m not sure sure about what I didn’t find. That is, I suspect false negatives, but not false positives, based on how I’ve done this so far.

Incoming Data Deluge

2009-11-24T16:11:00.001-08:00

I am waiting with anticipation... I am about to be drowning in data... >64 million quality filtered paired-end sequencing reads, constituting 5.4 Gigabases! Woo!!

The files are currently transferring to an FTP site that I will then download from. Tomorrow afternoon, I may legitimately be a genomicist!

Thanks to my friend and his tech for getting this done...

Recalculating

2009-11-20T10:19:00.000-08:00

On Wednesday, I put up some tables of data with bad calculations of molecules per cell. Rosie's concerns about the data were misplaced; they should've been concerns about my brain!

I'd been saying that 8 ng of 200 bp DNA / 200 ul of competent cells was ~ 1 molecule per cell. But I was off. Not just because I calculated it for pg instead of ng, but because I'd screwed up other parts of the math...

Here's a more correct way to approximate the number of molecules per cell starting with one of Rosie's Universal Constants:

(1e-18 g / one 1 kb DNA) * (one 200bp DNA / five 1kb DNA) =
(2e-19 g / one 200 bp DNA) * (1e9 ng / 1 g) =
2e-10 ng / one 200 bp DNA.
So that's about how much a single 200 bp DNA molecule weighs.

Then the inverse gives molecules per ng:
1 / (2e-10 ng / one 200 bp DNA) =
(5e+9 molecules DNA / ng DNA) * (8 ng) =
4e+10 molecules DNA / 8 ng DNA.

Adding 8 ng of DNA to 2e+8 cells (in 200 ul of competent cells), then gives:
200 molecules per cell...

That's better. And more in-line with the lab's previous work...
I have changed the Tables in Wednesday's post to reflect this correction.

Bacterial Genome Sequence Links

2009-11-19T17:08:00.000-08:00

I keep meaning to put up these links, so I can find them in the future. They get you to NCBI's genome FTP site. If you click on the genome you want, you get a listing of a bunch of file formats. For the complete genomes, the choice for the FastA of the complete sequence is the .fna file. But there’s a bunch of others, which makes getting the correctly formatted data for application X easier:

Bacterial Genomes FTP

For the draft assemblies, the way to get the data is at this link:
Whole Genome Shotgun FTP

But to find out which accession you want is a little irritating and the best way is via searching this list:
Whole Genome Shotgun HTML

Using the supplied links from this last list is not so useful for downloading the data as is the FTP link, since you would have to click a maddening amount to get all the contigs from a particular whole genome shotgun sequencing project.

Saturation and Chaser

2009-11-18T17:03:00.000-08:00

EDITED: Used bad mol per cell calculation earlier.

How much DNA does an average cell in a competent culture take up? I did two experiments today, which were aimed at trying to decide how much degeneracy we can tolerate in our degenerate USS experiment. There's several issues, but one simple one is just to know how much of an "optimal" USS is taken up in the presence of competing "suboptimal" USS...

In the first experiment, I added different amounts of hot USS-1 DNA (the new-fangled ones designed for Illumina sequencing... ~200 bp) to a fixed amount of competent cells to find out how much DNA was enough to saturate the cells with DNA. Here's the results for several amounts of USS-1 incubated with 200 ul of competent cells for 20 minutes:
Looks like less than 8 ng (or ~200 molecule per cell) is sub-saturating and greater than 8 ng is near saturation.

In the second experiment, I looked at how well USS-1 is taken up in the presence of a competitor mutant USS (USS-V6, which is taken up 10X worse than USS-1). Based on the results of the first experiment, I used two different amounts of hot USS-1, either sub-saturating (4ng) or saturating (16ng). For each of these, I had four samples in which different amounts of cold USS-V6 were added. Here's the results for subsaturating USS-1:

And now for saturating USS-1:

Though adding the competitor did reduce uptake, uptake decreased only by ~1/4, even when there was several-fold more competitor. Also, it appears that using saturating amounts of USS-1 yields a different decline in uptake with excess inhibitor than when USS-1 is sub-saturating.

So what's this all mean? I think it means that even when "optimal" USS are only 1/10 as common as "suboptimal" USS, uptake proceeds quite nicely. What does this mean practically for making our fancy oligo purchase? Not too much, except that I think we can err on the side of more degeneracy without too many problems. (This will be good for several reasons, probably the most important of which is that our preliminary data collection plans will work best with more degeneracy.)

USS uptake with Illumina adaptors

2009-11-17T17:39:00.000-08:00

(I'm somehow deeply amused by the above image. I think it has something to do with the arrow going from the DNA molecule to the sun.)

I've finally gotten around to doing some uptake experiments with the control sequences I got for our degenerate USS experiments. The important thing was to make sure that the constructs would behave as our older ones do with our new-fangled design that will allow for Illumina sequencing directly from the purified periplasmic DNA (i.e. no library construction steps)...

And it worked quite well! I also scaled these experiments down, so that I could more easily do some saturation curves tomorrow.

Results today for uptake of 8 ng of 200 bp DNA (~160 million molecules) by ~200 million competent cells (in 200 ul):

O.G.-USS1 -> 33.5%
New-USS1 -> 35.4%
New-USSV6 -> 3.7%
New-USSR -> 0.7%

Right on target!

Missing Posts

2009-11-17T17:36:00.000-08:00

Whew! I keep falling behind on my blogging.... Here are the posts that I haven’t written:

Using R:

(1) Simulation of degenerate USS experiment, assuming the USS motif defines uptake specificity.

(2) Scoring genomes for USS sites.

(3) Scoring alignments for USS sites. (This last bit is finally begun, but far from complete.)

I've got the figures for these posts; now I just need text!

Weening myself off of Excel

2009-11-06T16:59:00.001-08:00

In some sense, the computing I did today isn’t really useful, since I already worked out these things using Microsoft Excel. But I’ve been ordered by my bioinformatics consultants to stop with the Excel already. So as practice, I worked out some of the expected features of degenerate oligos again, but this time using R.

The main motivation for doing this besides practice is that I am fairly sure we should be ordering degenerate oligos with more degeneracy than we have previously considered. I won't make that argument here, but just repeat some analytical graphs I'd previously made.

It took a while (since I’m learning), but was still much more straight-forward than doing it in a spreadsheet. The exercise was extremely useful, as I learned a bunch of stuff (especially about plots in R), while doing the following:

Problem #1: Given a percentage of degeneracy per base, d, in an n length oligo, what is the proportion of oligos with k mismatches?
Answer #1: Use the binomial distribution. For a 32mer with different levels of degeneracy (shown in legend):
Problem #2: Given a million instances of such an oligo, how well would each possible oligo with k mismatches be observed?
Answer #2: Simply adjust each of the above values by dividing the number of classes within each of k mismatches (i.e. choose(n, k)):
Problem #3: If some number of bases, m, in the n-length oligo are “important”, what proportion of oligos with k mismatches will have x “hits”?
Answer #3: Use the hypergeometric distribution. The below plot is as for Problem #1 for 0.12 degeneracy, but with the # of hits broken down for each k:
I didn't try super-hard to make the perfect graphs, but it did take some effort to make a stacked bar plot...

Funding versus Science

2009-11-06T10:45:00.001-08:00

Yesterday, I ended up doing none of the things I intended to do, but cleaned them instead: (1) No bench-work, but a clean bench (and fridge and freezer spaces)!; (2) No emails, but a clean inbox!; (3) No thinking, but a clean mind! (Well, moreso than usual.)

My ears were also itching, because I realized that the NIH panel that's review my postdoc grant application was meeting, and was extra-worried, since now the Genome BC grant almost certainly will depend on the outcome of those scores.

Anyways, towards the end of the day, I switched over to browsing journals, which I rarely do these days...

I found a lot of good stuff I probably already should've known about, but I also found this opinion piece in PLoS about "what's wrong with funding of research". Pretty much sounds about right. I am not even in a position to have the stresses described in the article, but feel like if I don't get some kind of postdoc fellowship of my own, I'll be that much less likely to get hired as an independent researcher, since so much of what defines success is the ability to get money. But all of the PIs I know are overwhelmed by exactly these issues. It's especially daunting to realize that a full-sized NIH R01 grant can barely support a lab with 2-3 people.

On the other hand, my recent experiences with Rosie do bring home the fact that grant-writing is a good way to think rigorously about one's plans, and the Genome BC experience in particular (whether the grant is funded or not) helped our brains wrap around exactly what we'll be doing in the next couple of months.

Decompression

2009-11-04T13:30:00.001-08:00

Whew! Rosie and I just made what I consider a heroic effort to produce a grant application to Genome BC to use DNA sequencing to measure recombination biases during H. influenzae natural transformation.

It was heroic not only because we finally decided to apply only late last week, but because our co-funding support is tenuous at best. Genome BC requires that we match their funds at least equally with funds from another source. Rosie has funding, but it was applied for too long ago and the proposal only indirectly relates to our planned sequencing. We also applied for a CIHR grant recently, but will not have reviews until after the Genome BC committee meets in early January. The best hope of adequate co-funding comes from my NIH postdoctoral fellowship grant application (a resubmission) late this summer, for which I should have scores (or lack thereof) within a couple of weeks. If the application gets a good score, we can tell Genome BC that the major threat to the success of our application is ameliorated.

Almost immediately after submitting the grant application with Rosie last night, I had to turn to editing my buddy's manuscript (which I am an author on), which takes on the weighty topic of detecting rearrangements in complex mammalian genomes from limited sequencing data. I just turned my edits over to the corresponding author, and now, after letting the excess nitrogen out of my bloodstream, I need to decide what to do next.

Since aforementioned buddy also has taken several DNA samples off my hands for sequencing, I think I'd best turn to purifying uptake DNA from the periplasm of competent cells. I've already gotten things fairly well under way (see here, here, and here), but there's a bunch of uptake experiments waiting to be done. So tonight, I'll inoculate some cultures, so I can try a large-scale periplasmic DNA prep tomorrow, and tomorrow I'll also order some more radiolabel for doing more sensitive uptake experiments.

Sets of Snps

2009-10-27T13:13:00.000-07:00

So I’ve got pure DNA from my transformants, pretty much ready to send off for our first Illumina runs. I’m just doing a few simple checks with PCR and digestion to make sure everything is kosher.

But there is this little fear in my head about what I’ll do when I get the massive datasets. I got a hold of some example Illumina GA2 paired-end data taken from E. coli K12. Since I don’t have any data of my own yet from H. influenzae, this seemed like a good dataset to start learning how to do the initial data analysis.

I decided to go with the most widely-used reference-mapping software, called “Maq” for “Mapping and Alignment with Quality”. I can see why it’s widely-used; it is a breeze to install and use, which is the primary requirement for end-users to like a given piece of software. I’ve started dealing with just single-end reads from only one lane. The data is several years old, so there are “only” a few million reads of 35 base-pairs. Nevertheless, this represents nearly 20X sequence coverage of the E. coli K12 genome.

I’ll keep the Maq overview brief, since I went through it in lab meeting and Maq’s documentation is largely quite good (with the caveat that there are some challenges interpreting the meaning of all the columns in the different outputs). In short, Maq takes a “FastQ” file (which is the sequence data, including quality scores) and a “FastA” file (the sequence reference), converts them into a binary format (which is apparently helpful for making the mapping algorithm fast), and then maps individual reads to the reference, allowing up to 2 mismatches. The “pile-up” at each position is used to decide the “consensus base” at that position, based on the majority base and associated quality scores. The mappings can be extracted from the binary output using additional commands.

Here, I’ll focus on the .snps output (obtained by cns2snp... from the man page), since this will be the most straight-forward way for us to define recombinant segments in our transformants. I’m keeping things simple still, so there are a lot of other issues I could get tangled up in, other than the ones I’ll discuss here.

So I “Maqed” this one FastQ dataset from E. coli K12 (I’m calling it s11) against two different reference sequences:

K12 -> NC_000913.fna
O57 -> NC_002695.fna

The first is a sort of control for how the consensus base caller is working. Since the sequencing is from K12 DNA, most consensus bases called by Maq should match the reference. Occasionally some “SNPs” might appear. These could have several sources:

Spurious due to low coverage and high error
Actual mutations between the original completely sequenced K12 and the K12 used for the Illumina sequencing
Clonal variation where the consensus-calling favors the variant sequence.

The second is to call SNPs between the two strains. I would expect that the number of SNPs called against O57 would vastly exceed that of K12.

This is easily decided with a little UNIX command on the two files to count the lines (which each correspond to a SNP call):

> wc -l K12s11.snp O57s11.snp
6051 K12s11.snp
64217 O57s11.snp

Indeed, there are 10X more SNPs running s11 against O57 than against K12. The several thousand SNPs called against K12 are likely mostly errors. Maq doesn’t always assign the consensus base with A, C, G, or T, but with any of the other IUPAC nucleotide codes, so many of these “SNPs” are probably quite low confidence.

That’s all fine and good, but now what? How can I tell how well this s11 dataset did at calling all the SNPs between K12 and O57? I resorted to using Mummer’s Dnadiff program (previously discussed here) to compare the two reference sequences to each other and extract the SNPs. If the s11 dataset really did a good job sequencing the E. coli genome, there should be a strong overlap between the SNPs called by Mummer and those called by Maq.

(Here’s a GenomeMatcher DotPlot run with Mummer)
Again, UNIX came to the rescue, thanks to this page I found that provides UNIX commands for working with “sets”.

First, I had to get the appropriate three columns from the two SNP files:

the reference genome position (in O57)
the reference base (also from O57)
the query base (the SNP in K12; the consensus base call for Maq or the SNP call for Mummer).

Here's how I made those files:

> cut -f 2 -f 3 -f 4 O57s11.snp > maqSnp.txt
> cut -f 1 -f 2 -f 3 O57vsK12.snps > mumSnp.txt

This provided me with my two sets, one called maqSnp.txt and the other mumSnp.txt. Here’s their lengths:

> wc -l maqSnp.txt mumSnp.txt
64217 maqSnp.txt
76111 mumSnp.txt

Notably, Mummer called many more SNPs than Maq. I think this is largely because the SNP output from Mummer includes single-nucleotide indels, which Maq misses since is does an ungapped alignment. I’m not sure how to deal with this, but in our real experiments, they should still be discoverable, since we’ll map our reads to both donor and recipient genomes. Also, there are numerous “SNPs” in the Maq file that are non-A, C, G, T consensus bases, which will largely be absent from the Mummer comparison.

So, then it was a piece-of-cake to count the SNPs in the intersection of the two sets. Indeed there were several ways to do it, my favorite of which was:

> grep -xF -f maqSnp.txt mumSnp.txt | wc -l
57179

So:

That’s a whole lotta overlap! Happy! When I get the real transformant data, this will help me tremendously.

Transformants Produced!

2009-10-23T14:20:00.000-07:00

I'll set aside the hack bioinformatics posts for now and give an update on my transformation experiments... We’ve gotten access to a few lanes of Illumina GA2 sequencing for some preliminary studies, and right now I’m drying the genomic DNA samples that we plan to sequence.

The notion is to sequence several independent transformants of Haemophilus influenzae to get some idea of how much donor DNA taken up by cells finds its way into recipient chromosomes. This pilot study will go a long way in informing our planned large-scale experiments and give us a chance to learn how to handle the data.

Here’s what I did to produce the material...
...some transformations, of course!

First, I PCR amplified the gyrA and gyrB alleles from the MAP7 strain (which confer nalidixic acid and novobiocin resistance, respectively). MAP7 is a derivative of our recipient strain KW20 containing several point mutation that confer antibiotic resistances.

I used these PCR products to transform our donor strain 86-028NP to provide two selectable markers in the donor. I’ve been calling this strain 1350NN.

Then I extracted DNA from this strain and used it as the donor DNA to transform KW20 competent cells. By selecting for one or both markers, I can ensure that clones chosen for DNA extraction and sequencing were indeed derived from competent cells that got transformed.

Our baseline expectation is that there will be a large segment (10-50kb) of donor alleles in the transformants at selected sites and 2-3 additional large segments elsewhere in the genome.

Originally, we were going to do this transformation with only a single marker, but we realized that having two would allow us to measure the frequency of co-transformation.

Here’s what the transformation rates looked like:
I used MAP7 DNA as a donor as a control. Since MAP7 is more closely related to KW20 than 86-028NP, it is perhaps unsurprising that transformation rates were higher when using MAP7 as donor.

As for co-transformation, here’s the frequency of double transformants versus expected:
That corresponds to ~25-35% of the cells in the competent cell preparation actually being competent. I’ve been wracking my brain unsuccessfully trying to figure out how to do a back-of-theenvelope calculation as to how many independent molecules we expect to transform any given recipient. I just can’t figure out a concise or reasonable way to do it. Suffice it to say, I estimate a minimum of 20 kb of donor DNA in each transformant (1% of the genome), up to perhaps 100 kb (5% of the genome).

There’s only one way to find out…

Computing Bootcamp

2009-10-15T11:48:00.001-07:00

Whew, I’ve really fallen behind on my blogging... Last week, a good friend of mine came into town for a “northern retreat”, in which he hoped to get work done on a paper. Instead, he and I drank enormous amounts of beer and did an enormous amount of computing with the Haemophilus infuenzae genome (at least by my standards). While the beer probably didn’t help anything, the computing did.

I’ll go over some of what we did in future posts, but right here I just want to outline some of the computing lessons I learned looking over his shoulder over the week. Many of these lessons have been given to me before and are likely quite basic for the real computationalists out there, but somehow I’ve emerged from the computing emersion with a lot more competence and confidence than I had before...

Here's three useful things I'm getting better at:

(1) Avoid using the mouse. The more that can be accomplished from the command line and using keystrokes, the better. From the command line, tab-completion and cursor-control of the command history make issuing commands far more efficient. The coolest keystroke I’ve now picked up the habit of in the Mac Leopard OS is Cmd-Tab, which takes you to the last active open application (and repeated Cmd-Tabs cycle through the open applications in order of their previous usage). This is perfect for toggling between the command-line and a text-editor where one can keep track of what one is doing.

(2) Toggle between the command-line and a text-editor constantly. Rather than trying to write a whole script and then run it from the command-line, it was far easier and faster to simply try commands out, sending them to the standard output, and cobble together the script line-by-line, adding the working commands to a text document. This has three useful effects: (1) Bugs get worked out before they even go into a script, (2) It forces one to document one’s work, as in a lab notebook. This also ended up being quite useful for my lab meeting this week, in which I decided to illustrate some stuff directly from the terminal. (3) It is forcing me to work “properly”, that is sticking with UNIX commands as much as possible.

(3) Learn how the computer is indexing your data. This point is probably the most important, but also the one that is taking me the most effort. I’ll illustrate with an example (which I’ll get into in more scientific detail later):

The output of one of our scripts was a giant 3 column X 1.8 million row table. I wanted to look at a subset of this huge table, in which the values in some of the cells exceeded some threshold. At first I was doing this (in R) by writing fairly complicated loops, which would go through each line in the file, see if any cells fit my criteria, and then return a new file that only including those rows I was interested in. When I’d run the loop, it would take several minutes for finish. And writing the loop was somewhat cumbersome.

But the extremely valuable thing I learned was that R already had all the data in RAM indexed in a very specific way. Built-in functions (which are extremely fast) allowed me to access a subset of the data using a single simple line of code. Not only did this work dramatically faster, but was much more intuitive to write down. Furthermore, it made it possible for me to index the large dataset in several different ways and instantly call up whichever subset I wanted to plot or whatnot. I ended up with a much leaner and straightforward way of analyzing the giant table and I didn’t need to make a bunch of intermediary files or keep track of as many variables.

Next time, I’ll try and flesh out some of the details what I was doing...

Corrected Logos

2009-10-01T15:04:00.001-07:00

Yesterday, I attempted to make some logos out of the degenerate USS simulated data that Rosie sent me. Turns out, I was doing it wrong. After checking out the original logo paper, I was able to figure out how to make my own logos in Excel. I wasn't supposed to plot the "information content" of each base; I was supposed to take the total information content (in bits) of each position (as determined by the equation in the last post) and then, for each base, multiply that amount by the observed frequency of that base to get the height of that element in the logo. So, below I put the proper logos for the selected and unselected degenerate USS sets (background corrected):

Here's the selected set:
Here's the unselected set:
Woo!

Struggling with the Background

2009-09-30T17:59:00.000-07:00

Previously, I had shown some preliminary analysis of Rosie’s simulated uptake data of chromosomal DNA fragments. Rosie also sent me simulated uptake data of degenerate USS sequences (using a 12% degeneracy per position in this USS consensus sequence :

5’-AAAGTGCGGTCAAATTTCAGTCAATTTTT-3’).

So what can I do with this dataset? Well, first, since Rosie also provided me with the “scores” for each sequence, I could plot a histogram of the scores for the 100 selected and 100 unselected sequences, showing that the uptake algorithm seems to work pretty well.

Here's the histogram:
Notably, even the unselected sequences have rather high scores, when compared to the same analysis of genomic DNA fragments. This is unsurprising, since the sequences under selection in this degenerate USS simulation are all rather close to the consensus USS.

Here’s the histogram from the genomic uptake simulation again (just to compare):
(I think the reason for the difference in the “selected” distributions is due to a different level of stringency when Rosie produced the two simulated datasets.)

That’s all fine and good, but now what? With the genomic dataset, I could use the UCSC genome browser to plot the location of all the fragments I was sent, but this consists of two alignment blocks of sequence that look markedly similar.

The obvious thing to do was to make Weblogos of the two different datasets… The unselected set should have very little information in it, while the selected set should contain information. In doing this, I discovered a rather important issue… The on-line version of Weblogo does NOT, I repeat, does NOT account for the background distribution.

This is a problem. It means that whenever you make a Weblogo (on the webserver) from your alignment block, it is assuming that each base is equally likely to occur at a random position. This is why the y-axis in all Weblogos plots always has a maximum of 2 bits when using DNA sequence. Why is this a problem? First of all, if one is using an AT-rich genome, as we are, then the information content of any G or C is underestimated and any A or T is overestimated.

So how does Weblogo calculate the information content of each position in an alignment block? From the Weblogo paper (link found at the Weblogo website):
Rseq = the information at a particular position in the alignment
Smax = the maximum possible entropy
Sobs = the entropy of the observed distribution of symbols (bases)
N = the number of distinct symbols (4 for DNA)
pn = the observed frequency of symbol n

The log2 is there to put everything in terms of bits.

So for DNA (4 bases), the maximum entropy at a position is 2 bits. Makes perfect sense: 2 bits a base. However, this only makes sense if each base is equally probable for a randomly drawn sequence. Now for purposes of gaining an intuition for different motifs, this isn’t really a big deal, although it does complicate comparing motifs between genomes.

When this isn’t the case (probably much of the time), then different measures have been used, namely the “relative entropy” of a position. This is an odds ratio of the observed probability and the background probability. Apparently, the off-line version of Weblogo can account for non-uniform base composition, but I haven’t tried installing it yet, nor any other software out there that handles variation in GC content.

Why? Because what we need for our degenerate sequences is a different background distribution at each position! So, the first position in the core is 88% A, but the third position is 88% G!

To illustrate the problem, here is a Weblogo of the selected set:
Here’s the unselected set:
Looking closely, it is clear that there are differences in the amount of “information” at each position. So in the strong consensus positions of the USS, the selected set has higher “information” than the unselected set, while at weak consensus positions, that’s less true.

But the scaling of each base here is completely wrong. There isn’t nearly a bit of information at the first position in the unselected set. We expected 88% A. The fact that there are mostly As in the alignment block at the first position is NOT informative. In fact, if all was well, we’d get zero bits at all the unselected positions!

What to do? I tried to make my own logo, using the known true background distribution at each position. I won’t belabor the details too much at the moment, except to say that I had to figure out what a “pseudocount” was and how to incorporate it into the weight matrix, so as to not ever take the logarithm of zero.

Here’s the selected set of 100 sequences:
Here’s the unselected set of 100 sequences:
(Note that I somehow lost the first position when I did this, so the motif starts at the first position of the core.)

This actually looks quite a bit better, or at least more sensible. A few things worth noting:

I think if we did thousands of unselected sequences, we’d pretty much get zero information from that alignment, which is what we would want, since that’s just the background distribution.
Some values are negative. This is expected. Since these are scaled to log-odds ratios, when the frequency of seeing a certain base is less than the expected background frequency, a negative number emerges.
The scale is extremely reduced. Every position is worth less than 0.3 bits. This is also expected. One description I’ve seen about how information content can be thought about is how “surprised” one should be when making an observation (there’s even a unit of measure called a “suprisal”!). Since we are drawing from an extremely non-uniform distribution that actually favors the base that’s expected to be taken up by cells better, we are basically squashing our surprise way down. That is, getting an A at the first position of the core is highly favored, but it’s the most likely base to get anyways, even in the absence of selection.
The unimportant bases in the USS have the most information content in the selected set. At first this bothered me, but then I realized it was utterly expected for the same reason as above. For example, at position #18 above (sorry, it’s position 19 in the Weblogos), the selection algorithm doesn’t really care what base is there. That means that the selected set will let mutations at that position (from A to something else) come through, which will be surprising, when compared to the background distribution!

(ADDED LATER: Actually, this last point is wrong. The reason for so much information at the weak positions is related to the matrix that was used to select the sequences, not from surprise. I'll try and get a proper dataset later and redo this analysis. To some extent, the positions will still have some information, as partially explained in my erroneous explanation above, but not nearly so much.)

Whew! I’ve gotta quit now. There’s a lot more to think about here.

Mismatch repair versus Segregation

2009-09-28T12:45:00.000-07:00

Things have gone swimmingly with my strain construction plans, and indeed today I am extracting DNA that will presumably be sequenced. To recap, I made a couple of clinical isolates (86-028NP and PittGG) resistant to novobiocin (NovR) by transforming them with a bit of left-over NovR allele of the former postdoc. I then isolated the new strains’ DNA, and used these to transform the standard KW20 Rd strain. By selecting for NovR, we can be certain that the clones I pick took up DNA and recombined it into their genomes.

One technical issue arose, however, which required a little bit of thought: Should I have streaked for single colonies? I.e. once I had my transformants, it might be a good idea to streak out individual colonies to make sure I purified them away from any background or broke apart any doublet colonies. No big deal, but after talking it out with Rosie, we decided to skip it. Why? So that we might get lucky and distinguish recombination followed by mismatch repair versus recombination followed by segregation. In the following figures, I illustrate what I mean by this…

In this first one, the donor DNA is shown in red, and the recipient chromsome is shown in two colors, blue and green, to distinguish the strands. The lowercase letters indicate polymorphic sites in the donor genome. Little a is meant to be the selectable marker, in this case an allele of gyrB:
Donor DNA is incubated with competent recipient cells, and recombination of single-stranded DNA leaves patches of heteroduplex in the genome, shown as small red patches on either the blue or green strands.

After this, the cells have a chance to perform mismatch correction to fix any heteroduplex. I select for cells that have little a by plating to novobiocin plates, so only cells that end up a/a will survive an make colonies. (I am not going to show any examples of restoration repair, in which donor alleles are repaired back into recipient alleles… this will be invisible in our analysis.)

In the below example, I show the A/a and B/b heteroduplexes getting mismatch repaired into a/a and b/b, whereas C/c and D/d heteroduplexes remain unrepaired (they escape correction). What will happen in such as case is the generation of a sectored colony, in which (in principle) half the cells would have one genotype and the other half a different genotype:
In the above example, the original transformant segregates the c and d alleles into different cells, while a and b end up in all cells. If the whole resulting colony is grown up and sequenced, the a and b alleles will be the only ones observed, while at the other two loci, there will be a mix of C and c, along with a mix of D and d. We wouldn’t be able to tell “phase”, i.e. whether c and d were on the same or different chromosomes, unless we did streak for singles and the sequenced several clones. But as a first pass, this could be a really interesting analysis. It will also serve as excellent proof-of-principle for our more intense sequencing plans.

There is a caveat, however, which means we need to get a little bit lucky to be able to distinguish these phenomena (mismatch repair versus segregation). We won’t see two different genotypes, if the A/a heteroduplex isn’t mismatch corrected:
The issue isn’t that segregation didn’t happen; the problem is that one of the segregants dies under selection for little a.

Thus, if we see a pure genotype, then either all mismatches were corrected, or our selectable marker didn’t mismatch correct.

When I pre-screen my transformants to make sure they’re not spontaneous mutants, I might be able to pick a colony where I think segregation is occurring. If I get the standard sequencing traces back and see mixed bases in the chromatograms that corresponde to donor and recipient alleles, I’ll pick that kind of clone for sequencing…

One sort of sad note here, in terms of the more distant future, is that mismatch repair mutants, which should be quite useful for understanding transformation, will need to be transformed without selection if we hope to recover isolated segregants from individual transformants.

More simulated uptake

2009-09-25T17:33:00.001-07:00

Thanks to Rosie eliminating the int function from her Perl model, I got to take a look at some more simulated uptake data. Last time, there were several issues, which now seem solved. This time to model uptake, she used the real genomic USS position weight matrix to stochastically select 500 bp fragments from the first 50 kb of the Haemophilus genome. I got 200 from the forward strand and 200 from the reverse complement strand, along with a set of random fragments. This is 4X spanning coverage of the 50 kb…

Below is the way the data looked in the UCSC genome browser, added as custom tracks (click on the figures to enlarge).
From the top, the tracks are:
(1) Chromosome position
(2) 400 random fragments (shades of brown).
(3) 400 selected fragments (shades of blue).
(4) Positions of “perfect core” USS motifs (5’AAGTGCGGT-3’) on either strand
(5) RefSeq gene annotations.

Here’s a bit of the 50 kb zoomed in:
That looks pretty good for such low coverage! (In our real experiment, we expect to get several hundred times more data.) It’s starting to look like a real model of how uptake might look! The random fragments look roughly randomly distributed, and the selected fragments clearly show a punctate distribution around the “perfect core” USS sites.

Here’s a histogram of the scores of the best site on a given fragment for the random and selected datasets:
Indeed, the distributions are quite distinct, though notably the distribution of random fragments looks bimodal. This may simply be a feature of the genome, since there are so many USS sites… Worth thinking about though.

There are other details obscured in the browser figures: the shading indicates the relative score of the best site on the fragment (on a log scale), and each fragment also has orientation shown as a small arrowhead within the box. I’ve also associated each fragment with its score. So in the browser, I can easily check things out more carefully:
In this zoom, it’s clear that there is an excellent site to the left (just under 18,000), a weaker site to its right (~18,400; fragments are overlapping with the left-most site; no perfect core), and a pair of sites on different strands with different scores to the right (~19250). I can also retrieve the sequence associated with a given fragment to see if I can spot the USS site within it. And if I really zoom in, the DNA sequence is listed at the top.

It’s not really so easy to see what’s going on with all these overlapping fragments, so my next task will be to convert this data from fragment positions to spanning coverage per chromosomal position (though I’ll probably bin positions, perhaps every 100 bp to keep things reasonably small for now). I will take a stab at doing this properly (with a script) but may wimp out and do it in Excel. If I can then muscle this data into a WIG formatted file, then I’ll be able to plot the data in a way where “good” sites will look like peaks in coverage…

E-Z Strain Construction

2009-09-24T16:13:00.000-07:00

As preliminary data for our genome-wide recombination analysis (outlined in this post from Rosie), we want to sequence the whole genome of a single transformed clone in the next couple of months. The idea is to transform our standard KW20 Rd strain with DNA from one of the other completely sequenced strains (probably 86-028NP, possibly PittGG), select a single transformed colony, and sequence its genome.

This will provide us with all sorts of useful preliminary results:

Show that we can indeed handle the type (and amount) of data we’ll be obtaining.
Estimate the total amount of donor DNA a single recipient recombines (and fixes) into its genome.
Estimate the length of recombination tracts (gene conversions) / the strength of “linkage”.
Estimate mosaicism of donor and recipient sequences (mismatch repair).
Estimate the transformation rates for different classes of single-nucleotide differences (for example, the number of A->T transformation events observed versus the total A->T differences between the strains)

In particular, item (2) will be crucial for estimating the total amount of sequencing we would need to measure transformation rates per polymorphism across the genome. Simple transformation assays with DNA from the multi-antibiotic resistant MAP7 strain suggest that possibly 20-50kb of DNA may be replaced in a single transformant, but this type of analysis is restricted to only a few different sites in the genome and is very roughly calculated.

The analysis of a single transformed genome will still be preliminary with regards to (3)-(5), for which we will want genome sequences for several independent clones. In the future we are likely to barcode and pool independent transformants, since we expect that a single lane of Illumina sequencing will be overkill for a single Haemophilus genome of less than 2 Mb (250X sequence coverage).

Anyways, one issue with producing the material for this first sequencing experiment is that we need to make sure that the clone we select comes from a cell that was indeed competent and did indeed get transformed. Since only a fraction of cells in a competent culture are competent, we would be wasting a lot of time and money, if we accidentally just re-sequenced our recipient genome.

In order for this to work, we need our donor strain to carry an antibiotic resistance marker. By selecting for recipients that become resistant, we can be sure the clone we select took up DNA that got recombined into the genome. (This may also create a bias for donor alleles near the selected site, due to “linkage”.)

To this end, I am doing the following:

Made a couple strains (KW20, 86-028NP, and PittGG) resistant to novobiocin. I just did this. It worked like a charm thanks to the former postdoc having a well-organized lab notebook and a well-organized freezer box containing a tube with a NovR allele of gyrB already prepared for me. This was also my first time doing overnight transformations. I couldn’t believe how easy it was: Add a frozen aliquot of cells and some DNA to some sBHI media, let the cells grow overnight, and plate them the next day. There were plenty of resistant colonies this morning.
Prepare DNA from the newly produced 86-028NP NovR and PittGG NovR strains. I’ll do this tomorrow from the overnight cultures I just inoculated.
Transform KW20 with this DNA. I’ll use competent cells I already have tomorrow, after my DNA prep.
Saturday, assuming I have NovR transformants, I’ll pick and grow up some transformed colonies overnight.
Sunday, I can prepare this DNA, and that’ll be what we can send for sequencing!

So if all goes well, we should have our material in a few days! Then we wait. Then the real work begins…

(As a side note, 86-028NP indeed appears to already be resistant to another antibiotic, nalidixic acid. I will check to see if this resistance is transformable when I have the 86-028NP NovR DNA in hand.)

Fake Periplasmic Data

2009-09-23T11:48:00.000-07:00

UCSC Microbial Genomes Database was a nice find for me, since they host the Haemophilus influenzae KW20 genome. It has pretty much made me forget about my own plans to make a custom browser for the moment. Though that will change, as when we have our own data, we’ll absolutely need some off-line way to browse our datasets, since they will be so large…

As a first fake experiment to explore how our periplasmic DNA pools might look, Rosie sent me two sets of 200 sequences. One set was 200 randomly chosen 100mers from the first 10 kb of the Haemophilus genome, and the other set were 200 sequences (100mers) stochastically selected for the presence of a USS using her Perl scripts. All I had to do was turn her data into a BED formatted file, which only took a few minutes. As usual, I made the BED file using Microsoft Office, rather than a more savvy command-line way, which would've probably used Grep or something.

Here’s what her data looks like plotted as a custom track (squished) in the UCSC genome browser:

RANDOM
SELECTED
It looks like it sort of worked! There’s a prominent peak containing nearly half the sequences in the selected pool, while the random fragments look just like they ought to.

One issue here is that we know there are two other perfect matches to the core USS motif in the first 10 kb, and these weren’t captured by the selection algorithm. It’s slightly unclear why that is, but might have something to do with the USS position-weight matrix that was used. (Actually, there are six USS in the interval, but we were only searching one strand this time...)

A beginning!

The Last Straw

2009-09-17T23:06:00.000-07:00

Yesterday, Rosie kindly ran her Perl script over the USS construct I designed. The final thing I was worried about was whether or not my design had any USS or USS-like sequences in it, other than the one it's supposed to have. I'd checked the construct for any core USS motifs (5'-AAGTGCGGT-3'), but since we think that the motif is more complex than this, it was important to make sure that there were no extra sequences that got high scores using the USS position-weight matrix.
Fortunately, the construct looks good, so I can go ahead and order the control oligos and have high expectations that they'll work...

Here's how every 32 base pair window over the 199mer looks when scored with the USS PWM:
There's a single prominent high-scoring site right where it should be, and all of the surrounding area scores near background. The USS in the construct has a score (~10^-8) more than 10 orders of magnitude better than the next best sites. There's a slight increase for windows immediately adjacent to the USS, presumably because the AT-tracts in the USS are still contained in those windows. The rest of the construct only has scores at background.

Just to show that these other sites really do represent background levels of USS score, Rosie also ran a randomized version of the sequence:
Nothing better than 10^-18. Excellent.

Scale UP!

2009-09-17T15:40:00.000-07:00

How much periplasmic DNA can I hope to get using my current protocols, and how much DNA will I need? As promised, here are some rough calculations regarding the oligo purchases that we want to make.

I used a molecular weight calculator available on-line to determine the size of the dsDNA I described in the last post.
I alternatively could’ve used Rosie’s Universal Constants (660 g / mol of base pair and 10^-18 g / single 1 kb DNA molecule) to make this calculation, but since I’m dealing with a known sequence, I might as well get an exact molecular weight. (I also made a minor mistake in the last post, and the molecule I describe is actually only 199 bp).

So, for our USS molecule, MW = 122,828.6 g / mol. And the oligo synthesis service we’re planning on using will be at the 1 micromole scale. That means if we took all of the two oligos, annealed, extended, and purified, we’d end up with 0.123 grams of input DNA pool! That’s really a very large amount.

My previous concerns about needing to do PCR to maintain the pool are unfounded. This scale should be sufficient for hundreds (or even thousands) of experiments...

What follows are my preliminary assumptions about yields from the periplasmic DNA prep. They are based on several different experiments, though I am erring on the side of being conservative with my estimates and are guides for future experiments only. In this post, I will address the issues of scale-up at the end; before that I’ll just refer to the approximate total culture volume and amount of DNA that I’d need to get a target amount of DNA, assuming all else works perfectly.

So, I’ve now done several experiments using a PCR fragment bearing the consensus USS, called USS-1. If I add 20 ng DNA / 1 ml competent cells, ~50% is taken up. That is, in rec-2 cells, my theoretical yield of periplasmic DNA is 10 ng. My actual yield is considerably lower; as evaluated by my radiolabeling experiments, I estimate I get ~25% of my theoretical maximum.

This means,

1 ml cells + 20 ng DNA → 2.5 ng recovered.
20 ml cells + 400 ng DNA → 50 ng.
40 ml cells + 800 ng DNA → 100 ng

But this is only for the consensus sequence. Our real experiments will be a mix of molecules, some of which will be efficiently taken up and others that won’t. For a cursory estimate, we might assume that ~50% of fragments will be “good” USS and the other half will be “bad”. This would further reduce the yield.

That means, I am likely to need ~80 ml cultures and a starting input DNA amount of ~1600 ng, just to get back a mere 100 ng of DNA back!

Most Illumina sequencing centers seem to want ~1 ug of DNA to make libraries, but a lot of ChIP-seq experiments seem to call for only ~100 ng. In our case, there will be no downstream library construction, so we can likely get away with small amounts of DNA, as long as it is quite pure and accurately quantified.

Regardless, this is going to take fairly large cultures, fairly large amounts of DNA, and a good scaled-up periplasmic prep.

BUT, one important thing to note is that our degenerate oligo preparation will be more than sufficient for a large number of experiments, even at this large scale. For the controls, I can merely buy minimum-scale synthesis long oligos at ~$200 a pop. Since I can safely PCR amplify these, I will be able to make a replenishable stock for use in scale-up experiments.

More on this in the future, but while I’m doing this, I might as well estimate what it will take to get a microgram of chromosomal DNA fragments out of competent cell periplasms.

My previous experiments with sonicated DNA gave pretty consistent DNA uptake measurements:

~50% of 200 ng 1-10kb DNA / 1 ml cells → 100 ng max. yield.
~10% of 200 ng 0.2-0.4kb DNA / 1 ml cells → 20 ng max. yield.

Given a 25% recovery rate from the periplasm, this means that for a microgram of DNA, I will need:

1-10kb DNA: 8 micrograms in a 40 ml culture
0.2-0.4kb DNA: 40 micrograms in a 200 ml culture (!)

This last is really asking a lot. That size of scale-up will require special thought…

Appendix on Scale-up Issues:

Purity: I have not been adding RNase. I need to get all the RNA away, in order to accurately quantify the DNA. I am also concerned about salt. The CsCl in my DNA precipitates may not be getting washed out adequately by a single 80% ethanol wash.
Cell concentration: It would help for technical reasons, if I could concentrate the cells quite a bit before doing the organic extractions. I have used a ratio of 1:1, cells : organic solvents. So a 1 ml competent cell prep (~a billion cells) gets mixed with 1 ml solvent. But I might be able to resuspend 10 ml of cells in 1 ml and then use 1 ml solvent. I just don’t know.
DNA concentration: I want to make sure that I am saturating with DNA for my initial experiments, but I haven’t yet done a proper saturation curve to know what I should be using. This will decrease the total efficiency of DNA uptake, but my total yields will be higher, and I will be biasing things towards the best uptake sequences (which is a good place to start).
DNase: I have not been treating cells with DNase prior to isolation. From what I can tell, this is not a problem, and the free DNA is washed away. But if I use very high DNA concentrations, I will probably want to use DNase, just to be sure I’m eliminating free DNA completely.
Details, details: Scale-up is never quite as simple as just increasing the volume of everything. I will need to make sure that there are appropriate centrifuges, shakers, tubes, and everything else. Growth rates of cells and competence induction may be poor when going to larger volume cultures. I am also concerned about scaling up the organic extractions. It turns out that not all conicals are created equal; I’ve had disasters where the phenol has torn through the bottom of 50 ml conicals when doing large-scale organic extractions, depending on the brand of conical and rotor used. I’ll need to make sure that things like this don’t happen in advance before I mess up somebody else’s equipment!

Reverse engineering

2009-09-17T10:45:00.000-07:00

Okay, so now that we’ve exposited all the brilliant experiments we’re planning to do while writing proposals, the actual reality of doing the experiments is starting to sink in. We’ve also managed to put down some fairly concrete goals for the next several months.

One of our experiments involves measuring the specificity of DNA uptake by naturally competent H. influenzae for fragments containing “the genomic USS motif”. The H. influenzae genome contains an abundant sequence motif, and fragments bearing it are taken up better than fragments that don’t. This “uptake signal sequence” was originally defined by its functional role in DNA uptake, but has since been characterized mostly by bioinformatics, with no direct uptake specificity data. The limited data from previous lab members suggests only an imperfect correspondence between the properties of the genomic motif and the specificity of DNA uptake.

The idea, then, is to feed competent cells small DNA fragments bearing a degenerate (highly mutated) version of the USS consensus sequence, recover those that are preferentially taken up, and sequence the resulting pool. USSs are ~32 bases, well within the reach of single-end Illumina reads, if they are positioned properly next to a sequencing primer.

I’ve previously discussed the expected properties of a degenerate USS pool. And though I think we need to consider this more, I will focus this post on the design of other parts of the construct that will allow us to circumvent subsequent sequencing library construction steps. Illumina sequencing uses specific sequences added to the ends of molecules to capture and sequence DNA of interest...

Properties needed for a USS-containing construct, where Illumina sequencing can be directly performed to sequence the USS:

(1) SIZE: ≥200 base pairs, dsDNA. 200 base fragments with USS are efficiently taken up by cells, and the size is sufficient for efficient cluster synthesis and sequencing using Illumina’s Genetic Analyzer.
(2) CAPTURE SEQUENCES: One end of a strand of each fragment needs to be able to anneal to one of the two “Flow Cell Primers” (FP) in the Illumina flow cell, while the other end of the same molecule needs to contain the reverse complement of the other FP.
(3) SEQUENCING PRIMER BINDING SITE: The reverse complement of Illumina’s sequencing primer needs to be immediately downstream of the reverse complement of the USS. (This could work the other way, but getting the “sense” USS directly from the sequencing reads seems optimal).
(4) TAG SEQUENCE: The first few (four) bases of each read should be in non-degenerate fixed sequence to facilitate the alignment of the degenerate USS reads.
(5) CONSTRUCTION: After consulting several oligo makers, we learned that we wouldn’t be able to get our degenerate constructs built into an oligo longer than 130 nt. This means that I will need to anneal two oligos together and extend with polymerase to generate a full-length construct.

The first trick was to actually find out what the normal Illumina adapter and primer sequences were. They were available on-line, and I think I’ve mostly reverse-engineered what the different bits do. And think I have a reasonable design:

I’ll order two oligos, one 130 nt and the other 106 nt. (At the end of this post, I will list the exact sequences of each part and some notes.) They’ll have 36 bp of reverse complementarity at their 3’-ends, so that I can anneal them and extend to produce full-length construct.
To illustrate what all the different parts of the construct are for, here’s a color-coded version, for which I’ll schematically diagram the Illumina cluster synthesis and sequence priming.
The key features are that the flow cell primers (FP) are on opposite strands on opposite ends and the sequencing primer (SP) sits adjacent to the USS (with the 4-bp tag at the beginning). I am using plasmid sequence present in the lab’s other USS constructs for the Gaps (1 and 2).

To sequence the 200mer (either before or after recovery from competent cell periplasms), the DNA would be melted and annealed to an Illumina flow cell. Below are shown two different parts of a flow cell surface, where the two different strands of a single molecule might anneal.
DNA synthesis from FP1 or FP2 generates a covalently attached version of each strand.
The original molecule is melted off and washed out of the flow cell, and a special in situ PCR method generates clusters of single-strands covalently bound to the flow cell surface. In each cluster, the strands are oriented in both directions.
Sequencing then proceeds from SP binding sites. In this design, the SP binding site will then read the complement of the USS (with the first four fixed bases), so the actual sequence generated would be the USS contained in the construct.
There are several small details to go over to make sure that this design will work. Because the oligos are so expensive, and the degenerate oligo will be precious, I also plan to buy several non-degenerate oligos corresponding to perfect consensus, randomized, and mutant USSs. These will act as controls for the annealing/extension step that generates the uptake substrates and as controls for measuring saturation curves to optimize the appropriate DNA uptake conditions. I will also be able to do PCR to regenerate the control constructs, while I should probably avoid amplifying the degenerate USS construct for fears of strongly biasing the representation of different sequences.

NEXT UP: Uh oh… What about yields? Dimensional analysis…

APPENDIX:

The different parts of the two oligos:

Notes on my reverse engineering:

FP1 (25 nt): Composed of putative 20mer FP1 + first 5 bases of one adaptor (calling it A)
SP1 (33 nt): Sequencing primer for single-end Illumina runs. Includes the 13 bases of the normal adaptor that normally results in a 13 bp inverted repeat palindrome on either side of adapted DNA fragments.
USS (36 nt): Includes 4-base tag (ATGC) upstream of a 32-base genomic Gibbs consensus sequence with a set level of degeneracy at each position.
G1 (36 nt): Additional sequence from pGEM7f ,corresponding to the portion of the spacer region where the two oligos are intended to anneal.
G2 (46 nt): More sequence from pGEM7f, corresponding to the spacer region only on one of the two oligo.
FP2’ (23 nt): Composed of the complement to the 20mer FP2 + first 3 bases of the other adaptor (calling it B).
Total length after annealing and extension is 200 bases, where the USS is located from position 63 (after the spacer) to position 94. In the flow cell, the use of SP1 as a sequencing primer should read the complement of the USS sequence, so the actual sequence obtained will correspond to USS (with the first four bases always ATGC).

Eating chromosomal DNA fragments

2009-09-09T15:59:00.000-07:00

Haemophilus influenzae cells will take up closely related DNA from the environment quite efficiently, when they are made naturally competent by resource limitation.

Previously, I had done some experiments using sonicated chromosomal DNA of two different size distributions. The take-away lesson was that, for a fixed DNA concentration, larger fragments were taken up better than smaller fragments. This could be due to two non-exclusive reasons:

Larger fragments are more likely to contain an uptake signal sequence.
The uptake machinery is saturated when I used the smaller fragments, since there are more fragments per unit mass.

I am not certain of the best way to measure the relative contributions of these two factors to the observed disparity in uptake, though I’m pretty sure a saturation curve would be the way to start things off, that is measuring the amount of uptake over a wide range of DNA concentrations.

But first, and more to a practical concern for our sequencing plans...

I repeated this experiment, but also prepared total DNA and periplasmic DNA (by the slick method of Kahn et al) to make sure that I could cleanly recover chromosomal DNA fragments trapped in the periplasm from bulk chromosomes, as I previously showed for a small USS-containing PCR fragment.

Here are the results of that experiment (in which I provided ~0.5 billion competent cells with 200 ng of end-labeled DNA fragments of two different size distributions, either 1-10kb or 200-400 bp, for 30 minutes):
In (a), the % uptake is clearly better for the larger size distribution than the smaller size distribution. In (b) and (c), I show that I can purify periplasmic chromosomal fragments away from the cell’s chromosome. In (b), the results using the larger fragments is shown, while in (c) the results using the smaller fragments are shown. (I ran gels with two different agarose concentrations to optimize the separation for the two different input pools).

One thing to note is that the size-distribution of DNA between the input and periplasmic preparation were effectively indistinguishable. I looked at traces of these lanes in the Molecular Dynamics ImageQuant software, and they looked pretty much exactly the same. This is a little bit confusing, given the two models discussed above and the fact that fewer small fragments were taken up compared with larger fragments. I might have expected that there would be a bias towards the larger fragments in the periplasm compared to the input, but this was not the case.

Another thing to note is that, unlike when I previously did this experiment with USS-containing PCR fragments, there is still evidence of periplasmic DNA in wild type after 30 minutes. I don’t think this is due to poor washing of free DNA away from the cells, but rather reflects that there had been insufficient time to translocate all of the DNA in the periplasm into the cytosol. There is also the possibility that some of the non-chromosomal DNA in the wild-type samples are indeed cytosolic, which I can’t tell without some way to distinguish ssDNA and dsDNA.

Homoplasy versus Recombination

2009-08-28T16:07:00.000-07:00

Making up for lost blogging time!

Just to get started doing something with this alignment, I just grabbed a single 68,580 bp alignment block out of the XMFA file I generated (as described in the last post), which was itself a simple FastA formatted file. So I could use whatever programs could handle a simple FastA alignment and not worry about converting around to different formats yet, or coping with contig and rearrangement boundaries.

Here’s a picture from MAUVE of the aligned block I used (KW20 Rd reference coordinates 401,890 to 467,839). It’s the orange block shown below.
Based on some recent reading, I tried a couple different programs designed to detect signals of recombination from multiple alignment files. There is a serious caveat here, in that I do not understand exactly how these programs work, and since they were designed for use on sexual eukaryotes, it’s quite possible that the assumptions of the programs will not be accurate when applied to Haemophilus. Nevertheless, I sallied forth, just to see if they’d at least work with my sequences. I can come back and try to understand the meaning of the output more slowly...

LDhat

First, I tried LDhat, which is composed of a suite of programs aimed at detecting recombination in a population sample. I didn’t have any difficulty downloading and compiling the program, which was nice. I wish I could tell you more about it, but I never actually made it work, so there’s not much point yet. I got the “Convert” program to make the requisite files from an only slightly modified FastA, and I eventually got the “Pairwise” program to work. This program generates a “look-up table” which is required for the other programs to operate.

At first, the program kept choking on my data and was unable to produce the requisite look-up table for my data. So I fed it a pre-made file that can be found at LDhat’s website. The look-up table I used was first converted to a 15 taxa file using their “LKgen” utility. Then “Pairwise” worked and generated me a new look-up table that should theoretically have worked with the other programs. I next used the “Interval” program, but I could never get past here. It wouldn’t do the analysis, because the look-up table was somehow lacking. It was not “exhaustive”, whatever that means, and there was apparently a disparity between “npt” (which I assume stands for Nonparametric test) and the data.

No clue. Moving on. (I’ll try to return to this some other time, since the program “Rhomap” looks like what I’m looking for).

PhiTest

Next I turned to PhiTest, which was also easily downloaded and compiled. This program uses a principle that seems pretty sensible. It examines “incompatibilities” in phylogenetic signals. If two lineages of bacteria diverge and never recombine, then adjacent polymorphisms will most likely be “compatible”, that is both polymorphic sites will have the same phylogenetic signal (i.e. they will support the same tree topology). On the other hand, “incompatible” sites could have two possible histories, one in which there was recurrent mutation in different lineages and one in which there was recombination between lineages. This is illustrated nicely in a figure from the paper reporting PhiTest:

"Figure 1:
The dual nature of incompatibility. Two possible histories for a pair of incompatible sites are shown: (a) two incompatible sites explained by a recombination event and (b) two incompatible sites explained by a convergent mutation. Mutations in the first site are indicated by open circles and mutations in the second site are indicated by solid circles. To explain the incompatibility between the pair of sites either a recombination event must be invoked or a homoplasy must have occurred in the history of one of the sites."

Figure 1A:
Figure 1b:
So in (a) the black mutation only occurred once and was subsequently shared between the central lineages, whereas in (b) the black mutation occurred independently in the two lineages.

They more rigorously define “compatibility” thusly: “Two sites i and j are compatible if and only if there is a genealogical history that can be inferred parsimoniously that does not involve any recurrent or convergent mutations (known as homoplasies as in Figure 1b). If the two sites are not compatible, they are termed incompatible. Under an infinite-sites model (KIMURA 1969) of sequence evolution, the possibility of a homoplasy does not exist, and so incompatibility for a pair of sites implies that at least one recombination event must have occurred, as in Figure 1a.”

So the idea of recombination detection algorithms is to exploit the notion of incompatibility, but it is extremely important to note that there are really two possibilities for any incompatibility: (a) It is a true “homoplasy”, i.e. a recurrent mutation; (b) it indicates recombination. The basic idea of the NSS method (which is apparently what LDhat is using) is to look at neighboring pairs of sites and find regions where incompatibilities are clustered together, suggesting perhaps a hotspot.

The other, possibly more eukaryotic, aspect is that sites that are more distant will tend to be less compatible, since they’d more likely have had crossovers between them. Our provisional model for natural transformation would not necessarily have this feature, however, since we do not expect our recombination events to be crossover-like.

The authors of PhiTest produced a different measure than the NSS model that looked at the “Pairwise Homoplasy Index” for each pair of SNPs in an alignment. This gives a test statistic of the minimum number of homoplasies in the alignment, essentially based on a parsimony criterion. They could then apply permutation tests to measure the significance of a given PHI for an interval. Happily, the program will also output the significance under the NSS and MaxChi models.

For my alignment, PhiTest gave the probability of no recombination as 0 for all three models. So it looks like there’s quite a bit of homoplasy in my alignment, evidence of recombination. (It’s worth noting that PhiTest and LDhat both explicitly ignore gaps in the alignment.)

But that’s not good enough, what I really want is a plot of recombination rates across the alignment block. The figure I made from this post was the output of that analysis. It is only showing the ~9500 polymorphic sites and ignoring all the gaps and identical sequences. It’s showing, for each pair of SNPs, the probability of recombination (assuming that all homoplasies are due to recombination and not recurrent mutation).

It pretty much looks like there’s evidence of recombination almost everywhere! There’s just a couple areas that look like they have extensive compatibility, so it may actually be a lot easier for us to detect places where there’s been little recombination than places where there’s been lots.

I’ll need to go back to the raw data and figure out which set of SNPs are showing evidence of low recombination and see where those spots are in the alignment. They could be horizontally transferred segments with few USS?

I’m not sure what to make of this. It complicates things. To what extent can I trust these algorithms to find recombination and not recurrent mutation? The divergence between my sequences isn’t that high (only ~3%), so recurrent mutation doesn’t seem likely to be a huge problem, but still...

Anyways, I went ahead and used the other utility that comes with PhiTest, called Profile. This does PhiTest on sliding windows, providing a P-value for each segment. Here’s what it looked like, when I did a “Profile” of non-overlapping 500 bp chunks of my alignment:
The x-axis is the position in the alignment, and the y-axis is the p-value of the PhiTest for that particular 500 bp chunk. The line indicates a P-value threshold of 0.05, so any segment below that line has evidence of recombination, based on PHI. Segments that have higher P-values may either have compatible SNPs, or for the very high P-values, be sites with very little SNPs.

This seems promising, if we can trust that we’re measuring something meaningful. I have a strange feeling that using homoplasy to measure recombination is going to be fraught with difficulties, but it’s pretty much the best principle that we can work with to detect recombination in these population-level samples. Luckily, both PhiTest and the other methods are sophisticated enough that they are evaluating the homoplasy index for a given pair of SNPs relative to other surrounding SNPs, so when incompatibilities are detected, they are standing out from other possible pairings. This should help with the problem of true homoplasy, but then again, maybe it doesn't. I'll need to talk to my friend Corbin about this stuff.

Whew! Okay, enough on this for the moment...

More alignments!

2009-08-28T12:27:00.000-07:00

I realized something disturbing a couple days ago, when I crashed my computer by running out of hard-drive space... The Haemophilus influenzae genome is actually replicating inside my computer! Not just on my lab bench! Several of the programs I've been using (namely the alignment programs) have been spitting out files containing the complete genome, and I think I must've had several hundred before I got rid of a bunch of intermediary files. In silico infection! A scourge!

Anyways, I guess that last entry deserves more explanation, but it'll take me more than a post to do it. First off, the motivation for doing these kinds of analyses:

There’s 4 complete Haemophilus influenzae genomes, 11 with draft assemblies (in 4-50 contigs), and 17 for which sequencing is in progress. On top of that, hundreds of strains have been sequenced at seven house-keeping genes (MLST studies).

Within this “population genomic” data, we should be able to do some estimates of recombination rates between isolates, as well as some sophisticated analyses of genetic variation at USS motifs. But how? Here’s the basic idea:

(1) Align the complete and draft genomes.
(2) Use programs written to detect recombination signals in population genetic data.
(3) Examine genetic variation in and around USS motifs (scored for information content across each genome) and see if they correlate with the determined recombination signals.

There are several issues with this at each step, but I’ll take the first one first. I’ll start with a post describing how I produced the multiple alignment file that I fed to the PhiTest program...

The first trick was downloading the files...
I already had the four complete sequences, but hadn’t had the easiest time getting the contigs of the draft assemblies. Do they really want me to download each contig by clicking one at a time? My problem was solved when I found the whole-genome shotgun (WGS) FTP site at NCBI. Of course, these files were named by their WGS code, so I simply hand-search for the correct four letter code, based on this list and downloaded the corresponding file. If I’d been more clever, I would’ve done this from the command line with a script containing the list of WGS data I wanted.

I used the search term "txid727[orgn]" in the Genome Database or the Taxonomy Browser to call up Haemphilus influenzae genomes. Using the Genome database yielded some spare entries that were plasmid sequences, and using Taxonomy browser yielded many isolates where the sequencing is still in progress. Anyways, I could get to a Genome page for each sequenced isolate, but when I clocked to get the contigs, it provided separate accessions for each contig. Boo! There had to be a better way, which took some digging, but was here. This is all the WGS data at NCBI. Whew! Luckily, I found the requisite four-letter codes from my other searches and downloaded the appropriate *.gbff.gz files and expanded them with gunzip.

Anyway, this gave me a bunch of GenBank files with all the contigs in each one. The contigs themselves were simply ordered by size, with the largest first and the smallest last, and the coordinate system corresponded to this ordering. While this is fine, it would probably help me looking at alignments, if I ordered the contigs based on the reference KW20 Rd genome. It would be a meaningless gesture in some ways, but would certainly facilitate looking at pictures.

I did this by using MAUVE’s MoveContig function. I reordered each set of WGS contigs pairwise against the reference KW20 Rd. (Again, I should've been using the command-line, but found it easier to just click my way through... shameful.) Here’s an example of what MoveContigs did comparing PittII to Kw20 Rd:

Alignment #1: Keeps the original order of the contigs (from largest to smallest).
KW20 Rd is on top; PittII is on bottom. The colored blocks indicate contiguous aligned blocks. The red lines indicate contig boundaries (KW20 has only two, for the start and stop of the chromosome). The diagonal line between the two genomes show which blocks belong together. It’s quite obvious from all the crosses that the order of the contigs is off.

The program then reordered the contigs and did another alignment. Then it did this a third time (it keeps iterating until it’s done). It decides this is as good as it can do and stops.

Alignment #3: Reordered contigs maximizing “synteny” between the reference and PittII.
Clearly reordering the contigs has made the picture clearer to see. There are far fewer alignment blocks and fewer diagonally crossing central connector lines. Of course, some of the contigs did span some breakpoints, so there’s obviously rearrangements between the genomes. At each contig boundary, there is always the possibility of some rearrangement between the strains, so this contig order is at best provisional. Another thing I had to keep in mind while tooling around with these alignment programs is that the genome is circular, but the programs don’t treat the strings that way, so there’s often things at the right or left ends that obviously are still syntenic, but just don’t look like it with these linear visualizations.

Here’s the alignments again, without the colored blocks, just so that the reordering of the contigs is quite obvious:

Alignment #1:
Alignment #3:
Notably, the little tiny contigs that failed to find something in the reference to align with remained at the right end.

Alright, so going through 11 of these MoveContig steps provided me with 11 new files that had the contigs in a different order (based on KW20 Rd) than what I got from GenBank. The only unfortunate thing about this was that MoveContigs only outputted FastA, instead of GenBank, so I lost all the annotations. Oh well. I would still be able to use the KW20 Rd annotations. I think that the annotations can still be added back, but I haven’t figured that out yet.

To show that this exercise wasn’t a waste, below I show the stripped-down alignments when I did the big 15-way multiple alignments (I hid the 3 other completely sequenced strains from view, and only the top several contig reorderings are visible.

Before Contig Reordering:
After Contig Reordering:
Much easier to look at... But still sort of meaningless.

Alright! So that’s a 15-way multiple genome alignment. And indeed, it’s associated with a gigantic alignment file in the XMFA format: an extension of the FastA format that allows for a description of alignments that have rearrangements / contig breakponts / etc between them.

The most difficult challenges:
(1) Handling the breakpoints, contig boundaries and real rearrangement breakpoints alike.
(2) Handling indels. Those within an aligned block should be scorable, but how?

So really, the first thing I can work with without much difficulty are contiguous aligned blocks. I grabbed one of these from the XMFA file to produce a single FastA file, which I called test1.fasta, which I planned to use for some preliminary analyses...

TO BE CONTINUED

Measuring recombination?

2009-08-26T22:18:00.000-07:00

Whew! I've been slacking on the blog posts! Above: Probability of recombination between pairwise SNPs in a 15 strain alignment block of ~65 kb. Yellow indicates low probability of recombination...

The figure above was generated by PhiTest, a program that tries to identify regions of recombination between different taxa in a multiple alignment. Last night, I used MAUVE to generate a 15-way multiple alignment between the genome sequences of 15 different Haemophilus infleunzae isolates. (This took a while, but still was impressively quick.)

I grabbed a reasonably large alignment block (65,480 bp, including gaps) and stuck it into Phitest. It did everything reasonably quickly and gave me the unshocking result that there was a p = 0E+00 chance of no recombination occurring on the interval. So I decided to use the -g switch to provide a graphical output giving the chance of recombination between every pair-wise combination of SNPs (of which there were 9422; the program ignored indels).

This was fun, because it wrote such a gigantic file that my computer died due to lack of hard-drive space! That doesn't happen often! But it was mainly because I have several dozen gigabytes of stuff on my laptop that doesn't belong. I nudged a fraction of this aside and repeated the analysis to give the plot above.

Uptake and Transformation with "Biorupted" samples

2009-08-13T12:15:00.000-07:00

To get an idea of how uptake and transformation would work with different sized donor chromosome fragments, I took MAP7 DNA and sheared it in a “bioruptor”. I ended up with several samples with different size distributions, three of which I used for a pair of experiments:

LARGE: >40 kb (unsonicated)
MEDIUM: 1-10kb (1 X 10 min sonication)
SMALL: 100-400 bp (5 X 10 min sonication)

My naïve assumption was that % uptake would go down as the fragment size decreased, since fewer fragments would contain “uptake signal sequences” (USS), which have an average density in the genome of ~1kb.

I also thought that transformation rates would also go down for smaller fragments, but would not necessarily correlate that well with uptake, since additional steps of translocation and recombination could also potentially influence the efficiency of transformation. So for example, transformation might drop off more quickly than uptake, if degradation would affect smaller fragments more than larger fragments. (This seemed to be the case in Pifer and Smith, 1985.)

Keeping in mind that these are just one-off experiments and need to be repeated (like pretty much every experiment I’ve reported in this blog), the above predictions look like they’re true, but I’m not certain if my reasons are necessarily correct...

First, I’ll show the uptake data. I end-labeled the MEDIUM and SMALL donor fragments using Klenow and did a simple uptake experiment (comparing total radiolabel to that in cell pellets after 30 min of uptake) using wild type cells and as donors, either MEDIUM or SMALL chromosome fragments, and either saturating (500 ng) or sub-saturating (100 ng) amounts of input DNA per 0.5 ml of wild-type competent cells. Here’s the data:
So clearly, several-fold less chromosomal DNA is taken up when smaller fragments (100-400 bp) are used than when larger fragments (1kb-10kb) are used. As explained above, one reasonable explanation for this is that fewer fragments in the SMALL sample contain USS, so maybe only 1/10 to 1/2 will have a USS, whereas in the MEDIUM sample most fragments will contain at least 1 USS.

But there is an alternative explanation that I’d like to be able to distinguish (but am pretty sure I can’t with this one experiment). It could be that a given competent cell only takes up a fixed number of DNA fragments, independent of fragment size. So since the SMALL sample is composed of ~10-100X more fragments per unit mass, it could be that I’ve simply saturated the system with fragments in the case of SMALL, but not in the case of MEDIUM. This issue was addressed by Deich and Smith, 1980, and they concluded that indeed this was the case (that the number of molecules taken up was independent of fragment size), but while they do mention USS, they do not bring up USS density as a potential reason for their data.

I’d hoped that by doing a second DNA concentration (100 ng) that I might get hints as to which of the two above models is correct (or if they are both correct and both contribute to the observation), but I don’t really think I can say too much without repeating this several times and getting some error bars on that graph. Furthermore, I’m not really entirely sure what the expectations are for the two models. I’ll have to think on this some more....
-----
Okay, what about transformation rates using my “biorupted” fragments? Below are two graphs reporting the transformation rates of the KanR and NovR alleles from MAP7 to KW20 for the three different DNA pools (LARGE, MEDIUM, and SMALL):
(Note: In the case of the NovR/CFU SMALL sample, the number reported is actually the limit of detection, so NovR/CFU(small) is less than 4.4e-6. I didn't observe any NovR transformants for the SMALL fragments, despite having a decent limit of detection.)

First, I'll look at the difference between the LARGE and MEDIUM fragments. Both markers showed ~6-fold decrease in transformation rate in the briefly sonicated sample, compared to the large intact fragments. Possible explanations:

Less USS per fragment: I doubt this is a major reason for the difference. While a larger percentage of fragments are expected to contain no USS in the MEDIUM sample, it shouldn’t be that large of a difference, since the mean density of USS motifs is ~1kb.
Degradation by translocation or cytosolic nucleases: This seems to be a reasonable explanation. From an old set of experiments using a defined plasmid donor, Pifer and Smith, 1985 estimated that an average of ~1.5 kb of a leading 3’ end is degraded during translocation. Maybe the medium-sized fragments simply don’t survive translocation as well as large taken up fragments.
Recombination efficiency: Maybe both the LARGE and MEDIUM fragments make it into the cytosol, but homology search and recombination are much better for larger fragments.

Things look a little more interesting when looking at the change between MEDIUM and SMALL fragments: While the KanR rate only changed modestly (less than 2-fold), the NovR rate went down below my limit of detection. Possible explanations:

Distance to USS: The nearest USS to the NovR allele is more than 4oo bp away, but less than 400 bp for the KanR allele: I like this explanation. I really really need to figure out the identity of these antibiotic resistance alleles. We’re pretty sure it’s a mutation in the gyrB gene, but I don’t know the actual change. I looked at gyrB and it does contain a USS core motif and two other core motifs a few hundred bases before the start codon, but the gene is ~2.5 kb, so the actual gyrB mutation could easily be too far away from these USS. When I looked at the putative gene responsible for KanR (the ribosomal S7 gene), there was a single USS near the start, but none within. Again, without knowing the causative lesion, I can’t tell whether this is within the size distribution of the SMALL fragments.
Differences in degradation rates at the two loci: This is possible. The other thing these data suggest is that, despite the ~1.5 kb average degradation reported by Pifer and Smith, 1985, there’s still plenty of small fragments that can recombine, since the KanR rates between MEDIUM and SMALL are not really dramatically different.
Recombination signals: Also possible. And probably the hardest to tell, since the other effects need to be canceled out.

The main take-away (assuming that the results replicate) is that not only do different markers transform at different rates, the change in transformation rate for different sized-fragments also varies for different markers. The underlying reasons for this are probably interesting. So what next? I need to repeat this, but next time I’d like to:

Extend this to additional markers to see if this variability also applies to other loci.
Measure linkage between Kan and Nov. I’ve previously seen the known linkage between Kan and Nov using large fragments, but would expect linkage to vanish for small fragments when the KanR and NovR alleles never share the same fragment.

And again, I really need to know what the lesions are that are responsible for the MAP7 antibiotic resistances. I’ve looked around a fair amount, but it seems that many antibiotic resistances can be produced by mutations in more than one different gene, so narrowing it down isn’t that straightforward.

Back to the grind

2009-08-10T11:52:00.001-07:00

Whew! Getting that albatross (= grant application) off my neck feels good, but now I need to formulate a plan for the next several weeks of lab work. What are my priorities? There’s quite a lot of things I could be doing, so I might as well make a list...

(A) Generate purified DNA from 86-028NP, PittEE, and PittGG donor DNA: Should be pretty easy. I tried this a couple times in the past couple weeks, but had a few problems:

Got a goopy mess. I don’t know if it was my lysis or something about my extractions, but when I would try to pull the aqueous layers off of my phenol/chloroform extractions, this goopy junk at the interface kept interfering with my ability to complete the preps. I will try try again, but this next time I will do a few things differently: (i) start with fewer cells per volume, (ii) do a straight phenol extraction first, (iii) let lysis proceed longer.
PittEE didn’t grow. I’d gone to the lab stocks, streaked for single colonies, grown up overnight cultures from single colonies, and stored them in glycerol in the -80. But when I returned to PittEE frozen stocks, they don’t grow up! Weird. Maybe PittEE is particularly bad at surviving outside log-phase, so I had merely frozen down a pile of dead cells. I’ll have to return to the lab stocks and restreak for single colonies.
I never checked the putative nalidixic acid resistant phenotype of 86-028NP. I really hope that marker will work to select for transformants from 86-028NP to KW20, as this would really help several steps of the project. In our strain database, 86-028NP is listed as NalR, but this might only be a clinically relevant phenotype. I need to check it on plates.

(B) Check the uptake and transformation phenotypes of WT, rec-2, rec-1, and pilB: I’ve already gone through this exercise with wt and rec-2, but as a negative control, pilB mutants should neither take up DNA nor be transformed by it, and for my translocation experiments, rec-1 should take up DNA normally but fail to transform. I have all these strains frozen down as M-IV competent cells, so I just need to do the experiments to verify that my mutants have the expected phenotypes. The rec-1 mutant in particular is important, since I plan to use it for cytosolic DNA purifications, but it could be somewhat leaky, something I need to know. (The only rec-1 mutant available is the original isolated mutant from back in 1972. It’s older than me!

(C) Quantify uptake and transformation for different-sized donor DNA: A coupe weeks ago, I used a “Bioruptor” to sonicate MAP7 DNA to different size distributions, and I tested the transformation rate of antibiotic resistance into competent cells of some of these different sized DNA. I also want to check how well they’re taken up. Once I’ve done this, I’ll report back the full set of results. The expectation is that while even small fragments will be taken up well, only large fragments will transform well. There’s also an interesting expected relationship between the amount of uptake expected and the size distribution of the DNA. If the DNA fragments are quite small, then only a fraction are expected to contain uptake signals, but many fragments containing such signals will be taken up. On the other side, large DNA fragments will mostly all contain uptake signals, but only a few molecules will be taken up, but these will be larger... Hmmm....

(D) Design the degenerate USS construct. This will require its own post, but suffice it to say, we have a pretty good plan on how to design this construct so that we will be able to directly sequence the degenerate USS motif by short single-end Illumina sequencing without any processing steps upfront.

(E) Transformation experiment: If indeed 86-028NP can grow on nalidixic acid plates, then I will just go ahead and do half of the experimental-side of one of my specific aims. All I need to do is transform KW20 competent cells (which I have) with DNA taken from 86-028NP (which I’ll make), select for colonies that are NalR, pick several of them, grow them up, and extract their DNA. This would provide up with plenty of material to produce sequencing libraries and then send to sequencing. Then we can really get this project off the ground.

(F) Find somewhere that will do sequencing with a short turn-around: The local Illumina sequencers (at UBC’s Genome Center) offer what looks like excellent sequencing services for a reasonable cost, but their turn-around time is a bit too slow for us to hope to have any preliminary data for Rosie’s grant applications. If we want to get anything done more quickly, we will have to find someone else to get started.

(G) Design and produce a large DNA fragment from the 86-028NP genome to use for developing a cytosolic ssDNA prep: Since the most difficult experimental part of our plans is likely to be purifying (and cloning) ssDNA from the cytosol, I’d like to start with a defined construct I can use for purification experiments. I think that a good choice would be a large clone from 86-028NP, because this might be useful for other experiments as well. I will want something at least several kb on a plasmid. I will need to think about whether or not it matters what this fragment contains. As a first guess, I would want it to contain the putative NalR allele and several kb of flank.

What else should I be doing?

Almost...

2009-08-06T08:45:00.000-07:00

...there... However it turns out, I'll finally have this grant I've been working on out of my hands tomorrow. I'm actually thrilled to get back into the lab when it's all over...

Food-to-sex ratio?

2009-07-27T11:53:00.001-07:00

In my recent experiments using radiolabeled USS-1 donor DNAs, I was impressed by just how well naturally competent H. influenzae will slurp up USS-containing DNA. I’d read about it, but observing it myself was really something. (It reminds me a bit of the first time I “saw” Mendelian segregation when I dissected my first yeast tetrads.)

The other thing that struck me, reading the antecedents of my uptake experiments, was that a large majority of taken up DNA is simply degraded and the nucleotide subunits used for DNA replication. In my experiments, while donor DNA remained intact in rec-2 cells, there was no intact donor DNA after an hour in wild-type cells. All the radiolabel was found in the chromosomal DNA.

Using linearized USS-containing plasmid donor DNA, Barouki and Smith (1985) nicely show that this chromosomal labeling is NOT dependent on recombination, as rec-1 mutant competent cells obtain similar levels of chromosomal labeling (rec-1 is the recA homolog in H. influenzae). Using restriction digestion, they also nicely show that in wild-type cells, only some of the donor DNA manages to recombine into the recipient chromosome.

The most important lanes in the above autoradiogram from Figure 3 of their paper are Lane E and Lane F, which show restriction-digested DNA after uptake in wild-type and rec-1. In Lane E (wild-type), the two arrowheads indicate the chromosomal restriction fragments indicative of transformation by the linear plasmid donor. In Lane F (rec-1), there is no appreciable transformation of the donor DNA into the chromosome (those restriction fragments are gone).

These results raise an interesting prospect: Could we interpret amount of recombination-independent radiolabeling relative to the recombination-dependent radiolabeling as a “food-to-sex” ratio? Undeniably, DNA is taken up and used by competent cells, but it’s clearly used in two different ways: Subunit recycling (food) and recombination (sex). By my eye, it seems like a lot more of the donor DNA, even in wild-type cells, is used for food than for sex.

Of course the amount of taken up DNA a cell could use for “sex” would be highly dependent on what DNA was taken up. In the case of the experiment above if homologous portions of the donor plasmid are removed, then the “sex” fragments disappear (Lanes I and J for wt and rec-1, respectively). Furthermore, the length of homologous DNA fragments taken up by competent cells is likely to matter, due to degradation during translocation and in the cytosol.

In Maughan and Redfield (2009), they show extensive natural variation among H. influenzae strains in the amount of uptake and transformation that competent cells will undertake. Do strains that can take up DNA well but fail to transform have a high food to sex ratio?

The assay that Barouki and Smith use doesn’t seem like the best way to measure a food-to-sex ratio, since the “sex” signal is somewhat buried behind the “food” signal. I wonder if there’s an experimental scheme that would allow one to measure such a ratio more accurately. Is there a way I could feed one strain’s chromosomes to another recipient strain and figure out how much DNA incorporated into chromosomes is recombination-dependent and how much is recombination-independent?

Anyway, thinking about this has led me to having a slightly clearer idea about the food versus sex hypotheses for the maintenance of natural competence by natural selection. Things are rarely black-and-white, so perhaps both models have their merits, but it seems like it might be possible to experimentally measure how much naturally competent cells use DNA for food or sex. Would this help in understanding these arguments?

Pile-ups

2009-07-24T16:31:00.000-07:00

As I try and address the critiques of my original NIH postdoc fellowship application in my resubmission application, I started fleshing out ways I which our planned periplasmic donor DNA sequencing experiments will investigate the mechanism of DNA uptake.

I’m playing around with different figures that illustrate what I’ll do and what it might reveal. Here’s the basic idea:

(1) Incubate sheared chromosomal donor DNA with recipient competent cell preparations.
(2) Recover the DNA that is taken up in to the periplasm.
(3) Obtain paired-end sequence data for periplasmic and input DNA libraries.
(4) Compare the abundance of different sequences between the periplasmic and input DNA libraries to calculate the periplasmic uptake efficiency for sequences across the genome.

What would this data look like?

One mapped back to the genome, each paired-end read will define the span of an individual fragment from the sequenced pool.

Below, I illustrate what a few dozen of these spans would look like mapped to a short stretch of chromosome containing an uptake signal sequence in the center.

The first diagram shows what the input DNA pool would look like. The blue spans do not contain USS, while the red spans do.
The second diagram shows what a similar amount of sequencing of the periplasmic uptake DNA pool might look like. I assume that the presence of a USS motif is sufficient and necessary to strongly stimulate uptake (“all-or-nothing” model of uptake). Thus, spans containing USS would be much more abundant than spans that do not.
These sequence data could be plotted in several ways:

(1) Spanning coverage at each genomic position: The input DNA is expected to have roughly equal spanning coverage of eaach genomic position, but spanning coverage in the periplasmic uptake library is expected to be higher closer to USS motifs. As the distance between a genomic position and USS increases, fewer spans will contain both. Peaks will indicate USS, and peak height will indicate the effectiveness of individual USS loci. I estimate that one lane of sequencing the input will provide ~2500X spanning coverage per nucleotide for 500 bp donor DNA fragments.

(2) End coverage at each genomic position: If the position of USS in a fragment is irrelevant to the uptake mechanism, then plotting end coverage would have a different shape than the spanning coverage around. Since any fragment containing USS will be effectively taken up, I expect a more sawtooth-shaped distribution of end reads at USS determined by the spanning fragment length. I estimate that one lane of sequencing the input will provide 100X end coverage per nucleotide for 500 bp donor DNA fragments.

Below, I illustrate an idealized case of spanning and end coverage.
Different USS loci may behave in different fashions. Here are what two other scenarios might look like:

(1) If uptake is polarized by USS, such that the position of USS on a fragment is important, or if there was an uptake blocking sequence nearby a USS, the distribution might be skewed:
(2) If a fragment’s uptake is equally efficient with one or two USS motifs (a USS interference model), then coverage around two nearby USS might look like this:
I wonder what other mechanistic details might be found in the data...

Dose Response

2009-07-17T14:54:00.000-07:00

How hungry are competent cells for DNA? I know that about a billion cells will consume ~65% of 20 nanograms tasty USS-1 fragment, but what if I offer the cells different amounts of USS-1?

To get a better hands-on feel for the DNA uptake process in wild-type and rec-2 mutant competent cells, I did a dose response experiment, where I incubated competent cells with different amounts of USS-1 DNA.

For this first experiment, I used 0.5 ml of competent cell cultures for each sample and did 6 different amounts of USS-1 DNA (12 samples total for wt and rec-2). I didn’t have enough radiolabeled fragment for all of my desired concentration, so I mixed in some cold USS-1 DNA to make up the difference. I let the DNA and cells incubate for 30 mins, then I washed the cells several times and determined the total radioactive counts in the cell pellet and washes to determine the % uptake and total uptake.

Here’s the results:

Total Uptake:

Percent uptake:

Interestingly, rec-2 does better at low concentrations of DNA than wild-type, but worse with high concentrations. The latter could be due to the periplasm getting too clogged with DNA, such that the outer membrane uptake machinery has to work too hard to get more DNA through, while in wild-type translocation of DNA frees up space in the periplasm. But the former (higher uptake in rec-2 at low DNA concentrations) doesn't really make much sense to me. Maybe not all free nucleotides created during degradation at the inner membrane remain in the cell, so that at low concentrations, rec-2 simply holds more label?

Standing upon the shoulders of giants

2009-07-12T12:33:00.001-07:00

It really is gratifying to have things work the way they're supposed to. Some kind of bug bit me on Saturday and I came in to see if the periplasmic DNA preparation reported by Kahn et al 1983 would work in my hands. And sure enough it did!

The experiment was much the same as before. I added radiolabeled USS-1 fragments to either wild-type or rec-2 competent cell preps, incubated for 5 minutes, and then either prepared total DNA or did the periplasmic extraction (TE/1.5M CsCl + phenol/acetone, 1:1).

Since wild-type cells will take up the fragment, but also incorporate labeled subunits from degradation of taken up DNA, I can tell if the periplasmic DNA prep managed to exclude chromosomal DNA. But first, I counted the radiolabel present in the different cellular fractions...

This time, about a quarter of the USS-1 fragment added was taken up within five minutes (wild-type: 26%; rec-2: 28%). I suspect these numbers are lower than the last time I did it, because my five minutes was really five minutes (whereas the first time, I think I was 2-3 minutes late).

The extraction: When I collected the aqueous phase, I also collected the organic phase, and the interface between the phases (which should contain the cells minus their outer membranes). I counted the radiolabel in these different fractions as before:
Wild-type cells had label in both the aqueous extract, as well as in the interface containing the cells, while rec-2 had nearly all the label in the aqueous extract. The organic phase had less than 1% rec-2.

But here's the important bit:
Lane 1: Input (1/3, or 4 ng)
Lane 2: Total DNA, wild-type + USS-1 for 5 min.
Lane 3: Total DNA, rec-2 +USS-1 for 5 min.
Lane 4: Peri DNA, wild-type + USS-1 for 5 min.
Lane 5: Peri DNA, rec-2 +USS-1 for 5 min.

The important point here is that in the total DNA extract of wild-type, both intact donor USS-1 and chromosomal labeling are evident, while in the peri-extract of wild-type, there is no chromosomal label.

This means that the extraction I did successfully purified periplasmic DNA over chromosomal DNA. Fabulous!

Now I need to scale this protocol up, and get cleaner DNA (i.e. use RNase), so hopfully I can see this without using radiolabel. If I can really get clean periplasmic DNA with little or no chromosomal contamination, I will move onto doing the "real" experiment with donor DNA made up from sheared genomic DNA of another isolate.

Yay!

Building a periplasm prep...

2009-07-10T14:52:00.000-07:00

After my failed attempts at doing a large-scale periplasm prep right off the bat, I decided to spend this week going a bit more slowly. I repeated what others have already done successfully using radio-labeled DNA fragments as donors. This means that I can do smaller scale experiments and don't need particularly pure DNA.

And this time the experiments all worked. Here's what I did:

(1) I made competent cells of KW20 (RR722) and KW20 rec-2 (RR622). I confirmed that the wild-type strain transformed normally and the rec-2 strain not at all (or at least below my limit of detection). This confirmed that my competent cell preps were okay, and that the rec-2 strain seems to be correct.

(2) Following the DNA uptake assay protocol of Maughan and Redfield, 2009, I showed that USS-1 is taken up very well, but USS-R is only poorly taken up. To do this, I simply incubated ~12 ng of either radio-labeled USS-1 or USS-R with 0.5 ml of competent cells for 20 minutes, and then compared the radioactive counts in a washed cell pellet compared to the total counts:
Wild type and rec-2 both take up USS-1 well, but USS-R poorly, as expected. But rec-2 seems to take up USS-1 slightly better than wild type. This is also true in the next experiment. This may be significant but could also reflect slight differences in the competent cell prep of the two strains.

Possibly the coolest part of this for me was that I got numbers that were spot-on the former post-doc's numbers (found in her notebook) and older papers describing % uptake. That is: ~65% uptake for ~20 ng / ml of cells. This was very encouraging to me.

(3) I repeated the uptake assay described above using ~12 ng USS-1 donor DNA and incubated wild-type and rec-2 cells for either 5 or 60 minutes. This gave results similar to those shown above:
Most uptake was finished after only 5 mins, though additional incubation increased the level of uptake. The results were nearly identical for 60 min incubation as for 20 min incubation, so I don't need to do it for so long.. The rec-2 strain again showed slightly more uptake at all time points.

After this, I took it a step further: I also extracted the DNA from the cell pellets and ran them out on a gel. I also included the input donor DNA as a control. I dried down the gel and exposed it to a phosphor screen. This is what the gel looked like:

Lane 1: Donor DNA (50% of input; ~6 ng)
Lane 2: Wild-type + USS-1 for 5 min. Total DNA.
Lane 3: Wild-type + USS-1 for 60 min. Total DNA.
Lane 4: rec-2 + USS-1 for 5 min. Total DNA.
Lane 5: rec-2 + USS-1 for 60 min. Total DNA.

Alright! That's exactly what I hoped for! (Well, not quantitatively between lanes: this was a sloppy first experiment.) The gel shows that the natural competence phenotypes of the two strains: wild-type and rec-2.

Intact uptake DNA is the smaller (lower) band, while chromosomal DNA is the high molecular weight species. In wild type, donor DNA gets degraded and nucleotides can be incorporated into the chromosome over time. (Importantly, the labeling of the chromosome is NOT from transformation, but from incorporation of degraded nucleotides into the genome by DNA replication.) In rec-2, the radio-labeled donor DNA is trapped in the periplasm and isn't degraded. So there is no chromosomal labeling in this case.

This is effectively a repeat of an experiment from Barouki and Smith, 1985.

Next, I'll try exactly the same thing, but I'll also try the extraction from Kahn et al., 1983. If this successfully yields pure periplasmic DNA, then I expect that the extraction will not yield radiolabeled chromosome, even for the wild type sample. If that works, I can work on scaling the protocol up to do a real purification of uptake DNA.

Onward!

Imagine it exists, and maybe it does!

2009-07-07T12:43:00.000-07:00

Our proposed experiments involve capturing DNA molecules at the different stages of natural transformation. One of the technical challenges we face will be producing a library of DNA molecules that have been translocated into the cytosol. We have some schemes for how we’ll do the purification of donor DNA from the cytosol, but even assuming that this works wonderfully, we still need to turn these into double-stranded DNA. We can’t use a specific primer to the 3’ends of translocated ssDNAs, because (a) we don’t know the exact 3’ ends and (b) it will be a complex mixture.

What to do?

Until now the only thing that had occurred to me is to use random priming of our ssDNA to convert cytosolic ssDNA into dsDNA (shown schematically above), but this approach has several limitations. The biggest problem is that we would only be able to accurately identify the 5’-end of translocated DNA. The 3’-end of the final dsDNA we produce would not represent the 3’-ends of the original ssDNA molecules using random primers. Furthermore, we would not know which end of our dsDNA was the original ssDNA’s 5’ or 3’ end. And finally, we would end up with a highly heterogeneous size distribution, which might complicate sequencing.

How can I circumvent this, get both ends, and know which is which? I need a strategy like RACE. I decided to imagine that a certain enzyme existed that might help me in this endeavor and then see if it actually existed and was already was commercially available. This strategy has worked for me in the past: Once, I’d wanted to know if there were restriction enzymes that only nicked at their recognition sites, so I typed “nickase neb” into Google, and sure enough NEB carries nickases! Go biotechnologists!

This time, I want to tack some type of single-stranded adaptor sequence onto the 3’ ends of my putative cytosolic ssDNAs, so I typed “ssDNA ligase” into Google, and Presto!... Epicentre produces a single-stranded ligase that they call CircLigase. Sweet!

This doesn’t fully solve the problem, since the ligase will normally take an ssDNA and circularize it (since the intra-molecular ligation will usually be favored). This is useful to plenty of folks who are interested in doing rolling circle amplification and rolling circle transcription, but I would rather not circularize my ssDNAs, but would like to favor ligation of an ssDNA adaptor specifically to the 3’-end. This will require a couple of bells and whistles.

If we take our ssDNA and then treat it with a phosphatase, we can rid the 5’-end of its terminal phosphate and both block circular ligation, as well as ligation of our adaptor to the 5’end. If our adaptor oligonucleotide also has a protected 3’end (not sure how to do this... an oligo with a terminal dideoxy nucleotide?), then we’d block the ligation of the adaptors to each other and force ligation only in the orientation we want (5’ of the adaptor to 3’ of the target).

Then, using a primer complementary to the adaptor, we can convert full-length ssDNA into dsDNA. Furthermore, the adaptor marks the original 3’ end of the fragment, so we can give a polarity to our cytosolic fragments. Here’s the scheme:

Afterwards, of course, we’d need to either amplify this product or de-protect both ends, so that we could ligate sequencing adaptors to the mixture.

This plan just might work, and I could make sure it works using defined substrates, rather than precious (as well as non-existent) cytosolic DNA fractions. The main thing I can’t think of off the top of my head is getting a hold of an oligo with a protected 3’-end (preferably reversibly so).

UPDATE: Looks like at least some oligo companies can include dideoxy bases in oligos. Awesome. This is not quite as ideal as a reversible protection of the 3' end...

UPDATE 2: Uh oh. How to amplify the product? There's no primer sequence at one end... This could involve a phosphorylation step and a ligation of a normal adaptor to that end? That's an extra unfortunate step. Also, the adaptor sequences will eat into the sequence read length, but that should be acceptable, since we only need to get tag sequences.

Mutation versus Transformation: Structural Variation

2009-07-02T18:19:00.000-07:00

Here's an idea: Compare structural changes (indels and rearrangements) in transformed and untransformed cultures.

One major goal of our planned experiments is to measure the transformation rate across the genome. This is a fairly ambitious prospect, mainly because of the amount of sequencing it will require.

If we assume that the average donor allele transforms recipient chromosomes 1% of the time, then we would need to sequence the average locus 100 times to see the donor allele just one time. But to get a reliable measurement of transformation rate would require considerably more... perhaps 10,000 times. This would give us an average of 100 donor alleles / 10,000 alleles sequenced. Even using the Illumina platform would get fairly expensive to measure the transformation rate for every SNP.

One way we could do a preliminary sequencing experiment would be to ignore single-nucleotide differences between donor and recipient and focus solely on indel and rearrangement differences (a.k.a. structural variation). This data would speak to interesting hypotheses regarding the role of natural competence in maintaining the core genome and diversifying the accessory genome, and it would require considerably less sequencing.

So the goal of the preliminary experiment would be to measure structural variation due to mutation versus structural variation due to transformation.

How would this work, and why would it be cheap and easy?

First of all, the experiment-side is extremely easy. Naturally competent recipient cultures would be split in two. One part would be used to prepare untransformed recipient chromosomes; the other part would be incubated with donor DNA for a while, allowed to recover, and then transformed chromosomes would be purified. (We could also include a selection for a donor marker at this point to increase the relative transformation rate.) That’s it.

The sequencing side would also be comparatively easy. Untransformed and transformed chromosome preparations would be sheared, end-repaired, and size-fractionated by gel to 500 bp (as precisely as possible). This DNA would be ready for sequencing library construction and paired-end sequencing.

For a 500 bp library, I previously estimated ~2500X “spanning coverage” in one lane of Illumina sequencing using conservative estimates of the sequencing parameters. “Spanning coverage” is defined as how many times a particular genomic position is found between two mapped paired-end reads. So we’d get to 10,000X spanning coverage in ~1/2 a full run.

How would this help us measure the transformation of donor structural variants into the recipient?

Let’s use an example to illustrate. In the alignment shown above (using GenomeMatcher), the donor genome (86-028NP) is shown on top, while the recipient genome is shown on bottom (Kw20). As is pretty clear, genes Hi_0512 and Hi_0513 are absent from the donor genome, indicating a deletion of those genes from the donor (or possibly an insertion into the recipient). These genes happen to be the HindII restriction enzyme and methylase. The flanking genes are syntenic (so Hi_0511 = tchA, etc.).

Since we know the size-distribution of the library (500 bp), it would be quite simple to spot the deletion allele. Paired-end reads with one end in tchA and the other in rpoC would define the deletion. By contrast the insertion allele would always have tchA and rpoC mappings on different fragments (with the other ends in the HindII methylase or restriction enzyme).

Here’s a way of illustrating what paired-end reads of different kinds of alleles relative to the recipient would look like:
So for our deletion, we’d see paired-ends that mapped to positions further apart than they should be. By having extremely high spanning coverage, we could count how often these kinds of mappings occurred versus the recipient mappings. This would give us the rate of deletion.

In our untransformed library, if we saw the deletion allele, we’d be seeing mutation, while in the transformed library we’d be seeing mutation and/or transformation. Since we know the donor genome sequence, we can distinguish paired-end reads that look like the donor sequence versus de novo mutations that occurred when we grew out the cells. Better still, we could spot transformation-induced de novo mutations by comparing to the untransformed library.

Why else do the untransformed chromosomes at all? Well, structural mutations are likely to occur at a much higher rate than single-nucleotide mutations in many instances. Independent of our interest in natural transformation, doing the control experiment may reveal regions of the genome that are unstable, along with the mutation rate of different types of structural variation. This last part is non-trivial. If, for example, we see a particular deletion that occurred on 50% of the untransformed chromosomes we looked at, it could be that this is an extremely frequent mutation, but it could also have simply occurred early in the grow-out of the culture.

As a control for transformation rates, we don’t have to worry about that, but doing the untransformed control would make us confident that changes we saw were induced by transformation and not simply due to such kinds of mutation.

Happy (Belated) Canada Day!

2009-07-02T16:07:00.000-07:00

(Image: The Canadian-built robot arm attaching the space shuttle docked to the Hubble Space Telescope with the Earth in the background.)

Oh yeah, and my second attempt at preparing periplasm DNA was... inconclusive. But it was pretty interesting to try out. In particular, the TE/CsCl/phenol/acetone extraction was quite compelling visually, involving small bubbles breaking up and reforming. I need to start these experiments out at a smaller scale.

And luckily, we've received radiolabeled dATP, so I can do some more sensitive and controlled experiments next week. The use of radiolabeled uptake fragments will be significantly more sensitive and allow me to use small cultures and follow small amounts of uptake DNA.

Once I've got a functioning uptake assay, I can work out the best purification method and scale up from there.

Periplasm Prep Planning II

2009-06-30T19:13:00.001-07:00

I tried a modification of a periplasmic protein prep to try and purify uptake DNA, which didn't work. There are several possible reasons why the experiment might not have worked, but one simple reason could be that I failed to dissociate DNA from the membranes and cells when I did the chloroform extraction.

I know! Maybe I should try an extraction that has already been used for purifying uptake DNA...

Kahn, Barany, and Smith (1983) PNAS 80:6927. Rather than describe the paper at length here, I just want to show Table 1 and Figure 4b, which relate to my extraction plans:

The first two columns describe the extraction conditions (rows 1-5). Competent cell cultures were incubated with a radiolabeled plasmid and pelleted after DNA uptake (5 or 60 min). Cell pellets were then resuspended in the indicated aqueous and organic solutions (columns A and B) in a 1:1 mixture, gently mixed, centrifuged to separate the phases, and radioactive counts in each fraction were measured.

The remaining columns indicate the relative amount of uptake in the different fractions and the identity of the radiolabeled DNA, either transformed into the chromosome (C) or still a double-stranded donor DNA molecule (D).

In the fourth condition (row 4), the aqueous phase consists of mostly donor DNA! So chromosomal contamination is in the pellet, and the desired donor molecules are in the aqueous phase. Sounds like a scheme. That’s what I’ll proceed with tomorrow.

Why not condition 1, just TE and phenol? Looks good, right? Because of Figure 4B:

The bottom line is that the phenol condition degraded the donor molecules (lane D), whereas the phenol/acetone condition did not (lane I). Here’s the gory details:

Lanes A-D describe phenol extraction of intact donor DNA molecules (Table 1, row 1):
A: input plasmid donor molecule.
B: total DNA from cells 4 min into uptake.
C: DNA extracted and dialyzed out of the pellet.
D: DNA extracted into the aqueous layer (TE).

Lanes E-J describe the phenol/acetone extraction of intact donor DNA molecules:
E: input again.
F: input cut with HindIII.
G: total DNA after 8 min uptake.
H: same as G, but digested with HindIII.
I: phenol/acetone extracted DNA after 8 mins.
J: same as I, but digested with HindIII.

The point of the HindIII digestion was as an additional test for whether molecules were donor or chromosomal. Chromosomal DNA is resistant to HindIII (this being Haemophilus influenzae after all), while donor DNA is not. Also, the HindIII digests show that the recovered DNA is double-stranded, since ssDNA won’t get cut.

Condition 4 it is, then:
Aqueous: TE/1.5 M CsCl
Organic: phenol/acetone, 1:1

That's what I'll try next.

Hmm, I’ll have to remember how to clean the DNA of CsCl after the extraction...I seem to remember doing something in particular at one point...

Periplasm Prep

2009-06-30T08:32:00.000-07:00

UPDATED BELOW

Experimental plan for the day: Medium-scale periplasm prep test-- to purify double-stranded DNA in the “protected state” (the periplasm).

I have two PCR products: (1) a “good” uptake sequence USS-1, and (2) a “bad” uptake sequence, USS-R . I want to compare their uptake into rec-2 cells, which can bring DNA through the outer membrane, but not the inner membrane. This means that if I can specifically enrich USS-1, but not USS-R--and can see the difference in a gel-- then I’ve got a functioning periplasmic DNA prep. Naturally, it’ll probably take several attempts to get working...

I will use a modification of this paper and see if it nets some DNA where it should be. In outline, I’ll: Add chloroform to washed cell pellets. Soak. Extract periplasm with TE. Clean and concentrate. Run on a gel.

Based on other studies with radiolabeled USS-1 uptake, I expect that for 20 ng added to a 1 ml culture, ~50% will be taken up. To see uptake DNA without radiolabel on a gel and for reasonable controls, I will need larger cultures of competent cells than I aliquoted and froze last week.

Here’s my protocol so far:
1) Defrost two tubes of rec-2 (0.3 OD/ml aliquot) into fresh sBHI@37 (2X25 ml); wait ~2-2.5 hrs.
2) At OD600=0.3 / ml, transfer cells to M-IV by filtration.
3) Incubate 100 min @37 to induce natural competence. (negative control: frozen tube of rec-2 (0.1 OD/ml) into fresh sBHI (25 ml).)
4) Split cultures 2X and incubate with 20ng 222bp PCR fragments (USS-1, USS-R, none) / 1 ml M-IV culture (~10^9 cells) for 15-30 min @37, DNase I, EDTA to kill DNase I and other nucleases. Also add USS-1 to non-competents.
5) Spin, wash pellet 3X PBS, chloroform (20-40ul), incubate 20 min @RT (chloroform pellet DNA extraction?, save washes).
6) Extract with 100-200ul cold TE, proteinase? RNase?, p/c extraction, PCR clean-up column (or ppt?) to concentrate.
7) 1.2% agarose gel. Lanes:

Size standard
USS-1 input (2X dilution)
USS-R input (2X dilution)
rec-2 + USS-1 -> chloroform extract
rec-2 + USS-R -> chloroform extract
rec-2 + no dna -> chloroform extract
non-competent rec-2 + USS-1 -> chloroform extract

UPDATE:

Didn't work. A few little mishaps aside (mainly that the chloroform and cells really didn't mix well), I got no USS out of the prep, but did have a fair amount of chromosomal contamination. So clearly, I didn't really get the periplasm specifically, but since I didn't see any USS come through, it may also be that my competent cells weren't really.

I'll try to go through this again tomorrow, but instead of going straight for the periplasm preparation, I'll just lyse the cells, extract the DNA, run it over a mini-prep column, and run it on a gel. I'm not going to worry about the periplasm specifically, but simply that the cells are taking up DNA. When the radiolabel shows up, I can repeat the uptake assay the lab has typically done.

Plan for the next few weeks

2009-06-26T11:41:00.001-07:00

One proposal down, one to go... The next one isn’t due until August 8th, so I’ve got just over a month to get it done. This time, however, I am going to manage my time better, since I need to get some preliminary data and still keep learning how to use a computer.

So here’s my plan for the next several weeks:

(1) Work on the proposal for a limited time each day (~1-2 hrs). I’ll start by developing a detailed outline of what I want to say and the order I want to say it, rather than leaping straight into writing. Based on my experience with this last one and in the past, I find I am an extremely inefficient writer (both with my time and with my words), so hopefully I can improve by having more focused daily goals.

(2) Work on the computational stuff only 1-2 hrs per day. Still need to fix the browser to display the Hin genome. Still want to work out the best way to align the genomes and report differences... particularly enumerating structural variation (non-SNPs). I also need to keep a mind towards what file formats I expect to get from sequencing. It might be particularly useful to try and simulate the kinds of results I might expect from sequencing periplasmic uptake DNA, etc.

(3) The rest of the time will be dedicated to lab work. The priority is to use defined fragments (USS-1 and USS-R) to work out a periplasmic DNA purification protocol. I’ve got cleaned amplified USS-1 and USS-R fragments, and I’m making competent cells of wild-type and rec-2 today. I tried to grow up a pilA mutant to use as a no uptake control, but something was wonky with the strain. All I need now is label. And a calibrated Geiger counter. I’ll get these things done today.

(The image above was made using the Mac-specific application GenomeMatcher. It represents a BLAST alignment between KW20 and 086-28NP across an interval containing several inversions. It also has a bunch of useful-seeming utilities that I'd like to figure out. Now if I could just get it to use my MUMmer program, like it’s supposed to...)

But which would win in a fight?

2009-06-25T16:34:00.000-07:00

The budding and fission yeasts provide for several interesting comparisons. Some of these somewhat illustrate my considerable confusion about the “evolution of sex”. Here’s one:

While Saccharomyces cerevisiae (the budding yeast) prefers to spend time as a diploid in the G1 phase of the cell cycle (with two unreplicated genomes), Schizosaccharomyces pombe (the fission yeast) prefers being a G2 haploid (with one replicated genome). These “preferences” are inferred from two pieces of evidence from each of the yeasts:

(1) Budding yeast cells mate under rich conditions (become diploid), and sporulate under starvation (become haploid). Fission yeast cells remain haploids in rich conditions, while mating and immediately sporulating under starvation (zygotic meiosis).

(2) In cycling cells, budding yeast spends most of the cycle in G1. S-phase and mitosis are separated by only a short G2. By contrast, fission yeast cells reside in G2 most of the time, with only a short gap between mitosis and S-phase.

The interesting thing here is that, while the two yeasts prefer to maintain different ploidy levels (2N versus 1N), they both prefer to have two copies of the genome present (2C). Why might this be?

One suggestion is that this means that there will always be a template for recombination in the event of a DNA double-stranded break (DSB). DSBs are a particularly challenging form of DNA damage for cells. They must paste back the correct broken ends, which requires an intact homologous (identical-by-descent) template to happen with high accuracy.

Another possibility relates to newly introduced recessive mutations. In both cases, a recessive mutation will not normally affect the phenotype of the cell it occurs in, since there are two copies of the genome. But cell division has different consequences in diploid versus haploid cells. In a diploid, mitosis will maintain the heterozygosity of the locus sustaining the new mutation, while in a haploid the wild-type and mutant alleles will immediately segregate in mitosis. So natural selection will act differently on populations that are predominantly diploid or haploid.

An important point regarding the G2 preference of fission yeast: Even a new lethal recessive mutation will not usually kill a particular cell. A cell sustaining such a mutation in G2 segregates the lethal allele to only one of its progeny sister cells, maintaining the wild-type allele in the other.

Diploid cells can still segregate the wild-type and mutant alleles, even in the absence of meiosis. “Loss-of-heterozygosity” (LOH) can occur, if there is a crossover between the heterozygous locus and the centromere of homologous chromosomes. Such crossovers, though rare, occur at measurable rates due to recombinational repair of DSBs and collapsed replication forks. So even in the absence of proper sex, a diploid can still expose its new mutations to natural selection, just at a substantially lower rate than haploids (over many cell generations only some LOH will occur in the budding yeast, but will happen immediately for new fission yeast mutations).

This comparison illustrates contrasting life style choices. While being haploid or diploid has distinct population genetic consequences, the preference of both yeasts to exist with two genomes in a cell says something interesting... I’ll need to think a little more about what that something is before articulating it clearly...

Proposaling continues

2009-06-23T20:42:00.001-07:00

I continue on the saga of writing this soon-to-be-due proposal. I keep struggling with things that I think I understand in my head, but can't succinctly talk about in the text. Luckily, when I hit a wall, I can do some referencing and figure-making...

A taste of DNA

2009-06-18T18:07:00.001-07:00

I gave Rosie a draft of my grant application, so decided to dither with human USS motifs again.

First, for fun, I looked at where several of them were located... to see which genes taste best... and ran across USS in all sorts of random genes. (For example: a kinesin, an adductin, a phosphatidic acid phosphatase, a phosphatidic acid kinase(!), a cadherin-associated protein, a few hypothetical genes and transcriptions factors, etc. etc.) There were also several located outside genes and in gene-poor regions. It would be funny to do a GO annotation analysis, but most certainly a waste of time...

But to address Rosie’s comment: I found 762 10mer USS motifs in the human genome using a Tagscan search. This looks like it meets random expectations, since using Tagscan to count arbitrary 10mer motifs having similar base composition (and a CpG) gave similar numbers.

However, my analytical calculation of expected seemed way off. Where did I go wrong? (One reason why I might actually care about this is so that I could do a more precise analytical calculation of the number off USS motifs I’d expect in the Haemophilus genome.)

I had calculated the chance of getting a USS motif for a random 10mer as G* (25%)^10, for a 25% chance of drawing the correct base at each position (where G is genome size = 3.16e9 bases):

So G * (25%)^10 = G * 1 / 1,048,576 = G * 9.54 e -7 = 2,956 instances.

But since the human genome has only 41% GC content, I might adjust this measurement to be (20.5%)^5 * (29.5%)^5, for a 20.5% chance of drawing the correct base at GC positions and a 29.5% chance of drawing the correct base at AT positions:

So G * (20.5%)^5*(29.5%)^5 = G * 8.09 e -7 = 2,556 instances

These are for only a single strand of DNA. Since there’s also the reverse complement, we do have to multiply these values by 2. I can’t think why this wouldn’t be. So...

50% GC : 5922
41% GC : 5112

So even accounting for %GC failed to bring this calculation down to what Tagscan found.

What about dinucleotide composition?

It’s well known that mammalian genomes have a dearth of CpG dincucleotides, since these are used as sites of gene regulation by cytosine DNA methylation. Methylated cytosines tend to deaminate into thymidines, so there is a mutational pressure on CpG to go to TpG dinucleotides.

Anyways, I managed to find a nice table in this paper reporting dinucleotide composition in several genomes. The paper the table’s authors used for humans was dated from 1962, so these are probably not particularly precise numbers but sufficient for my purposes. It looks like the CpG dinucleotide is ~4-fold less than would be expected for a random genome with base composition like humans’.

This paucity of CpG dinucleotides in the human genome could account for the discrepency. So just for giggles, I ran a few additional Tagscans for 10mers with the correct base composition but lacking a CpG dinucleotide, giving the following numbers: 9143, 6698, and 4343 motifs. Those numbers look a lot more on-target (not terribly precise, but more accurate).

But how can I actually use the known distribution of dinucleotides in the human genome to arrive at an even more accurate estimate of how many USS motifs I’d expect to see? I thought up some crummy ways to sorta-kinda account for the CpG deficit, but would ideally like to be able to use arbitrary dinucleotide frequencies and arbitrary %GC to produce an expected value. I have a terrible feeling I need to use Markov or Ising models and just don’t have the heart for it right now...
----
As an aside, UNIX is pretty awesome. To figure out the number of motifs Tagscan was finding (their output is a list of matches), I simply typed:

> wc -l filename

And that gave me the number of lines in the file. The first row was the header, and the last row was blank, so subtracting 2 from the number gave me the total number of motifs.

If I had Tagscan (and the human genome) on my computer, I could fairly easily set up a script to iterate through a whole bunch of 10mers with specified parameters and draw up a distribution. If I was really cool, I'd exclusively use UNIX commands to do this. This would then allow me to ask what the significance of the number of USSs would be. (Obviously p > 0.05, but that would be one way to do a real statistical test, even if I never figured out how to work out the expected value by an analytical method.)

They're eating our DNA!!!

2009-06-16T21:58:00.000-07:00

I am in the midst of proposal-writing, which makes blogging tougher, but when I hit a writing-block at the end of the day, I decided to dally with a random thing I'd been meaning to figure out.

Several friends of mine, when I tell them about my plans with the naturally competent bacteria, have said, "Dude, you should feed them human DNA!"

Sound silly? Maybe it is, but I went ahead and used TAGSCAN to search the human genome for the 10-bp core USS motif: AAAGTGCGGT and found 762 instances of the USS core. BLAST and BLAT didn't want to deal with me due to how short the query USS motif is.

That means Haemophilus influenzae might find nearly 800 bits of our genomes tasty! I remember from my reading that there's nearly 200 micrograms of DNA per milliliter of our lung mucus, which seems like a heck of a lot. Much of this DNA is probably human...

(For context, on average Haemophilus have a USS less than every 2 kilobases, but humans have an average USS density lower than one in 20 megabases. So while humans have half as many USSs, there's a 10,000 times lower density.)

It isn't surprising the human genome contains at least some "USS". A random 10mer string would have a 1/(4^10), or a little less than 1 / 1,000,000, chance of being the USS motif. But the human genome is more than 3 billion base pairs. So randomly, we might expect to find more than 6000 USS (3 billion / 1 million X 2). (The 2 is to also count the reverse complement.)

If I did that right, then the observed number of USS motifs in the human genome is rather less than expected. That's interesting...

However, that was a crude estimate of the expected number of USS motifs. The USS core sequence has 5 GC bases and 5 AT bases, giving a GC base composition of 50%, whereas the human genome has only 41% GC content. There's also the issue of dinucleotide frequencies. For example, the CpG dinucleotide is underrepresented in the human genome, but happens to appear in the USS. I've been trying to figure out a rational way to incorporate this type of information into my estimate of expected, but so far have failed to do so properly. Regardless, my estimate of expected is certainly too high. The question is: how much so?

To produce control numbers, maybe tomorrow I'll run a few other TAGSCANs for arbitrary 10mers with the same GC content and a CpG and see how many come up. I haven't really looked at the distribution of USS, except that the number of USS per chromosome is highly correlated to chromosome size (R^2 = 0.91). I might also predict that USS will fall into more GC-rich regions of the genome.

But for fun, assuming that the result held, what might it mean if the USS motif is significantly underrepresented in the human genome? I can hardly imagine that Haemophilus could be responsible ("it ate them!!!"), but maybe it could work in the other direction? Perhaps the USS motif is only coincidentally somewhat rare in humans, and as such makes a good sequence to use for preferring conspecific DNA uptake? If the USS motif was extremely abundant in the human genome, then naturally competent Haemophilus might not take up conspecific DNA as easily? Hmmm...

Regardless, think about that when you fall asleep, folks... There's bacteria inside you, and THEY'RE EATING YOUR DNA...

UPDATE:

Putting my proposal together continues, but in a bout of procrastination I did go ahead and run a few random 10mers with the same base composition (and a CG dinucleotide) through TAGSCAN, and it looks like the observed number of USS motifs is about expected for 10mers of similar composition.

The USS motif:
AAAGTGCGGT 762

Randomized USS motifs:
GTACGTAAGG 414
CGAGAGAGTT 761
GAATACGTGG 1112
AGATAGCGTG 714
GTAACGAGTG 458
AGACGTTAGG 713

Nothing to see here... Move along...

Supragenomicisticexpialidocious

2009-06-09T22:19:00.000-07:00

Howzabout that “supragenome”?

Polymorphic gene content

Genome sequencing efforts have captured substantial variation in the gene content among closely related bacteria. For example, this nice study by Hogg et al. 2007 compared the genome sequences of thirteen Haemophilus influenzae: They found a “core genome” of ~1500 genes, along with another accessory (or contingency) genome of ~1300 genes. Any given isolate had a subset of the accessory genome, numbering around a few hundred extra genes beyond the core in each isolate.

So any two isolates have substantial amounts of DNA that is unshared. For example, Hogg et al. report that the Kw20/Rd and 86-028NP isolates differ by nearly 400,000 bp within only ~250 indels (Table 5-- mean size: ~1.5 kb, median size: ~300 bp). The genomes are < 2 Megabases long, so that’s about 20% of these chromosomes that is non-homologous.

(As an aside, the Hogg et al. paper’s methods section introduced me to MUMmer, which I previously discussed. The indels and other rearrangments are essentially defined by breaks in the alignment produced by the nucmer utility. More on this in the future...)

There are several possible arguments related to uptake specificity and the "supragenome":

Uptake specificity for variation?

Large variation in orthologous gene content suggests to Hogg et al. and others a “distributed genome hypothesis”, in which natural transformation can shuffle (or re-assort, or segregate) the accessory genes between isolates. This would then presumably allow for the rapid acquisition and loss of different genes (diversification) from within a given genetic background and thus perhaps rapid adaptation to environmental changes (or shifting host defenses). The “distributed genome hypothesis” is then implicitly related to the “sex hypothesis” for the maintenance of natural competence.

Uptake specificity for conservation?

On the other hand, natural transformation could also maintain the “core genome”. Thus, if there is plenty of conspecific DNA uptake, any bit of “core genome” taken up could be replaced in a cell that had lost it. This nice study by Treangen et al. 2008 using several neisserial genome sequences showed that DNA uptake sequences (DUS, the neisserial equivalent of USS) existed at a higher density in “core” regions of the genome than in the substantial alignment gaps between isolates (containing indel poymorphism).

Again, there is some indication of the “sex hypothesis” for the maintenance of natural competence, but it works in the opposite direction, maintaining the core rather than shuffling the accessory. I think the argument goes like this: (1) The core genome likely defines the more essential portions of the genome, since by definition any accessory genes are not required to live. (2) DUS could have been selected for within this partition of the genome, since it would help to maintain the more essential gene functions within a population. (3) Therefore the high number of DUS sequences could be a product of natural selection to maintain the integrity of the “core genome”.

Uptake specificity for no reason in particular?

However, there’s another possibility the authors partially explore that does not involve selection for DUS distributed throughout the genome, but represents almost the opposite model. Instead of selection, pehaps DUS accumulate due to happenstance intrinsic biases in the uptake and/or recombination machinery by a neutral molecular drive. So sequence variants that arise with a higher chance of being taken up later are more likely to spread through populations than variants with a lower chance of uptake. Thus the “core genome” could partially be that way, i.e. conserved across isolates--not exclusively because of essentiality or usefulness--but also by virtue of containing lots of DUS. So rather than DUS being selected for in order to maintain the core genome, segments of DNA containing DUS are simply mre easily replaced in lineages that lost them.

An affiliated idea suggests that if some accessory genes were from distant relatives and arrived by horizontal transfer by some mechanism besides natural competence, these sequences would not have had time to accumulate uptake sequences yet. Thus the paucity of DUS in the accessory genome might be in part due to the more recent arrival of that sequence in the genome, so the effects of drive have not yet become evident, rather than a specific selection pressure to maintain DUS in more important segments of the genome.

(The Treangen et al. paper introduced me to another genome alignment tool called M-GCAT. I’ve played with it a bit and managed to produce some figures effectively the same as what appears in their supplementary data-- the picture above (along with alignment files resembling multi-FASTA format) but have the unfortunate problem of being unable to re-load analyses I’ve performed later due to some kind of Python error. More in this later as well...)

How to analyze the core and accessory genomes myself?

I’ve clearly got a lot more thinking to do regarding these core and accessory genomes... Especially in light of the horizontal gene transfer issue.

But first I’d better figure out simply how to define the core and accessory genomes more specifically.

I’ve begun this by examining the gaps in the .rdiff and .qdiff output of dnadiff (a pairwise comparison of two genomes) to try and do some basic analysis myself. In a future post, I’ll report on my progress with this, but for now, I’ll just mention that most of the gaps are not strictly insertions or deletions, but are rather insertional deletions. Alignment gaps include both reference and query bases. But I still need to try and understand how dnadiff produced its .report output before I can get much further...

Joint molecules...

2009-06-05T15:20:00.001-07:00

There are several classes of transformation event that could be mediated by natural competence and recombination. The first thing that I had to wrap my head around (after having come from the double-stranded break world of recombination) was to realize that uptake DNA recombining into a host chromosome is single-stranded. I’m more used to drawing recombination models involving two broken ends of a double-stranded molecule.

So to kick things off, in the above figure, I’ve drawn what I think the joint molecule intermediates look like between polymorphic donor ssDNA invaders and recipient dsDNA substrates for the four major classes of transformation products I can envision.

The red (red/pink) lines indicate an ssDNA that’s found its way from outside the cell to the host chromosome (in blue/light blue). Aligned regions are base-paired. Unaligned regions are not...

The amount of transformation will depend not only on how many of a given joint molecule are formed, but also how they are resolved and the rate and directionality of “mismatch” correction.

One important aspect of the strand invasions depicted above will be the extent of heterology. The more polymorphisms in the donor sequence, the less stable the joint molecule will be.

Each class of polymorphism deserves special mention:

(1) Single nucleotide polymophisms: There are 12 possible heteroduplexes involving single-nucleotides, if we keep track of donor and recipient. So any of the 4 bases could go to any of the remaining 3... 3 X 4 = 12. In my preliminary analysis, I found nearly 39,029 SNPs between our reference Rd strain and the clinical isolate 86-028NP falling into each one of these different classes (with a paucity of G->C and C->G changes). Even without precisely measuring the frequency of transformation for each individual SNP, we may be able to estimate the relative transformation frequencies for these different classes of SNPs.

(2) Insertions: If the donor has an insertion relative to the recipient, recombination could yield an insertion. In the case of insertions, the uptake sequence can be anywhere on the fragment. base pairing of the substrate will require the formation of a loop. The size distribution of the input DNA will also dictate the size of possible insertional recombination. How the mismatch correction machinery handles insertional joint molecules like this is unclear to me. My preliminary analysis indicated 149 insertions in 86-028NP relative to Rd (mean size = 1162 bp, but the median only a few hundred).

(3) Deletions: I made this distinct from insertions for two reasons: (1) The uptake sequence (obviously) must flank the deleted segment. An absent piece of DNA can’t be used for uptake. This imposes directionality on deletions that’s distinct from insertions. (2) It’s unclear how the mismatch machinery would act on the two central joint molecules depicted... It restoration repair more likely in one case than in the other? My preliminary analysis indicates 137 deletions in 86-028NP relative to Rd (mean size 1904 bp, again with a much lower median).

(4) Rearrangements: I’ll give this very short shrift for now. If a particular piece of ssDNA spans a rearrangement breakpoint, then strand invasion of each end into two separate chromosomal positions could mediate an inversion or several other kinds of rearrangement. My preliminary analysis with MUMmer suggests ~9 inversions and 26 transpositions between the two strains.

There are many many ways in which rearrangements in the donor might produce fragments that could mediate a rearrangment of the recipient. It’s going to take some time to enumerate these. Indeed, even if the donor was perfectly colinear with the recipient, some fragments could still mediate rearrangements.

For example, looking at a close-up of the MAUVE output between Rd and 86-028NP shows that the rRNA gene cluster of 23S and 16S (in Red) often span the rearrangement breakpoints between the two strains:
This could arise in a traditional mutational way (i.e. recombination between inverted repeats would cause an inversion), but could also be mediated by transformation. For example, if a DNA fragment that contained a bit of rDNA and a bit of flank first invaded the wrong rDNA copy, the other end could then grab the original flank... causing a rearrangement.

That is, in the rearrangements-involving-repeats scenario, it would be quite possible to produce rearrangements distinct from either the donor or recipient DNA.

This is going to take a while to flesh out, as there are several possible outcomes to this sort of thing. It doesn't make things any easier that I also now have to consider that the chromosome is circular... sigh... I'd better try and remember what plectonemic and paranemic mean...

Nevertheless, the upside is that it will be possible to measure indel and rearrangement transformation rates with considerably more ease than the SNP class, due to the use of “spanning coverage”.