Workshop 6 - Training course in data analysis for genomic surveillance of African malaria vectors
Module 1 - Haplotypes#
Theme: Data
In previous workshops we have analysed data on individual genetic variations, such as single nucleotide polymorphisms. This module introduces a new type of genetic variation data called “haplotypes”. The benefit of haplotype data is that it conveys the combinations of mutations that occur together on the same segment of a DNA sequence. This in turn allows us to study how those DNA sequences are inherited and transmitted between generations. Many novel and powerful methods of analysis are possible with haplotype data, including analyses to detect genome regions under strong positive selection, which we will study in the subsequent modules of this workshop.
Learning objectives#
At the end of this module you will be able to:
Explain what haplotypes are.
Explain why haplotypes are useful for vector genomic surveillance.
Explain how haplotypes are inferred from genotype data (phasing).
Access haplotype data from MalariaGEN and explain how the data are stored.
Lecture#
English#
Français#
Please note that the code in the cells below might differ from that shown in the video. This can happen because Python packages and their dependencies change due to updates, necessitating tweaks to the code.
What are haplotypes?#
Please see the lecture videos for an explanation of what haplotypes are, why they are useful, and how we infer haplotypes from next-generation sequence data.
Accessing MalariaGEN haplotype data#
MalariaGEN routinely generates haplotype data from whole-genome sequencing of Anopheles mosquitoes. These data can be accessed as follows.
Setup#
%pip install -q --no-warn-conflicts malariagen_data
Note: you may need to restart the kernel to use updated packages.
import malariagen_data
Note that authentication is required to access data through the package, please follow the instructions here.
ag3 = malariagen_data.Ag3()
ag3
| MalariaGEN Ag3 API client | |
|---|---|
| Please note that data are subject to terms of use, for more information see the MalariaGEN website or contact support@malariagen.net. See also the Ag3 API docs. | |
| Storage URL | gs://vo_agam_release/ |
| Data releases available | 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9 |
| Results cache | None |
| Cohorts analysis | 20240418 |
| AIM analysis | 20220528 |
| Site filters analysis | dt_20200416 |
| Software version | malariagen_data 9.0.0 |
| Client location | unknown |
Haplotype data can be accessed via the haplotypes() function. Let’s look at the function documentation.
ag3.haplotypes?
Let’s demonstrate this function, accessing haplotype data for a small genome region on chromosome arm 2L from position 2,422,600 to 2,422,700. To keep things simple we’ll also only access haplotypes from 20 randomly chosen samples from Cameroon. This function returns an xarray dataset.
ds_haps = ag3.haplotypes(
region="2L:2,422,600-2,422,700",
analysis="gamb_colu_arab",
sample_sets = '3.0',
sample_query="country == 'Cameroon'",
cohort_size=20,
)
ds_haps
<xarray.Dataset> Size: 1kB
Dimensions: (variants: 26, alleles: 2, samples: 20, ploidy: 2)
Coordinates:
variant_position (variants) int32 104B dask.array<chunksize=(26,), meta=np.ndarray>
variant_contig (variants) uint8 26B dask.array<chunksize=(26,), meta=np.ndarray>
sample_id (samples) object 160B dask.array<chunksize=(11,), meta=np.ndarray>
Dimensions without coordinates: variants, alleles, samples, ploidy
Data variables:
variant_allele (variants, alleles) |S1 52B dask.array<chunksize=(26, 1), meta=np.ndarray>
call_genotype (variants, samples, ploidy) int8 1kB dask.array<chunksize=(26, 3, 2), meta=np.ndarray>
Attributes:
contigs: ('2R', '2L', '3R', '3L', 'X')
analysis: gamb_colu_arabThere are 20 samples and thus 40 haplotypes in this dataset. We also have data for 26 variants (SNPs).
n_samples = ds_haps.sizes["samples"]
n_samples
20
n_haps = n_samples * 2
n_haps
40
n_variants = ds_haps.sizes["variants"]
n_variants
26
The “variant_position” array has data on the genomic coordinates of the SNPs.
pos = ds_haps["variant_position"].values
pos
array([2422601, 2422606, 2422609, 2422610, 2422611, 2422613, 2422617,
2422621, 2422625, 2422629, 2422630, 2422634, 2422643, 2422645,
2422651, 2422652, 2422666, 2422678, 2422680, 2422683, 2422685,
2422687, 2422688, 2422691, 2422695, 2422697], dtype=int32)
The “variant_allele” array has the reference and alternate alleles for the SNPs.
alleles = ds_haps["variant_allele"].values
alleles
array([[b'T', b'A'],
[b'T', b'G'],
[b'C', b'T'],
[b'C', b'A'],
[b'T', b'A'],
[b'C', b'T'],
[b'C', b'A'],
[b'T', b'G'],
[b'C', b'A'],
[b'G', b'A'],
[b'C', b'A'],
[b'T', b'C'],
[b'A', b'G'],
[b'G', b'A'],
[b'T', b'C'],
[b'A', b'T'],
[b'C', b'T'],
[b'G', b'T'],
[b'C', b'G'],
[b'A', b'C'],
[b'A', b'G'],
[b'C', b'A'],
[b'G', b'T'],
[b'T', b'A'],
[b'C', b'T'],
[b'T', b'A']], dtype='|S1')
ref = alleles[:, 0]
ref
array([b'T', b'T', b'C', b'C', b'T', b'C', b'C', b'T', b'C', b'G', b'C',
b'T', b'A', b'G', b'T', b'A', b'C', b'G', b'C', b'A', b'A', b'C',
b'G', b'T', b'C', b'T'], dtype='|S1')
alt = alleles[:, 1]
alt
array([b'A', b'G', b'T', b'A', b'A', b'T', b'A', b'G', b'A', b'A', b'A',
b'C', b'G', b'A', b'C', b'T', b'T', b'T', b'G', b'C', b'G', b'A',
b'T', b'A', b'T', b'A'], dtype='|S1')
The haplotypes themselves are stored in the “call_genotype” array. The shape of this array is (number of variants, number of samples, number of haplotypes per sample).
gt = ds_haps["call_genotype"].values
gt.shape
(26, 20, 2)
This array uses the numeric encoding for the alleles. We can see that by inspecting its dtype and also inspecting the data values.
gt.dtype
dtype('int8')
gt
array([[[0, 0],
[0, 0],
[0, 0],
...,
[0, 0],
[0, 0],
[0, 0]],
[[0, 0],
[0, 0],
[0, 0],
...,
[0, 0],
[0, 0],
[0, 0]],
[[0, 0],
[0, 0],
[0, 0],
...,
[0, 0],
[0, 0],
[0, 0]],
...,
[[0, 0],
[0, 0],
[0, 0],
...,
[0, 0],
[0, 0],
[0, 0]],
[[0, 0],
[0, 0],
[0, 0],
...,
[0, 0],
[0, 0],
[0, 0]],
[[0, 0],
[0, 0],
[0, 0],
...,
[0, 0],
[0, 0],
[0, 0]]], dtype=int8)
The basic numpy representation of this array is a bit hard to understand, so to help build our intuition of this data structure let’s make a simple plotting function to visualise haplotypes.
def plot_haplotypes(
ds,
hide_text=False,
hide_non_seg=False,
):
import plotly.express as px
import allel
# convert to 2D array of haplotypes
pos = ds["variant_position"].values
gt = ds["call_genotype"].values
ht = gt.reshape(gt.shape[0], -1)
# size rows and columns
if hide_text:
# can make smaller because text not shown
row_height = 8
col_width = 5
else:
# make a bit bigger so we can see the text
row_height = 14
col_width = 12
# deal with non-segregating sites
if hide_non_seg:
ac = allel.HaplotypeArray(ht).count_alleles(max_allele=1)
loc_seg = ac.is_segregating()
ht = ht[loc_seg]
# plot a heatmap
fig = px.imshow(
ht.T,
text_auto=not hide_text,
labels=dict(
x="Variants",
y="Haplotypes",
),
width=ht.shape[0] * col_width + 200,
height=ht.shape[1] * row_height + 200,
zmin=0,
zmax=0,
color_continuous_scale="greys",
template="simple_white",
aspect="auto"
)
# visual styling
fig.update_layout(
coloraxis_showscale=False,
plot_bgcolor="#cccccc",
)
fig.update_traces(xgap=0, ygap=1)
fig.update_xaxes(tickvals=[])
fig.update_yaxes(tickvals=[])
return fig
This function will visualise haplotypes as a heatmap, with haplotypes as rows and variants as columns. A reference allele has a white background, and an alternate allele has a black background.
Let’s use the function to visualise the haplotypes we accessed above.
plot_haplotypes(ds_haps)
In this set of haplotype data we can see there are only two variants which are segregating in this set of samples. All the other variants are non-segregating, which means all haplotypes carry the same allele (in this case a “0” meaning the reference allele).
The non-segregating variants are not particularly informative, so let’s hide them.
plot_haplotypes(ds_haps, hide_non_seg=True)
We can see now that within this genome region, and within this set of 20 samples, we have two segregating sites and three unique haplotypes.
Well done!#
Hopefully this module has been a useful introduction to haplotype data. There are no practical exercises for this module, but you might like to launch this notebook for yourself, and try running the code examples above, or changing them to load data from a different genome region or a different set of samples.