Adar1 cloud data access

This notebook provides information about how to download data from the MalariaGEN Vector Observatory Anopheles darlingi Genomic Surveillance Project, for Anopheles darlingi via Google Cloud. This includes sample metadata, raw sequence reads, sequence read alignments, 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 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_adar_release_master_us_central1 bucket, which is a single-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 some Python packages to be installed:

%pip install -q malariagen_data

To make accessing these data more convenient, we’ve created the malariagen_data Python package. This is experimental so please let us know if you find any bugs or have any suggestions. See the Adar1.0 API docs for documentation of all functions available from this package.

Import other packages we’ll need to use here.

import numpy as np
import dask
import dask.array as da
from dask.diagnostics.progress import ProgressBar
# silence some warnings
dask.config.set(**{'array.slicing.split_large_chunks': False})
import allel
import malariagen_data

Adar1 data access from Google Cloud is set up with the following code:

adar1 = malariagen_data.Adar1()
adar1

MalariaGEN Adar1 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 Adir1 API docs.
Storage URL	gs://vo_adar_release_master_us_central1
Data releases available	1.0
Results cache	None
Cohorts analysis	20260303
Site filters analysis	dt_20260301
Software version	malariagen_data 15.6.0.post407+1430baf
Client location	Iowa, United States (Google Cloud us-central1)
Data filtered for unrestricted use only	False
Data filtered for surveillance use only	False
Relevant data releases	1.0

Sample sets

Data are organised into different releases. As an example, data in Adar1.0 are organised into 6 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 = adar1.sample_sets(release="1.0")
df_sample_sets

	sample_set	sample_count	study_id	study_url	terms_of_use_expiry_date	terms_of_use_url	release	unrestricted_use
0	1357-VO-BR-SALLUM-VMF00326	272	1357-VO-BR-SALLUM	https://www.malariagen.net/partner_study/1357-...	2027-11-30	https://www.malariagen.net/data/our-approach-s...	1.0	False
1	1358-VO-GF-GENDRIN-VMF00327	139	1358-VO-GF-GENDRIN	https://www.malariagen.net/partner_study/1358-...	2027-11-30	https://www.malariagen.net/data/our-approach-s...	1.0	False
2	1359-VO-GY-NILES-ROBIN-VMF00328	18	1359-VO-GY-NILES-ROBIN	https://www.malariagen.net/partner_study/1359-...	2027-11-30	https://www.malariagen.net/data/our-approach-s...	1.0	False
3	1360-VO-PE-GAMBOA-VMF00329	89	1360-VO-PE-GAMBOA	https://www.malariagen.net/partner_study/1360-...	2027-11-30	https://www.malariagen.net/data/our-approach-s...	1.0	False
4	1361-VO-VE-GRILLET-VMF00330	126	1361-VO-VE-GRILLET	https://www.malariagen.net/partner_study/1361-...	2027-11-30	https://www.malariagen.net/data/our-approach-s...	1.0	False
5	1362-VO-CO-QUINONES-VMF00331	449	1362-VO-CO-QUINONES	https://www.malariagen.net/partner_study/1362-...	2027-11-30	https://www.malariagen.net/data/our-approach-s...	1.0	False

For more information about these sample sets, you can read about each sample set from the URLs under the field study_url.

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 all samples in the Adar1.0 release into a pandas DataFrame:

df_samples = adar1.sample_metadata(sample_sets="1.0")
df_samples

	sample_id	partner_sample_id	contributor	country	location	year	month	latitude	longitude	sex_call	...	admin1_name	admin1_iso	admin2_name	taxon	cohort_admin1_year	cohort_admin1_month	cohort_admin1_quarter	cohort_admin2_year	cohort_admin2_month	cohort_admin2_quarter
0	VMF00326-0001	Coari_AM252_1	Maria Anice Mureb Sallum	Brazil	Coari	2022	11	-4.137	-63.150	UKN	...	Amazonas	BR-AM	Coari	darlingi	BR-AM_darl_2022	BR-AM_darl_2022_11	BR-AM_darl_2022_Q4	BR-AM_Coari_darl_2022	BR-AM_Coari_darl_2022_11	BR-AM_Coari_darl_2022_Q4
1	VMF00326-0002	Coari_AM252_2	Maria Anice Mureb Sallum	Brazil	Coari	2022	11	-4.137	-63.150	UKN	...	Amazonas	BR-AM	Coari	darlingi	BR-AM_darl_2022	BR-AM_darl_2022_11	BR-AM_darl_2022_Q4	BR-AM_Coari_darl_2022	BR-AM_Coari_darl_2022_11	BR-AM_Coari_darl_2022_Q4
2	VMF00326-0003	Coari_AM252_3	Maria Anice Mureb Sallum	Brazil	Coari	2022	11	-4.137	-63.150	UKN	...	Amazonas	BR-AM	Coari	darlingi	BR-AM_darl_2022	BR-AM_darl_2022_11	BR-AM_darl_2022_Q4	BR-AM_Coari_darl_2022	BR-AM_Coari_darl_2022_11	BR-AM_Coari_darl_2022_Q4
3	VMF00326-0004	Coari_AM252_4	Maria Anice Mureb Sallum	Brazil	Coari	2022	11	-4.137	-63.150	UKN	...	Amazonas	BR-AM	Coari	darlingi	BR-AM_darl_2022	BR-AM_darl_2022_11	BR-AM_darl_2022_Q4	BR-AM_Coari_darl_2022	BR-AM_Coari_darl_2022_11	BR-AM_Coari_darl_2022_Q4
4	VMF00326-0005	Coari_AM252_5	Maria Anice Mureb Sallum	Brazil	Coari	2022	11	-4.137	-63.150	UKN	...	Amazonas	BR-AM	Coari	darlingi	BR-AM_darl_2022	BR-AM_darl_2022_11	BR-AM_darl_2022_Q4	BR-AM_Coari_darl_2022	BR-AM_Coari_darl_2022_11	BR-AM_Coari_darl_2022_Q4
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
1089	VMF00331-0445	Tagachi_HER0112286021	Martha Quiñones	Colombia	Tagachí	2022	4	6.221	-76.726	UKN	...	Chocó	CO-CHO	Quibdó	darlingi	CO-CHO_darl_2022	CO-CHO_darl_2022_04	CO-CHO_darl_2022_Q2	CO-CHO_Quibdó_darl_2022	CO-CHO_Quibdó_darl_2022_04	CO-CHO_Quibdó_darl_2022_Q2
1090	VMF00331-0446	Tagachi_HER0112286022	Martha Quiñones	Colombia	Tagachí	2022	4	6.221	-76.726	UKN	...	Chocó	CO-CHO	Quibdó	darlingi	CO-CHO_darl_2022	CO-CHO_darl_2022_04	CO-CHO_darl_2022_Q2	CO-CHO_Quibdó_darl_2022	CO-CHO_Quibdó_darl_2022_04	CO-CHO_Quibdó_darl_2022_Q2
1091	VMF00331-0447	Tagachi_HER0112286023	Martha Quiñones	Colombia	Tagachí	2022	4	6.221	-76.726	UKN	...	Chocó	CO-CHO	Quibdó	darlingi	CO-CHO_darl_2022	CO-CHO_darl_2022_04	CO-CHO_darl_2022_Q2	CO-CHO_Quibdó_darl_2022	CO-CHO_Quibdó_darl_2022_04	CO-CHO_Quibdó_darl_2022_Q2
1092	VMF00331-0448	Tagachi_HER0112286024	Martha Quiñones	Colombia	Tagachí	2022	4	6.221	-76.726	UKN	...	Chocó	CO-CHO	Quibdó	darlingi	CO-CHO_darl_2022	CO-CHO_darl_2022_04	CO-CHO_darl_2022_Q2	CO-CHO_Quibdó_darl_2022	CO-CHO_Quibdó_darl_2022_04	CO-CHO_Quibdó_darl_2022_Q2
1093	VMF00331-0449	Tagachi_HER0112286037	Martha Quiñones	Colombia	Tagachí	2022	4	6.221	-76.726	UKN	...	Chocó	CO-CHO	Quibdó	darlingi	CO-CHO_darl_2022	CO-CHO_darl_2022_04	CO-CHO_darl_2022_Q2	CO-CHO_Quibdó_darl_2022	CO-CHO_Quibdó_darl_2022_04	CO-CHO_Quibdó_darl_2022_Q2

1094 rows × 46 columns

The sample_id column gives the sample identifier used throughout all Adir1.0 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.

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("taxon").size()

taxon
darlingi    1094
dtype: int64

SNP calls

Data on SNP calls, including the SNP positions, alleles, site filters, and genotypes, can be accessed as an xarray Dataset.

E.g., access SNP calls for contig 2 for all samples in Adar1.0.

ds_snps = adar1.snp_calls(region="2", sample_sets="1.0")
ds_snps

<xarray.Dataset> Size: 569GB
Dimensions:                       (variants: 30593416, alleles: 4,
                                   samples: 1094, ploidy: 2)
Coordinates:
    variant_position              (variants) int32 122MB dask.array<chunksize=(524288,), meta=np.ndarray>
    variant_contig                (variants) uint8 31MB dask.array<chunksize=(524288,), meta=np.ndarray>
    sample_id                     (samples) <U36 158kB dask.array<chunksize=(273,), meta=np.ndarray>
Dimensions without coordinates: variants, alleles, samples, ploidy
Data variables:
    variant_allele                (variants, alleles) |S1 122MB dask.array<chunksize=(524288, 4), meta=np.ndarray>
    variant_filter_pass_darlingi  (variants) bool 31MB dask.array<chunksize=(524288,), meta=np.ndarray>
    call_genotype                 (variants, samples, ploidy) int8 67GB dask.array<chunksize=(300000, 50, 2), meta=np.ndarray>
    call_GQ                       (variants, samples) int8 33GB dask.array<chunksize=(300000, 50), meta=np.ndarray>
    call_MQ                       (variants, samples) float32 134GB dask.array<chunksize=(300000, 50), meta=np.ndarray>
    call_AD                       (variants, samples, alleles) int16 268GB dask.array<chunksize=(300000, 50, 4), meta=np.ndarray>
    call_genotype_mask            (variants, samples, ploidy) bool 67GB dask.array<chunksize=(300000, 50, 2), meta=np.ndarray>
Attributes:
    contigs:  ('2', '3', 'X')

xarray.Dataset

• Dimensions:
  • variants: 30593416
  • alleles: 4
  • samples: 1094
  • ploidy: 2
• Coordinates: (3)
  • variant_position
    (variants)
    int32
    dask.array<chunksize=(524288,), meta=np.ndarray>

    Array Chunk Bytes 116.70 MiB 2.00 MiB Shape (30593416,) (524288,) Dask graph 59 chunks in 1 graph layer Data type int32 numpy.ndarray

  • variant_contig
    (variants)
    uint8
    dask.array<chunksize=(524288,), meta=np.ndarray>

    Array Chunk Bytes 29.18 MiB 512.00 kiB Shape (30593416,) (524288,) Dask graph 59 chunks in 1 graph layer Data type uint8 numpy.ndarray

  • sample_id
    (samples)
    <U36
    dask.array<chunksize=(273,), meta=np.ndarray>

    Array Chunk Bytes 153.84 kiB 63.14 kiB Shape (1094,) (449,) Dask graph 6 chunks in 13 graph layers Data type

• Data variables: (7)
  • variant_allele
    (variants, alleles)
    |S1
    dask.array<chunksize=(524288, 4), meta=np.ndarray>

    Array Chunk Bytes 116.70 MiB 2.00 MiB Shape (30593416, 4) (524288, 4) Dask graph 59 chunks in 5 graph layers Data type |S1 numpy.ndarray

  • variant_filter_pass_darlingi
    (variants)
    bool
    dask.array<chunksize=(524288,), meta=np.ndarray>

    Array Chunk Bytes 29.18 MiB 512.00 kiB Shape (30593416,) (524288,) Dask graph 59 chunks in 1 graph layer Data type bool numpy.ndarray

  • call_genotype
    (variants, samples, ploidy)
    int8
    dask.array<chunksize=(300000, 50, 2), meta=np.ndarray>

    Array Chunk Bytes 62.34 GiB 28.61 MiB Shape (30593416, 1094, 2) (300000, 50, 2) Dask graph 2448 chunks in 7 graph layers Data type int8 numpy.ndarray

  • call_GQ
    (variants, samples)
    int8
    dask.array<chunksize=(300000, 50), meta=np.ndarray>

    Array Chunk Bytes 31.17 GiB 14.31 MiB Shape (30593416, 1094) (300000, 50) Dask graph 2448 chunks in 7 graph layers Data type int8 numpy.ndarray

  • call_MQ
    (variants, samples)
    float32
    dask.array<chunksize=(300000, 50), meta=np.ndarray>

    Array Chunk Bytes 124.68 GiB 57.22 MiB Shape (30593416, 1094) (300000, 50) Dask graph 2448 chunks in 7 graph layers Data type float32 numpy.ndarray

  • call_AD
    (variants, samples, alleles)
    int16
    dask.array<chunksize=(300000, 50, 4), meta=np.ndarray>

    Array Chunk Bytes 249.36 GiB 114.44 MiB Shape (30593416, 1094, 4) (300000, 50, 4) Dask graph 2448 chunks in 7 graph layers Data type int16 numpy.ndarray

  • call_genotype_mask
    (variants, samples, ploidy)
    bool
    dask.array<chunksize=(300000, 50, 2), meta=np.ndarray>

    Array Chunk Bytes 62.34 GiB 28.61 MiB Shape (30593416, 1094, 2) (300000, 50, 2) Dask graph 2448 chunks in 8 graph layers Data type bool numpy.ndarray

• Indexes: (0)
• Attributes: (1)

  contigs :

  ('2', '3', 'X')

The arrays within this dataset are backed by Dask arrays, and can be accessed as shown below.

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 (e.g., 2RL) can be accessed as Dask arrays as follows.

pos = ds_snps["variant_position"].data
pos

Array Chunk Bytes 116.70 MiB 2.00 MiB Shape (30593416,) (524288,) Dask graph 59 chunks in 1 graph layer Data type int32 numpy.ndarray

alleles = ds_snps["variant_allele"].data
alleles

Array Chunk Bytes 116.70 MiB 2.00 MiB Shape (30593416, 4) (524288, 4) Dask graph 59 chunks in 5 graph layers Data type |S1 numpy.ndarray

Data can be loaded into memory as a NumPy array as shown in the following examples.

# read first 10 SNP positions into a numpy array
p = pos[:10].compute()
p

array([11, 12, 13, 28, 36, 43, 45, 53, 54, 55], dtype=int32)

# read first 10 SNP alleles into a numpy array
a = alleles[:10].compute()
a

array([[b'C', b'A', b'G', b'T'],
       [b'T', b'A', b'C', b'G'],
       [b'G', b'A', b'C', b'T'],
       [b'C', b'A', b'G', b'T'],
       [b'A', b'C', b'G', b'T'],
       [b'G', b'A', b'C', b'T'],
       [b'A', b'C', b'G', b'T'],
       [b'G', b'A', b'C', b'T'],
       [b'A', b'C', b'G', b'T'],
       [b'G', b'A', b'C', b'T']], dtype='|S1')

Here the first column contains the reference alleles, and the remaining columns contain the alternate alleles.

Note that a byte string data type is used here for efficiency. E.g., the Python code b'T' represents a byte string containing the letter “T”, which here stands for the nucleotide thymine.

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 Adar1 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.

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 as a dask array
filter_pass = ds_snps['variant_filter_pass_darlingi'].data
filter_pass

Array Chunk Bytes 29.18 MiB 512.00 kiB Shape (30593416,) (524288,) Dask graph 59 chunks in 1 graph layer Data type bool numpy.ndarray

# read filter values for first 10 SNPs (True means the site passes filters)
f = filter_pass[:10].compute()
f

array([ True,  True,  True,  True,  True,  True,  True,  True,  True,
        True])

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 are 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 = ds_snps["call_genotype"].data
gt

Array Chunk Bytes 62.34 GiB 28.61 MiB Shape (30593416, 1094, 2) (300000, 50, 2) Dask graph 2448 chunks in 7 graph layers Data type int8 numpy.ndarray

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.:

df_samples = adar1.sample_metadata(sample_sets="1.0")
gt = ds_snps["call_genotype"].data
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.:

loc_darlingi = df_samples.eval("taxon == 'darlingi'").values
print(f"found {np.count_nonzero(loc_darlingi)} darlingi samples")

found 1094 darlingi samples

ds_snps_darlingi = ds_snps.isel(samples=loc_darlingi)
ds_snps_darlingi

<xarray.Dataset> Size: 569GB
Dimensions:                       (variants: 30593416, alleles: 4,
                                   samples: 1094, ploidy: 2)
Coordinates:
    variant_position              (variants) int32 122MB dask.array<chunksize=(524288,), meta=np.ndarray>
    variant_contig                (variants) uint8 31MB dask.array<chunksize=(524288,), meta=np.ndarray>
    sample_id                     (samples) <U36 158kB dask.array<chunksize=(182,), meta=np.ndarray>
Dimensions without coordinates: variants, alleles, samples, ploidy
Data variables:
    variant_allele                (variants, alleles) |S1 122MB dask.array<chunksize=(524288, 4), meta=np.ndarray>
    variant_filter_pass_darlingi  (variants) bool 31MB dask.array<chunksize=(524288,), meta=np.ndarray>
    call_genotype                 (variants, samples, ploidy) int8 67GB dask.array<chunksize=(300000, 45, 2), meta=np.ndarray>
    call_GQ                       (variants, samples) int8 33GB dask.array<chunksize=(300000, 45), meta=np.ndarray>
    call_MQ                       (variants, samples) float32 134GB dask.array<chunksize=(300000, 45), meta=np.ndarray>
    call_AD                       (variants, samples, alleles) int16 268GB dask.array<chunksize=(300000, 45, 4), meta=np.ndarray>
    call_genotype_mask            (variants, samples, ploidy) bool 67GB dask.array<chunksize=(300000, 45, 2), meta=np.ndarray>
Attributes:
    contigs:  ('2', '3', 'X')

xarray.Dataset

• Dimensions:
  • variants: 30593416
  • alleles: 4
  • samples: 1094
  • ploidy: 2
• Coordinates: (3)
  • variant_position
    (variants)
    int32
    dask.array<chunksize=(524288,), meta=np.ndarray>

    Array Chunk Bytes 116.70 MiB 2.00 MiB Shape (30593416,) (524288,) Dask graph 59 chunks in 1 graph layer Data type int32 numpy.ndarray

  • variant_contig
    (variants)
    uint8
    dask.array<chunksize=(524288,), meta=np.ndarray>

    Array Chunk Bytes 29.18 MiB 512.00 kiB Shape (30593416,) (524288,) Dask graph 59 chunks in 1 graph layer Data type uint8 numpy.ndarray

  • sample_id
    (samples)
    <U36
    dask.array<chunksize=(182,), meta=np.ndarray>

    Array Chunk Bytes 153.84 kiB 25.59 kiB Shape (1094,) (182,) Dask graph 7 chunks in 14 graph layers Data type

• Data variables: (7)
  • variant_allele
    (variants, alleles)
    |S1
    dask.array<chunksize=(524288, 4), meta=np.ndarray>

    Array Chunk Bytes 116.70 MiB 2.00 MiB Shape (30593416, 4) (524288, 4) Dask graph 59 chunks in 5 graph layers Data type |S1 numpy.ndarray

  • variant_filter_pass_darlingi
    (variants)
    bool
    dask.array<chunksize=(524288,), meta=np.ndarray>

    Array Chunk Bytes 29.18 MiB 512.00 kiB Shape (30593416,) (524288,) Dask graph 59 chunks in 1 graph layer Data type bool numpy.ndarray

  • call_genotype
    (variants, samples, ploidy)
    int8
    dask.array<chunksize=(300000, 45, 2), meta=np.ndarray>

    Array Chunk Bytes 62.34 GiB 25.75 MiB Shape (30593416, 1094, 2) (300000, 45, 2) Dask graph 2550 chunks in 8 graph layers Data type int8 numpy.ndarray

  • call_GQ
    (variants, samples)
    int8
    dask.array<chunksize=(300000, 45), meta=np.ndarray>

    Array Chunk Bytes 31.17 GiB 12.87 MiB Shape (30593416, 1094) (300000, 45) Dask graph 2550 chunks in 8 graph layers Data type int8 numpy.ndarray

  • call_MQ
    (variants, samples)
    float32
    dask.array<chunksize=(300000, 45), meta=np.ndarray>

    Array Chunk Bytes 124.68 GiB 51.50 MiB Shape (30593416, 1094) (300000, 45) Dask graph 2550 chunks in 8 graph layers Data type float32 numpy.ndarray

  • call_AD
    (variants, samples, alleles)
    int16
    dask.array<chunksize=(300000, 45, 4), meta=np.ndarray>

    Array Chunk Bytes 249.36 GiB 103.00 MiB Shape (30593416, 1094, 4) (300000, 45, 4) Dask graph 2550 chunks in 8 graph layers Data type int16 numpy.ndarray

  • call_genotype_mask
    (variants, samples, ploidy)
    bool
    dask.array<chunksize=(300000, 45, 2), meta=np.ndarray>

    Array Chunk Bytes 62.34 GiB 25.75 MiB Shape (30593416, 1094, 2) (300000, 45, 2) Dask graph 2550 chunks in 9 graph layers Data type bool numpy.ndarray

• Indexes: (0)
• Attributes: (1)

  contigs :

  ('2', '3', 'X')

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. E.g., the code below sets up a genotype array.

# use the scikit-allel wrapper class for genotype calls
gt = allel.GenotypeDaskArray(ds_snps["call_genotype"].data)
gt

<GenotypeDaskArray shape=(30593416, 1094, 2) dtype=int8>

	0	1	2	3	4	...	1089	1090	1091	1092	1093
0	0/0	0/0	0/0	0/0	0/0	...	0/0	0/0	0/0	0/0	0/0
1	0/0	0/0	0/0	0/0	0/0	...	0/0	0/0	0/0	0/0	0/0
2	0/0	0/0	0/0	0/0	0/0	...	0/0	0/0	0/0	0/0	0/0
...	...
30593413	0/0	0/0	0/0	0/0	0/0	...	0/0	0/0	0/0	0/0	0/0
30593414	0/0	0/0	0/0	0/0	0/0	...	0/0	0/0	0/0	0/0	0/0
30593415	0/0	0/0	0/0	0/0	0/0	...	0/0	0/0	0/0	0/0	0/0

Example computation

Here’s an example computation to count the number of segregating SNPs on the longest contig (???) that also pass site filters. This may take a minute or two, because it is scanning genotype calls at millions of SNPs in hundreds of samples.

# choose contig (longest contig)
region = "2"
# choose site filter mask

# choose sample sets
sample_sets = ["1357-VO-BR-SALLUM-VMF00326"]

# access SNP calls
ds_snps = adar1.snp_calls(region=region, sample_sets=sample_sets)

# locate pass sites
loc_pass = ds_snps[f"variant_filter_pass_darlingi"].values

# perform an allele count over genotypes
gt = allel.GenotypeDaskArray(ds_snps["call_genotype"].data)
with ProgressBar():
    ac = gt.count_alleles(max_allele=3).compute()

# locate segregating sites
loc_seg = ac.is_segregating()

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

[########################################] | 100% Completed | 6.29 sms

22073253

Running larger computations

Please note that free cloud computing services such as Google Colab provide only limited computing resources. Thus although these services are able to efficiently read Adar1 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.