4.3 Analysing amplicon sequencing data with AmpSeeker¶

Theme: Tools & technology | Analysis

In Module 4.2 we saw how a multiplex amplicon panel is designed. Here we analyse data from one: Ag-vampIR, a panel of 82 amplicons that genotypes 90 SNPs in An. gambiae s.l., including insecticide resistance mutations, SNPs that tag known selective sweeps, and ancestry informative markers (AIMs) [1]. We analyse it with AmpSeeker, a Snakemake workflow that goes from reads to genotypes, quality control and an interactive results book [1].

The data come from a colony established in 2023 from Siaya, western Kenya, about 50 km from Busia. Mosquitoes were exposed to deltamethrin nets (PermaNet 2.0) in extended WHO cone tests and scored dead or alive. We also run the VK7 colony (An. coluzzii, Burkina Faso), sequenced on the same MiSeq run.

This is where the course storyline closes. On Days 1–3 we found genes overexpressed in the Busia colony and selective sweeps in the Busia 2016 atlas cohort at the Cyp6aa/p cluster (2R), 2L ~34 Mb, the Gste cluster (3R) and Cyp9k1 (X). Ag-vampIR has markers at all four loci. So: do the loci we found with RNA-seq and selection scans predict which mosquitoes survive? We follow our worked example, CYP6AA1, which sits in the 2R sweep, to the end.

Learning objectives¶

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

• Explain how amplicon sequencing data are processed: demultiplexing, trimming, alignment to a reference genome, and variant calling at targets.

• Judge the quality of an amplicon run from coverage, sample and target failure, heterozygosity and negative controls.

• Configure and run AmpSeeker from FASTQ files in a GitHub Codespace, and find your way around its results book.

• Load AmpSeeker genotypes in Python and test markers for association with a phenotype using a binomial GLM with false discovery rate control.

• Judge whether the Busia RNA-seq and selection results are supported by survival data from Siaya.

Setup¶

Colab already has pandas, plotly and statsmodels. We add scikit-allel to read VCF files.

From reads to genotypes: how amplicon data are analysed¶

Whole-genome sequencing spreads reads over the whole genome at ~30× coverage. Amplicon sequencing does the opposite: a multiplex PCR copies a few dozen short regions (~200 bp each), and hundreds of mosquitoes are sequenced together on one run. Each target ends up with hundreds or thousands of reads, at a cost of around US$10 per sample for Ag-vampIR [1]. Because every read starts and ends at a primer, the analysis has its own steps and its own failure modes.

1. Demultiplexing. In the second PCR each sample gets a unique pair of index sequences (Ag-vampIR uses 96 i7 × 16 i5 indexes, so up to 1,536 samples per run). The sequencer reads the indexes separately, and demultiplexing (bcl-convert) sorts reads into one pair of FASTQ files per sample using a SampleSheet.csv. A small fraction of reads always get the wrong sample, through index misassignment on the flow cell or cross-contamination in the lab. AmpSeeker can start from the raw MiSeq output folder (from-bcl: True) or, as today, from FASTQs that are already demultiplexed (from-bcl: False); the SRA stores demultiplexed reads.

2. Trimming. fastp removes adapters and low-quality ends. With ~200 bp amplicons and 2 × 150 bp reads, the two reads of a pair overlap, and short products (primer dimers) read into the adapter.

3. Alignment to the genome. Reads are aligned with bwa mem to the whole AgamP4 reference, not just to the amplicon sequences. That way positions are in genome coordinates, the same as Ag3 and the selection atlas, and reads from off-target products map elsewhere instead of being forced onto a target.

4. Variant calling. bcftools mpileup and bcftools call genotype each sample, in two ways:

• targets: only the panel’s SNPs, listed in config/ag-vampir.bed. A genotype is reported even when the sample carries the reference allele, so “not called” means “no data”.

• amplicons: every position covered by reads, which finds SNPs the panel was not designed for.

Calls are annotated with snpEff (gene, amino acid change). Be careful with any SNP under a primer: the primer sequence is copied into every read, so the read shows the primer’s base, not the mosquito’s. A SNP in a primer site can also stop one allele from amplifying (allele dropout), which makes heterozygotes look homozygous.

5. Quality control. Amplicons amplify with different efficiency, so coverage is never even. Look for:

• Sample failure: low total coverage (poor DNA, failed PCR) and many missing calls.

• Target failure: an amplicon with low coverage in every sample (a poor primer, primer dimers, or a SNP under a primer in this population).

• Contamination: a sample containing DNA from two mosquitoes has too many heterozygous calls. AmpSeeker removes samples with outlying heterozygosity or outlying PCA positions within their cohort.

• Negative controls (no-template wells): they should have almost no reads. Reads in a negative control measure background contamination and index misassignment for the whole plate. AmpSeeker drops samples whose names contain Negative or dil (dilution replicates) from the analysis.

A first QC check, before running anything¶

The number of reads per sample is public in the ENA metadata. Let us compare the negative controls with real samples.

Exercise 1¶

How many reads do the three negative controls have, compared with the median sample in the same colony? What might explain reads in a negative control, and what would you do if a negative control had as many reads as a real sample?

Practical: run AmpSeeker in your codespace¶

AmpSeeker needs conda and several GB of disk for its environments and the reference genome, so we run it in the same GitHub Codespace as RNA-Seq-Pop (Module 1.4). When the codespace was created, .devcontainer/setup.sh cloned AmpSeeker v0.7.0 to /workspaces/AmpSeeker.

Step 1: look at the workflow¶

cd /workspaces/AmpSeeker
git describe --tags          # v0.7.0
ls config/ workflow/rules/ workflow/notebooks/

• workflow/Snakefile includes rule files for each stage: qc.smk, map-call-illumina.smk, snpeff.smk, analysis.smk, ag-vampir.smk (species ID and kdr origins) and jupyter-book.smk.

• Most analyses are Jupyter notebooks in workflow/notebooks/, run with papermill, which passes parameters (file paths, cohort columns) into the notebook before running it. The run notebooks are then built into a results book with Jupyter Book.

• config/ag-vampir.bed lists the panel’s targets: contig, start, end, amplicon, marker name, and the reference and alternative allele of the marker.

Step 2: get the data, the reference and the config¶

data/prepare_siaya.sh (in the course repo) downloads the FASTQs from ENA and names them the way AmpSeeker expects when from-bcl: False: resources/reads/{sample_id}_1.fastq.gz and _2.fastq.gz. It also copies the course config and sample sheet into config/, and downloads the AgamP4 genome and annotation from Ensembl.

The full dataset is 360 samples (~0.8 GB) and makes ~4,000 Snakemake jobs. In class we run a subset of 12 samples per cohort (51 samples) so it finishes during the session; the Python sections below use the precomputed run on all samples.

cd /workspaces/anopheles-omics-training
bash data/prepare_siaya.sh -n 12         # drop "-n 12" to download all 360 samples

Look at what it made:

cd /workspaces/AmpSeeker
ls resources/reads | head
column -t config/siaya-metadata.tsv | head
cat config/siaya.yaml

Things to notice in config/siaya.yaml:

• from-bcl: False and no fq1/fq2 columns in the metadata, so AmpSeeker finds the reads by name.

• cohort-columns: [location, taxon, cohort]: results are summarised by each of these. taxon starts as unassigned and is filled in by the species-ID step.

• custom-snpeffdb: True builds a snpEff database from our reference instead of downloading one.

• sample-total-reads-threshold: samples below this total depth over all targets are removed.

Step 3: dry run and DAG¶

snakemake --cores 4 --use-conda --configfile config/siaya.yaml -n
snakemake --configfile config/siaya.yaml --rulegraph | dot -Tsvg > rulegraph.svg

The dry run lists every job and a count per rule. Open rulegraph.svg in the file explorer to see how the rules connect. (Use --dag instead of --rulegraph to draw every job; with hundreds of samples that picture is not readable.)

Step 4: run¶

nohup snakemake --cores 4 --use-conda --configfile config/siaya.yaml > ampseeker-run.log 2>&1 &
tail -f ampseeker-run.log          # Ctrl+C stops watching, not the run

The first run builds AmpSeeker’s conda environments (~15 minutes), then indexes the genome with bwa (~10 minutes). The per-sample jobs are quick. If the codespace stops, restart it and run the same command: Snakemake continues from where it stopped.

Step 5: open the results book¶

python -m http.server 8000 -d results/ampseeker-results/_build/html

Open port 8000 from the Ports tab. While your run finishes, use the book from the full dataset (Siaya, VK7, Ghana and The Gambia) published with the paper: https://sanjaynagi.github.io/agvampir002-results/

Exercise 2¶

From the dry run: which rules run once per sample, and which run once for the whole dataset? Which rules would be added if we had started from the MiSeq output folder with from-bcl: True?

Exploring the results book¶

The book has one page per notebook. Work through it in this order, using your own run or the published book.

Exercise 3¶

Use the results book to answer:

1. What taxon are the Siaya mosquitoes assigned, and why is that surprising for western Kenya? (Hint: compare the AIMs on chromosome 2 with those on chromosome 3.)

2. Which kdr alleles and origins are present in Siaya, and in VK7?

3. Is the Cyp6p4-I236M allele more common in the Siaya_alive or the Siaya_dead cohort?

Loading AmpSeeker’s outputs in Python¶

The results book is for browsing. To test hypotheses we load AmpSeeker’s output files ourselves. The ones we need, with their paths inside the AmpSeeker folder:

We download a bundle with these files from a full run, unless the folder already exists.

Now the coverage matrix. We plot log10 depth, amplicons as rows and samples (sorted by cohort) as columns. Horizontal stripes are weak amplicons; vertical stripes are weak samples.

Reading genotypes at the panel’s targets¶

The targets VCF has one record per target position. For each marker we count copies of the marker allele given in the BED file (0, 1 or 2; missing if the sample was not called). Note that the BED gives the allele the panel was designed to track: for a few tagging SNPs that is the reference base, and at Vgsc_V402L two alleles are tracked (C,T), so we make one row for each.

The Busia sweep markers¶

We label the markers that sit in the four Busia sweeps, and kdr for comparison.

Before testing anything, look at the plot. Several markers have a frequency of 0 (or 1) in Siaya: there is nothing to test there. We come back to this below.

Genotype:phenotype association¶

The model. For each marker we fit a binomial generalised linear model (logistic regression):

where $x$ is the number of copies (0, 1, 2) of the marker allele. This is an additive model: each copy multiplies the odds of surviving by the same amount, the odds ratio $e^{\beta_1}$. An odds ratio above 1 means the allele is associated with survival, below 1 with death. This is the test used in the Ag-vampIR paper [1].

Filters. We only test markers with at most 20% missing calls and a minor allele frequency of at least 2% in Siaya. A marker that is (almost) fixed cannot explain differences between mosquitoes.

Multiple testing. With dozens of tests, some will have $p < 0.05$ by chance. We control the false discovery rate (FDR) with the Benjamini–Hochberg method: among markers with FDR < 0.05, we expect fewer than 5% to be false positives.

Linkage. Markers in the same sweep sit on the same haplotype, so they give correlated results. Ten significant SNPs at the Cyp6 cluster are one signal, not ten.

Design. All the mosquitoes come from one colony and were tested in the same assay, which limits confounding by population structure. But a colony is not the wild population, and the result of a 40-minute cone test is not the same as surviving a night under a net.

Exercise 4¶

The single-marker tests find signals at the Cyp6aa/p cluster, at 2L ~34 Mb and at kdr. Are these independent, or could one locus be explaining another (for example through linkage, as 34 Mb and kdr are both on 2L)? Fit one GLM with Cyp6p4_I236M, 34mb_tag2 and Vgsc_L995F together and compare the odds ratios with the single-marker ones.

Exercise 5¶

Which Ag-vampIR target is closest to each of your candidate genes, and is it associated with survival in Siaya? The code below reads gene positions from the AgamP4 annotation.

Closing the loop: CYP6AA1 and the Busia sweeps¶

Our worked example, CYP6AA1. Through the course we have followed CYP6AA1 (AGAP002862):

• Day 3: it sits in the Busia 2016 sweep on 2R, and in Uganda and Kenya it is duplicated on the “triple mutant” haplotype: a Cyp6aa1 duplication, the Cyp6p4-I236M substitution and a ZZB transposable element insertion, which Njoroge et al. (2022) found spreading rapidly and linked to high-level metabolic resistance [2].

• Day 2: it is a known resistance gene, but it is not significantly overexpressed in the AnoExpress Busia survivors vs Kisumu contrast (log2FC −0.29), so the expression evidence comes from other studies, not from that comparison.

• Today: amplicons cannot see the duplication itself (a ~200 bp amplicon does not measure copy number). But Cyp6p4_I236M and Cyp6_tag8 sit on the triple-mutant haplotype, and in Siaya they predict survival after deltamethrin exposure (OR ≈ 3.4 per copy, ~63% frequency) [1].

So the haplotype that carries the CYP6AA1 duplication, under selection in Busia, predicts which mosquitoes 50 km away survive a deltamethrin net. That is strong field evidence for the locus. It does not tell us which of the three mutations (or CYP6AA1 versus its neighbours CYP6P3 and CYP6P4) causes the resistance: they are inherited together, so the association is with the whole haplotype. Separating them needs functional work, such as expressing the enzymes or genotyping populations in which the mutations occur apart. And it is cheap to monitor: I236M is a single SNP on the panel, and Njoroge et al. proposed the triple mutant as a molecular diagnostic [2].

The other Busia sweeps, as found in the published analysis [1] (check them against your output):

For Gste and Cyp9k1 we cannot say. Their tagging SNPs were chosen from sweeps in the first phase of Ag1000G, mostly in West and Central African samples (the Cyp9k1 tags track the sweep in Burkina Faso An. coluzzii, which is why VK7 is fixed for them). A sweep at the same gene in East Africa can be on a different haplotype, which these SNPs do not tag. No variation means no test, not no effect.

Two ways forward, both using tools from this course:

1. Rerun the association on all SNPs in the amplicons (Exercise 6). The paper found extra signals at Cyp4j5 and the TEP locus this way [1].

2. Design new tagging SNPs from the East African swept haplotypes (haplotype clustering, Module 3.3), and add them to the panel (AnoPrimer and multiply, Modules 4.1–4.2).

Other caveats before you write your dossier verdict: Siaya is a laboratory colony established in 2023, not a wild population, and it is not Busia; a 40-minute cone test is not a night under a net; an association is evidence, not proof of cause; and a marker is only as good as its linkage to the causal variant.

Which of your own candidate genes sit in a Busia sweep, and could the panel test them?

Exercise 6 (extension)¶

The targets VCF only covers the panel’s 90 SNPs. The amplicons VCF has every variable position on the amplicons. Repeat the association on all SNPs in the amplicons (split multi-allelic sites into one row per alternative allele). How many SNPs pass the filters, and which loci reach FDR < 0.05?

Quiz¶

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

Summary¶

• Amplicon sequencing genotypes a few dozen targets very deeply in hundreds of samples. The analysis is demultiplex → trim → align to the genome → call at the targets (and across the amplicons) → QC.

• QC is about coverage: failed samples, failed amplicons, contamination (excess heterozygosity) and negative controls, which in these data are not empty.

• AmpSeeker runs all of this with Snakemake from FASTQs (from-bcl: False, reads named {sample_id}_1.fastq.gz), and builds a results book with species ID, kdr origins, allele frequencies and population structure.

• A per-marker binomial GLM with FDR control links genotypes to survival. In Siaya, markers at the Cyp6aa/p cluster, the 2L 34 Mb sweep and kdr predict survival after deltamethrin exposure.

• Cyp6p4_I236M and Cyp6_tag8 tag the triple-mutant haplotype that carries the CYP6AA1 duplication, so the locus of our worked example predicts survival in the field. Which of its linked mutations causes resistance is still open.

• Markers must vary in the population you test. The Gste and Cyp9k1 tags are monomorphic in Siaya, so the panel cannot test those Busia sweeps; new East African tags are needed.

Well done!¶

You have followed the Busia loci from differential expression, through selection scans and assay design, to survival in the field. On Day 5 you bring this together in your gene dossier.

References¶

1. Nagi SC, Lucas ER, Ashraf F, et al. (2025). Targeted genomic surveillance of insecticide resistance in African malaria vectors. bioRxiv. Nagi et al. (2025)

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

3. Mölder F, Jablonski KP, Letcher B, et al. (2021). Sustainable data analysis with Snakemake. F1000Research 10:33. Mölder et al. (2021)

4. Danecek P, Bonfield JK, Liddle J, et al. (2021). Twelve years of SAMtools and BCFtools. GigaScience 10:giab008. Danecek et al. (2021)

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

AmpSeeker documentation: https://sanjaynagi.github.io/AmpSeeker/. Ag-vampIR results book: https://sanjaynagi.github.io/agvampir002-results/.