Author: kat

Bacterial genomics researcher in Melbourne, Australia

Update to Comparative Bacterial Genomics tutorial

by David Edwards

In 2013, Kat and I wrote what turned out to be a very popular Beginner’s guide for comparative bacterial genome analysis. After four years and 120,000+ downloads of the guide, we thought it might be time to update the hands-on tutorial that was included. 

As with any science, there have been advances in this time. We don’t have time to update all aspects, but felt it was important to update the recommended assembler from Velvet to SPAdes. The latter has become the ‘go-to’ assembler with our lab and many others over the last few years. Unfortunately, SPAdes does not work with Windows, but Windows users can use the original Velvet assembler if they wish to attempt their own assembly.

Also, Ryan Wick in our lab has developed a way to visualise the assembly graphs produced by SPAdes and other assemblers, in the form of a software program called Bandage. This allows us to examine and compare the properties of assembly graphs, useful if you are trying to assemble the same set of reads with different methods or parameter settings. 

The other changes in version 2 are mainly to fix broken links to the E. coli sequences that have now been archived by NCBI, kindly pointed out to us by Michael Hall and others via email.

We continue to recommend Artemis and ACT for visualising and comparing annotated bacterial genome sequences, and both tools are still actively maintained at the Sanger Institute. While BRIG is no longer actively maintained, we continue to recommend it as it appears to be stable across newer versions of Java and BLAST, and it remains incredibly useful.

Hands-on tutorial v2 (6 Mb PDF): ComparativeGenomicsTutorialV2

Original article: Beginner’s guide to comparative bacterial genome analysis using next-generation sequence data

Advertisements

Global genomic framework for typhoid

It’s been over a year since we published the first global whole-genome snapshot of nearly 2000 genomes of the typhoid bacterium, Salmonella Typhi in Nature Genetics.

That paper focused on the emergence and global dissemination of what we’ve been calling for years the “H58” clone (see this blog post). This clone accounted for nearly half of all the isolates sequenced, and is a big deal because it tends to be multidrug resistant (MDR), carrying a suite of resistance genes that render all the cheap, first-line drugs like chloramphenicol, ampicillin, and trimethoprim-sufamethoxazole useless for treatment. Detailed genomic epi studies show the local impact of the arrival of MDR H58 in countries as widespread as Malawi and Cambodia; and the emergence of fluoroquinolone resistant H58 sublineage in India and Nepal recently stopped a treatment trial because the current standard of care – ciprofloxacin – was resulting in frequent treatment failure.

While H58 is important, the global Typhi population contains a lot of genomic diversity outside the H58 clone, and we’ve turned our attention to the rest of the population now in a new paper in Nature Communications: “An extended genotyping framework for Salmonella enterica serovar Typhi, the cause of human typhoid

First, we decided that we needed to revisit the haplotyping scheme of Roumagnac et al (from which H58 gets its name), which was based on just ~80 genes, using the whole genome phylogeny. Here is the tree inferred from core genome SNPs in 1832 Typhi strains, with the old haplotypes indicated by the coloured ring around the outside. It’s pretty easy to see that some haplotypes (like H52 and H1) actually comprise multiple distinct phylogenetic lineages (low resolution), while others subdivide lineages (excessive resolution).

FigS1_GlobalTreeColouredByHgroups

Whole genome SNP tree for 1832 strains, outer ring indicates haplotypes based on mutations in 80 genes as defined in Roumagnac et al, Science 2006.

We used BAPS to define genetic clusters at various levels (thanks to Tom Connor for running this). We settled on 3 levels of hierarchical clustering, indicated in the tree below:

• 4 nested primary clusters (inner-most ring; yellow, green, blue, red). These have 100% bootstrap support and are each characterised by >20 SNPs

• Clusters are further divided into 16 clades (middle ring and labels). The median pairwise distance between isolates in the same clade is 109 SNPs, while the inter-clade SNP distance averages 243 SNPs.

• Clades are further divided into 49 subclades, indicated by alternating background shading colours. The median pairwise distance between isolates in the same subclade is 25 SNPs.

Fig1_noH58_clade_colours_subcladebg_190515_labelled_ED.png

Tree indicating new phylo-informed genotypes. Primary clusters 1-4 are indicated in the inner ring. Branch colours indicate clades, which are also labelled on the outside and coloured in the outer ring. Subclades are indicated by alternating background shading.

subcladerectOne of the key reasons we wanted to define the phylogenetic lineages in this way is to make them easier to identify and talk about. I’ve always been a fan of MLST for this reason, since it’s much easier to talk about K. pneumoniae ST258, ST11, ST15 etc than ‘that lineage that has reference strain X in it’. So we introduce a hierarchical nomenclature system, similar to the one currently in use for Mycobacterium tuberculosis, where the 4 primary clusters (1, 2, 3, 4) are subdivided into 16 clades (1.1, 1.2; 2.1, 2.2, etc) which in turn are subdivided into 49 subclades (1.1.1, 1.1.2, etc). This has the advantage of conveying hierarchical relationships between groups – e.g. 2.2.1 and 2.2.2 are sister subclades within clade 2.2, which is a sister clade of 2.1.

The subclades are easier to distinguish in the collapsed rectangular tree on the right, where each subclade is represented by just one strain.

Some BAPS clusters were polyphyletic and consisted of isolates belonging to rare phylogenetic lineages whose common ancestor in the tree coincided with the common ancestor of an entire clade (n=9) or primary cluster (n=2). These groups contain isolates that, given increased numbers, may emerge as distinct clusters that form sister taxa within the parent clade (or primary cluster), and were given the suffix ‘.0’ rather than a defined cluster number (e.g. 3.0 or 3.1.0) to indicate non-equivalence with the properly differentiated sister clades (n=16) or subclades (n=49). As more genomes are added, these are expected to be more clearly differented into distinct groups and given proper clade/subclade designations.

Next we defined a set of 68 SNPs that can be used to genotype isolates into these groups. We chose one SNP for each primary cluster, clade and subclade (preferentially choosing intragenic SNPs in well-conserved core genes). The SNPs are detailed in a supplementary spreadsheet, and we provide a script to assign strains to genotypes based on an input BAM or VCF file generated by mapping to the reference genome for Typhi strain CT18.

An isolate that belongs to a differentiated subclade such as 2.1.4 will be hierarchically identified by carrying the derived allele for primary cluster 2 (but not the nested clusters 3 and 4); the derived allele for clade 2.1 (but no other clades) and the derived allele for subclade 2.1.4 (but no other subclades). It is possible for an isolate to carry derived alleles for a primary cluster and clade with no further differentiation into subclade.

The clone formally known as H58
Under the new scheme, the infamous H58 clone is named subclade 4.3.1, which so far has no sister clades. I suspect those of us familiar with Typhi population genomics will keep referring to it informally as H58 for some time, since that name is now well known… but I will try to re-train myself to call it 4.3.1 (H58).

Now the fun part: exploring the geographical distribution of these lineages.

Fig1c_worldmap_pies_subclades

Figure 1c from the paper. Pie colours indicate clades found in each WHO region in the global data set (key is in the tree figure above).

In the paper we go on to show that:
• clades are widely geographically distributed, while subclades are geographically constrained (see heatmap below);
• genotyping can be used to predict the geographical origin of travel-associated typhoid in patients in London;
• even better predictions can be obtained based on genome-wide SNP distances to our reference panel of >1800 isolates… but of course that involves a lot more computationally intensive comparisons than a quick screen of a new isolate’s BAM file.

Screen Shot 2016-08-29 at 11.47.15 pm.png

Figure 2 of the paper, showing the geographical distribution of subclades, which shows most subclades are restricted to a single region. For this analysis, the effect of local outbreaks has been minimised by replacing groups of strains that share the same subclade and year and country of isolation with a single representative strain.

You can read the full details in the paper, but here I just want to highlight that you can now explore the global genomic framework for Typhi – including genotype designations as well as temporal and geographic data – interactively in MicroReact.

tree.png

 

How does genotyping help with studying local populations?

We have already begun using the new genotyping scheme in local typhoid studies. I find this a really helpful way to describe/summarise the local populations, and place them in the context of the global population without resorting to large trees.

For example in this recent Nigerian study, we described the population like this: “The majority of isolates (84/128, 66%) belonged to genotype 3.1.1 , which is relatively common across Africa, predominantly western and central countries. In the wider African collection genotype 3.1.1 was represented by isolates from neighbouring Cameroon and across West Africa (Benin, Togo, Ivory Coast, Burkina Faso, Mali, Guinea and Mauritania) suggesting long-term inter-country exchange within the region. Most of the remaining isolates belonged to four other genotypes (4.1, 2.2, 2.3.1 and 0.0.3).”

Of course genotype assignment is not the end of the story – we still want to build whole-genome trees to explore the relationships of local isolates with those from other countries. Importantly, working with genotypes means that we can achieve this without needing to build a megatree of all isolates in the local + global collections (n>2000). Instead, we can use the genotypes to identify which strains from the global collection are relatives of the Nigerian isolates, and build a much smaller tree that still captures all of the information about transmission/transfer between Nigeria and other countries:

journal-pntd-0004781-g001

The tree and map were made using MicroReact, you can recreate theme here: http://microreact.org/project/styphi_nigeria To get this colour scheme just click on the eye icon (bottom left) and select ‘country’; and to get the fan style tree, click the settings button (top right) and click the fan shape.

Another example is in our recent paper on isolates collected in Thailand before and after the introduction of their national vaccination program (pre-print here):

  • Genotype 3.2.1 was the most common (n=14, 32%), followed by genotype 2.1.7 (n=10, 23%)
  • Genotypes 2.0 (n=1, 2%) and 4.1 (n=3, 7%) were observed only in 1973 (pre-vaccine period)
  • Genotypes 2.1.7 (n=10, 23%), 2.3.4 (n=1, 2%), 3.4.0 (n=2, 5%), 3.0.0 (n=3, 7%), 3.1.2 (n=2, 5%), were observed only after 1981 (post-vaccine period)
  • Genotypes 3.2.1 and 2.4.0 were observed amongst both pre- and post-vaccine isolates, but the subclade phylogenies show that these more likely to represent re-introduction of strains from neighbouring countries than persistence within Thailand throughout the immunisation program.

Elizabethkingia anophelis

DATA: raw Illumina reads (CDC) under SRP072035

Assemblies & analyseshttps://github.com/katholt/elizabethkingia

SpeciesTree_OutbreakTree2

Left – core SNP tree, created from assembled genomes using Parsnp. Right – core SNP tree, created by mapping all outbreak genomes to our 9-contig assembly for SRR3240412, using our RedDog pipeline. Details below; all assemblies and mapping outputs are here.


March 19, 2016: I saw on twitter today that there was an outbreak of a weird bacteria I’ve never heard of before (Elizabethkingia anophelis) in Wisconsin, which had infected >50 people and killed almost 20.

The Wisconsin Health Services department has posted some information here (click the “For Health Professionals” tab to get info on the bacteria and antibiotic resistance).

CDC has deposited Illumina reads from 18 outbreak strains into SRA under project SRP072035 so I pulled the data and had a look. I managed to download the readsets in a few minutes (using bionode-ncbi) but it took a really long time to unpack these into fastq files using sra-toolkit.

As I have no idea about this species, I thought I’d start by looking for antibiotic resistance and plasmids in the first 6 read sets using our SRST2 software, while waiting for the rest of the reads to unpack… this ran quickly and showed me the same results for all 6 strains: GOB-10 (1 SNP from the closest allele in the ARG database) and B-2 (3 SNPs), at depths of 35-65x. For example:

SRR3240397 ARGannot.r1 GOB-1_Bla GOB-10_821 100.0 47.59 1snp 0.115 873 0.043 188 821 no;yes;GOB-10;Bla;AY647247;1-873;873
SRR3240397 ARGannot.r1 B-1_Bla B-2_1160 100.0 59.891 3snp 0.4 750 0.03 314 1160 no;no;B-2;Bla;AF189300;1-750;750

Because the matches were not identical, I pulled the consensus sequences (based on read mapping) using –report_all_consensus option in SRST2:

>314__B-1_Bla__B-2__1160 no;no;B-2;Bla;AF189300;1-750;750
ATGTTGAAAAAAATAAAAATAAGCTTGATTCTTGCTCTTGGGCTTACCAGTCTGCAGGCA
TTTGGACAGGAGAATCCTGACGTTAAAATTGATAAGCTAAAAGATAATCTGTATGTATAC
ACAACCTACAATACATTTAACGGGACTAAATATGCCGCTAATGCAGTATATCTGGTAACT
GATAAGGGTGTTGTGGTTATAGACTGTCCGTGGGGAGAAGACAAATTTAAAAGCTTTACG
GACGAGATTTATAAAAAACACGGAAAGAAAGTTATTATGAATATTGCAACACATTCTCAT
GATGATCGTGCCGGAGGTCTTGAATATTTTGGTAAAATAGGTGCAAAAACTTATTCTACT
AAAATGACAGATTCTATTTTAGCAAAAGAGAATAAGCCAAGAGCACAATATACTTTTGAC
AATAATAAATCTTTCAAAGTAGGAAAATCCGAGTTTCAGGTTTACTATCCCGGAAAAGGA
CATACAGCAGATAATGTGGTGGTATGGTTTCCAAAAGAAAAAGTATTGGTTGGAGGTTGT
ATTATAAAAAGCGCTGATTCAAAAGACCTGGGGTATATTGGAGAAGCATATGTAAACGAC
TGGACGCAGTCTGTACACAATATTCAACAAAAGTTTTCCGGTGCTCAGTACGTTGTTGCA
GGGCATGATGATTGGAAAGATCAAAGATCAATACAACGTACACTAGACTTAATCAATGAA
TATCAACAAAAACAAAAGGCTTCAAATTAA
>188__GOB-1_Bla__GOB-10__821 no;yes;GOB-10;Bla;AY647247;1-873;873
ATGAGAAATTTTGTTATACTGTTTTTCATGTTCATTTGCTTGGGCTTGAATGCTCAGGTA
GTAAAAGAACCTGAAAATATGCCCAAAGAATGGAACCAGACTTATGAACCCTTCAGAATT
GCAGGTAATTTATATTACGTAGGAACCTATGATTTGGCTTCTTACCTTATTGTGACAGAC
AAAGGCAATATTCTCATTAATACAGGAACGGCAGAATCGCTTCCAATAATAAAAGCAAAT
ATCCAAAAGCTCGGGTTTAATTATAAAGACATTAAGATCTTGCTGCTTACTCAGGCTCAC
TACGACCATACAGGTGCATTACAAGATCTTAAAACAGAAACCGGTGCAAAATTCTATGCC
GATAAAGAAGATGCTGATGTCCTGAGAACAGGGGGGAAGTCCGATTATGAAATGGGAAAA
TATGGGGTGACATTTAAACCTGTTACTCCGGATAAAACATTGAAAGATCAGGATAAAATA
ACACTGGGAAATACAATCCTGACTTTGCTTCATCATCCCGGACATACAAAAGGTTCCTGT
AGTTTTATTTTTGAAACAAAAGACGAGAAGAGAAAATATAGAGTTTTGATAGCTAATATG
CCCTCCGTTATTGTTGATAAGAAATTTTCTGAAGTTACCGCATATCCAAATATTCAGTCC
GATTATGCATATACTTTCAAAGCAATGAAGAATCTGGATTTTGATATTTGGGTGGCCTCC
CATGCAAGTCAGTTCGATCTCCATGAAAAACGTAAAGAAGGAGATCCGTACAATCCGCAA
TTGTTTATGGATAAGCAAAGCTATTTCCAAAACCTTAATGATTTGGAAAAAAGCTATCTC
GACAAAATAAAAAAAGATTCCCAAGATAAATAA

I had a quick look at the accessions for these genes (they are hiding in the fasta headers above) and found that they are carbapenemase genes reported from Elizabethkingia meningosepticum (previously called Chryseobacterium meningosepticum), reported in these papers: Bellais 2000, Antimicrob. Agents Chemother and Yum 2010, J Microbiol. These genes confer resistance to carbapenems like meropenem and imipenem, which probably contributes to these bacteria causing hospital-acquired infections as they will be selected for by carbapenem exposure.

SRST2 didn’t find any other acquired antibiotic resistance genes (from the ARG-Annot database) or known plasmid replicons (at least those in the PlasmidFinder database), which is consistent with the Wisconsin health services reports that these strains are susceptible to lots of readily accessible drugs including fluoroquinolones, rifampin and trimethoprim/sulfamethoxazole.

Running a quick NCBI BLAST search of the carbapenemase gene sequences shows that these new sequences, which are from outbreak strains identified definitively as Elizabethkingia anophelis by CDC, are closest to sequences annotated in NCBI as originating from Elizabethkingia meningosepticum. (The trees below are just straight out NCBI BLAST, obtained by clicking “Distance tree of results” and then downloading the newick tree files to view in FigTree.)

Screen Shot 2016-03-19 at 4.21.25 pm

The outbreak strain’s GOB-10 gene had 1 synonymous SNP compared to the reference sequence, while the B-2 gene had 1 synonymous and 2 non-synonymous SNPs (affecting codons 31 & 34, which is outside the beta-lactamase domain).

I am guessing that species assignations are pretty tricky for this genus, as few labs will have access to definitive tests to discern them, so we shouldn’t read much into this. However if it is true that the outbreak strains are Elizabethkingia anophelis and the close-matching genes in NCBI did come from Elizabethkingia meningosepticum, this would suggest that there is horizontal gene transfer between these species.

Note 1: while writing this, the fastqs finished extracting and I ran SRST2 and found the same antibiotic resistance gene results on all 18. I’m now running some SPAdes assemblies which I’ll post here later, to save others the trouble…

Note 2: the assemblies (SPAdes fasta and fastg; plus Prokka annotated in GenBank format), and various analyses including trees created using Parsnp (from assemblies) and our RedDog pipeline (mapping of reads to reference genome strain NUHP1 =CP007547) are here in github: https://github.com/katholt/elizabethkingia

The assemblies are a bit variable, but mostly ~3.9 Mbp (the reference is 4,369,828) but the best one was for SRR3240413 – 32 contigs with 3,911,053 bp total. Viewing the SPAdes assembly graph in our Bandage program shows that 3,910,660 bp are in a single linked graph, which corresponds to the chromosome. (The other little bits do not look like plasmids, just leftover bits of sequence and probably adapters, that SPAdes spit out in teeny bits of a few hundred bp each.)

SRR3240413_bandage_graph

 

The genomes look pretty similar at first glance, but interestingly 4 of them share a deletion of ~80 genes. That’s a great little epidemiological marker for the investigations.

Screen Shot 2016-03-19 at 7.24.21 pm

This was detected by our mapping pipeline RedDog, which I used to map the reads to reference genome NUHP1 CP007547 (this may not be the best reference, I just picked one randomly). The assemblies confirm it: genes BD94_0888 to BD94_0962, and the end of BD94_0963, are missing in these 4 strains (although reads do map to BD94_0948, because this is present in a second copy elsewhere in the genome).

Here’s that tree with a bit more detail (red = # core SNPs). The tree was made using FastTree, with NUHP1 reference genome as an outgroup (the outbreak strains are >40,000 SNPs away from this reference).

Screen Shot 2016-03-19 at 8.26.50 pm


UPDATE 20/3

David Edwards has been playing with Hybrid StriDe + SPAdes assembly recently, and tried this with SRR3240412. The SPAdes assembly (here) was 3,913,666 bp in 41 contigs. The hybrid assembly (here) is 3,917,367 bp in 9 contigs. This is what the graph looks like… I’m showing it coloured by BLAST matches to the NUPH1 reference genome so you can see what the likely path through the graph is (rainbow, red -> purple, indicates matches from start -> end coordinates of the reference genome).

scaffolds_resolved_blast2.png

March 22: Hybrid assemblies for all 18 outbreak genomes (9-15 contigs each; ranging 3,830,044 – 3,912,928 bp) are now in GitHub (thanks to David Edwards for this).

David re-ran the RedDog mapping pipeline using this genome assembly as the reference, and got a very similar tree (files in GitHub):

SRR3240412_RedDogTree

And here is a core genome SNP tree (made from genome assemblies using Harvest), which shows the outbreak strains are a novel lineage of E. anophelis, compared to currently available data (tree file in GitHub):

species_tree_with_NCTC10588

Note: Sylvain Brisse has shown the same thing using his core genome MLST scheme (see BioRxiv preprint posted March 19). I have used his nomenclature here (lineage A, lineage B). Note that Lineage B strains were originally identified as E. meningoseptica, but belong to E. anophelis. I have also included here the genome of E.meningoseptica NCTC10588, which was sequenced as part of the Sanger/PHE/PacBio type strain project, in this tree… it clusters within lineage A and is clearly E. anophelis.

UPDATE: This tweet-fest led to a collaborative project between myself, Sylvain Brisse (Institut Pasteur) and CDC, which was eventually published in Nature Communications a year after this blog post: “Evolutionary dynamics and genomic features of the Elizabethkingia anophelis 2015 to 2016 Wisconsin outbreak strain

 

Genomics of atypical enteropathogenic E. coli

Well our paper on “atypical” enteropathogenic E. coli is finally out, in the new journal Nature Microbiology. The hard work on this paper was done by Danielle Ingle, whom I co-supervise together with Roy Robins-Browne and will soon become the first of my students to complete a PhD. Unfortunately the paper is paywalled at the moment, so I’ll try to summarise the key elements here.

For those who aren’t familiar with the wacky world of E. coli pathotypes, the LEE genomic island is the darling of E. coli pathogenicity genetics. It encodes a type three secretion system, essentially a syringe structure that enables the bacteria to inject its own proteins (known as secreted effector proteins) directly into human cells and – ahem – mess them up (for more technical details see the impressive body of work from Melbourne University’s very own Jacqueline Pearson and Liz Hartland). E. coli that possess shiga toxin as well as the LEE cause nasty bloody diarrhoea and are defined as enterohaemorrhagic E. coli (EHEC), while those with LEE + the bfp adhesion factor cause diarrhoea and are defined as enteropathogenic E. coli (EPEC). Strains carrying the LEE but lacking both shiga toxin and bfp can also cause diarrhoea, but are also sometimes isolated from healthy individuals with no symptoms of diarrhoaea, and it is this big group of LEE + ??? strains that are defined as atypical EPEC (aEPEC), while those with bfp are dubbed typical EPEC.

No one has been able to find a common virulence factor that differentiates diarrhoeagenic aEPEC from asymptomatic strains… which is not that surprising considering that aEPEC is really just an umbrella term for a diverse group of organisms that happen to carry the LEE genomic island.

In this paper, we set out to define the population structure of aEPEC strains using genomics, and also to investigate variation in the LEE island itself. Our strain set is ~200 newly sequenced aEPEC that were isolated from children with diarrhoea, and from asymptomatic age-matched controls, as part of the Global Enteric Multicenter Study (GEMS) – a massive study into the etiology of childhood diarrhoea across seven sites in Africa and Asia, funded by the Gates Foundation. We also included lots of publicly available genomes from NCBI, which included ~60 additional aEPEC.

So what did we find? Well from a population structure point of view, we confirmed what we suspected from the beginning – that the ~200 aEPEC strains actually represent dozens of distinct lineages or clonal groups within the E. coli whole genome phylogeny. We tried making the core genome tree in a few different ways, including mapping reads to a reference genome vs using Roary to generate core gene alignments from assemblies; with and without removing recombinant regions identified using ClonalFrameML. The alternative trees all tell the same population structure story, as you can see below. An interactive version of the mapping-based, recombination-filtered tree (which appears in Figure 1 of the paper, panel a below) is available to play with in MicroReact.

 

Supplementary Figure 1

The 10 clades highlighted in the tree are those containing >5 aEPEC in our collection, which represent the most common aEPEC lineages. The figures below show that these 10 aEPEC groups are present across the Asian and African GEMS sites; most also appear in non-GEMS collections from Europe as well as North and South America, indicating they are globally distributed.

aEPEC_distribution

Next we looked in depth at the LEE genomic island itself. ClonalFrame analysis identified lots of recombination across the locus, particularly affecting the proteins that make direct contact with host cells:

Screen Shot 2016-01-19 at 9.46.45 am

Danielle also did a lot of work to examine selection pressures on the various LEE genes, and found evidence for co-evolutionary constraints on groups of LEE genes (see the paper for details).

We used the recombination-free ClonalFrame tree, which represents the underlying clonal or vertical evolutionary relationships between LEE sequences, to define LEE lineages and subtypes. We also made a MLST-style scheme for the LEE, to make it easy to identify types from sequence data without having to compare to reference trees. The sequence files are included in the SRST2 distribution and can be downloaded from https://github.com/katholt/srst2.

Screen Shot 2016-01-19 at 9.50.41 am

There were three clearly distinct LEE lineages, separated from one other by 4-7% divergence (similar to genus-level differences between bacteria) and with different preferences for chromosomal insertion site. One of the LEE lineages (lineage 1) was entirely novel and we found it only in aEPEC isolates from the GEMS study. Functionality of the LEE was confirmed in all GEMS strains, including these novel lineage 1 examples.

Importantly, different aEPEC lineages had distinct LEE variants, often integrated at different chromosomal sites, confirming that aEPEC lineages have evolved numerous times independently via distinct LEE acquisition events. Also, some aEPEC lineages contained strains carrying different LEE lineages and/or subtypes, sometimes integrated at different sites, indicating ongoing evolution within some aEPEC clades.

aEPEC_tree_LEE_subtypes

Unsurprisingly, some of the “aEPEC” lineages also include examples of EHEC and/or typical EPEC isolates, i.e. subclades within the LEE-containing lineage that have subsequently acquired shiga toxin (via phage) or a bfp plasmid. We also found extensive variation in the complement of known LEE-secreted effectors (which are scattered around the genome and known as non-LEE encoded effectors or nle genes), both within and between lineages.

There’s lots more work to do to unravel the genetic basis for pathogenicity in E. coli, including identifying more of the colonisation factors, regulators and effectors that are important for causing disease in humans. But I think this analysis is a great step towards clarifying the population structure of EPEC and the LEE more generally, which opens the door for much more informative diagnostics and epidemiological studies. This work should also help immensely in guiding future pathogenicity research… most research on EPEC or EHEC pathogenicity has so far focused on a handful of strains, but our genomic analysis shows there is a much wider diversity of strain backgrounds, LEE variants and effector gene content that also needs to be considered.

Microbial Genomics methods

MGen (Microbial Genomics) was created last year by the UK’s Microbiology Society, with the aim of becoming the go-to journal for microbial genomics research. As a Senior Editor, I was asked to help mark the occasion of MGen’s six month anniversary (on Jan 15, 2016), by reviewing the 24 papers published in the journal’s first six months since launch.
The full review is is available over on the MGen site, but I wanted to draw particular attention here to trends in the analysis tools used by the genomic epidemiology papers published so far in MGen, which largely reflects what’s been happening across the field generally.

tools word cloud

One third of the articles are cutting edge genomic epidemiologyin action, using genomics to investigate the evolution and transmission of a range of pathogens, from anthrax to dysentery and food poisoning.

Interestingly all of these studies used similar methodologies, reflecting the maturity of this field of research.

Common approaches include: (i) sequencing large numbers of isolates using high-throughput Illumina platforms; (ii) the identification of SNPs (single nucleotide polymorphisms) using read mapping approaches (with BWA, SMALT, SAMtools and GATK being popular tools); and (iii) uniform use of RAxML for generating maximum likelihood phylogenies.

Some also used BEAST to estimate mutation rates, divergence dates and phylogeographical patterns. Interestingly, half of these papers utilised MLST (multi-locus sequence typing) to identify clades, showing that this sub-genomic approach based on capillary sequencing of ~7 gene fragments is still considered useful by many genomic researchers. [NOTE: nice to see that most who were inferring MLST from Illumina data were using our SRST or SRST2 software.]

Almost all of the genomic epidemiology studies took steps to remove SNPs introduced via recombination, in order to capture the underlying signals of vertical inheritance that are so important for transmission studies. Popular tools were BratNextGen, Gubbins and ClonalFrameML, which were all published within the last 3 years.

For pan genome analysis, Velvet and SPAdes were the most popular tools for bacterial genome assembly, with Prokka and Prodigal for gene annotation, and LS-BSR and related approaches being commonly used to cluster orthologous groups of genes.

 

Most of this won’t be any surprise to people working in bacterial genomic epi, but I think it’s great to see that consensus is emerging on how best to do these sorts of analyses, and to at last have some reliable tools for detecting and accounting for recombination.

The area of least agreement remains SNP calling – which mainly comes down to which read mappers to use, and which SNP calling algorithms and filtering to go with? This is a complex area, as highlighted in the recent review “Best practices for evaluating single nucleotide variant calling methods for microbial genomics” in Frontiers in Genetics, which did a very thorough job of examining the issues that need to be considered, but (quite deliberately) doesn’t provide an answer to “which tool is best?”.

Although there is still no real consensus on exact methods for SNP calling, I think most of the tools people are using (ie a good, stable read mapper followed by SNP calling with an established tool like SAMtools or GATK, with some basic filtering to remove low-evidence or ambiguous calls) end up with very similar answers (as we saw with the NGS outbreak analysis challenge session held at the ASM NGS meeting in September 2015). All in all it seems to me that the use of genomics for public health & diagnostic microbiology is in far better shape in this respect than clinical human genomics, which is going through something of a crisis involving wide discrepancies in variant calling as well as uncertainty around data interpretation.