Is "Ecological Epistasis" a Good Term?

I've been inspired by a couple of recent twitter conversations I've had to write a post and basically lay out why using the phrase "ecological epistasis" triggers my population genetics spidey-sense.

The first conversation happened about a month ago while I was sitting through previews for the completely satisfying movie Pacific Rim.

Epistasis between the nuclear genome and #microbiome underscores the #hologenome Not a far fetched idea and soon to be non-controversial.
— Seth Bordenstein (@Symbionticism) July 14, 2013

The second happened yesterday (Ian is live tweeting a bunch of talks at BEACON and Maren Freisen (@symbiomics) was talking about her Medicago research)

Friesen: Although she phrased it as ecological epistasis. #2013beacon.
— Ian Dworkin (@IanDworkin) August 14, 2013

I'm using this space as a way to crystallize my thoughts and to try and solicit other opinions. I'm also going to try and be involved in Seth Bordenstein's G+ chat on the hologenome in a couple of weeks, so consider this a bit of a warmup.

Epistasis is a tricky word. Problems arise when different (yet highly related and somewhat overlapping) groups start to use the same word yet mean different things. One of the people responsible for making me the scientist I am today has written on this topic (here and here), so I won't go too much into it. To summarize though, you can define epistasis in the quantitative genetics sense (multiple loci interacting in a non-additive way), in the small genetic sense (two proteins actually interact or function in the same pathway as would be found by genetic screen) or you can define it in the larger genetic sense (interactions between multiple genes). Neither is wrong per se, but use of the same term can get confusing depending on your audience. I don't have a problem with any of these definitions, but I just think that making headway in biology always becomes more difficult when you have to start referencing quotes from Justice Potter Stewart.

So here's a couple of my problems with "ecological epistasis". This uses the latter definition of the term that I mention above, gene interactions writ large. If you have two co-evolving organisms, genes from one organism interact with genes from the other organism in the population genetics sense, fitness of one organism is dependent on the other organism, all well and good (if I'm misinterpreting, please let me know). We already have language to describe these circumstances though, in terms of Gene by Environment (GxE) interactions without the need to invoke the "e" word. In this case each organism is the co-evolving organisms environment variable. Why not modify this term a little bit and call it GxEG (environment/genetic) interactions? There is nuance in this specificity that isn't captured by using the term epistasis. This isn't really my area of expertise, but I'm guessing that dynamics that apply to interactions between genes residing in different genomes may be inherently different than those in linked together and vertically inherited.

The second idea that needs clarifying IMHO is what the limits on interactions between organisms are when speaking in population genetic terms. When I've seen the phrase "ecological epistasis" used, it's in reference to interactions between intimately co-evolving organisms. However, if you are going to define the term as interactions between genes in different organisms without specificity, you could extend the definition in absurdum. Much of my work focuses on plant pathogenic bacteria, and genomes of both the host and pathogen encode for proteins that mediate interactions between the two. Is this "ecological epistasis"? Stepping further back, this morning I killed a cricket that tormented my household last night (keeping my 8 month pregnant wife awake more than usual...I have no regrets about my actions). In result this is no different than a pathogen killing a host, just that co-evolutionary interactions are weaker between me and the cricket population at large than in typical pathogen/host dynamics. My foe and I both have genomes that encode for proteins that ultimately mediated our interactions this morning. Is this "ecological epistasis"? Ian Ziering managed to fight off a shark with a chainsaw:

                                                          Is this "ecological epistasis"?

I'll save my hologenome critiques (great term, needs limits on the definition) for a future blog/G+ chat. My point is simply that if you start defining interactions between organisms, that these interactions can take a wide variety of forms that you may not inherently consider. Specific wording could avoid me having to make sharknado references.

It's not that I think that using "epistasis" in the context of interacting organisms is improper. I think that the term is muddled enough as is that it doesn't make sense to use it for the sake of linking onto an already established (and muddled) term. Using "ecological epistasis" doesn't clarify things in the way that a more nuanced term could, at least to me, but maybe I'm just missing something?

Update: Maren Friesen has clarified what she was referencing in her talk:

ecological epistasis!=GxE; it's non-additive effects of spp on ecosystem @Symbionticism @surt_lab @IanDworkin @Graham_Coop @liveinsymbiosis
— ML Friesen (@symbiomics) August 16, 2013

So...4) non-additive interactions between species (not genes)

Update 2: Great response by Maren Friesen

What if Diet Soda Wasn't Diet?

Ideas are cheap, actually pulling off the experiments is the difficult part. Sometimes these experiments aren't even possible to do at the present time. I'm probably not the only one who has a running list of experiment ideas in a text document, many of which will never see the light of day. I'm going to start something new around here by posting about research/experiment ideas that I think would be interesting and informative, but which I have absolutely no time to carry out right now (however, if you're up for collaborating definitely shoot me an email!). I'm naturally curious, so it would give me great pleasure to see SOMEONE figure out the answers to these observations or actually carry out experiments. Hell, someone might have even already done the experiments (if so, please send me a link in the comments!). Use these posts for inspiration or even just to get a feel for how I think about science, especially if you're keen on being my grad student or postdoc in the future. Point is that ideas are cheap but my mind keeps grinding. So without further delay here's where it goes sometimes...

Since my undergraduate days I've had a thing for "diet" drinks. Soda, fruit juice, etc...I always go for the "light" version. First it was the deliciously aspartame-filled Diet Coke (I definitely don't have phenylketonuria) and I've since transitioned into deliciously sucralose-filled products. Supposedly, drinking diet products can help you shed weight (see here but also here). Diet soda et al. have no calories because they contain artificial sweeteners that can't be metabolized by your body. I've always believed this, I could be completely wrong but this seems right. Relevant to this story, it does seem as though drinking diet soda can actually make you gain weight and can increase the incidence of type II diabetes (Hmmmm...)

Here's the thing. Your body is also teeming with microbes, especially in your digestive tract, billions of them. Some of these can even aid digestion by breaking down products. If there is one thing I know that microbes are good at, it's adapting to use novel resources. Unexploited potential energy sources are just another niche that microbes can thrive in. I don't see why microbes can't break down, or easily evolve to break down, aspartame, sucralose, and Truvia.

So here's a couple of potential experiments. I'd like to take some gnotobiotic mice, as gut flora may influence their weight. In the lab I'd adapt a suite of common gut microbes to growing on one of the artificial sweeteners. Then I'd transplant these bacteria back into the gnotobiotic mice in one group, and "ancestral" bacteria that can't break down the sweetener into another group. Next I'd feed different groups of mice a diet supplemented with one of the three sweeteners (as well as a regular control diet). The null hypothesis in this case would be that there will be no different in weight gain attributable to evolved vs. un-evolved microbes. A second experiment is really just a converse of the first. Basically I'd feed mice with "normal" gut flora a diet supplemented with one of the three sweeteners or the control diet with none. Then I'd measure if the ability of gut microbes to digest the artificial sweeteners changes over time.  Null hypothesis here is that there would be no change in the microbe's abilities to break down artificial sweeteners over time.

So that's the outline. Thoughts? Has this been done? If someone does this will the artificial sweetener industry put a hit out on them?

UPDATE: Thanks for the input folks! Definitely understand now that there is much less artificial sweetener in diet soda than regular. Maybe not the best example, but, doesn't change the thought experiment. I know people who replace regular sugar with sucralose or Truvia in coffee and baking. They use the exact same amounts so, plus or minus differences in the molecular formulas, there's roughly the same potential mass going in.

Yes Mom, I do study GMOs

Reading Amy Harmon's great piece on GMOs and citrus greening inspired me to write this post. What follows is a slightly fictionalized account of a conversation I had with my mom. Don't worry, these conversations actually happened pretty much how I describe. I recently stopped home for a couple of days (my favorite conference to attend is but 2 hours away from my parent's house in VT) and eventually found myself arguing with her about the benefits of genetically modified organisms (GMOs). Her main comment was something along the lines of "How do you know what happens when you stick a lemon gene into corn. There could be horrible side effects". I found myself making the case that substantial scientific evidence exists concerning on the safety of GMOs and human health as well as describing how corn was completely different from it's non-domesticated (and hence non-genetically modified) ancestor teosinte. Standard stuff really, and the conversation ended in rhetorical standstill as is par for the course when I disagree with my parents.

A few hours later my mom asked me about my own research program. I started to tell her about horizontal gene transfer (HGT) in microbes, how the transfer of such genes is a driving force for microbial evolution, and finished by describing how we know very little about the side effects of HGT. Then it hit me, my research links up perfectly with the discussion about the side effects of GMOs. HGT is a natural process that is effectively indistinguishable from the creation of GMOs. At a forest through the trees level specific genes start out in species A and are transferred to species B. In the case of HGT, the vector for transfer can be a plasmid/phage/transposon/etc whereas for GMOs the vector can be a plasmid/phage/transposon/etc. In the former, random chance (and many other factors such as environmental proximity) determine which HGT events occur, whereas in the latter it's humans that determine which occur. The only (arguably subtle) difference between HGT and GMOs is what structures selection pressures. In the case of HGT, natural selection culls out unproductive combinations of genes and backgrounds whereas with GMOs humans directly select and screen for the most "productive" combinations. You could even argue, thinking about the selection pressures on the movement of antibiotic resistance genes in microbial pathogens, that there is substantial overlap even in selection pressures. If you just focus on the movement of genes and don't worry about the how, the natural process of HGT and artificial process of GMO creation are exactly the same. What we learn about the side effects of HGT will be directly applicable to understanding the side effects of GMOs, i.e. for figuring out how badly a single lemon gene would screw up your tasty corn. My research can actually be able to address my mom's original question.

"Ahh...but Dave", you might say, "microbes are different than corn". Well, it turns out that HGT occurs much more frequently in multicellular eukaryotes (like corn) than we previously thought. Aphids come in different colors because they have acquired carotenoids from fungus. A substantial portion of the genome that codes for your steak is potentially derived from snakes (arguing about the precise percentage can get a little hand-wavy since this may only be one HGT event). Michael Douglas may have gotten oral cancer because of viral HGT. Perhaps most relevant to this discussion, there is a gene in sorghum and rice that is has been acquired from a parasitic plant. The list goes on and on and will only grow as more genomes are sequenced. Yes mom, even though I study the transfer of microbial genes, I'm still studying nature's GMOs.

Just because you can do something doesn't mean you should...on using dN/dS with horizontally transferred regions

(Disclaimer: I'm writing this post because I've recently seen a couple of papers that measure dN/dS within obviously horizontally transferred genes and I view these analyses as sketchy for the reasons described below. If I'm making obvious mistakes or there is a body of literature that goes over how these analyses are OK, please point me in that direction!)

One morning during my senior year in college, I was driving my car around town and had the sneaky suspicion that my car's lighter wasn't working. I pushed the lighter in for a couple of minutes and then examined the end. Since the coils weren't red it only fueled my suspicion. What I did next stands firmly in the pantheon of stupid things I've done in my life...I placed my thumb over the coils to see if they were hot. I honestly can't tell you why I did this, but every sensation for the next couple of seconds has been burned into my memory. It took about four months for my thumbprint to grow back. Long story short, just because you can do something doesn't mean you should.

Molecular evolutionary analyses have gotten extremely easy to perform in the last couple of decades. There exists a variety of plug and play and freely accessible computational tools and programs that will chew up any DNA sequences you input and spit out numerical results related to evolutionary scenarios. These programs will work even if you input completely fake sequences and made up taxa. Like many programs (looking at you t-test in Excel) they are agnostic to the underlying assumptions. It's up to the user to make sure comparisons are valid.

dN/dS and it's cousin Ka/Ks are often used to measure whether positive selection has acted on a set of nucleotide sequences. Explained very briefly, for a given nucleotide sequence, dN is the number of non-synonymous nucleotide changes in the focal sequence compared to a reference sequence. dS is the number of synonymous changes in those same sequences. For these analyses, synonymous mutations are assumed to be selectively neutral and thus represent the background rate of evolution. If there is no selection acting on a protein sequence, it is assumed that non-synonymous changes should occur as frequently as synonymous. Therefore, all else equal (in a very unrealistic world), proteins under no selective pressure should have a dN/dS ratio around 1. A dN/dS ratio > 1 indicates positive selection because in this case non-synonymous changes occur more frequently than the baseline synonymous rate and selection must therefore be increasing dN. A dN/dS ratio < 1 indicates purifying selection because non-synonymous changes are being selected against and are therefore found in lower numbers than expected. I'm glossing over a lot of nuance in these analyses but that's the 30,000 foot view.

While you can calculate dN/dS for a bunch of sequence blind to their evolutionary relationships, this number doesn't really give you much insight into the evolutionary process. Alternatively, if you have a set of sequences for comparison and a phylogenetic tree showing evolutionary relationships between these sequences you can orient where and at what rates these nucleotide changes have occurred. If you use a program uses statistics to infer evolutionary models (like PAML), you will be asked to input a phylogeny for this very reason. To avoid confounding yourself with circular calculations, it's important to use phylogenies estimated from loci outside the scope of the analysis, even though in many cases phylogenies built from your loci of interest will (for the most part) match those built from the rest of the genome.

Genes undergoing horizontal gene transfer (by definition) have evolutionary histories and phylogenies that differ from the rest of the genome. You can certainly build a phylogeny using these sequences, and you can use this as your tree for calculation of evolutionary rates,  and PAML will accept your input files without questioning them. You will be given a numerical answer when calculating dN/dS (or more specifically omega). In the best case scenario you will then go on and publish papers, generate press releases, and maybe even sit in on a podcast to describe your awesome results. The problem that arises though, is you have no clue as to the true evolutionary relationships between the genes of interest. The phylogeny just tells you how similar a given set of sequences are, but it doesn't necessarily tell you the evolutionary history of horizontally transferred genes. Take for instance a specific gene present on a highly transmissible plasmid (Fig. 1).

Figure 1: If Your Gene of Interest Undergoes Horizontal Gene Transfer, Good Luck Guessing It's Evolutionary History

This plasmid can be passed around from strain to strain, and may experience a number of different environments. In some environments/backgrounds your gene of interest may be under positive selection, in some cases negative, in others no selection at all. Unless you know the complete transmission history of this plasmid you are stuck sampling gene sequences from the last strain it's found in. In the example in Fig. 1, you might assume that your gene of interest is under positive selection in Species C compared to A but in reality it's neutral in both cases. Moreover, if you were to compare this gene from species B and C you would think that there was no difference in selection between the two. When you build a phylogeny with your genes of interest, and sit down to calculate dN/dS (or even worse omega), all you can see is that non-synonymous changes have occurred at some point since the divergence point between two sequences.

Alternatively, bacterial loci may recombine with divergent orthologues from closely or distantly related strains through processes like natural transformation. In this case only a small fraction of the codons in a gene may be replaced with divergent sequences, but such events can totally skew dN/dS ratios (Fig. 2).

Figure 2: Site Specific Recombination Disrupts Interpretation of dN/dS

Although the neutral theory provides a decent null hypothesis for understanding selection on mutational events in natural populations, there is no comparable model for bacterial recombination. Are a majority of homologous recombination events neutral? Are a majority of recombination events subject to positive selection? We have very little understanding of how selection and homologous recombination interact in natural bacterial populations, but we do know that recombination is very frequent within some lineages even across "housekeeping" genes. Along these lines, a recent study investigating Staphlococcus aureus and Clostridium difficile suggests that frequent recombination events bring with them many more synonymous polymorphisms than non-synonymous. Offhand, I can think of many scenarios where dN/dS ratios wouldn't match the underlying evolutionary dynamics because of horizontal transfer but the fact of the matter is that we as a community have no clue as to what to expect.

All this brings me to the ultimate take home message of this post. If there's a hint that horizontal gene transfer may be affecting your analyses of selection, probably best to avoid having me as a reviewer because it's going to take a heck of a lot to convince me that the methods are sound.

Update: Matt Barber (@MattFBarber) reminded me of a great review paper from Jesse Shapiro relevant to this post, but which goes much more in depth into a variety of analyses. Also came across an awesome recent paper from Sheppard et al. (h/t Mark Pallen, who referrs to the "et al." as the Sheppard, Didelot, Falush posse) on how to use bacterial recombination to your advantage for GWAS studies.

The story behind "Exploring the costs of horizontal transfer"

My new review in Trends in Ecology and Evolution went live last week. This paper, while not really experimental, took a little bit of a circuitous route and a bit of luck. For all of you out there sitting on ideas for reviews/opinions but not knowing how to get these published, here's how it happened (with a bit of philosophy of how to be a researcher thrown in).

One of the important skills I first learned in grad school was how to delve headfirst into a completely new topic and figure out the salient and relevant points. Sure, you might get sent a paper of interest by a colleague or have one picked for journal club, and these are great for getting a cliff's notes version of an area of research, but to understand a topic you need to know the background. One of the best ways for me to figure out the background (i.e why is the question interesting, what has been done in the past, what are directions for future research) was to find a review article and use that to point towards new references. These references can be found in the articles themselves, but its also quite helpful to search for other articles that cite the review. It doesn't even have to be a new article, because the rabbit hole of references eventually leads to the present day.

Somewhere along the way I developed a soft spot for the Trends family of journals (I know...Elsevier is evil...fully acknowledged, Frontiers is turning into a great place for reviews in the future). Trends articles were clear, concise, opinionated, and a great foothold for jumping into new areas. While there are other equally great resources for reviews like the Annual Reviews family, Bioessays, MMBR, etc...I set it as a goal early on in grad school to publish a first author article in Trends.

There are two different paths to write a Trends article: 1) you can receive an invitation from one of the editors or 2) you can submit a short pre-proposal and see if an Editor likes the idea. I was part of an invited review once, but this new article was the product of submitting and resubmitting pre-proposals.

When I started my lab I wanted to break away from what my postdoctoral advisor (Jeff Dangl) is known for and get back to my roots in microbial evolution. Sure, I still do a lot of phytopathology work in my lab, but I try to do this in the context of understanding how microbial pathogens evolve and adapt to new environments. This strategy has its positives and negatives, but long story short I found myself writing grant applications where the first page was devoted to explaining why the questions I was asking were interesting and important. This was too much space devoted to justifying my questions. Part of this is my lack of grant writing experiences, but part of this was that I felt I had to explain why I was asking the questions I was asking because, while there were many other articles preceding my ideas, there was no article (*as far as I can tell, point me towards these if they exist and I missed them please!) that laid these questions out in the context I was thinking about. I was simply using too many words to say something that could be easier said in a long review article and then cited.

My first stab at getting this review published was a pre-proposal to Trends in Microbiology. Obviously, this didn't work. I took a little bit of time, re-jiggered the ideas and foci, and submitted another pre-proposal to Trends in Ecology and Evolution (TREE). This time a very kind editor thought enough of what I wrote to give me a chance at a full article. I hadn't actually written the article yet, but the ideas were circulating in my head and most of the references I ended up using were gleaned from iterations of grants I was writing. I had about 3 months to write the full article, which is both plenty of time and not close to enough time, but I was able to pull together a draft and circulate it amongst my colleagues before submitting for formal peer review. This piece actually started out as an opinion, simply because I wasn't really sure if people were thinking about things the way I was. There is a feeling in science where you are both terrified and excited in the same moment. Either you are an idiot who is seeing things that other people have seen before and simply don't recognize this, or you are actually seeing things in a new way. The same feeling occurs when dealing with great new experimental results, except there is an extra option...either your experiment is 1) awesome sauce! 2) a trivial result that other people have seen before but which you haven't realized for one reason or another, or 3) a lab mistake.

 I got the first reviews back and realized that I was onto something that other people were thinking about already (so...no controversial opinion needed), but that there was definitely a place for what I was writing. It's always difficult to read critique, but the reviewers actually did a great job pointing me towards other papers and helping me discover/emphasize other interesting research findings and directions. Specifically, the first iteration was too heavy on history and specifics and too light on the evolutionary implications of horizontal gene transfer. I made these changes, added a figure (in retrospect, should have changed the font and cut down on whitespace but I'm very happy with the information that it conveys), and sent it back in for a second review. This time it went through without a  hitch, and I celebrated accordingly. Side note: I'm a slightly large person and have a tendency to break things in the house and the lab (My wife calls me Shrek). The sequence of the disruptive protein in figure 1 is an homage to my tendencies. A slightly less than subtle Easter egg, but I like trying to put those in my papers (there are plenty more subtle ones in other papers...). I'm also not alone in doing this.

The take home message is not to worry about whether something is good or not, just submit and see what happens. For a long time I was scared/worried about writing pre-proposals to the Trends editors, but it was a fairly seamless process once I went through it. You don't have to wait for the magical email invitation, just take the initiative and see if your idea flies. I'm writing my first grant after the review is out, and it's definitely much easier to write the background now. The last couple of paragraphs of the review also nicely set up some other manuscripts I'll be submitting this summer.

Follow the biology (pt. 5) Season Finale

Yesterday I described assembling a reference genome for Pseudomonas stutzeri strain 28a24, in order to identify causative mutations for a brown colony/culture phenotype which spontaneously popped up while I was playing around with this strain. It's a pretty striking phenotype:

In addition to Illumina data for the reference parent strain, I also received sequencing reads for this brown colony phenotype strain as well as a couple of other independent derivatives of strain 28a24 (which will also be useful in this mutation hunt). There are a variety of programs to align reads to the draft genome assembly, like Bowtie2 and Maq. However, because it's a bit tricky to align reads back to draft genomes and because I'd like to quickly and visually be able to inspect the alignments, I'm going to use a program called Geneious for this post.

The first step is to import the draft assembly into Geneious, and in this case I've collapsed all of the contigs into one fasta read so that they are separated by 100 N's (that way I can tell what contigs I'm artificially collapsing from real scaffolded contigs). Next I import both trimmed Illumina paired read files, and perform a reference assembly vs. the draft genome. In Geneious it's pretty easy to extract the relevant variant information, like places where coverage is 0 (indicating a deletion) or small variants like single nucleotides or insertions/deletions. Here's what the output looks like:

For the reads from the brown phenotype strain there are basically 167 regions where there are no reads mapped back to the draft genome. A quick bit of further inspection shows that these are all just places where I inserted the 100 N's to link together contigs. There are also 327 smaller variants that are backed up with sufficient coverage levels. My arbitrary threshold here was 10 reads per variant, but my coverage levels are way over that across the board (between 70 and 100x).

Here is where those other independent derivatives of the parent strain come into play. There are inevitably going to be assembly errors in the draft genome, and there are going to be places where reads are improperly mapped back to the draft genome. Aligning the independent non-brown isolate reads back against the draft genome and comparing lists of variants allows me to cull the list of variants I need to look more deeply at by disregarding variants shared by both. After this step I'm only left with 3  changes in the brown strain vs. the reference strain.

The first change is at position 2,184,445 in the draft genome, but remember that the number here is arbitrary because I've linked everything together. This variant is a deletion of a G in the brown genome.

Next step is to extract ~1000bp from around the variant and use blastx to give me an idea of what this protein is.

Basically it's a chemotaxis signaling gene. Not the best candidate for the brown phenotype, but an indication that the brown strain probably isn't as motile as 28a24. Next variant up is at position 3,555,964. It's a G->T transversion potentially involved in choline transport...still not the best candidate for the brown phenotype.

The last variant is the most interesting. It's a T->G transversion in a gene that codes for homogentisate 1,2-dioxygenase (hgmA)

To illustrate the function of this gene, I'm going to pull up the tyrosine metabolism pathway from the KEGG server (hgmA is highlighted in red).

The function of HgmA is to convert homogentisate to 4-Maleyl acetoacetate. Innocuous enough and I'm no biochemist so in pre-Internet world I'd be somewhat lost right now. Luckily I have the power of google and knowledge of the brown phenotype so voila (LMGTFY). Apparently brown pigment accumulation in a wide variety of bacteria is due to a build up of homogentisic acid. There is even this paper in Pseudomonas putida. Oxidation of these compounds leads to quinoid derivatives, which spontaneously polymerize to yield melanin like things. Without doing any more genetics I'm pretty sure this is what I've been looking for. It's a glutamate to an aspartate change at position 338 in the protein sequence, a fairly innocuous change but which is (if this truly is the causal variant) in a very important part of the protein sequence. From here, if I were interested, I would clone the wild type version of hgmA and naturally transform it back into the brown variant of 28a24 to complement the mutation and demonstrate direct causality. I might also try and set up media without tyrosine to see if this strain is auxotrophic, which is a quality of other brown variants (see the P. putida paper above). For now I'll just leave it at that and move on to another interesting project that I can blog about.

Follow the biology (pt. 4). Now we are getting somewhere...

First off, apologies on the complete lack of updates. In all honesty, there hasn't been much going on with this project since last summer but now I'm finally at the point of figuring out what kinds of genes are behind this mysterious brown colony phenotype in P. stutzeri. 

Picking up where I left off , I've been through many unsuccessful attempts to disrupt the brown phenotype with transposon mutagenesis. The other straightforward option for figuring out the genetic basis for this effect is to sequence the whole genome of the mutant strain and find differences between the mutant and the wild type. Luckily this is 2013 and is easily possible. Confession time, while I know my way around genome scale data and can handle these type of analyses, I am by no means a full fledged bioinformaticist. I'm a microbial geneticist that can run command line unix programs and program a little bit in perl and python, out of necessity. If you have advice or ideas or a better way to carry out the analyses I'm going to describe below, please leave a comment or contact me outside of the blog. I'm always up for learning new and better ways to analyze data!

Long story short, I decided to sequence a variety of genomes using Illumina 100bp PE libraries on a HiSeq. For a single bacterial genome one lane of Illumina HiSeq is complete and utter overkill, when you can multiplex and sequence 24 at a time it's still overkill but less so. In case you are wondering, my upper limit it 24 for this run because of money not because of lack of want. It's a blunt decision, but one of many cost-benefit types of decisions you have to weigh when running a lab on a time dependent budget.

Before I even start to look into the brown mutant genome (next post), the first thing that needs to be done is sequencing and assembly of the reference, "wild type" strain. The bacteria I'm working with here is Pseudomonas stutzeri, specifically strain 28a24 from this paper by Johannes Sikorski. I acquired a murder (not the right collective noun, but any microbiologist can sympathize) of strains from Johannes a few years ago because one of my interests is on the evolutionary effects of natural transformation in natural populations. At the time I had just started a postdoc in Jeff Dangl's lab working with P. syringae, saw the power of Pseudomonads as a system, and wanted to get my hands on naturally transformable strains. Didn't quite know what I was going to do with them at the time, but since I started my lab in Tucson this 28a23 strain has come in very handy as an evolutionary model system.

I received my Illumina read files last week from our sequencing center, and began the slog that can be assembly. While bacterial genome assembly is going to get much easier and exact in the next 5 years (see here) right now working with Illumina reads is kind of like cooking in that you start with a base dish that is seasoned to flavor. My dish of choice for working with Illumina reads is SOAPdenovo. Why you may ask? Frankly, most of the short read assemblers perform equally on bacterial genomes for Illumina PE. Way back when I started using Velvet as an assembler, but over time became frustrated with the implementation. I can't quite recall when it happened, but it might have been when Illumina came out with their GAII platform and the computer cluster I was working with at the time didn't have enough memory to actually assemble genomes with Velvet. Regardless...I mean no slight to Daniel Zerbino and the EMBL team, but at that point I went with SOAPdenovo and have been working with this since due to what I'd like to euphemistically call "research momentum".

Every time I've performed assemblies with 100bp PE reads, the best kmer size for the assemblies always falls around 47 and so I always start around that number and modify as necessary (In case you are curious, helpful tutorials for how De Bruijn graphs work can be found here). One more thing to mention, these reads have already been trimmed for quality. The first thing I noticed with this recent batch of data was that there was A HUGE AMOUNT of it, even when split into 24 different samples. The thing about bacterial genome assemblies is that they start to crap out somewhere above 70 or 80x coverage, basically because depth and random errors confuse the assemblers. If I run the assembler including all of the data, this is the result:


2251 Contigs of 2251 are over  bp -> 4674087 total bp
218 Contigs above 5000bp -> 1612110 total bp
30 Contigs above 10000 bp -> 352858 total bp
0 Contigs above 50000 bp -> 0 total bp
0 Contigs above 100,000 bp -> 0 total bp
0 Contigs above 200,000 bp -> 0 total bp
Mean Contig size = 2076.44913371835 bp
Largest Contig = 20517

The output is from a program Josie Reinhardt wrote when she was a grad student in Corbin Jones' lab (used again because of research momentum). The first line gives the total amount of contigs and scaffolds in the assembly, as well as the total assembly size. In this case, including all of the data yields 2251 total contigs, a total genome size of 4.67Mb, and a mean/largest contig of 2076/25,517bp. Not great, but I warned you about including everything.

Next, I'm going to lower the coverage. Judging by this assembly and other P. stutzeri genome sizes, 4.5Mb is right where this genome should be. For 70x or so coverage I am going to only want to include   approximately 1.6 million reads from each PE file ( (4,500,000*70) / 200). When I run the assembly this time, the output gets significantly better:


467 Contigs of 467 are over  bp -> 4695960 total bp
224 Contigs above 5000bp -> 4375593 total bp
158 Contigs above 10000 bp -> 3890697 total bp
11 Contigs above 50000 bp -> 709014 total bp
0 Contigs above 100,000 bp -> 0 total bp
0 Contigs above 200,000 bp -> 0 total bp
Mean Contig size = 10055.5888650964 bp
Largest Contig = 81785

I can show you more of this type of data with different amounts of coverage, but it's not going to change the outcome. The only other parameter to change a bit is kmer size. The figures above are for 47, but what happens if I use 49 with this reduced coverage?


382 Contigs of 382 are over  bp -> 4694856 total bp
178 Contigs above 5000bp -> 4429753 total bp
130 Contigs above 10000 bp -> 4065927 total bp
18 Contigs above 50000 bp -> 1317727 total bp
2 Contigs above 100,000 bp -> 211768 total bp
0 Contigs above 200,000 bp -> 0 total bp
Mean Contig size = 12290.1989528796 bp
Largest Contig = 110167

Even better (and the best assemblies I've gotten so far). Increasing kmer size to 51 makes it incrementally worse:


479 Contigs of 479 are over  bp -> 4633694 total bp
164 Contigs above 5000bp -> 4361629 total bp
124 Contigs above 10000 bp -> 4057109 total bp
19 Contigs above 50000 bp -> 1389519 total bp
1 Contigs above 100,000 bp -> 120161 total bp
0 Contigs above 200,000 bp -> 0 total bp
Mean Contig size = 9673.68267223382 bp
Largest Contig = 120161

Last question you might have, since there are other P. stutzeri genomes publicly available, is why not use those and carry out reference guided assembly?  I haven't really looked into this much yet, but with some quick analyses it doesn't seem like these genomes are really that similar to one another (~85% nucleotide identity, aside from many presence/absence polymorphisms) so it's not that easy to actually line them up by nucleotide sequence. Maybe better when I've got protein sequences, but that's for another post.

So now I'm ready to compare the brown phenotype genome vs. this assembled draft of 28a24. I know what the answer is, but I'm going to leave that until the next post.