3.3 (Optional) Following up the Busia Cyp6 sweep¶

Theme: Analysis

The selection atlas told us (Module 3.2) that Busia 2016 has a strong sweep centred on Cyp6aa1, and in Module 3.1 we found the Cyp6aap_Dup1a duplication and CYP6P4 I236M at high frequency there. These are two of the three mutations of the East African “triple mutant” [1]. But a signal is a summary. To design an assay on Day 4 we need to know which haplotypes were swept and which SNPs mark them.

In this module we zoom in on the sweep: we recompute H12 around the locus, cluster haplotypes to find the swept haplotype(s), test whether they carry I236M and the duplication, and pick SNPs that tag them. These SNPs are candidate targets for the AnoPrimer assays in Module 4.1. The same code works for the Cyp9k1 sweep on X: change one setting.

This module builds on haplotype clustering from the PAMCA course (W7-M3) and the selection atlas follow-up guide.

Learning objectives¶

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

Recompute H12 in a small region around an atlas signal and compare it with the atlas intervals.
Use haplotype clustering to identify swept haplotypes and follow them across cohorts and years.
Test whether a swept haplotype carries known resistance mutations (SNPs and CNVs).
Find SNPs that tag a swept haplotype and judge whether they would make good assay targets.

Setup¶

This module is Colab-friendly because it only loads haplotypes for small regions (a few Mb for H12, 30 kb for clustering). It takes a few minutes to run.

%pip install -q --no-warn-conflicts malariagen_data==15.9.0

import sys
import os
os.environ["MGEN_SHOW_PROGRESS"] = "0"
import warnings
warnings.filterwarnings("ignore")

from pathlib import Path
import allel
import numpy as np
import pandas as pd
import plotly.express as px
import plotly.graph_objects as go
import plotly.io as pio
from scipy.cluster.hierarchy import fcluster, linkage
from scipy.spatial.distance import pdist
import malariagen_data
pio.templates.default = "simple_white"  # course style for every Plotly figure

# Colab needs its own renderer; elsewhere, save figures as Plotly JSON so the book can show them
pio.renderers.default = "colab" if "google.colab" in sys.modules else "plotly_mimetype+notebook_connected"
pd.set_option("display.width", 200)

ag3 = malariagen_data.Ag3()

# TODO(data): the atlas signal table is not yet published at a public URL (see Module 3.2 and data/README.md).
SIGNALS_URL = "https://raw.githubusercontent.com/anopheles-genomic-surveillance/selection-atlas/main/h12-signal-detection-all.csv"

# Ag-vampIR amplicon panel targets (pinned AmpSeeker commit), used on Day 4
AGVAMPIR_BED = "https://raw.githubusercontent.com/sanjaynagi/AmpSeeker/864ba69/config/ag-vampir.bed"

# Replace with the genes you were assigned on Day 2 (see data/candidate-genes.tsv)
CANDIDATE_GENES = ["AGAP002862", "AGAP009193"]  # CYP6AA1, GSTE4

The worked example follows the sweep at the first default gene, Cyp6aa1. Set LOCUS = "Cyp9k1" to repeat everything for the X-chromosome sweep instead (Exercise 1).

LOCUS = "Cyp6"  # or "Cyp9k1"

LOCI = {
    "Cyp6": dict(
        atlas_contig="2RL",
        gene="AGAP002862",                        # CYP6AA1
        h12_region="2R:26,500,000-30,500,000",    # a few Mb around the signal
        cluster_region="2R:28,480,000-28,510,000",  # 30 kb over Cyp6aa1-Cyp6p2
        snp_transcript="AGAP002867-RA",           # CYP6P4, shows I236M under the dendrogram
        cnv_contig="2R",
        cnv_allele="Cyp6aap_Dup1a",
        marker=("CYP6P4 I236M", "2R:28497967-28497967", "C"),  # (name, position, allele)
    ),
    "Cyp9k1": dict(
        atlas_contig="X",
        gene="AGAP000818",                        # CYP9K1
        h12_region="X:13,500,000-17,500,000",
        cluster_region="X:15,230,000-15,260,000",
        snp_transcript=None,
        cnv_contig="X",
        cnv_allele="Cyp9k1_Dup8",
        marker=None,
    ),
}
cfg = LOCI[LOCUS]

SAMPLE_SETS = [
    "AG1000G-UG",                     # Tororo and Kanungu, 2012
    "1178-VO-UG-LAWNICZAK-VMF00025",  # Busia 2016-17, Mayuge 2017
    "1288-VO-UG-DONNELLY-VMF00168",   # Namayingo 2017-19
    "1288-VO-UG-DONNELLY-VMF00219",   # Busia 2017 and 2019
    "1274-VO-KE-KAMAU-VMF00246",      # Teso North (Busia County, Kenya) 2019
]
BUSIA_2016 = "cohort_admin2_quarter == 'UG-E_Busia_gamb_2016_Q2'"  # the atlas cohort
CLUSTER_COHORTS = ["UG-E_Busia_gamb_2016", "UG-E_Busia_gamb_2019", "KE-04_Teso-North_gamb_2019"]
CLUSTER_QUERY = f"taxon == 'gambiae' and cohort_admin2_year in {CLUSTER_COHORTS}"
DISTRICTS = ["Busia", "Namayingo", "Mayuge", "Tororo", "Teso North", "Kanungu"]
REGION_QUERY = f"taxon == 'gambiae' and admin2_name in {DISTRICTS}"

Step 1: the atlas signal¶

We take the Busia 2016 signal at the locus from the atlas table. Our loci are on 2R and X, where atlas and arm coordinates are the same (see Module 3.2 for 2L and 3L).

signals = pd.read_csv(SIGNALS_URL)
gene = ag3.genome_features(attributes=["ID", "Name"]).query("ID == @cfg['gene']").iloc[0]
sig = (
    signals.query("cohort_id == 'UG-E_Busia_gamb_2016_Q2' and contig == @cfg['atlas_contig']")
    .query("span2_pstart <= @gene.end and span2_pstop >= @gene.start")
    .iloc[0]
)
sig[["cohort_id", "contig", "pcenter", "stat_max", "delta_i",
     "focus_pstart", "focus_pstop", "span1_pstart", "span1_pstop"]]

Step 2: an H12 zoom¶

The atlas computed H12 over whole chromosomes. We can reproduce the local picture from haplotypes in a few Mb around the signal, using the same window size the atlas calibrated for Busia 2016 (2,000 SNPs, shown on the cohort page). H12 is the sum of the squared frequencies of the two most common haplotypes plus the others squared (Garud’s H12 [2]); allel.garud_h() returns it as its second value. For whole chromosomes use ag3.h12_gwss() as in W6-M3.

H12_WINDOW = 2000  # atlas-calibrated H12 window size for UG-E_Busia_gamb_2016_Q2


def h12_region(region, sample_query, window_size=H12_WINDOW, **kwargs):
    # Windowed H12 over a region, as in ag3.h12_gwss but for a small region only.
    ds = ag3.haplotypes(region=region, analysis="gamb_colu", sample_sets=SAMPLE_SETS,
                        sample_query=sample_query, **kwargs)
    ht = allel.GenotypeArray(ds["call_genotype"].values).to_haplotypes()
    pos = ds["variant_position"].values
    h12 = allel.moving_statistic(ht, statistic=lambda h: allel.garud_h(h)[1], size=window_size)
    x = allel.moving_statistic(pos, statistic=np.mean, size=window_size)
    return pd.DataFrame({"position": x, "h12": h12}), ds.sizes["samples"]


df_h12, n = h12_region(cfg["h12_region"], BUSIA_2016)
print(f"{n} mosquitoes, {len(df_h12)} windows")

24 mosquitoes, 693 windows

def plot_h12(df_h12, sig, gene, title):
    fig = go.Figure()
    for name, alpha in [("span2", 0.08), ("span1", 0.15), ("focus", 0.3)]:
        fig.add_vrect(x0=sig[f"{name}_pstart"] / 1e6, x1=sig[f"{name}_pstop"] / 1e6,
                      fillcolor="firebrick", opacity=alpha, line_width=0,
                      annotation_text=name, annotation_position="top left")
    fig.add_trace(go.Scatter(x=df_h12["position"] / 1e6, y=df_h12["h12"], mode="markers",
                             marker=dict(size=5, color="black"), name="H12"))
    fig.add_vline(x=(gene.start + gene.end) / 2e6, line_dash="dot", annotation_text=gene.Name)
    fig.update_layout(title=title, xaxis_title=f"{gene.contig} position (Mb)", yaxis_title="H12",
                      yaxis_range=[0, 1], width=850, height=380, showlegend=False)
    return fig


plot_h12(df_h12, sig, gene, f"H12 around {gene.Name}, Busia 2016 (window {H12_WINDOW} SNPs)")

For Cyp6, H12 rises to a peak inside the atlas focus interval: the highest window (H12 ≈ 0.80 at ~28.49 Mb) matches the atlas stat_max for this signal. Each window of 2,000 phased SNPs spans about 6 kb here, and H12 stays high over the whole gene cluster, so the scan alone cannot say which gene or variant was selected. For that we need the haplotypes themselves.

Step 3: haplotype clustering¶

We cluster haplotypes in a 30 kb window over the locus, for Busia in 2016 and 2019 and Teso North (Kenya) in 2019. In a sweep, many haplotypes are identical or nearly so and form a large cluster with branch lengths close to zero. Colour shows the cohort.

fig, leaf_data = ag3.plot_haplotype_clustering_advanced(
    region=cfg["cluster_region"],
    analysis="gamb_colu",
    sample_sets=SAMPLE_SETS,
    sample_query=CLUSTER_QUERY,
    snp_transcript=cfg["snp_transcript"],
    snp_filter_min_maf=0.05,
    color="cohort_admin2_year",
    cluster_threshold=5,
    min_cluster_size=5,
    width=1000,
)

The row(s) under the dendrogram show the cluster assignment (cutting the tree at 5 SNP differences) and, for Cyp6, the non-synonymous SNPs in CYP6P4 above 5 % frequency. Hover over the leaves to see which mosquito each haplotype comes from.

To work with the alleles on each haplotype we redo the clustering ourselves: load the haplotypes, compute pairwise distances (number of SNP differences) and cut the tree.

ds_haps = ag3.haplotypes(region=cfg["cluster_region"], analysis="gamb_colu",
                         sample_sets=SAMPLE_SETS, sample_query=CLUSTER_QUERY)
ht = allel.GenotypeArray(ds_haps["call_genotype"].values).to_haplotypes()
pos = ds_haps["variant_position"].values
alleles = ds_haps["variant_allele"].values.astype(str)
sample_ids = ds_haps["sample_id"].values

seg = ht.count_alleles(max_allele=1).is_segregating()
ht_seg = ht.compress(seg, axis=0)
dist = pdist(ht_seg.T, metric="hamming") * ht_seg.shape[0]  # number of SNP differences
clusters = fcluster(linkage(dist, method="complete"), t=5, criterion="distance")

swept_id = np.bincount(clusters).argmax()
is_swept = clusters == swept_id
print(f"{ht.shape[1]} haplotypes, {ht_seg.shape[0]} segregating SNPs; "
      f"largest cluster: {is_swept.sum()} haplotypes ({is_swept.mean():.0%})")

258 haplotypes, 890 segregating SNPs; largest cluster: 215 haplotypes (83%)

df_meta = ag3.sample_metadata(sample_sets=SAMPLE_SETS, sample_query=CLUSTER_QUERY).set_index("sample_id")
hap_cohort = np.repeat(df_meta.loc[sample_ids, "cohort_admin2_year"].values, 2)  # two haplotypes per mosquito
pd.crosstab(hap_cohort, np.where(is_swept, "swept", "other"), normalize="index").round(2)

The swept haplotype is at a similar frequency in Busia 2016, Busia 2019 and Teso North 2019: the sweep crosses the border.

Step 4: what does the swept haplotype carry?¶

First, known SNPs. We load the Ag-vampIR panel targets and, for those in our window, compare the alternate allele frequency on swept and other haplotypes.

panel = pd.read_csv(AGVAMPIR_BED, sep="\t", header=None,
                    names=["contig", "start", "end", "amplicon", "target", "ref", "alt"])
panel["position"] = panel["end"]  # BED is 0-based; end is the 1-based position
contig = cfg["cluster_region"].split(":")[0]
in_window = panel.query("contig == @contig and @pos.min() <= position <= @pos.max()")


def haplotype_allele_frequencies(position, alt):
    # Frequency of `alt` at `position` on swept and other haplotypes (NaN if not a phased site).
    ix = np.flatnonzero(pos == position)
    if len(ix) == 0 or alt not in alleles[ix[0]]:
        return np.nan, np.nan
    carries = ht[ix[0]] == list(alleles[ix[0]]).index(alt)
    return carries[is_swept].mean(), carries[~is_swept].mean()


rows = []
for t in in_window.itertuples():
    f_swept, f_other = haplotype_allele_frequencies(t.position, t.alt)
    rows.append(dict(target=t.target, position=t.position, alt=t.alt,
                     frq_swept=f_swept, frq_other=f_other))
pd.DataFrame(rows).round(2)

For Cyp6, every swept haplotype carries the Cyp6_tag8 (and Cyp6_tag1) allele, but so do some other haplotypes, so these are good but imperfect tags. Most other Ag-vampIR Cyp6 targets are either not variable here or not among the phased sites (NaN); the panel’s tags were designed across many African populations, so only some of them are informative in Busia.

CYP6P4 I236M is NaN because that SNP is not in the phased haplotype data. We can still test it: compare each mosquito’s number of swept haplotypes (0, 1 or 2) with its number of 236M alleles in the unphased SNP calls. If the swept haplotype carries 236M, the two counts should match. (For Cyp9k1 there is no known functional SNP, so this cell does nothing.)

n_swept = pd.Series(is_swept.reshape(-1, 2).sum(axis=1), index=sample_ids, name="n swept haplotypes")

if cfg["marker"] is not None:
    marker_name, marker_region, marker_allele = cfg["marker"]
    ds_snp = ag3.snp_calls(region=marker_region, sample_sets=SAMPLE_SETS, sample_query=CLUSTER_QUERY)
    snp_alleles = list(ds_snp["variant_allele"].values[0].astype(str))
    gt = ds_snp["call_genotype"].values[0]
    n_marker = pd.Series((gt == snp_alleles.index(marker_allele)).sum(axis=1),
                         index=ds_snp["sample_id"].values, name=f"n {marker_name} alleles")
    display(pd.crosstab(n_swept, n_marker.reindex(n_swept.index)))

For Cyp6: all mosquitoes with two swept haplotypes are homozygous for 236M, and most with one are heterozygous. The swept haplotype carries CYP6P4 I236M. (The few mismatches are expected from genotyping or phasing errors in an amplified region.) On Day 4 you will see that Cyp6p4_I236M and Cyp6_tag8 are the Ag-vampIR markers associated with surviving deltamethrin in Siaya.

Next, the duplication. CNVs are not in the haplotype data, so we compare each mosquito’s number of swept haplotypes with whether it carries the duplication allele.

ds_cnv = ag3.cnv_discordant_read_calls(contig=cfg["cnv_contig"], sample_sets=SAMPLE_SETS,
                                        sample_query=CLUSTER_QUERY)
ix = list(ds_cnv["variant_id"].values).index(cfg["cnv_allele"])
dup = pd.Series(ds_cnv["call_genotype"].values[ix], index=ds_cnv["sample_id"].values,
                name=f"{cfg['cnv_allele']} carrier")

df = pd.concat([n_swept, dup], axis=1, join="inner")
df["cohort"] = df_meta.loc[df.index, "cohort_admin2_year"]
pd.crosstab([df["cohort"], df["n swept haplotypes"]], df[dup.name])

Mosquitoes carrying at least one swept haplotype carry the duplication. Because the duplication is carried by nearly everyone around Busia, this table cannot tell us much about mosquitoes without it. Exercise 1 repeats the test for Cyp9k1, where the amplification is at intermediate frequency and the test is much more powerful.

Step 5: SNPs that tag the swept haplotype¶

A good tag SNP has one allele on (almost) every swept haplotype and (almost) never on other haplotypes. We scan all phased SNPs in the window, keep those whose alternate allele is on at least 95 % of swept haplotypes, and rank them by the difference in frequency between swept and other haplotypes.

alt_is_1 = ht_seg == 1  # allele 1 = the alternate allele in the phased data
tags = pd.DataFrame({
    "position": pos[seg],
    "ref": alleles[seg][:, 0],
    "alt": alleles[seg][:, 1],
    "frq_swept": alt_is_1[:, is_swept].mean(axis=1),
    "frq_other": alt_is_1[:, ~is_swept].mean(axis=1),
})
tags["difference"] = tags["frq_swept"] - tags["frq_other"]
tags = tags.query("frq_swept >= 0.95").sort_values("difference", ascending=False)
print(len(tags), "candidate tag SNPs")
tags.head(15).round(2)

162 candidate tag SNPs

At Cyp6 no SNP separates the swept haplotype perfectly: the best tags are still found on about a sixth of the “other” haplotypes. The dendrogram explains why: a small group of haplotypes sits just a few SNPs away from the swept cluster and shares most of its alleles. Those relatives may carry the duplication too, which the next check helps to judge.

Phased haplotypes only include SNPs that pass the phasing filters, so a tag SNP should be checked in the unphased SNP calls too, and in more places. Below we compute the frequency of the top tag alleles in every district and year around Busia with snp_allele_frequencies(), alongside the duplication (from Module 3.1). A tag that tracks the duplication everywhere, including in Tororo and Kanungu where the duplication is less common, is more useful than one that only works in Busia.

top = tags.head(10)
frq = ag3.snp_allele_frequencies(
    region=cfg["cluster_region"],
    cohorts="admin2_year",
    sample_sets=SAMPLE_SETS,
    sample_query=REGION_QUERY,
).reset_index()
frq_tags = frq.merge(top[["position", "alt"]], left_on=["position", "alt_allele"], right_on=["position", "alt"])
ag3.plot_frequencies_heatmap(frq_tags.set_index("label"), index=None,
                             title=f"Candidate tag SNPs at the {LOCUS} sweep");

Compare these frequencies with Module 3.1. In every cohort the top tags have almost exactly the frequency of CYP6P4 I236M, including Kanungu (about 80 %), where only a quarter of mosquitoes carry Cyp6aap_Dup1a. So these SNPs tag the I236M haplotype background on which the duplication arose, not the duplication itself. Around Busia, where nearly all I236M haplotypes also carry the duplication, that makes little difference; elsewhere you would need a direct CNV assay (copy number qPCR or a PCR across the duplication breakpoint). (Remember that CNV frequencies in Module 3.1 are carrier frequencies: if a fraction p of chromosomes carry a duplication, about 1 − (1 − p)² of mosquitoes do.)

Finally we save the candidate tags for Day 4.

out = Path("course-work")
out.mkdir(exist_ok=True)
tags.assign(contig=contig, locus=LOCUS).to_csv(out / f"tag-snps-{LOCUS}.csv", index=False)
print("saved", out / f"tag-snps-{LOCUS}.csv")

saved course-work/tag-snps-Cyp6.csv

What makes a good assay target?¶

When you choose SNPs for primers and probes on Day 4, think about:

Tagging accuracy: close to 100 % on the swept haplotype and close to 0 % elsewhere, in the populations you will monitor (here Busia, Teso North and Siaya).
Position: SNPs inside an amplified region can give distorted genotype calls, because a mosquito with extra copies is not diploid there (PAMCA W7-M1). A tag just outside the duplication can be easier to genotype.
Neighbouring variation: other common SNPs within ~20 bp would sit under primers or probes. AnoPrimer checks this against Ag3 (Module 4.1).
Function vs linkage: a tag only predicts resistance while it stays linked to the causal mutation(s). A functional variant (like I236M) or the CNV itself is more robust over time.

Exercise 1¶

Set LOCUS = "Cyp9k1" and rerun from Step 1. How strong is the H12 peak? How common is the swept haplotype? Does it carry Cyp9k1_Dup8? Do any Ag-vampIR Cyp9k1 tags work in Busia?

Exercise 2¶

Recompute the H12 zoom for Busia 2019 (74 mosquitoes). H12 depends on sample size, so downsample to 24 mosquitoes with cohort_size=24 (passed through to ag3.haplotypes). Has the peak changed since 2016?

busia_2019 = "cohort_admin2_quarter == 'UG-E_Busia_gamb_2019_Q2'"
df_h12_2019, n = h12_region(cfg["h12_region"], busia_2019, cohort_size=24, random_seed=42)
plot_h12(df_h12_2019, sig, gene, "H12, Busia 2019 (24 mosquitoes)")

The peak is still centred on the Cyp6 cluster; compare its height with 2016. Remember that the atlas intervals shaded on this plot were fitted to the 2016 data.

Exercise 3¶

Change the clustering window from 30 kb to 10 kb and to 100 kb (cfg["cluster_region"]). What happens to the size of the largest cluster and to the number of tag SNPs? Why?

Exercise 4¶

Pick two tag SNPs you would take forward to Day 4 and justify your choice using the criteria above. Record them in your gene dossier.

Quiz¶

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

%pip install -q jupyterquiz==2.9.6.4
from jupyterquiz import display_quiz
display_quiz("https://raw.githubusercontent.com/sanjaynagi/anopheles-omics-training/main/quizzes/3-3-sweep-follow-up.json")

Summary¶

An atlas signal can be followed up by recomputing H12 locally and clustering haplotypes inside its focus.
Around Busia, most haplotypes at the Cyp6 cluster belong to one swept haplotype that carries CYP6P4 I236M and the Cyp6aap_Dup1a duplication, the East African triple mutant; the same haplotype is found in Teso North, Kenya.
The Cyp9k1 sweep is the Cyp9k1_Dup8 amplification, carried by about a quarter of haplotypes.
Tag SNPs make CNVs and haplotypes easy to genotype, but must be checked in the populations you will monitor: the best Cyp6 tags follow the I236M background rather than the duplication, and the Ag-vampIR Cyp9k1 tags are not variable around Busia.

Well done!¶

You now have candidate SNPs for the assays you will design in Module 4.1, and a hypothesis to test with the Siaya amplicon data in Module 4.3.

References¶

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. Njoroge et al. (2022)
Garud NR, Messer PW, Buzbas EO, Petrov DA (2015). Recent selective sweeps in North American Drosophila melanogaster show signatures of soft sweeps. PLoS Genetics 11:e1005004. Garud et al. (2015)
Lucas ER et al. (2019). Whole-genome sequencing reveals high complexity of copy number variation at insecticide resistance loci in malaria mosquitoes. Genome Research 29:1250–1261. Lucas et al. (2019)
Nagi SC, Harding NJ, Lawniczak MKN, Donnelly MJ, Miles A. An atlas of positive selection in the genomes of major malaria vectors (in preparation). https://anopheles-genomic-surveillance.github.io/selection-atlas/
AmpSeeker and the Ag-vampIR panel: https://github.com/sanjaynagi/AmpSeeker

References¶

Njoroge, H., van’t Hof, A., Oruni, A., Pipini, D., Nagi, S. C., Lynd, A., Lucas, E. R., Tomlinson, S., Grau‐Bove, X., McDermott, D., Wat’senga, F. T., Manzambi, E. Z., Agossa, F. R., Mokuba, A., Irish, S., Kabula, B., Mbogo, C., Bargul, J., Paine, M. J. I., … Donnelly, M. J. (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(16), 4307–4318. 10.1111/mec.16591
Garud, N. R., Messer, P. W., Buzbas, E. O., & Petrov, D. A. (2015). Recent Selective Sweeps in North American Drosophila melanogaster Show Signatures of Soft Sweeps. PLOS Genetics, 11(2), e1005004. 10.1371/journal.pgen.1005004
Lucas, E. R., Miles, A., Harding, N. J., Clarkson, C. S., Lawniczak, M. K. N., Kwiatkowski, D. P., Weetman, D., & Donnelly, M. J. (2019). Whole-genome sequencing reveals high complexity of copy number variation at insecticide resistance loci in malaria mosquitoes. Genome Research, 29(8), 1250–1261. 10.1101/gr.245795.118