Pf7 Data Access

This page provides information about how to access data from Plasmodium falciparum version 7 (Pf7) project via Google Cloud. This includes sample metadata and single nucleotide polymorphism (SNP) calls. This release spans multiple MalariaGEN projects including Pf community project, GenRe-Mekong and SpotMalaria, and a collaboration between 82 studies spread around the globe.

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 (shaped like a rocket) at the top of the page and select one of the cloud computing services available.

For a quick overview on the spatial and geographical distribution of the samples available, please visit our Pf7 web-app.

Setup

Running this notebook requires some Python packages to be installed. These packages can be installed via pip or conda. E.g.:

!pip install -q --no-warn-conflicts malariagen_data

To make accessing these data more convenient, we’ve created the 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.

Now import the packages we’ll need to use here.

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

To access the pf7 data stored on google cloud use the following code:

pf7 = malariagen_data.Pf7()

Metadata

Data on the samples that were sequenced as part of this resource are available. It includes the time and place of collection, quality metrics, and accession numbers.

To see all the information available, load sample metadata into a pandas dataframe:

pf7_metadata = pf7.sample_metadata()

pf7_metadata.head()
Sample Study Country Admin level 1 Country latitude Country longitude Admin level 1 latitude Admin level 1 longitude Year ENA All samples same case Population % callable QC pass Exclusion reason Sample type Sample was in Pf6
0 FP0008-C 1147-PF-MR-CONWAY Mauritania Hodh el Gharbi 20.265149 -10.337093 16.565426 -9.832345 2014.0 ERR1081237 FP0008-C AF-W 82.16 True Analysis_set gDNA True
1 FP0009-C 1147-PF-MR-CONWAY Mauritania Hodh el Gharbi 20.265149 -10.337093 16.565426 -9.832345 2014.0 ERR1081238 FP0009-C AF-W 88.85 True Analysis_set gDNA True
2 FP0010-CW 1147-PF-MR-CONWAY Mauritania Hodh el Gharbi 20.265149 -10.337093 16.565426 -9.832345 2014.0 ERR2889621 FP0010-CW AF-W 86.46 True Analysis_set sWGA False
3 FP0011-CW 1147-PF-MR-CONWAY Mauritania Hodh el Gharbi 20.265149 -10.337093 16.565426 -9.832345 2014.0 ERR2889624 FP0011-CW AF-W 86.35 True Analysis_set sWGA False
4 FP0012-CW 1147-PF-MR-CONWAY Mauritania Hodh el Gharbi 20.265149 -10.337093 16.565426 -9.832345 2014.0 ERR2889627 FP0012-CW AF-W 89.74 True Analysis_set sWGA False
print("The data set has {} samples and {} fields".format(pf7_metadata.shape[0],pf7_metadata.shape[1]))
The data set has 20864 samples and 17 fields

We can explore each of the fields:

  • The Sample column gives the unique sample identifier used throughout all Pf7 analyses.

  • The Study refers to the partner study which collected the sample.

  • The Country & Admin level 1 describe the location where the sample was collected from.

  • The Country latitude, Country longitude, Admin level 1 latitude and Admin 1 longitude contain the GADM coordinates for each country & administrative level 1.

  • The Year column gives the time of sample collection.

  • The ENA column gives the run accession(s) for the sequencing read data for each sample.

  • The All samples same case column identifies samples set collected from the same individual.

  • The Population column gives the population to which the sample has been assigned. The possible values are: Africa - West (AF-W), Africa-Central (AF-C), Africa - East (AF-E), Africa - Northeast (AF-NE), Asia - South - East (AS-S-E), Asia - South – Far East (AS-S-FE), Asia - Southeast - West (AS-SE-W), Asia - Southeast - East (AS-SE-E), Oceania - New Guinea (OC-NG), South America (SA).

  • The % callable column refers to the % of the genome with coverage of at least 5 reads and less than 10% of reads with mapping quality 0.

  • The QC pass column defines whether the sample passed (True) or failed (False) QC.

  • The Exclusion reason describes the reason why the particular sample was excluded from the main analysis.

  • The Sample type column gives details on the DNA preparation method used

  • The Sample was in Pf6 column defines whether the sample was included in the previous version of the data release (Pf6) or if it is new to Pf7.

The python package Pandas can be used to explore and query the sample metadata in different ways. For example, here is a summary of the numbers of samples grouped by the country they were collected in:

pf7_metadata.groupby("Country").size()
Country
Bangladesh                          1658
Benin                                334
Burkina Faso                          58
Cambodia                            1723
Cameroon                             294
Colombia                             159
Côte d'Ivoire                         71
Democratic Republic of the Congo     573
Ethiopia                              34
Gabon                                 59
Gambia                              1247
Ghana                               4090
Guinea                               199
India                                316
Indonesia                            133
Kenya                                726
Laos                                1052
Madagascar                            25
Malawi                               371
Mali                                1804
Mauritania                           104
Mozambique                            91
Myanmar                             1260
Nigeria                              140
Papua New Guinea                     251
Peru                                  21
Senegal                              155
Sudan                                203
Tanzania                             697
Thailand                            1106
Uganda                                15
Venezuela                              2
Vietnam                             1733
dtype: int64

Variant Calls

These files contain details of 10,145,661 discovered variant genome positions. These variants were discovered amongst all samples from the release.

4,397,801 of these variant positions are SNPs, with the remainder being either short insertion/deletions (indels), or a combination of SNPs and indels.

Data on variant calls, including the genomic positions, alleles, and genotypes, can be accessed as an xarray Dataset:

variant_dataset = pf7.variant_calls()
variant_dataset
<xarray.Dataset>
Dimensions:              (variants: 10145661, alleles: 7, samples: 20864,
                          ploidy: 2)
Coordinates:
    variant_position     (variants) int32 dask.array<chunksize=(8388608,), meta=np.ndarray>
    variant_chrom        (variants) object dask.array<chunksize=(4194304,), meta=np.ndarray>
    sample_id            (samples) object dask.array<chunksize=(20864,), meta=np.ndarray>
Dimensions without coordinates: variants, alleles, samples, ploidy
Data variables:
    variant_allele       (variants, alleles) object dask.array<chunksize=(699051, 1), meta=np.ndarray>
    variant_filter_pass  (variants) bool dask.array<chunksize=(10145661,), meta=np.ndarray>
    variant_is_snp       (variants) bool dask.array<chunksize=(10145661,), meta=np.ndarray>
    variant_numalt       (variants) int32 dask.array<chunksize=(8388608,), meta=np.ndarray>
    variant_CDS          (variants) bool dask.array<chunksize=(10145661,), meta=np.ndarray>
    call_genotype        (variants, samples, ploidy) int8 dask.array<chunksize=(167773, 100, 2), meta=np.ndarray>
    call_AD              (variants, samples, alleles) int16 dask.array<chunksize=(23968, 100, 7), meta=np.ndarray>

The default returns a basic set of data most commonly used for data analysis. However, for more complex analysis the full range of variables available in the zarr can be accessed by setting the extended flag to True, as shown below:

extended_variant_dataset = pf7.variant_calls(extended=True)
extended_variant_dataset
<xarray.Dataset>
Dimensions:                                        (variants: 10145661,
                                                    alleles: 7, samples: 20864,
                                                    ploidy: 2, genotypes: 3,
                                                    sb_statistics: 4,
                                                    alt_alleles: 6)
Coordinates:
    variant_position                               (variants) int32 dask.array<chunksize=(8388608,), meta=np.ndarray>
    variant_chrom                                  (variants) object dask.array<chunksize=(4194304,), meta=np.ndarray>
    sample_id                                      (samples) object dask.array<chunksize=(20864,), meta=np.ndarray>
Dimensions without coordinates: variants, alleles, samples, ploidy, genotypes,
                                sb_statistics, alt_alleles
Data variables: (12/92)
    variant_allele                                 (variants, alleles) object dask.array<chunksize=(699051, 1), meta=np.ndarray>
    variant_filter_pass                            (variants) bool dask.array<chunksize=(10145661,), meta=np.ndarray>
    variant_is_snp                                 (variants) bool dask.array<chunksize=(10145661,), meta=np.ndarray>
    variant_numalt                                 (variants) int32 dask.array<chunksize=(8388608,), meta=np.ndarray>
    variant_CDS                                    (variants) bool dask.array<chunksize=(10145661,), meta=np.ndarray>
    call_genotype                                  (variants, samples, ploidy) int8 dask.array<chunksize=(167773, 100, 2), meta=np.ndarray>
    ...                                             ...
    variant_RegionType                             (variants) object dask.array<chunksize=(4194304,), meta=np.ndarray>
    variant_SOR                                    (variants) float32 dask.array<chunksize=(8388608,), meta=np.ndarray>
    variant_VQSLOD                                 (variants) float32 dask.array<chunksize=(8388608,), meta=np.ndarray>
    variant_altlen                                 (variants, alt_alleles) int32 dask.array<chunksize=(1398102, 6), meta=np.ndarray>
    variant_culprit                                (variants) object dask.array<chunksize=(4194304,), meta=np.ndarray>
    variant_set                                    (variants) object dask.array<chunksize=(4194304,), meta=np.ndarray>

Each of the elements in this xarray dataset is a dask array. The individual dask arrays can be accessed as follows, replacing the string with the variable you are looking for:

pos = variant_dataset["variant_position"].data
pos
Array Chunk
Bytes 38.70 MiB 32.00 MiB
Shape (10145661,) (8388608,)
Dask graph 2 chunks in 1 graph layer
Data type int32 numpy.ndarray
10145661 1

Genotypes

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,

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

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

gt = variant_dataset["call_genotype"].data
gt
Array Chunk
Bytes 394.28 GiB 32.00 MiB
Shape (10145661, 20864, 2) (167773, 100, 2)
Dask graph 12749 chunks in 1 graph layer
Data type int8 numpy.ndarray
2 20864 10145661

Note that the columns of this array (second dimension) match the rows in the sample metadata. You can use this correspondance to apply further subsetting operations to the genotypes by querying the sample metadata. E.g.:

loc_colombia = pf7_metadata.eval("Country == 'Colombia'").values
print(f"found {np.count_nonzero(loc_colombia)} samples from Colombia")
variant_dataset_colombia = variant_dataset.isel(samples=loc_colombia)
variant_dataset_colombia
found 159 samples from Colombia
<xarray.Dataset>
Dimensions:              (variants: 10145661, alleles: 7, samples: 159,
                          ploidy: 2)
Coordinates:
    variant_position     (variants) int32 dask.array<chunksize=(8388608,), meta=np.ndarray>
    variant_chrom        (variants) object dask.array<chunksize=(4194304,), meta=np.ndarray>
    sample_id            (samples) object dask.array<chunksize=(159,), meta=np.ndarray>
Dimensions without coordinates: variants, alleles, samples, ploidy
Data variables:
    variant_allele       (variants, alleles) object dask.array<chunksize=(699051, 1), meta=np.ndarray>
    variant_filter_pass  (variants) bool dask.array<chunksize=(10145661,), meta=np.ndarray>
    variant_is_snp       (variants) bool dask.array<chunksize=(10145661,), meta=np.ndarray>
    variant_numalt       (variants) int32 dask.array<chunksize=(8388608,), meta=np.ndarray>
    variant_CDS          (variants) bool dask.array<chunksize=(10145661,), meta=np.ndarray>
    call_genotype        (variants, samples, ploidy) int8 dask.array<chunksize=(167773, 32, 2), meta=np.ndarray>
    call_AD              (variants, samples, alleles) int16 dask.array<chunksize=(23968, 32, 7), meta=np.ndarray>

The data on genomic variants can be loaded into memory as numpy arrays as shown in the following example, where we read genotypes for the first 5 SNPs and the first 3 samples:

g = gt[:5, :3, :].compute()
g
array([[[-1, -1],
        [ 0,  0],
        [-1, -1]],

       [[-1, -1],
        [ 0,  0],
        [-1, -1]],

       [[-1, -1],
        [ 0,  0],
        [-1, -1]],

       [[-1, -1],
        [ 0,  0],
        [-1, -1]],

       [[-1, -1],
        [ 0,  0],
        [-1, -1]]], 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 using the Colombian samples subset we created above.

# use the scikit-allel wrapper class for genotype calls
gt = allel.GenotypeDaskArray(variant_dataset_colombia["call_genotype"].data)
gt
<GenotypeDaskArray shape=(10145661, 159, 2) dtype=int8>
01234...154155156157158
00/00/0./.0/00/0...0/00/00/0./.0/0
10/00/0./.0/00/0...0/00/00/0./.0/0
20/00/00/00/00/0...0/00/00/0./.0/0
......
10145658./../../../../...../../../../../.
10145659./../../../../...../../../../../.
10145660./../../../../...../../../../../.

Genome Annotations

Gene annotations provide information on which regions of the genome contain DNA sequences that encode genes, which are transcribed and spliced into messenger RNA (mRNA) and then translated to make proteins.

For convenience, we’ve added some functionality to the malariagen_data package for loading these gene annotations into a pandas data frame as shown below:

genome_features = pf7.genome_features()
genome_features
contig source type start end score strand phase ID Parent Name alias
0 Pf3D7_01_v3 chado repeat_region 1 360 NaN + NaN Pfalciparum_REP_20 NaN NaN NaN
1 Pf3D7_01_v3 chado repeat_region 361 1418 NaN + NaN Pfalciparum_REP_15 NaN NaN NaN
2 Pf3D7_01_v3 chado repeat_region 2160 3858 NaN + NaN Pfalciparum_REP_35 NaN NaN NaN
3 Pf3D7_01_v3 chado repeat_region 8856 9021 NaN + NaN Pfalciparum_REP_5 NaN NaN NaN
4 Pf3D7_01_v3 chado repeat_region 9313 9529 NaN + NaN Pfalciparum_REP_25 NaN NaN NaN
... ... ... ... ... ... ... ... ... ... ... ... ...
40708 Pf_M76611 chado rRNA 5772 5854 NaN - NaN PF3D7_MIT04100.1 PF3D7_MIT04100 NaN NaN
40709 Pf_M76611 chado CDS 5772 5854 NaN - 0.0 PF3D7_MIT04100.1:exon:1 PF3D7_MIT04100.1 NaN NaN
40710 Pf_M76611 chado gene 5861 5954 NaN - NaN PF3D7_MIT04200 NaN RNA8 mal_rna_19
40711 Pf_M76611 chado rRNA 5861 5954 NaN - NaN PF3D7_MIT04200.1 PF3D7_MIT04200 NaN NaN
40712 Pf_M76611 chado CDS 5861 5954 NaN - 0.0 PF3D7_MIT04200.1:exon:1 PF3D7_MIT04200.1 NaN NaN

40713 rows × 12 columns

The above loads a default set of attributes "ID", "Parent", "Name", "alias". To access all features set attributes to "*".

pf7.genome_features(attributes="*")
contig source type start end score strand phase Dbxref Derives_from ... paralogous_to polypeptide_domain previous_systematic_id product product_synonym signal_peptide stop_codon_redefined_as_selenocysteine synonym translation transmembrane_polypeptide_region
0 Pf3D7_01_v3 chado repeat_region 1 360 NaN + NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
1 Pf3D7_01_v3 chado repeat_region 361 1418 NaN + NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2 Pf3D7_01_v3 chado repeat_region 2160 3858 NaN + NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
3 Pf3D7_01_v3 chado repeat_region 8856 9021 NaN + NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
4 Pf3D7_01_v3 chado repeat_region 9313 9529 NaN + NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
40708 Pf_M76611 chado rRNA 5772 5854 NaN - NaN NaN NaN ... NaN NaN NaN term=large subunit ribosomal RNA fragment D;db... NaN NaN NaN NaN NaN NaN
40709 Pf_M76611 chado CDS 5772 5854 NaN - 0.0 NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
40710 Pf_M76611 chado gene 5861 5954 NaN - NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
40711 Pf_M76611 chado rRNA 5861 5954 NaN - NaN NaN NaN ... NaN NaN NaN term=ribosomal RNA fragment RNA8;db_xref=PMID:... NaN NaN NaN NaN NaN NaN
40712 Pf_M76611 chado CDS 5861 5954 NaN - 0.0 NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

40713 rows × 35 columns

Or to get a specific set of attributes specify them in a list

pf7.genome_features(attributes=['alias','comment','product'])
contig source type start end score strand phase alias comment product
0 Pf3D7_01_v3 chado repeat_region 1 360 NaN + NaN NaN telomeric repeat region NaN
1 Pf3D7_01_v3 chado repeat_region 361 1418 NaN + NaN NaN 14bp repeat NaN
2 Pf3D7_01_v3 chado repeat_region 2160 3858 NaN + NaN NaN 65bp repeat NaN
3 Pf3D7_01_v3 chado repeat_region 8856 9021 NaN + NaN NaN 25bp repeat NaN
4 Pf3D7_01_v3 chado repeat_region 9313 9529 NaN + NaN NaN 26bp repeat NaN
... ... ... ... ... ... ... ... ... ... ... ...
40708 Pf_M76611 chado rRNA 5772 5854 NaN - NaN NaN fragment number 7 transcript contains a 3' oli... term=large subunit ribosomal RNA fragment D;db...
40709 Pf_M76611 chado CDS 5772 5854 NaN - 0.0 NaN NaN NaN
40710 Pf_M76611 chado gene 5861 5954 NaN - NaN mal_rna_19 NaN NaN
40711 Pf_M76611 chado rRNA 5861 5954 NaN - NaN NaN transcript contains a 3' oligo(A) structure fr... term=ribosomal RNA fragment RNA8;db_xref=PMID:...
40712 Pf_M76611 chado CDS 5861 5954 NaN - 0.0 NaN NaN NaN

40713 rows × 11 columns

Genome Reference

We mapped sequence reads for all samples against the P. falciparum 3D7 v3 reference genome.

For convenience, the reference genome sequence can be loaded as a dask array, e.g.:

ref = pf7.genome_sequence()
ref
Array Chunk
Bytes 22.25 MiB 412.03 kiB
Shape (23332839,) (421914,)
Dask graph 72 chunks in 17 graph layers
Data type |S1 numpy.ndarray
23332839 1

This can be loaded as a numpy array using the following

ref.compute()
array([b't', b'g', b'a', ..., b'a', b't', b'a'], dtype='|S1')

The reference can also be subset by contig.

The set of contigs used can be accessed as follows:

pf7.contigs
['Pf3D7_01_v3',
 'Pf3D7_02_v3',
 'Pf3D7_03_v3',
 'Pf3D7_04_v3',
 'Pf3D7_05_v3',
 'Pf3D7_06_v3',
 'Pf3D7_07_v3',
 'Pf3D7_08_v3',
 'Pf3D7_09_v3',
 'Pf3D7_10_v3',
 'Pf3D7_11_v3',
 'Pf3D7_12_v3',
 'Pf3D7_13_v3',
 'Pf3D7_14_v3',
 'Pf3D7_API_v3',
 'Pf_M76611']

To load a single contig

pf7.genome_sequence(region='Pf3D7_01_v3')
Array Chunk
Bytes 625.83 kiB 312.92 kiB
Shape (640851,) (320426,)
Dask graph 2 chunks in 1 graph layer
Data type |S1 numpy.ndarray
640851 1

To load multiple contigs specify them in a list. The data will be concatenated.

pf7.genome_sequence(region=['Pf3D7_07_v3','Pf3D7_02_v3','Pf3D7_03_v3'])
Array Chunk
Bytes 3.30 MiB 352.83 kiB
Shape (3460280,) (361302,)
Dask graph 12 chunks in 4 graph layers
Data type |S1 numpy.ndarray
3460280 1

You can also specify a specific region of the contig.

pf7.genome_sequence(region=['Pf3D7_07_v3','Pf3D7_02_v3:15-20','Pf3D7_03_v3:40-50'])
Array Chunk
Bytes 1.38 MiB 352.83 kiB
Shape (1445224,) (361302,)
Dask graph 6 chunks in 6 graph layers
Data type |S1 numpy.ndarray
1445224 1

If you know the gene name you would like to access, but aren’t sure what the ID would be you can access this through the annotations. Below is an example for CRT.

gene_name = str(genome_features.loc[genome_features.Name == 'CRT'].ID.values)
print(gene_name)
['PF3D7_0709000']

You can then enter this as the region

pf7.genome_sequence(region='PF3D7_0709000')
Array Chunk
Bytes 3.86 kiB 3.86 kiB
Shape (3957,) (3957,)
Dask graph 1 chunks in 2 graph layers
Data type |S1 numpy.ndarray
3957 1