Ag3 cloud data access

This page provides information about how to access data from Anopheles gambiae 1000 Genomes project (Ag1000G) phase 3 via Google Cloud. This includes sample metadata and single nucleotide polymorphism (SNP) calls.

This notebook illustrates how to read data directly from the cloud, without having to first download any data locally. This notebook can be run from any computer, but will work best when run from a compute node within Google Cloud, because it will be physically closer to the data and so data transfer is faster. For example, this notebook can be run via MyBinder or Google Colab which are free interactive computing service running in the cloud.

To launch this notebook in the cloud and run it for yourself, click the launch icon () at the top of the page and select one of the cloud computing services available.

Data hosting

All data required for this notebook is hosted on Google Cloud Storage (GCS). Data are hosted in the vo_agam_release bucket, which is a multi-region bucket located in the United States. All data hosted in GCS are publicly accessible and do not require any authentication to access.

Setup

Running this notebook requires the following Python packages to be installed: numpy, dask, zarr, gcsfs, fsspec. These packages can be installed via pip or conda. E.g.:

!pip install -q numpy dask[array] zarr gcsfs fsspec

To make accessing these data more convenient, we’ve also created a malariagen_data Python package, which is available from PyPI. This is experimental so please let us know if you find any bugs or have any suggestions. The malariagen_data package can be installed via pip, e.g.:

!pip install -q malariagen-data

Once installed, data access from Google Cloud is set up with the following code:

import malariagen_data
ag3 = malariagen_data.Ag3("gs://vo_agam_release/")

Sample sets

Data in this release are organised into 28 sample sets. Each of these sample sets corresponds to a set of mosquito specimens contributed by a collaborating study. Depending on your objectives, you may want to access data from only specific sample sets, or all sample sets.

To see which sample sets are available, load the sample set manifest into a pandas dataframe:

df_sample_sets = ag3.sample_sets()
df_sample_sets
sample_set sample_count release
0 AG1000G-AO 81 v3
1 AG1000G-BF-A 181 v3
2 AG1000G-BF-B 102 v3
3 AG1000G-BF-C 13 v3
4 AG1000G-CD 76 v3
5 AG1000G-CF 73 v3
6 AG1000G-CI 80 v3
7 AG1000G-CM-A 303 v3
8 AG1000G-CM-B 97 v3
9 AG1000G-CM-C 44 v3
10 AG1000G-FR 23 v3
11 AG1000G-GA-A 69 v3
12 AG1000G-GH 100 v3
13 AG1000G-GM-A 74 v3
14 AG1000G-GM-B 31 v3
15 AG1000G-GM-C 174 v3
16 AG1000G-GN-A 45 v3
17 AG1000G-GN-B 185 v3
18 AG1000G-GQ 10 v3
19 AG1000G-GW 101 v3
20 AG1000G-KE 86 v3
21 AG1000G-ML-A 60 v3
22 AG1000G-ML-B 71 v3
23 AG1000G-MW 41 v3
24 AG1000G-MZ 74 v3
25 AG1000G-TZ 300 v3
26 AG1000G-UG 290 v3
27 AG1000G-X 297 v3

For more information about these sample sets, see the section on sample sets in the introduction to Ag1000G phase 3.

Sample metadata

Data about the samples that were sequenced to generate this data resource are available, including the time and place of collection, the gender of the specimen, and our call regarding the species of the specimen. These are organised by sample set.

E.g., load sample metadata for the AG1000G-BF-A and AG1000G-BF-B sample sets into a pandas dataframe:

df_samples = ag3.sample_metadata(sample_sets=["AG1000G-BF-A", "AG1000G-BF-B"])
df_samples
sample_id partner_sample_id contributor country location year month latitude longitude sex_call sample_set release aim_fraction_colu aim_fraction_arab species_gambcolu_arabiensis species_gambiae_coluzzii species
0 AB0085-Cx BF2-4 Austin Burt Burkina Faso Pala 2012 7 11.150 -4.235 F AG1000G-BF-A v3 0.024 0.002 gamb_colu gambiae gambiae
1 AB0086-Cx BF2-6 Austin Burt Burkina Faso Pala 2012 7 11.150 -4.235 F AG1000G-BF-A v3 0.038 0.002 gamb_colu gambiae gambiae
2 AB0087-C BF3-3 Austin Burt Burkina Faso Bana 2012 7 11.233 -4.472 F AG1000G-BF-A v3 0.982 0.002 gamb_colu coluzzii coluzzii
3 AB0088-C BF3-5 Austin Burt Burkina Faso Bana 2012 7 11.233 -4.472 F AG1000G-BF-A v3 0.990 0.002 gamb_colu coluzzii coluzzii
4 AB0089-Cx BF3-8 Austin Burt Burkina Faso Bana 2012 7 11.233 -4.472 F AG1000G-BF-A v3 0.975 0.002 gamb_colu coluzzii coluzzii
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
278 AB0533-C BF13-18 Austin Burt Burkina Faso Souroukoudinga 2014 7 11.235 -4.535 F AG1000G-BF-B v3 0.021 0.002 gamb_colu gambiae gambiae
279 AB0536-C BF13-31 Austin Burt Burkina Faso Souroukoudinga 2014 7 11.235 -4.535 F AG1000G-BF-B v3 0.025 0.002 gamb_colu gambiae gambiae
280 AB0537-C BF13-32 Austin Burt Burkina Faso Souroukoudinga 2014 7 11.235 -4.535 F AG1000G-BF-B v3 0.029 0.002 gamb_colu gambiae gambiae
281 AB0538-C BF13-33 Austin Burt Burkina Faso Souroukoudinga 2014 7 11.235 -4.535 F AG1000G-BF-B v3 0.018 0.002 gamb_colu gambiae gambiae
282 AB0408-C BF14-20 Austin Burt Burkina Faso Bana 2014 7 11.233 -4.472 F AG1000G-BF-B v3 0.986 0.002 gamb_colu coluzzii coluzzii

283 rows × 17 columns

The sample_id column gives the sample identifier used throughout all Ag1000G analyses.

The country, location, latitude and longitude columns give the location where the specimen was collected.

The year and month columns give the approximate date when the specimen was collected.

The sex_call column gives the gender as determined from the sequence data.

To load metadata for all wild-caught samples, you can use the shortcut “v3_wild”, e.g.:

df_samples = ag3.sample_metadata(sample_sets="v3_wild")
df_samples
sample_id partner_sample_id contributor country location year month latitude longitude sex_call sample_set release aim_fraction_colu aim_fraction_arab species_gambcolu_arabiensis species_gambiae_coluzzii species
0 AR0047-C LUA047 Joao Pinto Angola Luanda 2009 4 -8.884 13.302 F AG1000G-AO v3 0.945 0.001 gamb_colu coluzzii coluzzii
1 AR0049-C LUA049 Joao Pinto Angola Luanda 2009 4 -8.884 13.302 F AG1000G-AO v3 0.933 0.001 gamb_colu coluzzii coluzzii
2 AR0051-C LUA051 Joao Pinto Angola Luanda 2009 4 -8.884 13.302 F AG1000G-AO v3 0.937 0.002 gamb_colu coluzzii coluzzii
3 AR0061-C LUA061 Joao Pinto Angola Luanda 2009 4 -8.884 13.302 F AG1000G-AO v3 0.938 0.002 gamb_colu coluzzii coluzzii
4 AR0078-C LUA078 Joao Pinto Angola Luanda 2009 4 -8.884 13.302 F AG1000G-AO v3 0.926 0.001 gamb_colu coluzzii coluzzii
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2779 AC0295-C K92 Martin Donnelly Uganda Kihihi 2012 11 -0.751 29.701 F AG1000G-UG v3 0.026 0.002 gamb_colu gambiae gambiae
2780 AC0296-C K93 Martin Donnelly Uganda Kihihi 2012 11 -0.751 29.701 F AG1000G-UG v3 0.029 0.003 gamb_colu gambiae gambiae
2781 AC0297-C K94 Martin Donnelly Uganda Kihihi 2012 11 -0.751 29.701 F AG1000G-UG v3 0.026 0.002 gamb_colu gambiae gambiae
2782 AC0298-C K95 Martin Donnelly Uganda Kihihi 2012 11 -0.751 29.701 F AG1000G-UG v3 0.029 0.002 gamb_colu gambiae gambiae
2783 AC0299-C K96 Martin Donnelly Uganda Kihihi 2012 11 -0.751 29.701 F AG1000G-UG v3 0.029 0.002 gamb_colu gambiae gambiae

2784 rows × 17 columns

Pandas can be used to explore and query the sample metadata in various ways. E.g., here is a summary of the numbers of samples by species:

df_samples.groupby("species").size()
species
arabiensis                          368
coluzzii                            675
gambiae                            1571
intermediate_arabiensis_gambiae       1
intermediate_gambiae_coluzzii       169
dtype: int64

Note that samples within a sample set may belong to different species. For convenience, we have made a species call for all samples in Ag1000G phase 3, using the genomic data. Calling species is not always straightforward, and we have used two different methods for species calling, ancestry informative markers (AIM) and principal components analysis (PCA). When loading the sample metadata, the AIM species calls are included by default. The results of these two different methods generally agree, although there are some populations where results are different, particularly in Guinea-Bissau, The Gambia, Kenya and Tanzania. If you have any questions about how to interpret these species calls, please get in touch.

SNP sites and alleles

We have called SNP genotypes in all samples at all positions in the genome where the reference allele is not “N”. Data on this set of genomic positions and alleles for a given chromosome arm (e.g., 3R) can be accessed as dask arrays as follows:

pos, ref, alt = ag3.snp_sites("3R")
pos
Array Chunk
Bytes 208.91 MB 2.10 MB
Shape (52226568,) (524288,)
Count 101 Tasks 100 Chunks
Type int32 numpy.ndarray
52226568 1
ref
Array Chunk
Bytes 52.23 MB 524.29 kB
Shape (52226568,) (524288,)
Count 101 Tasks 100 Chunks
Type |S1 numpy.ndarray
52226568 1
alt
Array Chunk
Bytes 156.68 MB 1.57 MB
Shape (52226568, 3) (524288, 3)
Count 101 Tasks 100 Chunks
Type |S1 numpy.ndarray
3 52226568

Data can be loaded into memory as numpy arrays as shown in the following examples.

# read first 10 SNP positions into a numpy array
p = pos[:10].compute()
p
array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10], dtype=int32)
# read first 10 SNP reference alleles
r = ref[:10].compute()
r
array([b'C', b'C', b'T', b'C', b'T', b'A', b'C', b'G', b'T', b'T'],
      dtype='|S1')
# read first 10 SNP alternate alleles
a = alt[:10].compute()
a
array([[b'A', b'T', b'G'],
       [b'A', b'T', b'G'],
       [b'A', b'C', b'G'],
       [b'A', b'T', b'G'],
       [b'A', b'C', b'G'],
       [b'C', b'T', b'G'],
       [b'A', b'T', b'G'],
       [b'A', b'C', b'T'],
       [b'A', b'C', b'G'],
       [b'A', b'C', b'G']], dtype='|S1')

Note that we have chosen to genotype all samples at all sites in the genome, assuming all possible SNP alleles. Not all of these alternate alleles will actually have been observed in the Ag3 samples. To determine which sites and alleles are segregating, an allele count can be performed over the samples you are interested in. See the example below.

Site filters

SNP calling is not always reliable, and we have created some site filters to allow excluding low quality SNPs.

Because different species may have different genome accessibility issues, we have created three separate site filters:

  • The “gamb_colu” filter is design for working only with An. gambiae and/or An. coluzzii samples.

  • The “arab” filter is designed for working with An. arabiensis samples.

  • The “gamb_colu_arab” filter is suitable for when analysing samples of any species together.

Each set of site filters provides a “filter_pass” Boolean mask for each chromosome arm, where True indicates that the site passed the filter and is accessible to high quality SNP calling.

The site filters data can be accessed as dask arrays as shown in the examples below.

# access gamb_colu_arab site filters for chromosome arm 3R as a dask array
filter_pass = ag3.site_filters("3R", mask="gamb_colu_arab")
filter_pass
Array Chunk
Bytes 52.23 MB 300.00 kB
Shape (52226568,) (300000,)
Count 176 Tasks 175 Chunks
Type bool numpy.ndarray
52226568 1
# read filter values for first 10 SNPs (True means the site passes filters)
f = filter_pass[:10].compute()
f
array([False, False, False, False, False, False, False, False, False,
       False])

SNP genotypes

SNP genotypes for individual samples are available. Genotypes are stored as a three-dimensional array, where the first dimension corresponds to genomic positions, the second dimension is samples, and the third dimension is ploidy (2). Values coded as integers, where -1 represents a missing value, 0 represents the reference allele, and 1, 2, and 3 represent alternate alleles.

SNP genotypes can be accessed as dask arrays as shown below.

gt = ag3.snp_genotypes("3R", sample_sets="v3_wild")
gt
Array Chunk
Bytes 290.80 GB 30.00 MB
Shape (52226568, 2784, 2) (300000, 50, 2)
Count 23827 Tasks 11900 Chunks
Type int8 numpy.ndarray
2 2784 52226568

Note that the columns of this array (second dimension) match the rows in the sample metadata, if the same sample sets were loaded. I.e.:

contig = '3R'
df_samples = ag3.sample_metadata(sample_sets="v3_wild")
gt = ag3.snp_genotypes(contig=contig, sample_sets="v3_wild")
len(df_samples) == gt.shape[1]
True

You can use this correspondance to apply further subsetting operations to the genotypes by querying the sample metadata. E.g.:

import dask.array as da
loc_gambiae = df_samples.query("species == 'gambiae'").index.values
print(f"found {len(loc_gambiae)} gambiae samples")
gt_gambiae = da.take(gt, loc_gambiae, axis=1)
gt_gambiae
found 1571 gambiae samples
Array Chunk
Bytes 164.10 GB 30.00 MB
Shape (52226568, 1571, 2) (300000, 50, 2)
Count 33627 Tasks 9800 Chunks
Type int8 numpy.ndarray
2 1571 52226568

Data can be read into memory as numpy arrays, e.g., read genotypes for the first 5 SNPs and the first 3 samples:

g = gt[:5, :3, :].compute()
g
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]]], dtype=int8)

If you want to work with the genotype calls, you may find it convenient to use scikit-allel:

 !pip install -q scikit-allel

E.g., this code sets up a genotype array:

# use the scikit-allel wrapper class for genotype calls
import allel
gt = ag3.snp_genotypes("3R", sample_sets=["AG1000G-BF-A", "AG1000G-BF-B"])
gt = allel.GenotypeDaskArray(gt)
gt
<GenotypeDaskArray shape=(52226568, 283, 2) dtype=int8>
01234...278279280281282
00/00/00/00/00/0...0/00/00/00/00/0
10/00/00/00/00/0...0/00/00/00/00/0
20/00/00/00/00/0...0/00/00/00/00/0
......
52226565./../../../../....0/0./../../../.
52226566./../../../../...../../../../../.
52226567./../../../../...../../../../../.

Here’s an example computation to count the number of segregating SNPs on chromosome arm 3R that also pass gamb_colu_arab site filters. This may take a little while, because it is scanning genotype calls at millions of SNPs in hundreds of samples:

# import dask progress bar
from dask.diagnostics.progress import ProgressBar

# import numpy
import numpy as np

# choose chromosome arm
contig = "3R"

# choose site filter mask
mask = "gamb_colu_arab"

# choose sample sets
sample_sets = ["AG1000G-BF-A", "AG1000G-BF-B"]

# locate pass sites
loc_pass = ag3.site_filters(contig=contig, mask=mask).compute()

# perform an allele count over genotypes
gt = ag3.snp_genotypes(contig=contig, sample_sets=sample_sets)
gt = allel.GenotypeDaskArray(gt)
ac = gt.count_alleles(max_allele=3)

# locate segregating sites
loc_seg = ac.is_segregating()

# count segregating and pass sites
n_pass_seg = da.count_nonzero(loc_pass & loc_seg)

# run the computation
with ProgressBar():
    n_pass_seg = n_pass_seg.compute()

n_pass_seg
[########################################] | 100% Completed |  3min  6.7s
11226723

Running larger computations

Please note that free cloud computing services such as Google Colab and MyBinder provide only limited computing resources. Thus although these services are able to efficiently read Ag3 data stored on Google Cloud, you may find that you run out of memory, or computations take a long time running on a single core. If you would like any suggestions regarding how to set up more powerful compute resources in the cloud, please feel free to get in touch via the malariagen/vector-data GitHub discussion board.

Feedback and suggestions

If there are particular analyses you would like to run, or if you have other suggestions for useful documentation we could add to this site, we would love to know, please get in touch via the malariagen/vector-data GitHub discussion board.