2.1 RNA-seq: differential expression in Busia¶

Theme: Analysis

Yesterday you launched RNA-Seq-Pop on a subsampled copy of the Busia data. Today we work with the results of the same workflow (v2.3.0, same config.yaml) run on the full dataset: 16 libraries, ~50 million read pairs each.

We start with the theory you need to trust an RNA-seq result: how reads become counts, why counts must be normalised, how DESeq2 tests for differential expression and why we correct for multiple testing. Then we explore the Busia results: quality control, sample clustering, volcano plots for both contrasts, and gene set enrichment. At the end you are assigned the candidate genes that you will follow for the rest of the week.

Learning objectives¶

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

• Explain how RNA-seq reads become gene counts (library preparation, kallisto pseudoalignment, counts vs TPM).

• Describe how DESeq2 normalises counts and tests for differential expression, and interpret log2 fold change, adjusted p-values and shrinkage.

• Check RNA-seq quality control and sample clustering, and make volcano plots from RNA-Seq-Pop output.

• Interpret a hypergeometric gene set enrichment test.

• Find your candidate genes in the Busia results.

Setup¶

From mosquitoes to reads¶

An RNA-seq experiment measures which genes are being transcribed, and how much, at one moment. For the Busia libraries [1]:

1. RNA extraction. Total RNA from a pool of five 4-day-old females per library. Pooling averages out individual variation and gives enough RNA; replicates (4–6 pools per group) capture the variation between pools.

2. mRNA enrichment. Poly-A selection keeps mRNA and removes most ribosomal RNA.

3. Fragmentation and cDNA synthesis. The mRNA is broken into short pieces and reverse-transcribed into cDNA.

4. Adapters and PCR. Sequencing adapters (with a sample index) are ligated and the library amplified.

5. Sequencing. 150 bp paired-end reads on an Illumina NovaSeq 6000.

Every step can add bias (PCR duplicates, 3′ bias from degraded RNA, GC content). That is why we look at QC before trusting any result. And a key point: a gene with more mRNA gives more fragments, but so does a longer gene, and so does a library that was simply sequenced more deeply.

From reads to counts: pseudoalignment with kallisto¶

To count reads per transcript we could align every read to the genome base by base. Kallisto [2] takes a faster route. It builds an index of all k-mers (31 bp words) in the transcriptome. For each read it looks up its k-mers and finds the set of transcripts that contain all of them: the read’s compatibility class. It never works out the exact alignment position, hence pseudoalignment. An expectation–maximisation (EM) algorithm then shares reads that are compatible with several transcripts (e.g. isoforms of one gene) among them.

Kallisto reports for each transcript:

• est_counts: the estimated number of reads from that transcript.

• tpm: transcripts per million, which corrects for transcript length and sequencing depth.

where $c_i$ is the count and $\ell_i$ the effective length of transcript $i$. TPMs are good for comparing genes within a sample. For testing differences between groups, DESeq2 wants the raw counts, because its statistical model is built on counts. RNA-Seq-Pop sums est_counts over each gene’s transcripts, rounds them, and passes them to DESeq2.

Kallisto also runs 100 bootstraps (-b 100), re-estimating counts on resampled reads. Sleuth uses these to test individual isoforms (results/isoformdiff/).

The short gene has four times the TPM of the long gene, even though both have 400 reads: per unit of length, it is more highly expressed.

Quality control¶

Before any biology, check the data. RNA-Seq-Pop collects QC from fastp, kallisto, HISAT2, samtools, mosdepth and bcftools into one MultiQC report, results/qc/multiQC.html. Download it (in Colab: the folder icon on the left, then busia-results/results/qc/) and open it in your browser. Look at the fastp read quality and adapter content, and the HISAT2 alignment rates.

Here we plot the key numbers directly. First, how many reads kallisto pseudoaligned to the transcriptome:

Mean genome-wide depth looks low because most of the genome is not transcribed. Coverage is very high in expressed exons and zero elsewhere, which matters for variant calling (Module 2.2).

Normalisation¶

Raw counts are not comparable between samples. Two things differ:

• Sequencing depth. A library with twice as many reads has roughly twice the count for every gene.

• Library composition. If a few genes are hugely expressed in one sample, they use up reads and every other gene looks lower, even if its expression has not changed.

Dividing by the total count (counts per million) handles depth but not composition. DESeq2 uses the median-of-ratios method [3]:

1. For each gene, compute the geometric mean of its counts across all samples (a pseudo-reference sample).

2. For each sample, divide every gene’s count by that reference.

3. The sample’s size factor is the median of those ratios. Most genes are not differentially expressed, so the median is robust to the few that are.

4. Normalised counts = raw counts / size factor.

RNA-Seq-Pop writes both results/counts/rawcounts.tsv and results/counts/normCounts.tsv. Let’s compute the size factors ourselves and check them against DESeq2’s.

Do the samples cluster by group?¶

Before testing genes, check that replicates look alike and groups look different. RNA-Seq-Pop plots a PCA of variance-stabilised counts (results/counts/PCA.pdf). We do a similar PCA on log-transformed normalised counts of the 2,000 most variable genes.

In the paper, replicates of each group clustered together in both the PCA and the correlation heatmap [1]. Which group is most different? Does that fit with Kisumu being colonised in 1975 and the Busia colony in 2018?

The DESeq2 model¶

Why a negative binomial?¶

If sequencing were the only source of noise, the count for a gene would follow a Poisson distribution, whose variance equals its mean. Real biological replicates vary more than that: different pools of mosquitoes really do express genes at slightly different levels. DESeq2 [3] models counts with a negative binomial distribution:

where $\mu_{ij} = s_j q_{ij}$ is the size factor times the true expression, and $\alpha_i$ is the dispersion of gene $i$. With only 4–6 replicates, each gene’s dispersion is hard to estimate, so DESeq2 shares information across genes: it fits a trend of dispersion against mean expression and shrinks each gene’s estimate towards it.

Let’s check that the Busia replicates are over-dispersed:

For well-expressed genes the variance is far above the Poisson line: that extra variance is what the dispersion parameter captures.

Design, contrasts and the Wald test¶

RNA-Seq-Pop fits one model to all 16 samples with the design ~ treatment, then extracts each contrast. For BusiaParental_BusiaSelected, the coefficient is

so positive log2FC means higher expression in the second (case) group. A Wald test divides the log2FC by its standard error (lfcSE) to give the stat column and a p-value.

Multiple testing and the false discovery rate¶

We test ~13,000 genes. At p < 0.05, about 650 genes would be “significant” by chance even if nothing had changed. DESeq2 reports padj, p-values adjusted with the Benjamini–Hochberg method, which controls the false discovery rate (FDR): among genes with padj < 0.05, we expect about 5% to be false positives. DESeq2 also applies independent filtering: genes with very low counts are given padj = NA because they have almost no power to be significant, and dropping them makes the correction less severe for the rest.

Log2 fold change shrinkage¶

A gene with 2 reads in one group and 10 in the other has a log2FC of ~2.3, but that estimate is extremely noisy. DESeq2’s lfcShrink (e.g. with the apeglm method) pulls such estimates towards zero in proportion to their uncertainty, which makes ranking genes by fold change more reliable. RNA-Seq-Pop v2.3.0 reports unshrunken log2FC, so be careful with genes that have a low baseMean (mean normalised count) or a large lfcSE.

The differential expression tables¶

RNA-Seq-Pop writes one table per contrast, results/genediff/{contrast}.csv, and all of them as sheets of results/genediff/Ag_Busia_diffexp.xlsx. The columns are:

Worked example: CYP6AA1¶

Throughout the course we show each method first on CYP6AA1 (AGAP002862), then you repeat it on your own genes. CYP6AA1 is a cytochrome P450 in the Cyp6aa/Cyp6p cluster on chromosome arm 2R (~28.48 Mb), a region under strong selection in East Africa. In Uganda and Kenya a “triple mutant” haplotype has spread that combines a duplication of Cyp6aa1, the missense mutation Cyp6p4-I236M and a ZZB transposable element insertion, and it is associated with pyrethroid resistance [9]. This haplotype sits in the Busia 2R selective sweep you will meet on Day 3, and Ag-vampIR markers for it predict survival on deltamethrin nets in Siaya (Day 4).

A gene duplication matters for RNA-seq. With two copies of the gene, there is more template to transcribe, so reads over a duplicated gene inflate its apparent expression: a higher count may reflect copy number rather than stronger regulation. The same duplication shows up in the genome alignments as elevated coverage over the gene. RNA-seq alone cannot tell the two apart; on Day 3 you check CNV frequencies in wild Busia mosquitoes with WGS data.

Let’s look at CYP6AA1 in both contrasts: first parental vs selected (did selection change it?), then Kisumu vs parental (is it different between the strains?).

And the normalised counts in each sample. We define a small plotting function here and reuse it for other genes below.

Interpret the result with care. In the AnoExpress meta-analysis of the same Busia libraries, CYP6AA1 is not significantly different between selected survivors and Kisumu (log2FC −0.29), so we cannot claim it is overexpressed relative to Kisumu from that comparison. Compare the two RNA-Seq-Pop contrasts: a gene can respond to selection within the colony while showing no difference from a lab strain, especially one like Kisumu whose own history is unusual (Module 2.2). Write down the log2FC and padj for both contrasts in your dossier.

Volcano plots¶

A volcano plot shows effect size (log2FC) against evidence (−log10 padj). Genes top-left and top-right are strongly and significantly changed.

In the paper [1], all three sensory appendage protein genes (Sap1–3) were overexpressed in the selected line compared with its parent, with Sap2 up ~10.7-fold. SAP2 binds pyrethroids and was first linked to resistance in West African An. coluzzii [4]. Many cytochrome P450s, carboxylesterases and ABC transporters were also differentially expressed.

Noisy fold changes at low expression¶

Plot log2FC against mean expression (an MA plot) and colour by the standard error. The widest spread of fold changes is at low baseMean, where the estimates are least reliable. This is what shrinkage would correct.

Comparing the two contrasts¶

Some genes change in the same direction in both comparisons: higher in the Busia parental line than in Kisumu, and higher again after selection. The paper highlighted a carboxylesterase cluster on 2L (Coebe2c, Coebe3c, Coebe4c) with this pattern.

Look at the counts¶

Always go back to the normalised counts for a gene you care about. A significant result driven by one odd replicate looks very different from a consistent shift.

Gene set enrichment¶

A list of 5,000 DE genes is hard to interpret one by one. Enrichment analysis asks whether genes with a particular function are over-represented among them. Functions come from Gene Ontology (GO) terms, which describe molecular functions, biological processes and cellular components.

RNA-Seq-Pop uses over-representation analysis with the hypergeometric test. For each GO term:

• $N$ = genes with any GO annotation (the universe)

• $K$ = genes annotated with this term

• $n$ = genes in our list (RNA-Seq-Pop uses genes with padj < 0.05 and fold change > 2)

• $k$ = genes in our list that have this term

If the list were a random draw of $n$ genes, $k$ would follow a hypergeometric distribution. The p-value is the probability of seeing $k$ or more: $P(X \geq k)$. With thousands of GO terms tested, p-values are again corrected with Benjamini–Hochberg.

This differs from rank-based GSEA (as in the GSEA and fgsea tools), which uses the whole ranked list of genes rather than a significance cut-off. Over-representation depends on the thresholds you choose; rank-based methods do not, but ask a subtly different question.

A worked example: suppose 8,000 genes have GO annotations, 100 of them are “monooxygenase activity” (mostly P450s), and 400 genes are strongly upregulated after selection, 15 of which are monooxygenases.

RNA-Seq-Pop writes the results to results/gsea/genediff/{contrast}.de.tsv, with columns annotation (GO ID), pval, padj and descriptions. For the course run, GO terms come from the eggNOG-mapper annotations that ship with RNA-Seq-Pop (the paper used VectorBase’s GO file).

RNA-Seq-Pop also runs the same test on the genes with the top 5% of Fst values (results/gsea/fst/). We come back to Fst in Module 2.2.

Your candidate genes¶

Over the rest of the week you follow two or three genes from Busia to a field assay. The course pool is in data/candidate-genes.tsv. Each gene is overexpressed in Busia; the pool mixes known insecticide resistance genes (category = known_ir: P450s, GSTs, carboxylesterases) with less studied genes (category = novel). CYP6AA1 is the focal gene that everyone works through as the worked example. Other columns tell you whether the gene is overexpressed in other AnoExpress experiments, lies in a Busia selective sweep, or is near an Ag-vampIR amplicon; you will check these yourself on Days 2–4.

The busia_log2fc and busia_padj columns come from the AnoExpress meta-analysis of the same Busia data, where the contrast is selected survivors vs Kisumu. RNA-Seq-Pop gives us the two steps of that comparison separately. Because DESeq2 fits one model to all samples, the two log2FCs add up: log2(Selected/Kisumu) = log2(Parental/Kisumu) + log2(Selected/Parental).

How genes are assigned. Your facilitator gives you a participant number. You get one known IR gene and two novel genes, chosen so that the class covers the whole pool. (Self-paced learners: pick any number.) Copy your genes into the CANDIDATE_GENES cell at the top of every notebook from now on.

Start your gene dossier (assessment/gene-dossier-template.md) now: for each gene record its log2FC and padj in both contrasts, and whether its counts look consistent across replicates.

Exercise 1¶

In the BusiaParental_BusiaSelected contrast, how many genes are significantly upregulated (padj < 0.05, log2FC > 1) that are P450s (GeneName starting with CYP), carboxylesterases (COE), GSTs (GST) or ABC transporters (ABC)?

Exercise 2¶

Is Vgsc (AGAP004707), the pyrethroid target, differentially expressed between the parental and selected lines? What does this tell you about the limits of an expression-only analysis?

Exercise 3¶

Check that log2FCs are additive: plot log2FC_Sel_vs_Kis (computed above) against the AnoExpress busia_log2fc for the candidate genes. Are they similar? Why might they not match exactly?

Exercise 4¶

For each of your candidate genes, is the change mainly between Kisumu and the parental line, or during selection? Look at the volcano plots, the table above and plot_genes. Note your answer in your dossier.

Quiz¶

Check your understanding. The quiz runs in Colab and in the course book.

Summary¶

• Kallisto pseudoaligns reads to the transcriptome; RNA-Seq-Pop sums transcript counts per gene for DESeq2. TPMs compare genes within a sample; raw counts are used for testing.

• DESeq2 normalises with size factors (median of ratios), models counts with a negative binomial with gene-wise dispersion, and tests contrasts with a Wald test. log2FC > 0 means higher in the second group of the contrast.

• Use adjusted p-values (FDR); treat large fold changes of lowly expressed genes with caution, because RNA-Seq-Pop does not shrink log2FC.

• In Busia, thousands of genes differ between groups; SAP genes, P450s and a carboxylesterase cluster are upregulated after selection, while Vgsc is not differentially expressed.

• Over-representation analysis uses a hypergeometric test on a thresholded gene list.

• A duplicated gene, such as Cyp6aa1 on the Uganda/Kenya triple-mutant haplotype, can look overexpressed simply because it has more copies.

Well done!¶

You now have your candidate genes. In Module 2.2 we look at what the sequence in the same reads tells us about them, and in Module 2.3 whether they are overexpressed in other studies.

References¶

1. 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)

2. Bray NL, Pimentel H, Melsted P, Pachter L (2016). Near-optimal probabilistic RNA-seq quantification. Nature Biotechnology 34:525–527. Bray et al. (2016)

3. 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)

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

5. Zhu A, Ibrahim JG, Love MI (2019). Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences. Bioinformatics 35:2084–2092. Zhu et al. (2018)

6. Benjamini Y, Hochberg Y (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society B 57:289–300. Benjamini & Hochberg (1995)

7. Pimentel H, Bray NL, Puente S, Melsted P, Pachter L (2017). Differential analysis of RNA-seq incorporating quantification uncertainty. Nature Methods 14:687–690. Pimentel et al. (2017)

8. Subramanian A et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. PNAS 102:15545–15550. Subramanian et al. (2005)

9. Njoroge H 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)