Pf7 Data Access
Contents
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
|
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
|
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
0 | 1 | 2 | 3 | 4 | ... | 154 | 155 | 156 | 157 | 158 | ||
---|---|---|---|---|---|---|---|---|---|---|---|---|
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 | |
2 | 0/0 | 0/0 | 0/0 | 0/0 | 0/0 | ... | 0/0 | 0/0 | 0/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
|
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')
|
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'])
|
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'])
|
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')
|