alimanfoo

A vintage year for African malaria vector genomics

Mon, 06 Jan 2025 00:00:00 +0000

2024 was a great year for African Anopheles genomics. Many long labours all converged to produce a bumper crop of new research papers which have redefined our understanding of major malaria vector species and how they are evolving. Here is a selection of papers published last year, all of which made use of or generously contributed data to the Malaria Vector Genome Observatory.

Clara Rehmann et al. published Evaluating evidence for co-geography in the Anopheles–Plasmodium host–parasite system.

Luna Kamua et al. published The Anopheles coluzzii range extends into Kenya: detection, insecticide resistance profiles and population genetic structure in relation to conspecific populations in West and Central Africa.

Andrzej Pastusiak et al. published A metagenomic analysis of the phase 2 Anopheles gambiae 1000 genomes dataset reveals a wide diversity of cobionts associated with field collected mosquitoes.

Sanjay Nagi et al. published Parallel evolution in mosquito vectors - a duplicated esterase locus is associated with resistance to pirimiphos-methyl in Anopheles gambiae.

Jorge Amaya Romero et al. published Mitochondrial variation in Anopheles gambiae and Anopheles coluzzii: phylogeographic legacy and mitonuclear associations with metabolic resistance to pathogens and insecticides.

Rob McCann et al. posted the preprint Variation in spatial population structure in the Anopheles gambiae species complex.

Beniamino Caputo et al. published Population genomic evidence of a putative ‘far-west’ African cryptic taxon in the Anopheles gambiae complex.

Mahamadi Kientega et al. published Whole-genome sequencing of major malaria vectors reveals the evolution of new insecticide resistance variants in a longitudinal study in Burkina Faso.

Tristan Dennis et al. published Cryptic population structure and insecticide resistance in Anopheles gambiae from the southern Democratic Republic of Congo.

Naomi Dyer et al. published Mechanisms of transcriptional regulation in Anopheles gambiae revealed by allele-specific expression.

Joel Odero et al. published Discovery of knock-down resistance in the major African malaria vector Anopheles funestus.

Sophia Mwinyi et al. posted the preprint Genomic analysis reveals a new cryptic taxon of malaria vectors with a distinct insecticide resistance profile in the coast of East Africa.

Brian Polo et al. posted the preprint Genomic surveillance reveals geographical heterogeneity and differences in known and novel insecticide resistance mechanisms in Anopheles arabiensis across Kenya.

Eric Lucas et al. published Copy number variants underlie major selective sweeps in insecticide resistance genes in Anopheles arabiensis.

Rounding off the year in style, Marilou Bodde and Joachim Nwezeobi posted the preprint Genomic diversity of the African malaria vector Anopheles funestus.

There are a load more great papers in the works from these and other groups, so I’m looking forward to 2025, Happy New Year!

A coding assistant for genomic data analysis?

Sat, 16 Nov 2024 00:00:00 +0000

Large genomic datasets like the vector observatory can be difficult to analyse. To make analysing these data easier, we've created analytical software and an online training course, and these have definitely helped. But many potential users of these data are coming from a background without much experience of genomics or programming, and so coming up with the right code can still be a challenge.

AI coding assistants are getting very good. E.g., I recently saw a talk and paper by Josh Proctor on using an LLM to generate disease modelling code which performed scarily well. Would it be possible to create a specialised AI assistant to support generation of code to analyse genomic data on malaria mosquitoes?

Use case: investigating insecticide resistance¶

Consider a use case where an analyst would like to investigate evidence for insecticide resistance in mosquito populations within a country of interest. Assume the analyst has good knowledge of vector biology but patchy knowledge of genomics and the molecular biology of insecticide resistance. Assume the analyst has completed basic training in vector genomics data analysis but is not an experienced programmer and finds writing data analysis code challenging.

Typically an analyst might start with a relatively broad question like, "Is there resistance to pyrethroid insecticides in Uganda?" But there are multiple possible mechanisms of pyrethroid resistance that need to be considered, and there are also multiple species of malaria vector present in Uganda. So this question is likely to need some refinement. If I were assisting this user, I would ask a series of questions to narrow down the focus, and to help the user become aware of the relevant data that are available. I would then be able to identify a smaller set of analysis function calls and parameter values that are relevant.

Would it be possible to create a chatbot assistant to help with this refinement? E.g., here is a hypothetical conversation between a user and an assistant that would refine the scope of the initial question to reach a point where code generation would be easier:

User	Assistant
Is there resistance to pyrethroid insecticides in Uganda?
	Which mosquito species are you interested in? In Uganda there are data available for Anopheles gambiae, Anopheles arabiensis and Anopheles funestus.
Anopheles gambiae.
	Which molecular mechanism of pyrethroid resistance are you interested in? E.g., target-site resistance, Cyp-mediated metabolic resistance, or other resistance mechanisms?
Cyp-mediated metabolic resistance.

Code generation¶

If we can get to this point, we probably have just enough information to start suggesting some code. E.g., we know that copy number amplification of Cyp genes has been associated with pyrethroid resistance. If I were assisting the analyst, I might suggest to try analysing CNV frequencies at a selection of Cyp genes which have been previously linked to pyrethroid resistance. Here is some code which sets up an API to access data on the Anopheles gambiae complex, then uses two function calls to compute and then visualise gene CNV frequencies in Uganda, grouping mosquitoes by year and level one administrative units.

# Set up the API to access data from the Anopheles gambiae complex.
import malariagen_data
ag3 = malariagen_data.Ag3()

# Define genes of interest.
interesting_cyp_genes = [
    "AGAP002862",  # Cyp6aa1
    "AGAP013128",  # Cyp6aa2
    "AGAP002865",  # Cyp6p3
    "AGAP000818",  # Cyp9k1
    "AGAP008212",  # Cyp6m2
    "AGAP008218",  # Cyp6z2    
]

# Compute gene CNV frequencies.
df_cyp_cnv_frq = ag3.gene_cnv_frequencies(
    region=interesting_cyp_genes,
    sample_query="country == 'Uganda' and taxon == 'gambiae'",
    cohorts="admin1_year",
)

# Visualise CNV frequencies as a heatmap.
ag3.plot_frequencies_heatmap(df_cyp_cnv_frq)

This is a useful starting point for exploratory data analysis, because we can see for example that amplifications of Cyp6aa1 and Cyp69k1 are very high or fixed in almost all regions and years. Both of these genes are associated with pyrethroid resistance, and so this combination implies very strong metabolic resistance.

For Uganda there are data for multiple years and regions, and the user might also be interested in changes over time. I might therefore be tempted to suggest some code to visualise frequency time series, e.g.:

# Compute gene CNV frequencies.
ds_cyp_cnv_frq = ag3.gene_cnv_frequencies_advanced(
    region=interesting_cyp_genes,
    sample_query="country == 'Uganda' and taxon == 'gambiae'",
    area_by="admin1_iso",
    period_by="year",
)

# Plot frequency time series.
ag3.plot_frequencies_time_series(
    ds_cyp_cnv_frq,
    height=900,  # need enough height for three sub-plots
)

These time series give us some more information, e.g., they show that frequencies of Cyp6aa1 and Cyp9k1 have generally increased after 2012, although there are some geographical differences between the three regions of Uganda for which we have data.

Training data / examples¶

What training data or existing examples exist that could be used to train an assistant?

Probably the best source would be the learning materials for the data analysis training course we've developed. These materials comprise a collection of Jupyter notebooks with code examples weaved in with explanatory text.

We also create many Jupyter notebooks as part of previous and ongoing studies where we're collaborating with scientific colleagues to analyse real datasets. These tend not to have nearly as much explanatory text, rather they are focused on just analysing data to produce tables and figures for reports and papers. However, if we knew they were to be used to train an assistant, we could potentially go through and annotate some with more explanatory text.

In general, I wonder what the best way of using Jupyter notebooks to train an assistant would be.

For the API itself, there are docstrings and type annotations for all public methods. The docstrings generally don't include examples, but will have some brief text explaining what the methods are used for. These are also used to generate an API documentation website. What would be the best way to make use of this?

Finally, I wonder how an assistant could learn more about the content of the actual data. E.g., how to teach an assistant something about the samples that have been sequenced from different places, time points and mosquito species? This information is very useful when trying to craft parameter values for API calls. I guess this might be roughly analogous to problems where you want to train an assistant to help generate SQL queries for a specific database.

Learning methods¶

What would be the best way to train a generative model?

If we have a library of Jupyter notebook with text and code examples, could this be used for some form of retrieval-augmented generation? I.e., the user's question is used to retrieve the most relevant snippets of notebooks, then these are included in the prompt?

What about few-shot learning or fine-tuning? For one specific use case like the one described above, perhapse few-shot learning would be sufficient? But to cover a breadth of different use cases, perhaps fine-tuning is needed?

Integration with Google colab or Jupyter lab?¶

Colab has code generation integrated within the notebook UI, which is exactly where you'd want this assistant to be:

...is there a way to achieve this, ideally for colab but alternatively for any Jupyter lab interface?

Challenges and subtleties¶

If it was possible to generate code examples like the ones above, then this would definitely be getting us towards the original question, but it's worth mentioning that there are a lot of variations and other analyses we might suggest, as well as some subtleties.

E.g., if we wanted to analyse Anopheles funestus instead of Anopheles gambiae, we would need to provide a different set of gene identifiers for the interesting Cyp genes to include.

Also, there are actually 107 cytochrome P450 (Cyp) genes in the Anopheles gambiae genome in total, and although only some of these have been associated with resistance previously, it's possible that prior knowledge is incomplete. Should we analyse all 107 genes, or just the smaller set of validated genes? Furthermore, some studies are also starting to find SNPs in Cyp genes are markers of resistance. Should we analyse SNPs as well as CNVs? These are tricky decisions, because we could come up some analysis code that is perfectly sensible from a biological point of view, but would overwhelm the user with data. E.g., there are often dozens of SNPs within each Cyp gene, even after filtering out low frequency variants, and 107 genes each with dozens of SNPs is hard to visualise.

These are just a couple of potential challenges involved in considering what would be the most appropriate code examples to generate to address the original question, there are others.

How does Gemini do now?¶

For fun, how does vanilla Gemini do now without any additional context or specialisation?

Starting with something simple:

Good job. Try something a little more specific:

Nailed that one too. OK, now for a challenge:

The first couple of lines are not far off.

In the sample_metadata() call the parameter where should be sample_query, although to be fair the parameter value is correct.

In the line filtering down to Anopheles gambiae samples, the "species" field should be "taxon", but otherwise that code would work.

From here on it's making stuff up. E.g., there is no function gene_variants(). Although Gemini does seem to know that there is a gene called "Vgsc" which you would want to analyse, which is correct as Vgsc is the target site of pyrethroid insecticides.

How about a challenging but more specific prompt:

Lots of mistakes here too. E.g., there's no function gene_metadata(). There are functions sample_metadata() and gene_cnv() but these haven't been called correctly. The attempt to calculate CNV frequencies is entertaining!

Clearly a lot of room for improvement, but also a lot of base knowledge that could be built on.

Last word¶

Hopefully this is a useful use case to start thinking about what sort of assistant would be useful and if/how it could be implemented. Suggestions and thoughts very welcome.

Faster neighbour-joining

Fri, 18 Oct 2024 00:00:00 +0000

A new Python package called anjl provides fast implementations of the neighbour-joining algorithm and interactive plotting of unrooted trees.

Neighbour-joining trees (NJT) are a useful tool for exploratory analysis of genetic population structure. Previously I've used a Python package called biotite to compute NJT, but recently I was trying to build some bigger trees and things started to take a long time to run. I decided to see if I could create a faster implementation, and so anjl was born...

The canonical algorithm¶

Neighbour-joining is an iterative algorithm where each iteration involves searching a distance matrix for the pair of nearest neighbours which minimise a value $Q$, which are then joined to create a new internal node in the tree. The canonical algorithm performs a full search of the distance matrix in each iteration, and so expected performance is $O(n^3)$, meaning that speed of computation scales with the number of samples cubed. Thus even a very efficient implementation of the canonical algorithm will hit scalability limitations pretty quickly. Nevertheless, I thought I'd start with the canonical algorithm to have a reference implementation to compare to. After a bit of prototyping and tuning I had an implementation using numba that was around 10x faster than biotite. I decided this was worth packaging and so created a new Python packaged called anjl which stands for "A neighbour-joining library".

The rapid algorithm¶

I also found a paper by Martin Simonsen et al. from 2008 describing a "rapid" neighbour-joining algorithm. This algorithm maintains a copy of the distance matrix where each row is sorted, and uses this to avoid doing a full search of the distance matrix in each iteration. I added a numba implementation of neighbour-joining based on Simonsen's algorithm to anjl. Benchmarking with some mosquito genomic data showed performance was more than 2x better the canonical algorithm for larger distance matrices. I was pretty happy with this, until...

The dynamic algorithm¶

I probably should've searched harder before I started implementing the rapid algorithm, because I then found a more recent paper by Philip Clausen describing a dynamic algorithm with some impressive benchmark results. The paper is pretty terse and it took a while to wrap my head around it, but a bit like staring at a magic eye picture, eventually it started to make sense. The basic insight is that the neighbour-joining criterion $Q$ is weakened at each iteration. This means that, if you start by finding the minimum value of $Q$ within each row of the distance matrix, then in the next iteration all values of $Q$ within the same row will be greater. In other words, the minimum value of $Q$ within each row provides a lower bound for all values of $Q$ in the next (and subsequent) iterations. You can use this information to skip over searching a lot of rows of the distance matrix at each iteration. The only exceptions to this rule are new nodes created during the neighbour-joining process, so you have to check if they create a better join. It's an amazing insight which is also beautifully simple to implement (when you get it!). And it's ridiculously fast, around 10x better than the rapid algorithm on the data I was using to benchmark with.

Interactive plotting¶

I also decided to include in anjl a function to visualise a NJT as an unrooted tree using plotly. Here plotly is useful because you can interact with the plot, hovering over leaf nodes to get more information, zooming and panning, and showing or hiding data via the legend. Here's a simple example.

import anjl
D, leaf_data = anjl.data.mosquitoes()
Z = anjl.dynamic_nj(D)
anjl.plot(
    Z, 
    leaf_data=leaf_data, 
    color="taxon",
    hover_name="sample_id",
)

Some numba lessons¶

I've used numba before but learned a couple of new things from working on this project.

Boundscheck¶

By default, numba does not perform any checks to see if you are trying to index an array with a value that is too large (out of bounds). This is for performance reasons, otherwise the bounds checks would take some time. While I was working on the rapid algorithm I inadvertantly introduced an out of bounds indexing bug. I would expect this to cause program crashes, but actually the bug only caused a very occasional crash. It also didn't actually affect the correctness of the output. Long story short, it took me a while to figure out something was wrong. Eventually I realised I might have an indexing bug, and I set the boundscheck parameter to True in the number njit() decorator. Immediately this revealed that, yes, I had an index out of bounds error.

The moral of the story is, it's probably a good idea to turn boundscheck on during development and when running tests. I ended up putting a conditional block in my code like this...

# Detect whether we are running via pytest, if so run with boundscheck enabled to catch
# any out of bounds errors.
# https://docs.pytest.org/en/stable/example/simple.html#detect-if-running-from-within-a-pytest-run
if os.environ.get("PYTEST_VERSION") is not None:
    BOUNDSCHECK = True
else:
    BOUNDSCHECK = False

I then used this in all njit decorators, e.g.:

@njit(
    boundscheck=BOUNDSCHECK,
)
def dynamic_search(...):
    ...

Wraparound¶

When you are indexing an array, numba also checks to see if you are indexing with a negative value, in case it needs to handle wraparound logic (e.g., where a value like -1 means the last value in the array). But this wraparound logic also takes some time, and if you're doing a lot of array access, this can affect performance. You can circumvent this by using unsigned integer data types for variables that you are using to index arrays. However, if you are using for loops in numba, the induction variable is a signed integer by default (although this is different for parallel for loops). You can get around this by casting the induction variable to an unsigned integer data type like np.uintp. E.g., all my for loops ended up looking like this...

for _i in range(n):
    i = np.uintp(_i)
    # Use i for indexing arrays...

Profiling¶

What do you do if you think you've done everything you can to tune a numba function, but it's still running slower than you think it should? Ideally you'd do some profiling to find which section of code the program is spending most time in. While there are great tools for profiling normal Python code, there isn't yet a good tool for profiling code within a numba function. The only option I could find is a tool called profila which I did get to work some of the time, but which also crashed a lot. There is an open feature request on the numba github repo, it seems a better solution is possible but still some way off.

Algorithm is everything¶

Of course the ultimately lesson here is that finding the right algorithm is potentially much more important than tuning tricks! The dynamic algorithm makes a massive difference, which is entirely due to an insight into the structure of the problem.

Anopheles coluzzii in Kenya - neighbour-joining tree analysis

Thu, 05 Sep 2024 00:00:00 +0000

This post replicates the neighbour-joining tree analysis from the paper by Kamau et al. (2024) on Anopheles coluzzii in northwestern Kenya, and demonstrates the plot_njt() function in the malariagen_data Python package.

We recently worked with Luna Kamau's group at KEMRI to analyse genomic data from a survey of Anopheles mosquitoes across Kenya, which turned up the unexpected discovery of Anopheles coluzzii in northwestern Kenya. The finding was published earlier this year in Kamau et al. 2024, and Kelly Bennett from my team who led the analysis presented the results at a recent MalariaGEN journal club.

In this post I thought I'd replicate one of the analyses and figures from the paper, which is a neighbour-joining tree investigating how the Kenyan An. coluzzii were genetically related to An. coluzzii from countries in West and Central Africa. It's also a good opportunity to demonstrate the plot_njt() function in the malariagen_data Python API, which was recently added and used for this analysis.

Setup¶

To begin with, set up access to data from the Malaria Vector Genome Observatory.

import malariagen_data
ag3 = malariagen_data.Ag3(
    results_cache="results_cache",
)

The data are stored in Google Cloud Storage (GCS) in the us-central1 region. You can run this analysis on Google Colab, although you might need a bit of patience as Colab VMs can be in a different US region so may not be super close to the data, and also don't have a huge amount of compute power. I ran this analysis on a Vertex AI Workbench VM with 32 vCPU and 32 GB of RAM (e2-highcpu-32) and it took a couple of minutes. If you'd like to run this analysis yourself, see Vector Observatory Data Access for more info.

Sample selection¶

The samples we're interested in are in sample set 1274-VO-KE-KAMAU-VMF00246 which comes from Luna Kamau's study 1274-VO-KE-KAMAU. This sample set is included in the Ag3.9 release. I'm going to also include samples from the Anopheles gambiae 1000 Genomes Project (included in the Ag3.0 release) and from Fontaine et al. (2015) (included in Ag3.10) for comparison.

sample_sets = [
    # Samples from Kenya from Luna Kamau at KEMRI.
    "1274-VO-KE-KAMAU-VMF00246",
    # Ag1000G sample sets for comparison.
    "AG1000G-BF-A",
    "AG1000G-ML-A",
    "AG1000G-GH",
    "AG1000G-CI",
    "AG1000G-CF",
    "AG1000G-CM-B",
    "AG1000G-CM-C",
    "AG1000G-AO",
    # Samples from Fontaine et al. (2015).
    "fontaine-2015-rebuild",
]

For this analysis we are only interested in Anopheles coluzzii, so define a sample metadata query to select this taxon.

sample_query = "taxon == 'coluzzii'"

Take a look at how many samples will be included in the analysis.

df_count = ag3.count_samples(
    sample_sets=sample_sets,
    sample_query=sample_query,
    index=["sample_set", "country", "admin1_name"],
)
df_count

		taxon	coluzzii
sample_set	country	admin1_name
1274-VO-KE-KAMAU-VMF00246	Kenya	Turkana	26
AG1000G-AO	Angola	Luanda	81
AG1000G-BF-A	Burkina Faso	Hauts-Bassins	82
AG1000G-CF	Central African Republic	Bangui	18
AG1000G-CI	Cote d'Ivoire	Lagunes	80
AG1000G-CM-B	Cameroon	Far North	2
AG1000G-CM-B	Cameroon	North	5
AG1000G-CM-C	Cameroon	Adamaoua	2
		Centre	12
		Littoral	2
		North	1
		South	2
AG1000G-GH	Ghana	Central Region	25
		Eastern Region	1
		Greater Accra Region	14
		Western Region	24
AG1000G-ML-A	Mali	Koulikouro	27
fontaine-2015-rebuild	Burkina Faso	Hauts-Bassins	3
	Cameroon	Centre	4
	Cameroon	South	4

How many samples in total?

df_count["coluzzii"].sum()

np.int64(415)

Neighbour-joining tree analysis¶

For this analysis, I'm going to use data from Chromosome 3 (contig="3RL"). I'll request 100,000 SNPs, requiring SNPs to have no missing genotype calls (max_missing_an=0) and a minimum minor allele count of 10 (min_minor_ac=10) because rare variants are less informative. All SNPs meeting these criteria will be found and then thinned to approximately the requested number.

# Define the order that countries should appear in the legend.
category_orders = {
    "country": [
        "Kenya",
        "Burkina Faso",
        "Mali",
        "Cote d'Ivoire",
        "Ghana",
        "Central African Republic",
        "Cameroon", 
        "Angola",
    ],
}

# Compute and plot the neighbour-joining tree.
fig = ag3.plot_njt(
    metric="cityblock",
    title="Anopheles coluzzii",
    width=800,
    height=800,
    show=False,
    color="country",
    marker_size=5,
    line_width=1,
    category_orders=category_orders,
    region="3RL",
    n_snps=100_000,
    sample_sets=sample_sets,
    sample_query=sample_query,
    site_mask="gamb_colu",
    max_missing_an=0,
    min_minor_ac=10,
)

# Tidy up the plot.
fig.update_layout(
    legend=dict(
        title="Country",
    ),
    title_font=dict(
        size=20,
    ),
)

# Show the plot.
fig

This plot is interactive, you can hover over the leaf nodes (circle markers) which represent the individual mosquitoes and get more information about them such as location and date of sampling. The legend is also interactive, you can click on legend entries to show or hide different countries.

The Kenyan mosquitoes are plotted in blue, and are most closely connected to mosquitoes from Burkina Faso, Mali and northern Cameroon.

Topography and climate of northwestern Kenya¶

The Kenyan An. coluzzii came from Turkana which is in the northwestern part of the country, bordering South Sudan. More work is needed to understand how much connectivity and migration of An. coluzzii mosquitoes between Kenya and other countries, but I found it interesting to look at some maps of topography and climate and see how Turkana relates to surrounding areas.

Topographically, northwestern Kenya is more connected with South Sudan than with other regions of Kenya. Here's a map of topography:

Climate is also interesting, because northwestern Kenya is very arid for at least part of the year, similar to Mali, Burkina Faso, northern Ghana and northern Cameroon. E.g., here's a map of IPCC climate data from January which is the middle of the dry season north of the equator:

An. coluzzii is known for its ability to survive in regions with a longer dry season.

Neighbour-joining tree implementation notes¶

The implementation of the plot_njt() function in the malariagen_data package uses biotite under the hood to compute the neighbour-joining tree from a distance matrix. There are also implementation of the neighbour-joining algorithm in BioPython and scikit-bio, but I did some benchmarking and biotite was the fastest.

To create an interactive plot of the resulting tree, I implemented the equal angles layout, which is my favourite layout for unrooted trees, and plotting the results using Plotly. It would be amazing if there was a Python library that made interactive plots of different types of trees and layouts, like ggtree but using Plotly or Bokeh, but I couldn't find anything yet. Interactivity is so useful for exploratory work, being able to hover over samples and zoom in and out really helps when you're trying to make sense of complicated datasets.

Blogging with Jekyll, Jupyter notebooks and interactive plots (Plotly, Bokeh)

Tue, 20 Aug 2024 00:00:00 +0000

This post explores how to create blog posts from Jupyter notebooks which include interactive plots made with Plotly or Bokeh.

This blog is hosted on GitHub Pages and uses the Jekyll static site generator to create the site's HTML from the content I create. If I'm writing simple blog posts then authoring a markdown file is fine, but more often I'd like to create posts with code examples and outputs including plots of data, and so would like to author the post as a Jupyter notebook. I'd particularly like to be able to include interactive plots using Plotly and/or Bokeh.

I've previously figured out how to do this with Bokeh plots by using nbconvert to convert a notebook to a markdown file. However, this involved a bit of manual post-processing faff. So I thought I'd revisit this and see if I can find a simpler approach which works for recent versions of Bokeh and Plotly. It seems there's several different ways you can do this, but here is one possible approach.

Prerequisites¶

I have the GitHub repo containing my blog content cloned to my local filesystem. I also have things set up so I can generate and preview the site locally by doing:

$ bundle exec jekyll serve

I also have jupyterlab, nbconvert, plotly and bokeh installed into a conda environment on my local system. I'll explore later exactly how I'm doing this to ensure I have a locked environment associated with each blog post.

Creating the source notebook for a new blog post¶

To author a new post I've run Jupyter Lab and created a notebook. Ultimately the content for the post will need to end up in the _posts folder, but I've put the source notebook in a separate folder. E.g., this notebook is at src/2024-08-20-jupyter-blogging/post.ipynb.

Front matter¶

Every blog post needs to have "front matter" which provides the title and tells Jekyll to display the content as a blog post. To add this front matter, I've put the following into the very first cell of this notebook, which I've selected to be a "Raw" cell type:

---
layout: post
title: Blogging with Jekyll, Jupyter notebooks and interactive plots (Plotly, Bokeh)
---

I also add an empty line after the front matter, which ensures the front matter won't get mangled together with whatever cell comes after.

Disabling liquid tags¶

A problem that can occur when blog posts contain code or interactive plots is that Jekyll mistakenly tries to process a piece of code as a liquid tag, which is a system that Jekyll uses to inject variables into posts. I don't need to use liquid tags at all, so I'll disable them for the entire post. This can be done by adding a "Raw" type cell at the beginning of the notebook containing:

{% raw %}

Then add another "Raw" type cell right at the end of the notebook containing:

{% endraw %}

Plotly interactive plots¶

Check what version of Plotly is installed.

import plotly
plotly.__version__

'5.24.0'

I've found that Plotly depends on RequireJS, and although this will automatically be loaded when using Jupyter Lab, this won't be automatically included in the HTML generated by nbconvert. So this cell manually ensures RequireJS is present.

%%html
<script
src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.7/require.min.js"
integrity="sha512-J5ha2LF4Le+PBQnI5+xAVJDR+sZG9uSgroy4n/A6TLjNkvYQbqZA8WHZdaOvJ0HiKkBC9Frmvs10rFDSHKmveQ=="
crossorigin="anonymous"
referrerpolicy="no-referrer"></script>

If I don't want to see this cell in the blog post, I could add a cell tag to remove the input.

I'd like to create Plotly plots and have them display interactively, both within the Jupyter notebook as I'm authoring it, and in the final blog post displayed on the Jekyll site. To get this to work, I found I needed to add the following cell to the notebook before creating the plot.

import plotly.io as pio
pio.renderers.default = "notebook_connected+plotly_mimetype"

Apparently this tells Plotly to render any plots in a way that is suitable for multiple environments. The plotly_mimetype renderer is suitable for viewing the source notebook within Jupyter Lab. The notebook_connected renderer is necessary for the plot to appear in the Jekyll site.

Now try to create a plot.

import plotly.express as px

df = px.data.iris()
df

	sepal_length	sepal_width	petal_length	petal_width	species	species_id
0	5.1	3.5	1.4	0.2	setosa	1
1	4.9	3.0	1.4	0.2	setosa	1
2	4.7	3.2	1.3	0.2	setosa	1
3	4.6	3.1	1.5	0.2	setosa	1
4	5.0	3.6	1.4	0.2	setosa	1
...	...	...	...	...	...	...
145	6.7	3.0	5.2	2.3	virginica	3
146	6.3	2.5	5.0	1.9	virginica	3
147	6.5	3.0	5.2	2.0	virginica	3
148	6.2	3.4	5.4	2.3	virginica	3
149	5.9	3.0	5.1	1.8	virginica	3

150 rows × 6 columns

px.scatter(
    data_frame=df, 
    x="sepal_length", 
    y="sepal_width", 
    color="species",
    width=600,
    height=500,
    title="Plotly demo",
)

A scatter plot should be visible in the cell output above, both in the source notebook and the blog post. It should also be interactive, e.g., pan, zoom and hover tooltips should work.

Bokeh interactive plots¶

Check what version of Bokeh is installed.

import bokeh
import bokeh.plotting
import bokeh.palettes
import bokeh.transform
bokeh.__version__

'3.5.2'

I've found that Bokeh depends on jQuery, and although this will automatically be loaded when using Jupyter Lab, this won't be automatically included in the HTML generated by nbconvert. So this cell manually ensures jQuery is present.

%%html
<script
  src="https://code.jquery.com/jquery-3.7.1.slim.min.js"
  integrity="sha256-kmHvs0B+OpCW5GVHUNjv9rOmY0IvSIRcf7zGUDTDQM8="
  crossorigin="anonymous"></script>

If I don't want to see this cell in the blog post, I could add a cell tag to remove the input.

Now call output_notebook() to get plots showing as output within the notebook.

bokeh.plotting.output_notebook()

Loading BokehJS ...

Above should say "BokehJS ... successfully loaded." both in Jupyter Lab and the final post.

Now create a Bokeh plot.

# Define variables to plot.
x_variable = "sepal_length"
y_variable = "sepal_width"
color_variable = "species"

# Set up mapping from values to colors.
color_factors = df[color_variable].unique()
palette = bokeh.palettes.d3['Category10'][len(color_factors)]
color_map = bokeh.transform.factor_cmap(
    field_name=color_variable,
    palette=palette,
    factors=color_factors,
)

# Create a figure.
fig = bokeh.plotting.figure(
    width=500, 
    height=400,
    tools="pan,box_zoom,wheel_zoom,save,reset,hover,help",
    tooltips=[
        ("species", "@species"),
        (x_variable, f"@{x_variable}"),
        (y_variable, f"@{y_variable}"),
    ]
)

# Add a scatter trace.
fig.scatter(
    source=df, 
    x=x_variable, 
    y=y_variable, 
    color=color_map,
    legend_field=color_variable,
    size=10,
    alpha=0.7,
)

# Tidy up.
fig.title = "Bokeh demo"
fig.xaxis.axis_label = x_variable
fig.yaxis.axis_label = y_variable

bokeh.plotting.show(fig)

A scatter plot should be visible above, and should be interactive, e.g., pan, zoom and hover tooltips should work.

Converting the notebook to HTML¶

Once finished editing the notebook, I now need to convert it to something that Jekyll can use. nbconvert can convert a notebook to either markdown or HTML, and previously I have converted to markdown. However, currently I'm finding it slightly easier to convert to HTML. Jekyll is quite happy to be given an HTML file.

This from the command line:

$ jupyter nbconvert post.ipynb --to html --template basic --no-prompt

...will generate a file post.html alongside post.ipynb.

Previewing the new post¶

Above I mentioned that I put the source notebook at src/2024-08-20-jupyter-blogging/post.ipynb. This means the converted HTML is now at src/2024-08-20-jupyter-blogging/post.html. Because it's outside the _posts folder, it won't yet be visible on the front page of the blog. To get the new post listed on the home page I did:

$ cd _posts
$ ln -s ../src/2024-08-20-jupyter-blogging/post.html 2024-08-20-jupyter-blogging.html

This creates a symlink within the _posts folder. The post is now visible at http://127.0.0.1:4000/2024/08/20/jupyter-blogging.html.

Customising CSS¶

I added a bit of CSS to style the blog post. Here's some rules I added to _sass/_base.css:

div.input_area pre {
    padding: 8px 12px;
    background-color: #eef;
    border: 1px solid #e8e8e8;
    border-radius: 3px;
}

div.input, div.output {
    font-size: 12px;
}

hr {
    margin-block-start: 1em;
    margin-block-end: 1em;
}

Managing environments¶

I mentioned above that I'm creating a separate Python environment for each blog post. The rationale for this is that software changes over time, and if I ever need to go back and edit or update a blog post created from a Jupyter notebook, I'd like a record of which packages and versions I used. When I move on to the next blog post, I am free to use a different set of packages and versions, without disrupting any previous blog posts.

I'm doing this using pixi. E.g., when I first created the directory for the source notebook for this post, I also initialised a pixi environment...

$ mkdir -pv src/2024-08-20-jupyter-blogging
$ cd src/2024-08-20-jupyter-blogging
$ pixi init .
$ pixi add python=3.11 jupyterlab nbconvert plotly pandas bokeh

To launch Jupyter Lab I then did:

$ pixi run jupyter lab

...and to convert the notebook to HTML:

$ pixi run jupyter nbconvert post.ipynb --to html --template basic --no-prompt

To get this conversion running automatically on a change to the source notebook:

$ while inotifywait -e close_write post.ipynb; do pixi run jupyter nbconvert post.ipynb --to html --template basic --no-prompt; done

Here's a bash script I use to run the automatic conversion, adding in options to remove cell input or output or whole cells based on cell tags:

#!/bin/bash
while inotifywait -e close_write post.ipynb
do
    pixi run jupyter nbconvert post.ipynb \
	 --to html \
	 --template basic \
	 --no-prompt \
	 --TagRemovePreprocessor.enabled=True \
	 --TagRemovePreprocessor.remove_input_tags remove_input \
	 --TagRemovePreprocessor.remove_all_outputs_tags remove_output \
	 --TagRemovePreprocessor.remove_cell_tags remove_cell
done

Using pixi this way, everything I need for each post is then self-contained and versions are locked.

Summary of dev setup¶

When developing a new blog post...

Open a terminal, cd to the root directory of my blog repo, and run Jekyll:

$ bundle exec jekyll serve

Open a second terminal, create a source directory for the new blog post, initialise pixi, add packages, run Jupyter Lab:

$ mkdir -pv src/YYYY-MM-DD-new-blog-post
$ cd src/YYYY-MM-DD-new-blog-post
$ pixi init .
$ pixi add python=3.11 jupyterlab nbconvert plotly pandas bokeh
$ pixi run jupyter lab

...and create a notebook called "post.ipynb".

Open a third terminal, start a watcher to convert the notebook to HTML on save:

$ cd src/YYYY-MM-DD-new-blog-post
$ while inotifywait -e close_write post.ipynb; do pixi run jupyter nbconvert post.ipynb --to html --template basic --no-prompt --TagRemovePreprocessor.enabled=True --TagRemovePreprocessor.remove_input_tags remove_input --TagRemovePreprocessor.remove_all_output_tags remove_output --TagRemovePreprocessor.remove_cell_tags remove_cell; done

Symlink the converted HTML into the _posts folder:

$ cd _posts
$ ln -s ../src/YYYY-MM-DD-new-blog-post/post.html YYYY-MM-DD-new-blog-post.html

Troubleshooting¶

On one occasion I found the post was not rendering at all, due to the fact that the front matter was not followed by a newline. I had placed another raw cell immediately after the first raw cell with the front matter, and the two raw cells were output in the HTML without any newline between them. Adding in a newline after the front matter solved it.

I hope that might be useful. Any suggestions or pointers to alternative approaches welcome in comments below :)

Diplotype clustering

Fri, 16 Aug 2024 00:00:00 +0000

This post introduces advanced diplotype clustering, a new function in the malariagen_data Python package for investigating insecticide resistance mutations. Cross-posted from Sanjay Nagi's blog, Sanjay did most of the work on this but it was a pleasure to collaborate.

Vectors of disease evolve rapidly in response to the interventions we throw at them. By monitoring genetic changes in these populations over time, we can detect emerging insecticide resistance mechanisms and monitor the spread of known mechanisms, with the aim of informing vector control strategies. In a world of limited active ingredients, this surveillance will be crucial for maintaining the effectiveness of front-line interventions like long-lasting insecticide-treated bed nets (LLINs) and indoor residual spraying (IRS).

To facilitate this, we have been developing systems and tools for genomic surveillance. The MalariaGEN Vector Observatory have developed the Python package malariagen_data, which provides tools for accessing and analysing genomic datasets from major malaria vectors. By developing innovative software, we can help to build capacity in genomic research, allowing more people to perform robust, complex genomic analyses that would otherwise be limited to a select few.

In this blog post, we wanted to share a new function that we have recently added to malariagen_data: plot_diplotype_clustering_advanced(). This method allows us to rapidly zoom in on a genome region of interest and identify selective sweeps, assess their size, detect potential gene flow events between countries or species, and investigate whether sweeps are driven by copy number variants (CNVs), amino acid mutations or both.

But what exactly are diplotypes, and why are they useful? A diplotype, sometimes referred to as a multi-locus genotype, is essentially the combination of two haplotypes from a single mosquito - one from each chromosome - at a particular genomic region. By analysing diplotypes rather than haplotypes, we can better capture the full genetic variation present in an individual, including complex structural variants like CNVs that can be difficult to phase onto haplotypes. Often, CNVs and multiallelic SNPs are ignored when analysing haplotype data. The more mosquitoes we sequence, the worse this problem gets - An. gambiae s.l. is so genetically diverse, eventually, a significant proportion of all SNPs become multiallelic.

The new diplotype clustering functionality in malariagen_data performs hierarchical clustering on diplotypes from a specified genomic region. It then visualises the results, displaying:

The clustering dendrogram
Sample metadata (e.g., species, collection location)
Heterozygosity of each sample (within this genomic region)
Copy number at genes of interest
Amino acid variants in a specified transcript

A case study¶

To illustrate the power of this approach, let's look at a case study of the Gste2 gene from some recent whole-genome data of An. gambiae s.l from Obuasi, central Ghana. Anopheles mosquitoes from this area are highly resistant to multiple classes of insecticides (Mugenzi et al. 2022). The Gste2 gene is known to be involved in resistance to DDT (and potentially other insecticides), through either copy number variation (Lucas et al. 2019), amino acid mutations (Mitchell et al. 2014), or both. Gste2-I114T and Gste2-L119V are the major amino acid mutations at this locus known to confer resistance.

Here is an example of a plot generated by the plot_diplotype_clustering_advanced() function:

We've annotated this figure with some diplotype clusters which are particularly interesting. For example, cluster B contains a large number of diplotypes which are all genetically identical and have very low heterozygosity - this is what you expect when a selective sweep has occurred, and you find many individuals that are homozygous for the haplotype under selection. All individuals in cluster B also carry the Gste2-I114T substitution. If we did not already know something about this mutation, this figure would give us a clue that the mutation is a potential driver of selection and insecticide resistance. In fact, we already know this mutation causes insecticide resistance, so it is no surprise to find it linked to a selective sweep in this dataset.

We can also see another large cluster (cluster C) which does not harbour either I114T or L119V, but instead, the Gste2-F120L mutation. This cluster is homozygous for F120L and shows low heterozygosity, again indicative of diplotypes which have two copies of the same swept haplotype. Two things about the Gste2-F120L mutation are convincing as a potential driver of resistance. Firstly, it is in very close physical proximity to known resistance mutations in codons 114 and 119. According to Riveron et al. (2014), the 120 codon is located at the active site of the enzyme and is therefore likely to interact with the insecticide. Secondly, there are no CNVs associated with this sweep, and no other amino acid variants except N3K, which is less likely to be causative due to its physical location away from the active site.

We also observe a small cluster of individuals (cluster A) which harbour a copy number variant (CNV) spanning Gste2, Gste1, Gste3 and Gste7. CNVs could be driving insecticide resistance by increasing the expression of the genes they encompass, allowing the mosquito to detoxify more of the insecticide as a result. This CNV, an amplification, seems to exist at a variable copy number. In total, we can see six or seven distinct CNVs in these samples.

This case study demonstrates how diplotype clustering can provide insights into the mutations causing insecticide resistance; in a single snapshot, we can explore amino acid and CNV data and really understand the nature of selection at a genomic region.

Here is the code which was used to generate the figure above. The function generates a Plotly figure, which is interactive and can be panned and zoomed as well as hovering to see some additional data.

# setup
import malariagen_data
ag3 = malariagen_data.Ag3(
    show_progress=False,  # hide progress messages for the blog post
)

# example function call
ag3.plot_diplotype_clustering_advanced(
    region="3R:28,597,000-28,600,000",         # The genomic region for clustering
    cnv_region="3R:28,594,000-28,605,000",     # The genomic region for CNV data
    snp_transcript="AGAP009194-RA",            # The transcript for amino acid variants
    sample_sets="1244-VO-GH-YAWSON-VMF00149",  # The sample set
    sample_query=None,                         # A query to filter samples
    site_mask="gamb_colu",                     # The site mask to use
    linkage_method="complete",                 # The linkage method to use
    color="taxon",                             # The metadata column to determine color
)

We hope that others will find it useful for their own research. Please feel free to get in touch if you have any questions or feedback :)

Arthropod Genomics Symposium 2024

Wed, 03 Jul 2024 00:00:00 +0000

I felt very lucky to attend the Arthropod Genomics Symposium in Nairobi last week. The scientific talks were excellent, and there aren’t many other meetings where you get the same crossover of folks working on disease vectors, agricultural pests and species in need of conservation - three different fields of applied genomics but lots of common challenges and new ideas. The last AGS I attended was in 2017 and so it was great to reconnect and see how much things have progressed.

I presented on accelerating genomic surveillance of malaria vectors in Africa. More than 23,000 African Anopheles genomes have now been sequenced through concerted efforts by MalariaGEN, PAMCA, GAARD, Target Malaria, PAMGEN and many others. The main question is, wow can we now build on this to address the challenges of resistance to dual AI LLINs and other evolutionary threats to vector control? My slides are here.

I also compiled some highlights below from other talks with links to papers where I could find them.

Drew Hammond (@DrewMarcHammond) showed that viable resistance to gene drives in Anopheles is possible even in ultra-conserved regions like doublesex. Resistance can be mitigated by multiplexing targets though. Paper coming soon from Mariaou et al.

Jonathan Kayondo gave an update on the @TargetMalaria project developing gene drives for Anopheles, and work happening on non gene drive male bias at @UVRIug in Uganda.

George Christophides (@ChristophidesG) presented transmission-blocking gene drives for malaria in the Transmission Zero project (@transm0) and work in Africa to develop and test a two-component gene drive system. See Hoermann et al. 2022.

Zhijian Jake Tu (@Zhijian_Tu) shared work on the sex determination system in Aedes. Female mosquitoes can be converted into fertile males by inserting just two genes on an autosome. Could be valuable for sex separation or gene drives to bias sex ratio. See Aryan et al. 2020 and Biedler et al. 2024.

Jeremy Herren (@JeremyHerren) shared work on microsporidia MB, a mosquito symbiont with potential to block malaria transmission. The @Symbio_Vector project has new data on how microsporidia infection affects mosquito gene expression and immunity. See Makhulu et al. 2023.

Nsa Dada spoke about the Anopheles mosquito microbiome and it’s potential role in insecticide resistance. Evidence that the microbiome differs between resistant and susceptible mosquitoes, replicated in two mosquito species. Is it causal? See Omoke et al. 2021 and Pelloquin et al. 2021.

Mary Ann McDowell spoke about population genetics and metagenomics of sandflies which transmit leishmania. Evidence for cryptic species in Brazil, associated with differences in male courtship song and microbiome. See Labbé et al. 2023.

Michael Branstetter (@bramic21) presented the Beenome100 project which is assembling genomes for bee species important for agriculture and conservation. See beenome100.org.

Matilda Gikonyo used trio binning to improve genome assemblies of small herbivorous insects with high genetic diversity, enabling resolution of maternal and paternal haplotypes.

Sarah Dyer presented the new rapid.ensembl.org site which releases new genomes every 2 weeks, keeping up with all the new genomes from tree of life and earth biogenome projects. Also the main @ensembl site has has a complete makeover, see beta.ensembl.org.

Joachim Nwezeobi presented on the population genomics of Anopheles funestus from across Africa. New results on population structure, chromosomal inversions, recent selection, copy number variation and adaptive gene flow. Paper coming soon from Joachim and @mariloubodde.

Tamar Carter is using population genomics to track the invasive malaria vector Anopheles stephensi. The spread of stephensi in Ethiopia seems to align with transport networks. Also stephensi in Yemen are structured into two distinct populations. See Carter et al. 2021, Samake et al. 2023 and Assada et al. 2024.

Fredrick Kebaso (@fredrickkebaso) assembled a genome for the camel ked, a biting fly that transmits various pathogens and affects the health of camels, making life difficult for marginalised communities in the Horn of Africa coping with climate change.

Billiah Bwana improved the annotation of genes in the tsetse fly expressed in antennae and involved in olfaction. Important because odour-baited traps are used for tsetse control. See Bwana et al. 2022.

Luna Kamau presented the discovery of Anopheles coluzzii in Kenya, which is genetically related to sahelian An. coluzzii from Mali and Burkina Faso and shares insecticide resistance alleles with these West African populations. See Kamau et al. 2024.

Megan Fritz (@MosquitoDoc) showed how the corn earworm has evolved resistance to transgenic (Bt) corn via copy number variation at a cluster of trypsin genes. Genomic monitoring could be used to detect resistance emerging in as little as 1 generation. See Pezzini et al. 2024 and Taylor et al. 2024.

Isaiah Debrah (@ISAIAHDEBRAH1) has found non-coding RNAs which are upregulated in insecticide resistant Anopheles funestus in western Kenya. The role of ncRNAs in regulating expression of detoxification genes needs more study. See Debrah et al. 2024.

Diana Omoke investigated gene expression profiles of Anopheles arabiensis populations from Western Kenya resistant to multiple insecticide classes. Resistant mosquitoes overexpressed salivary and cuticular proteins in addition to multiple detoxification genes. See Omoke et al. 2024.

Mahamat Gadji (@2017Kafir) reported genome-wide association analysis of Anopheles funestus from Mibellon in Cameroon highly resistant to pyrethroid insecticides. Strong selection at Cyp6 (rp1) and Cyp9 loci. See Gadji et al. 2024.

That’s it for the science. Going to be hard to beat this for a conference excursion, half a day at the Nairobi National Park…

Huge thanks to @SamRund, @Tung_Jowi, Mary Ann for organising a fantastic meeting, and to @pamcafrica, @KEMRI_Kenya, @ndeckinstitute, @Arthropod_i5k, @IGBIllinois, @kstatebio, @NIH, for making it possible.

Looking forward to the next AGS, wherever in the world it may be!

The Malaria Vector Genome Observatory is online

Fri, 07 Jun 2024 00:00:00 +0000

This post introduces malariagen.net/vobs, a new web page for the Malaria Vector Genome Observatory providing an entry point to data, training, software, research and more.

The main focus of my team and I over the last few years has been building the Malaria Vector Genome Observatory, which is a collaborative effort to create large-scale open data resources for research and surveillance of mosquitoes carrying malaria. The observatory is more than data, however, because we've also built software tools and training resources to help scientists in our community access and analyse these data. On top of that we've collaborated with a good number of other groups to analyse and publish research using these data.

Because the observatory is a large and relatively complex project, it's been difficult to communicate an overview of the whole project, and it's also been difficult for users and interested folks to find and navigate relevant information. Part of the problem was that the observatory didn't have a single web page giving an overview and entry point to all of these different parts, so we decided to make one, working together with the MalariaGEN communications team and the Sanger Institute web team. The result is now live at malariagen.net/vobs. This post gives a breakdown of the different parts and why they're there.

Here's the page header and navigation bar. Data is the main focus of the observatory, so this is front and centre, but there are five other components I wanted to highlight. The Training and API are key resources to help folks use the data. All data from the observatory is in Google Cloud, so I also wanted to include something about Compute resources that can be used to access data. Of course, the observatory is the result of contributions by many different Partner groups, and that needs to be represented. And there's an increasing amount of new Research coming out.

Problem statement¶

Next is a high-level statement of the big questions the observatory is aiming to help solve. This is mainly for people who've never heard of the observatory before, to give a sense of why it's relevant. But it's super hard to boil so much complexity down into a few short sentences, so I struggled with this. The different groups contributing to the observatory actually have a very broad range of interests and questions they're trying to answer, and there's also a range from basic biology and blue sky research through to very applied and short term disease control issues. That said, these three questions hopefully do a decent job of capturing the major focus of the vector observatory community and why we're interested in genomics.

Call to action¶

This is the call to action, a single statement trying to capture the central goal of the observatory. This was even harder to write than the problem statement! But the fun part was creating an image to go behind it. This was the first time I tried to use generative AI to do something serious, more on that below.

At a glance¶

Next is an overview of the six main components of the observatory, which correspond to the six parts of the page below. This is partly repeating information from the navigation bar, but I wanted to reinforce this structure, and provide a sentence for each to explain what it is. This also serves as another menu, people can click through to the section they're most interested in.

Data¶

Now we get into the main content. The first section covers data, giving some numbers to convey the size of the data resource, and providing links through to the data user guide which is a documentation site that has more technical information about the data and how to access it.

Training¶

Training comes next, this is a key entry point for many first time users of the data resource. The main link is through to the data analysis training course website.

API¶

We've built a Python API to access the data, but it's become much more than a simple access tool, there's a growing range of analytical and data visualisation functions, so I made a gallery of screenshots from various plots you can make with it to try and convey this capability. The bottom of the section has a link through to the API docs website.

Compute¶

The Compute section highlights three services, all of which can be used to access compute resources inside Google Cloud which then have good locality with the data. I put Colab first because it's free and so provides the best starting point for most users. But I also wanted to highlight some other options for folks who need more compute power. There are more options too like Vertex AI Workbench, but I thought I wouldn't overload.

Partners¶

Lots of people have contributed to the observatory, it's been a big community effort. I'd ultimately like to put an interactive map here, which would also be a bit easier to update as new contributors join.

Research¶

A key output from the observatory is new research, and so I wanted to highlight this. We're trying to keep this up to date, so this can be a reference point for papers using data from the observatory.

Implementation notes¶

Prototyping with Jekyll¶

I built a first prototype of the landing page using GitHub pages, starting from this one page creative Jekyll theme. It's a pretty old theme now and uses an older version of Bootstrap, but it's free and worked well for prototyping purposes. I was actually away from my computer when I built this, but I really had a bee in my bonnet and wanted to make something, so I built it on my phone with the GitHub file editor - not the easiest thing I've ever done, but sometimes you get an idea and it's hard to put it down! Working on my phone at least it forced me to learn something about the Bootstrap grid system and how to do mobile-first responsive design. The production landing page was then built by the MalariaGEN comms team (Jon and Sree) and the Sanger web team, using Wordpress which hosts the MalariaGEN website.

Generative AI images¶

I mentioned above I used generative AI to create one of the images. A colleague had pointed me at MidJourney previously and I'd had a tinker, but this seemed like a good excuse to try it out more extensively. I wanted something cool that gave a sense of scope and scale. The two themes I tried to prompt around where "A galaxy made of DNA" (riffing on the analogy with space observatories) and "A sea of data" (conveying that, hey, there's a shed-load of data here). Here are my two favourites.

A galaxy made from DNA¶

A vast swirling colourful sea of data¶

Bias¶

On the slightly less fun side, I also learned just how bad the gender and racial biases are in MidJourney. E.g., when I prompted with "computer scientist" almost all images came out with a bearded white man, unless I explicitly added information about diversity in the prompt. There's a lot written about this problem, e.g., here's one article talking about gender bias in MidJourney I found. Bottom line, if you're using generative AI, bear this in mind.

Future work¶

There's more to be done here, but hopefully this is a useful page for now. A couple of things in particular I'd like to add...

Contributing to the observatory¶

A key question that comes up often is how people can contribute to the observatory, and so I think we need some documentation about that, and to add a section to the landing page so it's easy to find.

Technical and scientific updates¶

We also need some way to surface and communicate updates, particularly about new data releases, new data features, new API releases, new research, etc., so there's an easy way to see what's new and get notified about important changes. We might need to make a separate site for that, and link through from the landing page, so it doesn't get cluttered.

20,000 Anopheles mosquito genomes

Thu, 16 Nov 2023 00:00:00 +0000

Thanks to the concerted efforts of many different contributing studies and institutions, the Malaria Vector Genome Observatory now contains analysis-ready data from whole-genome sequencing of 22,015 Anopheles mosquitoes. Within the last 12 months, more than 10,000 new genomes were added, doubling the total size of the data resource. The first data were also added for mosquitoes from the Anopheles funestus species subgroup, now accounting for 3,038 genomes.

To mark this milestone, here are three new visualisations of the data.

The first is a network depicting genetic relationships between mosquitoes from the Anopheles gambiae species complex. Each point represents a mosquito genome coloured by species. Each mosquito is connected to two other mosquitoes which are the closest in terms of genetic distance. The video shows an animation of this network being laid out using a force-directed graph.

I also created an alternative video showing the network using a different layout algorithm.

The second visualisation is a timeline showing when data releases where made available to contributing partner studies and how many mosquitoes each release contains. Data generation has accelerated particularly since 2020, thanks to support from the Bill and Melinda Gates Foundation and Wellcome.

The third visualisation is a map showing the geographical coverage of sampling locations in Africa. Each location is represented by a pie chart showing the number of samples by species. Data are available from 29 African countries, and a growing number of countries are now represented by longitudinal and subnational sampling.

A huge thank you to everyone who has supported and contributed to creating this data resource. So many people have gone to great lengths, overcome significant challenges and shown extraordinary patience, dedication and generosity to make this possible.

Training course in data analysis for genomic surveillance of African malaria vectors

Fri, 10 Mar 2023 00:00:00 +0000

This post introduces a new training course and explains some of the rationale and technology choices we made to create, publish and deliver the course online.

Part of our mission is to make malaria genomics as accessible as possible. Over the past year my team and I have been working hard to develop a training course that introduces scientists to the elements of mosquito genomic data analysis. The goal was to create a course that can help someone get as quickly as possible from little or no prior experience of genomic data analysis to a point where they can independently plan and run a suite of common analysis methods required for surveillance of mosquito populations. This post gives an introduction to the scope and structure of the course, along with some discussion of choices we made to develop, publish and deliver the course online.

Scope and structure¶

The course is structured into 8 workshops, each of which requires around 6-8 hours of learning time and so can be delivered in a day. Each workshop is focused around a single topic, such as population structure or detecting genes under recent positive selection. The final workshop then focuses on how to plan, execute and present a complete analysis. Here's the full workshop programme:

The topics we cover include some basic population genetics, such as analysing population structure, identifying crypting species and quantifying genetic diversity, and evolutionary threats to malaria vector control, such as analysing insecticide resistance mutations and how they spread between mosquito populations. There is obviously a very applied focus for these topics, as most members of our community are interested in using genomics to help monitor and control malaria vectors better.

Each workshop is divided into 4 modules which cover the topic from different angles. The modules have different themes. Everything builds towards the Analysis module, where we learn how to run one or more commonly needed analysis methods on some real data from the Malaria Vector Genome Observatory. However, I wanted to also include modules on the underlying Tools & Technology, Biology and Data so that trainees get some insight into everything that goes into that analysis. E.g., here are the modules from Workshop 2:

Some workshops also have a Journal Club module, where an existing paper relevant to the topic is presented by one of the authors, so trainees see how the ideas and methods have been applied.

Each module comprises a lecture video in French or English plus a lecture notebook which includes fully worked and executable code examples. After watching the lecture video, trainees launch the notebook in Google Colab and execute all the code examples for themselves using real genomic data, attempting some additional exercises along the way.

Executable lecture notebooks¶

The lecture notebooks from all workshops and modules have been used to build the training course website. This website was built using Jupyter Book, a relatively new technology that fortuitously arrived just before we started this project, and has been great. The source code for the training website is hosted on GitHub and the website is built and deployed to GitHub Pages. If you browse into the docs folder you'll see a file called _toc.yml which defines the table of contents for the site, and _config.yml which provides the site configuration. There is then one folder for each workshop (e.g., workshop-2), and inside each folder are some Jupyter notebooks, one for each workshop module. Jupyter Book renders each notebook to an HTML file and organises all the files into a static website according to the table of contents you define. As course developers, all we needed to do was create the notebooks and update the table of contents, and everything would be rendered and published automatically via a GitHub action.

A super nice feature of Jupyter Book is that it adds a button to each page to launch the notebook in Google Colab:

In one click, anyone following the training course can have the code examples opened and ready to run in Google Colab. All examples use real data which is hosted in Google Cloud Storage, and use the malariagen_data Python API which accesses data directly from GCS, so there is no need to move any data around before starting an analysis, and a single pip install command brings in all the software needed. This makes the transition from passive learning to active learning via practical, hands-on analysis about as easy as you can make it.

English and French¶

A substantial fraction of our community is based in francophone Africa, and we really wanted to do something to make the course accessible to different language groups. I'm extremely fortunate to have a native French speaker, Jon Brenas, in my team, as well as support from Eric Lucas at LSTM who is also bilingual. So we decided to translate and rerecord all of the lecture videos from English into French. There are 32 modules in total, each with a lecture between 30-50 minutes in length, often with very technical language. Jon was also often working to tight deadlines, as we typically only managed to finish the English content for each workshop about a week before we were due to deliver the workshop to the first cohort of trainees. It was a heroic effort, Jon is a legend.

We decided not to translate all of the text in the lecture notebooks as well, as this would have meant maintaining two different versions of each notebook, which would have been too much. However, we did translate all of the practical exercises into French within each notebook. At some point it would be amazing if Jupyter Book had some kind of support for internationalisation, so you could have text for different languages within the same notebook/page.

Running the course online¶

All of the lecture videos and notebooks are available via the course website, so anyone can follow the course in their own time via self-directed learning. But we also wanted to run the course as a series of online workshops, to provide a more structured experience with support from experienced teaching assistants. To make this happen we worked together with Vikki Simpson's team at GSU, particularly the wonderful Paballo Chauke, and with the PAMCA Anopheles genomics team led by Elijah Juma.

We decided to deliver the workshops via Zoom, because it was the platform we were most familiar with, and supports breakout rooms and video sharing. We also decided to prerecord all of the lecture videos, primarily so we could translate them into French, but also so that we could easily rerun the course multiple times in future, and publish them on YouTube to provide a self-directed learning option. At the beginning of each workshop, all trainees would log into Zoom, and we started with a plenary opening session of around 30 minutes where we welcomed everyone, gave a high-level introduction to the main topic of the workshop, and gave some logistical information about how the workshop would be run. For each module, we then split the 50 trainees into 5 pre-assigned breakout rooms, with around 10 trainees and 2 teaching assistants per room. One of the teaching assistants would play the lecture video, and then we would allow 30-50 minutes for trainees to run the lecture notebook for themselves and attempt the practical exercises, asking questions or support from TAs as needed. We ran two modules back-to-back in the morning, had a break for lunch, then ran two more modules in the afternoon, before bringing everyone back to the main room to give a short wrap-up talk and provide some final information to trainees.

Overall the format worked well, and feedback from the course about the way it was delivered was generally positive. The main feedback we received was that trainees would like more live interaction with lecturers, as watching prerecorded videos is not as stimulating. We were all conscious of this, but had to make a compromise here, as it would have been logistically challenging to deliver all lectures live. During the practical sessions, trainees were often quiet and perhaps tentative to ask questions, and it certainly was not the same as all being in the same computer lab and being able to walk around and chat or see what people are attempting or stuck on. But again we had to compromise, with trainees spread across many countries there is obviously a huge accessibility benefit to offering the course online. Any suggestions for how to improve interactivity online and help break the ice would be very welcome.

We're planning to run the course again later this year, please keep an eye out on Twitter or LinkedIn for the announcement if you know someone who might be interested.