2.2 The sequence in RNA-seq: variants, selection and ancestry¶

Theme: Analysis

In Module 2.1 we treated RNA-seq reads as counts. But every read is also a piece of sequence. The Busia libraries cover almost every expressed gene at high depth, so we can call SNPs, measure allele frequencies of known resistance mutations, scan for genetic differentiation between groups, estimate genetic diversity, and even infer chromosomal inversions and species ancestry.

In Busia this matters: the strongest single change during deltamethrin selection, the rise of the Vgsc-995S kdr mutation, cannot be seen in the expression data at all. This module shows how RNA-Seq-Pop extracts population genetic signals from RNA-seq, what can go wrong, and what the signals say about the Busia colony and your candidate genes.

Learning objectives¶

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

• Explain how variants are called from pooled RNA-seq, and the main caveats: uneven coverage, allele-specific expression and pooling.

• Read allele frequencies of known resistance mutations (e.g. Vgsc-995S) from RNA-Seq-Pop output.

• Interpret Fst and PBS per gene and in genomic windows, and genetic diversity statistics.

• Use karyotype and ancestry-informative marker results to check and interpret a study.

• Summarise the sequence evidence for your candidate genes.

Setup¶

Calling variants from RNA-seq¶

RNA-Seq-Pop’s variant branch runs in parallel with the expression branch:

1. HISAT2 aligns reads to the genome. It is splice-aware, so reads that span an intron are split across the two exons. Duplicate reads are marked with samblaster.

2. Freebayes calls SNPs jointly in all 16 samples, treating each group as a population, with --ploidy 10 and --pooled-discrete: each library is a pool of 5 diploid females, so it is modelled as one “individual” with 10 genome copies.

3. Indels are removed (they are unreliable from RNA-seq), and snpEff annotates each SNP (synonymous, missense, intronic ...).

4. For each analysis, SNPs are filtered on quality (QUAL ≥ 30) and missingness. In our config missingness: 1 keeps only SNPs with genotype calls in every sample.

A genotype call like 0/0/0/0/0/0/0/1/1/1 means the pool is estimated to carry 3 alternative alleles out of 10, an allele frequency of 0.3.

Caveats¶

RNA-seq is not whole-genome sequencing. Keep these in mind whenever you interpret variant results:

• Only expressed regions are covered. Depth is proportional to expression, so a highly expressed gene may have 1,000× coverage and a lowly expressed gene almost none. Genes not expressed in adult females (or in whole bodies) contribute no SNPs. Most intergenic and intronic sequence is missing.

• Allele-specific expression (ASE). If one allele is transcribed more than the other (for example because a cis-regulatory variant is linked to it), its frequency in the RNA is higher than in the DNA. For a gene overexpressed through a linked promoter variant or duplication, this bias is exactly where you are looking. Read-based allele frequencies are “expressed allele frequencies”.

• Pooling. Pools of 5 give at most 10 chromosomes per library, and individuals may not contribute equal amounts of RNA. You cannot recover individual genotypes or haplotypes.

• Mapping artefacts. Reads near splice junctions, from paralogous genes (e.g. P450 clusters) or from RNA editing sites can create false SNPs.

• Replicates are not independent populations. In a colony, replicates are drawn from the same population, so between-replicate variation reflects sampling of mosquitoes, not populations.

Despite these, allele frequencies from well-covered RNA-seq correlate well with DNA-based estimates [1]. The paper recommends 25–60 million reads per library and ≥ 10 individuals per pool for variant work.

First, where are the SNPs? RNA-Seq-Pop compares the proportion of SNPs in each genomic feature with the proportion of the reference genome made up by that feature:

Variants of interest¶

For known mutations we do not need freebayes at all. RNA-Seq-Pop runs samtools mpileup at each position listed in resources/exampleMutations.tsv and counts the reads carrying each base in each BAM. The allele balance (proportion of reads with the mutant base) estimates the expressed allele frequency. The list includes the kdr mutations in Vgsc (the voltage-gated sodium channel, target of pyrethroids and DDT), Cyp4j5-L43F, Cyp6p4-I236M, Gste2-I114T and L119V, Ace1-G280S and Rdl-A296G/S.

Vgsc mutations are named by their position in the Musca domestica protein in older papers; in An. gambiae numbering, the East African kdr mutation L1014S is L995S and the West African L1014F is L995F.

In the paper [2], Vgsc-995S rose from 25% in the parental line to 100% in the selected survivors, and Cyp4j5-L43F from 45% to 95%. Gste2-L119V rose from 11% to 52%, while Cyp6p4-I236M did not change (39% vs 38%). Gste2-I114T was absent from Busia but at 92% in Kisumu, a surprise for a “susceptible” strain. No secondary kdr mutations were found alongside 995S.

Genetic differentiation: Fst¶

Fst measures how much of the genetic variation lies between populations rather than within them. RNA-Seq-Pop uses Hudson’s Fst [3], computed from allele counts with scikit-allel:

where $H_w$ is the mean diversity within the two populations and $H_b$ the diversity between them. Fst is 0 if the two groups have identical allele frequencies and 1 if they are fixed for different alleles.

Selection on a locus raises Fst there, because the selected allele (and the haplotype around it) changes frequency more than the genomic background. But genetic drift also raises Fst, especially in small lab colonies, so the genome-wide background matters: we look for genes and regions that stand out.

Fst per gene¶

results/variantAnalysis/selection/FstPerGene.tsv has, for each gene with at least 2 SNPs, the average Fst for each contrast ({contrast}_zFst, despite the name it is not a z-score), its standard error, the number of SNPs and the gene midpoint.

Fst in windows¶

Per-gene Fst is noisy for genes with few SNPs. A genome-wide selection scan averages Fst over windows of consecutive SNPs instead. RNA-Seq-Pop writes windows of 1,000, 2,000 and 5,000 SNPs to results/variantAnalysis/selection/fst/{window}/{contrast}.Fst.{contig}.tsv (columns midpoint, Fst). Because SNPs from RNA-seq are clustered in genes, a 1,000-SNP window can span anything from a few kb to several Mb.

The paper found high Fst between the parental and selected lines across the genome (a small colony under strong selection drifts fast), but especially on 2L, with peaks around Vgsc and across the 2La inversion.

The population branch statistic (PBS)¶

Fst between two groups cannot tell you which group changed. PBS [4] uses three groups. Each pairwise Fst is turned into a branch length $T = -\log(1 - F_{ST})$, and the branch leading to population A is

With BusiaSelected_BusiaParental_Kisumu, PBS is the amount of change specific to the selected line since it split from the parental line, with Kisumu as the outgroup. RNA-Seq-Pop uses scikit-allel’s normalised PBS, which reduces the influence of loci where all branches are long. Per-gene values are in results/variantAnalysis/selection/PbsPerGene.tsv, windows in results/variantAnalysis/selection/pbs/.

Worked example: the Cyp6aa/Cyp6p cluster¶

In Module 2.1 we looked at the expression of CYP6AA1 (AGAP002862). Now we look at the sequence around it. Cyp6aa1 sits in a cluster of P450 genes on 2R (~28.45–28.55 Mb) that includes Cyp6p3 and Cyp6p4. In Uganda and Kenya a triple-mutant haplotype (a Cyp6aa1 duplication, the Cyp6p4-I236M missense mutation and a ZZB transposable element insertion) has risen rapidly and is associated with pyrethroid resistance [9]. If selection in the colony acted on this haplotype, we would expect raised Fst and PBS in the cluster and a change in Cyp6p4-I236M frequency.

Recall from the variants-of-interest heatmap that the paper found no change in Cyp6p4-I236M between the parental and selected lines (39% vs 38%), and that Cyp6p4 was only at the 89th percentile of per-gene Fst. Let’s check the whole cluster.

Duplications and RNA-seq. A duplicated gene has more copies to transcribe and more genomic template for reads to map to. In RNA-seq this inflates the gene’s expression estimate, and in the HISAT2 alignments it shows up as elevated coverage over the duplicated region compared with its neighbours. But RNA-seq coverage also depends on how strongly each gene is expressed, so from RNA-seq alone you cannot separate “more copies” from “more transcription per copy”. Duplications are best detected from whole-genome coverage, which is what you will do for wild Busia mosquitoes on Day 3.

If you ran RNA-Seq-Pop in your Codespace, you can look at the coverage yourself (in the Codespace terminal, from /workspaces/rna-seq-pop):

# mean depth per sample over Cyp6aa1 (2R:28,480,576-28,482,637) and its neighbour Cyp6p4
for bam in results/alignments/*.hisat2.bam; do
  echo -n "$(basename $bam .hisat2.bam) "
  samtools depth -a -r 2R:28480576-28482637 $bam | awk '{s+=$3} END {printf "Cyp6aa1 %.0f  ", s/NR}'
  samtools depth -a -r 2R:28497000-28499000 $bam | awk '{s+=$3} END {printf "Cyp6p4 %.0f\n", s/NR}'
done

(samtools is in RNA-Seq-Pop’s variants conda environment under .snakemake/conda/.)

Genetic diversity¶

Nucleotide diversity (π) is the average number of differences per site between two randomly chosen sequences. Watterson’s θ is based on the number of segregating sites. Both fall after a population bottleneck (e.g. colonisation) and around a selective sweep. RNA-Seq-Pop averages them over 20 kb windows (results/variantAnalysis/diversity/SequenceDiversity.tsv and WattersonsTheta.tsv), and also gives per-gene π and Tajima’s D.

Because SNPs are called only in expressed regions and windows with no data count as zero, the absolute values are not comparable with WGS estimates. Compare groups with each other.

In the paper, the parental Busia line had the highest diversity, the selected line less, and Kisumu the least (π ≈ 1.0 × 10⁻³, 0.7 × 10⁻³ and 0.6 × 10⁻³). Note that the paper down-sampled the Busia groups to four replicates to match Kisumu, whereas our run uses all six, which can shift the estimates slightly.

Karyotyping chromosomal inversions¶

An. gambiae carries large polymorphic inversions on chromosome 2. 2La (~21 Mb) is associated with aridity tolerance, Plasmodium infection and insecticide resistance. Inside an inversion, recombination between the inverted and standard arrangements is suppressed, so some SNPs are almost perfectly associated with each arrangement. RNA-Seq-Pop uses these tag SNPs (from compkaryo [5]) and reports, for each sample, the mean number of inverted alleles at the tag SNPs it could genotype (_mean, out of the ploidy of 10), the karyotype frequency (_freq) and the number of tag SNPs used.

The paper found 2La at 33% in the parental line and 86% in the selected survivors, supporting an earlier association between 2La and pyrethroid resistance in Uganda. 2Rb differed between Kisumu and both Busia lines.

This has a direct consequence for Fst: the whole 2La region (20.5–42.2 Mb on 2L) changed frequency as one block, so every gene inside it shows raised Fst. Cyp4j5 sits inside 2La. Is Cyp4j5-L43F itself selected, or is it hitchhiking on the inversion?

Ancestry-informative markers¶

The An. gambiae complex contains closely related species that can hybridise. Ancestry-informative markers (AIMs) are SNPs with (nearly) fixed differences between species. RNA-Seq-Pop counts, for each group and chromosome arm, the fraction of alleles at AIMs that match An. gambiae or An. coluzzii.

The paper found that the Busia colony is almost entirely An. gambiae, but Kisumu shows a large amount of putative An. coluzzii ancestry on the autosomes and none on X. Kisumu was colonised in 1975 from an area without An. coluzzii, so the most likely explanation is contamination by hybridisation in an insectary. The X chromosome resists introgression, which fits this. This is a good example of RNA-seq variants acting as a quality check on the experiment itself; it is also a reason to be cautious about Kisumu as a “reference” in the Kisumu_BusiaParental contrast.

Your candidate genes¶

Collect the sequence evidence for your genes: number of SNPs, Fst in each contrast (and its percentile), PBS, and diversity and Tajima’s D in each group. Is the gene inside 2La or 2Rb?

If a gene is missing from the table, it had fewer than 2 SNPs passing the filters: usually because it is lowly expressed (little coverage) or not variable. That in itself is worth noting in your dossier.

Exercise 1¶

Which genes have the highest PBS on the BusiaSelected branch? Take the top 20 genes with at least 10 SNPs and look them up: are any known resistance genes, and how many are inside 2La?

Exercise 2¶

Compare Vgsc (AGAP004707) in the two kinds of analysis: its differential expression result (Module 2.1: results/genediff/BusiaParental_BusiaSelected.csv) and its Fst percentile here. Now do the same for SAP2 (AGAP008052). What does each gene’s pattern suggest about the mechanism?

Exercise 3¶

The Cyp4j5-L43F frequency in the RNA of a pool may not equal its frequency in the DNA. Give two reasons, and say which direction each would bias the estimate.

Exercise 4¶

Fill in the “Genomic variation” part of your dossier for each candidate gene: number of SNPs, Fst percentile in both contrasts, PBS percentile, inversion, and any variant of interest in or near the gene. On Day 3 you will check whether the same genes carry CNVs or missense SNPs in wild Ugandan and Kenyan mosquitoes.

Quiz¶

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

Summary¶

• RNA-seq reads give SNPs in expressed genes. RNA-Seq-Pop aligns with HISAT2 and calls pooled variants with freebayes (ploidy 10 for pools of 5 diploids).

• Expressed allele frequencies can be biased by allele-specific expression, uneven coverage and pooling, but are useful when coverage is high.

• In Busia, selection raised Vgsc-995S to 100% and Cyp4j5-L43F to 95%, and the 2La inversion from 33% to 86%. Fst was highest on 2L.

• Fst measures differentiation between two groups; PBS isolates change on one branch using an outgroup. Inversions and drift in small colonies raise both, so compare against the background.

• Karyotypes and AIMs are checks on the experiment: Kisumu carries An. coluzzii ancestry on its autosomes.

Well done!¶

You have now combined expression and sequence evidence for your candidate genes. In Module 2.3 you check whether they are overexpressed in other studies with AnoExpress.

References¶

1. Jehl F et al. (2021). RNA-Seq data for reliable SNP detection and genotype calling: interest for coding variant characterization and cis-regulation analysis by allele-specific expression in livestock species. Frontiers in Genetics 12:655707. Jehl et al. (2021)

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. Bhatia G, Patterson N, Sankararaman S, Price AL (2013). Estimating and interpreting FST: the impact of rare variants. Genome Research 23:1514–1521. https://doi.org/10.1101/gr.154831.113

4. Yi X et al. (2010). Sequencing of 50 human exomes reveals adaptation to high altitude. Science 329:75–78. Yi et al. (2010)

5. Love RR, Redmond SN, Pombi M, Caputo B, Petrarca V, della Torre A, Besansky NJ (2019). In silico karyotyping of chromosomally polymorphic malaria mosquitoes in the Anopheles gambiae complex. G3 9:3249–3262. Love et al. (2019)

6. Clarkson CS et al. (2021). The genetic architecture of target-site resistance to pyrethroid insecticides in the African malaria vectors Anopheles gambiae and Anopheles coluzzii. Molecular Ecology 30:5303–5317. Clarkson et al. (2021)

7. Weetman D et al. (2018). Candidate-gene based GWAS identifies reproducible DNA markers for metabolic pyrethroid resistance from standing genetic variation in East African Anopheles gambiae. Scientific Reports 8:2920. Weetman et al. (2018)

8. Miles A et al. scikit-allel: explore and analyse genetic variation. Alistair Miles et al. (2024)

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)