2.3 AnoExpress: is your gene overexpressed elsewhere?¶

Theme: Analysis

In module 2.1 you found thousands of genes differentially expressed between the Busia colony and the susceptible Kisumu strain, and you were given two or three candidate genes. One experiment cannot tell you which of those genes matter for resistance. Many will differ because Busia and Kisumu are different colonies, kept in different insectaries for different numbers of generations.

A gene that is overexpressed in many independent resistant populations, compared with many different susceptible strains, is a stronger candidate. In this module you use AnoExpress, a meta-analysis of RNA-seq studies of insecticide resistance in Anopheles, to check whether your genes replicate. The Busia colony is one of the experiments in AnoExpress, so you can see where Busia sits among the others.

Our worked example is CYP6AA1 (AGAP002862). It is duplicated on the “triple mutant” haplotype that is spreading through Uganda and Kenya (a Cyp6aa1 duplication, Cyp6p4 I236M and a transposable element insertion) [7], and that haplotype sits in the Busia 2R selective sweep. We run each analysis on CYP6AA1 first; you then repeat it on your own genes.

Learning objectives¶

At the end of this module you will be able to:

• explain why replication across studies helps separate resistance genes from colony effects, and what limits a meta-analysis of existing RNA-seq studies;

• plot the expression of your candidate genes across all AnoExpress experiments and find the Busia contrast among them;

• run a genome-wide expression scan and spot clusters of overexpressed genes around your candidates;

• look at the gene family and the co-expressed partners of a candidate, and run a simple enrichment test;

• record the Replication section of your gene dossier.

Setup¶

We also load the full candidate gene pool, so we can compare everyone’s genes. The Busia contrast in AnoExpress is called BusiaSurvivors_v_Kisumu.

Why a meta-analysis?¶

A single RNA-seq experiment compares one resistant colony with one susceptible strain. Any difference between the two, not just resistance, can make a gene differentially expressed: genetic background, time in colony, diet, age, rearing conditions and batch. In the Busia RNA-seq, over 5,000 genes were significantly different between the parental colony and Kisumu [2]. Most of them are not resistance genes.

Replication is the simplest filter. If a gene is overexpressed in resistant mosquitoes from Burkina Faso, Côte d’Ivoire, Cameroon and Uganda, each compared with a different susceptible strain, colony effects become a less likely explanation. This logic was used to find the sensory appendage protein SAP2 [4] and other “meta-signature” genes from microarray data [3].

AnoExpress [1] collects RNA-seq studies of insecticide resistance in An. gambiae, An. coluzzii, An. arabiensis and An. funestus. The count data from each study was re-analysed in the same way: each resistant vs susceptible comparison was run separately in DESeq2 [5], and the log2 fold changes and adjusted p-values were collected into one table. Older microarray studies from IR-TEx [3] can be added with microarray=True.

It helps to be clear about what this is and is not:

• It is a collection of comparable fold changes, not a pooled statistical test. You judge replication by counting and looking, not from a single meta-analysis p-value.

• Fold changes are relative. A log2FC of 3 means 8-fold higher than that susceptible strain. If the reference strain has unusually high or low expression, the fold change changes with it.

• Counts matter too. A gene that goes from 0 to 10 reads has a huge fold change but may be biologically trivial. Always look at counts alongside fold changes.

• Analyses with more species contain fewer genes, because only genes with orthologs in all included species are kept. For Busia (An. gambiae) the gamb_colu analysis keeps the most genes.

The experiments¶

Let’s look at what is in AnoExpress. Each row is one resistant vs susceptible comparison.

Notice the heterogeneity in study design:

• Populations: resistant samples range from long-established lab colonies to recently colonised field populations (like Busia). Some were exposed to insecticide before sampling (the Busia survivors), others were unexposed (the Bouaké samples).

• References: each species uses one or two old, susceptible lab strains (Kisumu, Ngousso, Mali, Moz, Dongola, Fang). These strains have their own history and quirks [6].

• Methods: read counting differed between studies (HISAT2, STAR, kallisto, featureCounts). Tissue, age and sex (whole adult females, usually) are broadly similar, but not identical.

• Geography: West and Central Africa dominate. Busia is the only An. gambiae comparison from East Africa.

Keep this in mind: replication across West African An. coluzzii studies is useful evidence for a Ugandan An. gambiae gene, but it is not the same as replication in East Africa.

Your genes across all experiments¶

plot_gene_expression() draws two panels. On the left, each point is the log2 fold change in one comparison, coloured by species. On the right, each point is the read count in one resistant sample. Hover over the points to see which experiment they come from.

Plots are good for looking; tables are better for your dossier. The function below counts, for each gene, how many gamb_colu comparisons show significant overexpression (padj < 0.05 and log2FC > 0) or underexpression, and pulls out the Busia contrast.

CYP6AA1 is a surprise. It is significantly overexpressed in 10 of the 18 gamb_colu comparisons, and strongly so in the six Côte d’Ivoire An. coluzzii comparisons (log2FC about +2.9 to +3.8, so 7- to 14-fold). Yet in BusiaSurvivors_v_Kisumu it is not significant (log2FC −0.29). Its neighbour CYP6P3 is even lower in Busia (log2FC −3.4). GSTE4, by contrast, is modestly but significantly up in Busia.

So we cannot claim from this contrast that CYP6AA1 is overexpressed in Busia. Why might that be? Look at the counts in the Busia experiment alone (CYP6P3 is added for comparison).

CYP6AA1 is highly expressed in both groups (log2 counts about 9.5 to 10). CYP6P3 is expressed in Busia, but this batch of Kisumu expresses it more. Several explanations are possible, and they are worth writing down in your dossier:

1. The reference. The fold change is only as good as the comparison. A susceptible strain with high baseline expression of a gene hides overexpression in the resistant colony. The West African comparisons mostly use Ngousso or Mali as the reference, not Kisumu.

2. The haplotype may be rare in the colony. If the Cyp6aa1 duplication was rare or lost when Busia was colonised, the colony would not overexpress it even though wild Busia mosquitoes carry it. You can check its frequency in wild Busia mosquitoes on Day 3, and in the colony from the RNA-Seq-Pop results.

3. The mechanism may not be (only) expression. A sweep can be driven by an amino acid change (such as Cyp6p4 I236M, which sits on the same haplotype), or by expression in a specific tissue or after exposure only.

4. Colony history. Busia had been in colony for about 24 generations; resistance alleles can be lost or fixed during colonisation, and the survivors were selected for only a few generations.

The RNA-Seq-Pop run from Day 1 also gives selected vs parental, which removes the Kisumu reference altogether. Compare your genes in both contrasts.

Where does Busia sit?¶

A useful overview is to plot every gene’s fold change in Busia against its median fold change in all the other experiments. Genes in the top-right are overexpressed in Busia and elsewhere; genes high on the y-axis only are Busia-specific. Candidate pool genes are highlighted.

Exercise 1¶

For each of your candidate genes, answer the Replication questions from your dossier:

1. Is the gene significantly overexpressed in Busia vs Kisumu? By how much (fold change, not log2)?

2. In how many of the other gamb_colu comparisons is it significantly over- and underexpressed?

3. Which countries and species do those comparisons come from? (Hint: use fc.loc[gene], padj.loc[gene] and meta.)

4. Does the gene exist in the gamb_colu_arab_fun analysis? If yes, is it also overexpressed in An. arabiensis or An. funestus?

Genome-wide expression scans¶

Genes that sit next to each other are often overexpressed together. There are two common reasons:

• Copy number variants (CNVs). A duplication that includes several genes raises the expression of all of them. The Cyp6aa/Cyp6p cluster on 2R and the Gste cluster on 3R both carry many CNVs in African An. gambiae [7].

• Shared regulation. Neighbouring genes can share enhancers or chromatin state.

A genome-wide expression scan plots the fold change of every gene along a chromosome arm, with a moving median across windows of neighbouring genes. Peaks in the line show regions where many neighbouring genes are overexpressed across experiments.

xpress.contig_expression() returns the fold change of every gene on an arm in every experiment, plus the moving median. The function below plots it with the Busia contrast in red, the Busia selection signals from the selection atlas (Day 3) in gold, and your candidate genes as dashed lines.

The moving median peaks at the Cyp6aa/Cyp6p cluster (about 28.5 Mb), right in the middle of the Busia selection signal. The peak is driven by the other experiments (grey). The Busia points (red) at the cluster sit below zero, while a few genes nearby, at about 28.3 Mb, are strongly up in Busia.

The scan mixes all experiments. To see what Busia does in a region, pull out the fold changes for every gene in a span and draw a heatmap. Genes are sorted by position.

AnoExpress uses compound contigs for regions: 2RL is 2R followed by 2L, and 3RL is 3R followed by 3L. Positions on 2R and 3R are the same as usual; positions on 2L and 3L are shifted by the length of 2R (61,545,105 bp) or 3R (53,200,684 bp). region_around() below handles this for you.

The Cyp6aa and Cyp6p genes are marked with ▶. CYP6AA1, CYP6AA2, CYP6P3 and CYP6P4 are red in most West African experiments, which fits a regulatory or CNV mechanism spanning several genes, like the Cyp6aa1 duplication. In the Busia column (the last one) they are white or blue. Yet this region is the centre of a strong selection signal in Busia (Day 3). A little further out, but still inside that signal, several genes are strongly up in Busia: AGAP002850, AGAP002852 and CYP6Z4 (AGAP002894) are all in the candidate pool.

Look closely at the six Côte d’Ivoire columns (Agboville, Dabou, Tiassale): CYP6P1 to CYP6P5 have identical fold changes. That suggests reads from these very similar genes were counted together in that study. Processing choices in the original studies carry through into any meta-analysis.

Exercise 2¶

1. Run a genome-wide expression scan for the chromosome arm that carries your first candidate gene. Is the gene in a region where many neighbours are also overexpressed?

2. Draw a heatmap for 100 kb either side of each of your candidates using region_around() and region_heatmap(). Are neighbouring genes up in Busia too?

3. What would you look for on Day 3 to test whether a CNV explains the pattern?

Gene families¶

Most metabolic resistance genes belong to a few large families: cytochrome P450s, glutathione S-transferases (GSTs), carboxylesterases (COEs), UDP-glucuronosyltransferases (UGTs) and ABC transporters. Looking at the whole family shows whether your gene stands out from its relatives, or whether the whole family moves together.

plot_gene_family_expression() takes a Pfam domain (short name) or a GO term, and plots every gene that carries it. Useful Pfam names from the AnoExpress notebooks:

Co-expression¶

Genes that go up and down together across many experiments may share a regulator, sit on the same CNV, or work in the same pathway. A simple way to find co-expressed genes is to correlate fold changes across experiments. With only 18 gamb_colu comparisons the correlations are noisy, so treat the top hits as leads, not conclusions.

The AnoExpress gene regulatory network notebook goes further, using GENIE3 to link transcription factors to their likely targets.

The same_region column flags hits within 500 kb of your gene. Neighbouring genes at the top of the list point towards a shared CNV or regulatory region.

Exercise 3¶

1. List the genes most co-expressed with each of your candidates. Are any of them known resistance genes (P450s, GSTs, COEs, Sap genes)?

2. Are any of the co-expressed genes physical neighbours? What does that suggest?

3. Plot the fold changes of your gene and its top three partners with plot_gene_expression().

Enrichment¶

Finally, zoom out. Which kinds of genes are overexpressed in Busia, and which are overexpressed consistently across experiments? A hypergeometric test asks whether a gene list contains more members of a GO term or Pfam domain than expected by chance, given how many genes carry it genome-wide. AnoExpress corrects the p-values for multiple testing (Benjamini-Hochberg).

Here we test the genes that are significantly overexpressed in at least 10 of the 18 gamb_colu comparisons.

Exercise 4¶

1. Repeat the Pfam enrichment for genes that are significantly overexpressed in Busia with log2FC > 2. Which gene families are enriched? Are they the same as in the consistent set?

2. Run a GO enrichment on the same list with xpress.go_hypergeometric(). Which terms come up?

3. Is your gene in an enriched family? Does that make it a stronger or a weaker candidate?

Exercise 5: update your dossier¶

Write the Replication section of your gene dossier for each candidate. Include:

• the Busia log2FC and adjusted p-value;

• the number of gamb_colu comparisons with significant over- and underexpression, and where they are from;

• whether neighbouring genes are co-overexpressed (possible CNV), and the gene’s top co-expressed partners;

• your judgement: is the Busia result replicated, Busia-specific, or contradicted elsewhere?

Quiz¶

Summary¶

• One RNA-seq experiment cannot separate resistance genes from colony differences. Replication across independent studies is a cheap, powerful filter.

• AnoExpress gives comparable fold changes and adjusted p-values for many resistant vs susceptible comparisons, including Busia. It is a collection, not a pooled test: count, look and think.

• Fold changes depend on the reference strain. CYP6AA1, our worked example, is up 7- to 14-fold in Côte d’Ivoire but unchanged in Busia vs Kisumu, even though it lies in the Busia 2R sweep and is duplicated on the East African triple-mutant haplotype.

• Genome-wide expression scans and region heatmaps reveal clusters of co-overexpressed neighbours, which often point to CNVs, the topic of Day 3.

• Gene families, co-expression and enrichment put a candidate in context, but they add prior plausibility rather than direct evidence.

Well done!¶

You now know whether your genes replicate. Tomorrow you ask why they are overexpressed, and whether they are under selection in wild Ugandan and Kenyan mosquitoes.

References¶

1. Nagi SC, Ingham VA. AnoExpress: Anopheles gene expression in resistance studies. Python package and Colab notebooks, https://github.com/sanjaynagi/AnoExpress (version 0.3.3).

2. Nagi SC, Oruni A, Weetman D, Donnelly MJ (2023). RNA-Seq-Pop: Exploiting the sequence in RNA sequencing, a Snakemake workflow reveals patterns of insecticide resistance in the malaria vector Anopheles gambiae. Molecular Ecology Resources 23: 946–961. Nagi et al. (2023)

3. Ingham VA, Wagstaff S, Ranson H (2018). Transcriptomic meta-signatures identified in Anopheles gambiae populations reveal previously undetected insecticide resistance mechanisms. Nature Communications 9: 5282. Ingham et al. (2018)

4. Ingham VA, Anthousi A, Douris V, et al. (2020). A sensory appendage protein protects malaria vectors from pyrethroids. Nature 577: 376–380. Ingham et al. (2019)

5. Love MI, Huber W, Anders S (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15: 550. Love et al. (2014)

6. Ingham VA, Morgan J, Lees RS, Ranson H (2019). Characterisation of Anopheles strains used for laboratory screening of new vector control products. Parasites & Vectors 12. Williams et al. (2019)

7. Njoroge H, van’t Hof A, Oruni A, et al. (2022). Identification of a rapidly-spreading triple mutant for high-level metabolic insecticide resistance in Anopheles gambiae provides a real-time molecular diagnostic for antimalarial intervention deployment. Molecular Ecology 31: 4307–4318. Njoroge et al. (2022)