Author Archives: John Davey

It is an exciting time to be working on Heliconius butterflies, as new sequencing technologies allow us to get closer and closer to the genes underlying colour pattern variation. We’re going to start featuring new Heliconius papers on the blog, to highlight new findings and put them into context. In this post, we’re going to look at a new paper in Genome Research that describes the region of the genome that produces red patterns in Heliconius erato and reveals how this region varies between H. erato and its mimic, Heliconius melpomene [1].

Red bands and rays are one of the key features of Heliconius wing patterns and occur repeatedly in mimics across South America. For example, Heliconius melpomene and Heliconius erato mimic each other in many different countries, often sharing either bands or rays (see picture above [2]). For over a century, researchers have been trying to understand how this mimicry could evolve. One of the major parts of this puzzle is to discover the genes that produce the red patterns and how they vary between species.

It has been known for many years that these features are controlled by one large region of the genome less than one megabase long, known as the B locus in H. melpomene and the D locus in H. erato [3]. In recent years, we have been making great progress in describing and understanding this region. In 2010, the whole region was sequenced using bacterial artificial chromosomes (BACs) in H. melpomene and H. erato [4,5]; the region is around 700 kilobases long in H. melpomene and 1 megabase long in H. erato, and contains around 20 genes, which are not only the same genes in both species but are also conserved in the same order.

We now know that the optix gene, which is found within the B/D region, is expressed on the wing wherever red pattern elements are expressed across many different Heliconius species, and that genetic variations in optix are strongly associated with variations in red patterns in multiple species [6]. We also know that this gene has evolved to produce red patterns separately in H. melpomene and H. erato, converging on the same phenotype but differing genetically [2].

We now really want to know what’s happening across the whole B/D region. Are there any other genes strongly associated with red band patterns? What genetic variations in H. melpomene and H. erato exist in the wild? To answer these questions, we need to sequence the whole B/D region in many butterflies from each species. Unfortunately, until recently it has not been possible to do this. BAC sequencing of one copy of the whole B/D region was very laborious and expensive; it would be very difficult to sequence tens of butterflies this way. A short 800 bp fragment of optix was sequenced in many butterflies using polymerase chain reactions (PCRs) to associate optix with red patterns, but over 1000 similar experiments would be required to cover the whole region.

With next generation sequencing, we can now sequence whole genomes of many butterflies and use these sequences to study particular regions of the genome like B/D. The new paper in Genome Research by Megan Supple and colleagues reports on the sequencing of short reads from the whole genomes of 45 H. erato butterflies from four hybrid zones across South America using next generation sequencing. By taking many butterflies with rayed patterns and red banded (‘postman’) patterns, and looking for variations in the D locus sequences, we can see what’s going on across the whole region.

Next generation sequencing data is still fairly new, and as with many scientific projects the scientists usually learn how to analyse their data during the project. As Megan says, “I was just handed a large pile of data. At that point, I had no idea how to go about analyzing next-gen sequences. It was quite a learning process.”

It is not yet possible to sequence whole chromosomes in one go; current technologies can only produce short reads a few hundred bases long. The raw sequence data for the 45 H. erato genomes is hundreds of millions of 100 base pair (bp) sequence reads. The H. melpomene genome has been assembled, but the H. erato genome has not, so the first step of the analysis was to align the 100bp reads to the existing BAC sequence of the erato D region. By comparing reads between different butterflies at the same locations, variations in the D locus could be identified. This is done with tools like the Genome Analysis Tool Kit, a suite of software written to analyse human genome data sets such as the 1000 Genomes Project data, but as DNA is the same in butterflies and humans, we can use the same software to study the evolution of mimicry.

Megan found the same strong association with optix shown in previous papers, but also found a very strong association in a 65 kb long region upstream of optix, which does not appear to contain any genes. This strongly suggests that there are sequences that regulate other genes in this region, which are responsible for turning the genes with red patterns on and off. In the H. melpomene genome paper, this same region was shown to vary considerably between H. melpomene subspecies [7].

“What surprised me the most was that analyzing H. erato and H. melpomene separately identified almost the same exact region. I did not expect that the boundaries identified in H. melpomene would almost perfectly coincide with the boundaries in H. erato. That indicates that we have hit the boundary of something important – we believe that those boundaries are close to important functional regions that are under strong divergent selection.”

So are these regions actually the same in both species? Apparently not. In H. erato, there are 76 single base variations with one variant perfectly associated with the postman pattern and the other associated with rays; in H. melpomene, there are 430 such variations. When the sequences of the whole regions are compared, the sequences from H. melpomene individuals with rayed patterns do not group together with sequences from H. erato rayed individuals, and H. melpomene postman individuals do not group with H. erato postman individuals.

This strongly indicates that, although the H. melpomene and H. erato butterflies look very similar, the genetic sequences producing the patterns are not the same, and H. melpomene has evolved the pattern independently from H. erato. Megan Supple: “It was really compelling to see how perfectly the sequences in this region cluster by pattern within each species, but always keep the two mimics separate.”

Now the architecture of the B/D region is clear, and the full sequences and variations are available, we can search for the functional variants that produce the different red patterns. This will take some time, as there are 65,000 bases to search through, and confirming the function of any candidate sequences will require some difficult experimental work. But this new paper brings us very close to finally identifying these variants, and to providing a mechanical explanation for the evolution of red pattern mimicry in Heliconius.

1. Supple M.A., Hines H.M., Dasmahapatra K.K., Lewis J.J., Nielsen D.M., Lavoie C., Ray D.A., Salazar C., McMillan W.O., Counterman B.A. (2013). Genomic architecture of adaptive colour pattern divergence and convergence in Heliconius butterflies. Genome Research, doi:10.1101/gr.150615.112.

2. Hines H.M., Counterman B.A., Papa R., de Moura P.A., Cardoso M.Z., Linares M., Mallet J., Reed R.D., Jiggins C.D., Kronforst M.R., McMillan W.O. (2011). Wing patterning gene redefines the mimetic history of Heliconius butterflies. Proceedings of the National Academy of Sciences 108(49):19666-19671.

3. Sheppard P.M., Turner J.R.G., Brown K.S., Benson W.W., Singer M.C. (1985). Genetics and the Evolution of Muellerian Mimicry in Heliconius Butterflies. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 308(1137):443-610.

4. Baxter S.W., Nadeau N.J., Maroja L.S., Wilkinson P., Counterman B.A., Dawson A., Beltran M., Perez-Espona S., Chamberlain N., Ferguson L., Clark R., Davidson C., Glithero R., Mallet J., McMillan W.O., Kronforst M., Joron M., ffrench-Constant R.H., Jiggins C.D. (2010). Genomic Hotspots for Adaptation: The Population Genetics of Muellerian Mimicry in the Heliconius melpomene Clade. PLoS Genetics 6(2):e1000794.

5. Counterman B.A., Araujo-Perez F., Hines H.M., Baxter S.W., Morrison C.M., Lindstrom D.P., Papa R., Ferguson L., Joron M., ffrench-Constant R.H., Smith C.P., Nielsen D.M., Chen R., Jiggins C.D., Reed R.D., Halder G., Mallet J., McMillan W.O. (2010). Genomic Hotspots for Adaptation: The Population Genetics of Muellerian Mimicry in Heliconius erato. PLoS Genetics 6(2):e1000796.

6. Reed R.D., Papa R., Martin A., Hines H.M., Counterman B.A., Pardo-Diaz C., Jiggins C.D., Chamberlain N.L., Kronforst M.R., Chen R., Halder G., Nijhout H.F., McMillan W.O. (2011). optix Drives the Repeated Convergent Evolution of Butterfly Wing Pattern Mimicry. Science 333:1137-1141.

7. Heliconius Genome Consortium (2012). Butterfly genome reveals promiscuous exchange of mimicry adaptations among species. Nature 487:94-98.

In the last few years, the Heliconius community has generated a large quantity of next generation sequencing data, mostly from Illumina sequencers. One of the first steps in the analysis of this data is the alignment of short reads to reference sequences, whether BAC sequences, transcriptomes or whole genomes. There are over 70 read alignment tools available (Fonseca 2012), which has been thought to be too many. Which one(s) should we use for Heliconius butterfly data?

Several popular, fast aligners such as BWA and Bowtie are known to work well on human genomes. But Heliconius genomes, radiating over around 10 million years, are considerably more divergent than human genomes. In July 2011, I ran some aligner tests to see what the impact of divergence was on aligner performance. In those tests, Stampy performed considerably better than other aligners. Since then, we have mostly used Stampy for read alignment, notably for (I think) all of the alignments in the Heliconius genome + introgression paper.

Stampy has a serious drawback, which is that it is extremely slow compared to almost all other aligners – four or five times slower than Bowtie. This has been very inconvenient in the past and is only getting worse as data set size increases. I am about to analyse 1.8 billion reads from 70 whole genomes. Aligning this amount of data with Stampy will take many weeks even over many tens of compute cores. There have been several new aligner tools released in the time since I did the last test, so I thought it was time for another comparison.

Be warned, this is a very basic comparison done over a few days, with many flaws. If it was good enough to publish properly I would do so. It really only gets the ball rolling again, but I don’t really have the data to come to any firm conclusions, just to outline what the problems are.

Sequence Data

Three whole genome sequences were used, one from a Heliconius melpomene melpomene x Heliconius melpomene rosina cross, one from Heliconius ethilla (09-49) and one from Heliconius erato erato (NCS_2556). These three genomes should be increasingly divergent from the Heliconius melpomene melopmene reference genome (see tree at Tree of Life). Each sample has 4m paired end 100bp HiSeq reads or just over. This is only about 3 times coverage of the genome, so we don’t expect to see every base covered or very accurate SNPs called. But using such limited data means the aligners can be tested relatively quickly. (I leave for another day the serious issue of what we might expect from an alignment of H. erato and H. melpomene, and how much sense it makes to align erato or silvaniform clade genomes to the H. melpomene reference, given their somewhat large divergence.)

Analysis

1. Align with mostly default options.

2. Convert to sorted BAM with readgroups using Picard AddOrReplaceReadGroups.

3. Run samtools flagstat to get basic mapping statistics.

4. Run GATK UnifiedGenotyper to call reference and SNP bases. Use -out_mode EMIT_ALL_CONFIDENT_SITES, -stand_call_conf=1 and -stand_emit_conf=1; the samples used are very low coverage so good SNPs are likely to have qualities lower than the default threshold of 30. 1 is too low, but it may reveal weaknesses in the aligners.

5. Hack coverage and SNPs at various different qualities using one line Perl scripts.

This analysis is OK as far as it goes, but there are several problems with it.

1. It should be possible to improve aligner performance by tuning parameters. I did try a handful of alternative running modes (Bowtie 2’s end-to-end and local, for example) but have not searched numerical parameter spaces at all (see comment below on BWA).

2. Reporting what GATK makes of an alignment does not necessarily reflect the quality of the alignment (see Mihaela Pertea and Steven Salzberg’s takedown of a similar comparison).

3. Mapping quality scores and things like ‘properly paired’ counts are assigned by the aligners themselves and may not reflect actual quality; it would be possible to write an aligner that dumped reads randomly but assigned maximum mapping qualities, which would come out quite well on many of the following metrics (although it might be expected that GATK would call a large quantity of spurious SNPs, of which more later).

4. Performance times are unlikely to be accurate. All jobs were run on a 48-core server (4x 12-core 2.2GHz AMD Opteron 6174s) with 514 Mb RAM. Alignment jobs were run in parallel on 36 cores (for all jobs except bwa sampe, which can only be run on a single core). The server is quiet, but it is quite possible that other jobs or IO bottlenecks slowed down the alignment jobs, which were only run once. Or input data may have been loaded into RAM for later alignments. Still, I expect the orders of magnitude to be correct.

5. There are many other metrics one can think of to assess aligners, many of which have been used in other aligner comparisons, but I don’t like simulated data, and we do not have a well-validated set of real reads with known positions to test against at present.

Aligners

With over 70 aligners available, which to choose? I have, of course, been far from comprehensive. I really just wanted to try Bowtie 2, but I did attempt a mildly systematic search for other potentially useful aligners. I used the HTS Mappers page and limited myself to DNA aligners that process Illumina data, are quality aware, handle paired ends and run multithreaded. Note that I probably ignored several good mappers by doing this; for example, I would have excluded Stampy because it is listed as not being multithreaded, as this feature has only recently been added. I then excluded several aligners because previously published data showed they were likely to be inferior to other aligners (for example, the SeqAlto paper shows SNAP to be very fast but perform poorly against most other aligners, particularly on divergent data, so it was unlikely to be useful for Heliconius data). I also ignored a few because they appeared to be abandoned or were commercial (CLC, Novoalign – Novoalign tests in July 2011 showed it to perform worse than Stampy for our purposes, as also shown in the Stampy paper). If I have missed your favourite aligner and you think it is worth testing on divergent data, please do let me know.

BWA is used throughout as a kind of control, as it is what we used to use and what is still widely used by many others (and recommended by GATK for example). It does not perform very well in any test, but perhaps it could be tuned to do much better, as it uses a seed of 32 bases which is very large for divergent data. But as it is often tempting to simply run using defaults, it is worth seeing just how bad the performance can be with divergent genomes.

Results

H. m. melpomene x H. m. rosina read mappings

Most of the aligners have been tested on the H. m. melpomene x H. m. rosina cross individual, as that’s the dataset I’m most interested in now. I took only the best performing aligners on to the H. ethilla and H. erato comparisons, assuming that aligners that perform badly on nearby genomes will do worse on divergent genomes.

[table id=1 /]

Timing notes: the BWA timings show times for the aln step (one for each read) plus the sampe step. Early versions of Stampy could use BWA in a hybrid mode, rapidly mapping easy reads with BWA and then mapping the remaining unmapped reads with Stampy. Recent versions can now take any BAM file as input for this hybrid mapping, so any aligner can be used in this way. The X + Stampy rows show this hybrid mapping, with the time for aligner X given first and the subsequent Stampy run second.

Looking at numbers of reads mapped is not particularly informative, because exact mappers like BWA and Bowtie 2 with the default end-to-end mapping setting (where the whole read has to map) will inevitably do worse than the partial mappers like Stampy and Bowtie 2 with the local mapping setting (where only part of a read has to map). But it does indicate that Bowtie 2 with local and very-sensitive settings and SMALT are both competitive with Stampy, whereas SeqAlto is not, on this data at least. It also shows that small improvements can be achieved by running Stampy in hybrid mode with any other aligner, but the performance improvement is not that great.  Sadly, running Stampy in fast mode achieved nothing like the 50% speed up the manual claims, with slightly degraded results. Both Stampy runs lag way, way behind Bowtie 2, SMALT and SeqAlto; it will make life much easier if we can convince ourselves that one of these aligners is more accurate.

H. ethilla and H. erato read mappings

Only BWA (as a control), Bowtie 2, SMALT and Stampy were tested on the more divergent genomes. Bowtie 2 was run with the local and verysensitive options as these options did not degrade performance very much in the H. melpomene tests.

[table id=3 /]

[table id=4 /]

Clearly BWA fails horribly on H. erato and does badly on H. ethilla. Bowtie 2, SMALT and Stampy all appear to be competitive, and cannot be separated on these metrics. While SMALT maps more reads, Bowtie 2 has more properly paired. As expected, Stampy’s performance is poor; Bowtie 2 is fastest, but SMALT is acceptable.

How many bases mapped?

As mentioned above, counting number of reads mapped is not very informative when partial mapping is allowed. It is much better to count number of bases mapped. The following charts show the total numbers of bases mapped for H. melpomene, H. ethilla and H. erato, broken down by mapping quality as recorded by GATK in the MQ field. Each bin value is a lower bound, so ’30’ means all bases with mapping qualities from 30 to 39.

This chart shows that, with data from similar genomes, all aligners do reasonably well. Actually, quite a lot of bases are not covered across the genome (around 200 million out of 273 million), but that’s presumably because only 4 million reads have been aligned. SMALT actually performs slightly better than Stampy (an extra half a megabase aligned), but any of the aligners would do a good job. Using Stampy in hybrid mode does not beat the SMALT result, except when paired with SMALT, when an additional megabase can be aligned, at the cost of quadrupling the running time.

Also, intriguingly, the chart shows how varied mapping qualities can be. BWA and SMALT don’t call anything over 60, Bowtie 2 doesn’t call anything over 40, but Stampy hands out a large amount of qualities in the 90s. Clearly there is some degree of bullshitting here, if we interpret mapping qualities as true phred scores – a mapping quality of 99 would mean there is a 1 in 10 billion chance of it being wrong. However, I don’t necessarily see a problem with this as long as thresholds are chosen appropriately. For example, in the Heliconius genome paper, we used a threshold of 67 for mapping qualities. This would have excluded all base calls from BWA, but it worked well for Stampy.

This effect also means that it is possible comparisons such as Heng Li’s ROCs are not entirely appropriate; for example, perhaps these curves retain a lot of Stampy bases at Qs of 20 or 30, which Stampy considers to be poor quality bases that should be rejected, and if a higher threshold was used it would produce a much better curve. But then, perhaps it shouldn’t output such high mapping qualities.

For H. ethilla and H. erato, there is a reduction in number of bases aligned as diversity increases, as expected, but Stampy maps several million more bases than any other aligner. Both Bowtie 2 and Smalt do reasonably well, however, and, depending on the analysis, it may be preferable to trade coverage for performance.

How many SNPs mapped?

The variations in a genome are usually our useful data points, rather than the number of bases. The following charts show the number of SNPs called by GATK and their mapping qualities.

Interpreting this chart is not straightforward because we have no way of deciding whether the SNPs are real or not, so it is not clear what is better. The Stampy alignment produces the most SNPs, but if it places a lot of reads incorrectly with highly exaggerated SNP calls, then GATK will call a lot of high quality SNPs from that alignment.

It may help to compare the number of bases mapped to the number of SNPs called. We might expect that the number of SNPs called should shrink or grow proportionally to the number of bases called. However, although Stampy aligns half a million more bases than SMALT, it also produces half a million more SNPs than SMALT, which is considerably more than we might expect proportionally. This makes me suspect that SMALT might be a sensible, conservative alternative here, rather than simply going for the aligner with the most SNPs. But then, maybe BWA would be an even more sensible option, as it maps most of the bases but apparently creates many fewer SNPs, with lower mapping qualities overall.

The following slide is taken from a talk on SMALT by the authors, and indicates that Stampy may indeed be overexuberant when aligning:

This indicates that the extra few percent of mapping that Stampy does is actually mostly erroneous, at least for nearby genomes. So to reduce false positives, SMALT appears to be a better option. But this is a test using simulated data by the authors of the tool, so may not be true in the wild (and may be falling foul of the mapping quality exaggeration issue discussed above – perhaps we could easily exclude the erroneous mappings by just raising our mapping quality threshold).

For H. ethilla and H. erato, the number of SNPs called is more in line with the number of bases aligned, and Stampy has a considerable edge over the other aligners, as it did for number of bases mapped. (As we would hope, H. ethilla has more SNPs than H. melpomene. However, H. erato actually has less, presumably because the most divergent regions simply don’t align.)

Discussion

At present, the data above don’t give us enough information to make a choice once and for all. I think I would still prefer to use Stampy over other aligners if I was working with silvaniform clade or erato clade genomes and I had to align them to the H. melpomene genome. But for aligning melpomene clade genomes, I will probably use SMALT for now, and am using SMALT for the current whole genome alignments. It provides the best tradeoff between performance and coverage (and apparently accuracy, but it is too early to confirm that). It means I can get this set of alignments done in a week, rather than two to three weeks with Stampy.

There is obviously a lot more that could be done here, but I wonder if it is worth doing. Obviously if erroneous data is getting into our analyses, that’s a problem, especially if tools are labelling erroneous data with inflated quality scores. But I think I would prefer to figure out better biological filters to apply to the data after alignment than taking several weeks of work to assess aligners in more detail.

However, I would be keen to at least pool opinions on this – should we validate our aligners in more detail? Have we seen any evidence of systematic errors getting into final data sets because of aligner problems? Are there other simple tests we could or should do that would settle the issue?