Pf7: an open dataset of Plasmodium falciparum genome variation in 20,000 worldwide samples

1Despite global malaria eradication efforts in the mid-20th century and more recent advances in malaria control, Plasmodium falciparum remains endemic throughout Africa, Asia, South America and Oceania. According to the most recent World Malaria Report, each year over 200 million people suffer from malaria due to P. falciparum and over 600,000 die as a result. Most of the disease burden falls on Africa, and particularly African children. There is international commitment to control malaria more effectively and many countries are working towards the long-term goal of malaria elimination. However the parasites are continually evolving to resist antimalarial drugs and to evade host immunity, and this is a major challenge to sustainable malaria control and elimination.

2Our understanding of the evolutionary biology and population genomics of malaria parasites has advanced considerably over the past decade. There is now a substantial body of literature on the genomic diversity and global population structure of malaria parasites, on within-host genetic variation and what this tells us about superinfection and cotransmission, on the identification and monitoring of parasite drug resistance loci, on genetic variation in malaria vaccine antigens, and on methods of analysing genetic relatedness to understand patterns of malaria transmission. A useful summary of the current state of the field can be found in the recent review by Neafsey et al..

3This rapidly growing area of research is underpinned by open data on genome sequence variation in natural parasite populations. Here we report a new release of curated open data on Plasmodium falciparum genome variation from the MalariaGEN network. It includes samples featured in previous MalariaGEN data releases including the Pf3k Project, the Plasmodium falciparum Community Project and the GenRe Mekong Project. To avoid confusion between different MalariaGEN datasets we now identify each by a version number. In this new nomenclature, the previous version is called Pf6, and the version described here is called Pf7.

4Whole genome sequencing of all the samples in the Pf7 dataset was performed at the Wellcome Sanger Institute and a standardised analysis pipeline was used for variant discovery and genotyping. The Pf7 analysis pipeline was broadly similar to that used for the Pf6 dataset with some improvements that are described in more detail below. Sequence data and genotype calls were returned to partners for use in their own analyses and publications in line with MalariaGEN’s guiding principles on equitable data sharing.

5The Pf7 dataset comprises 20,864 samples of P. falciparum collected by 82 partner studies from 33 countries in Africa, Asia, South America and Oceania between 1984 and 2018 ( Table 1, Supplementary Tables 1 and 2, Supplementary Figure 1). Compared to the Pf6 dataset, this includes 13,752 new samples, 33 new partner studies and 5 additional countries. The majority of new samples (12,146) were collected since 2014, but there were also 379 samples collected prior to 2000 (Supplementary Figure 2).

6The most significant technical advance in the Pf7 dataset is that most of the new samples (12,891/13,752, 94%) came from the MalariaGEN SpotMalaria Project. SpotMalaria was designed to simplify and standardise the process of collecting dried blood spot (DBS) samples for parasite genetic analysis. The SpotMalaria protocol ensures that the vast majority of DBS samples are suitable for targeted genotyping of drug resistance loci, e.g. by amplicon sequencing, and that a significant proportion of samples are also suitable for whole genome sequencing. This requires an intermediate step known as selective whole genome amplification (sWGA), which makes it possible to obtain parasite genome sequence data from a very small sample, but at the cost of introducing considerable variability in sequencing coverage across the genome. This is the first study to analyse a large number of sWGA P. falciparum genomes and therefore it was important to establish that sWGA was not introducing significant biases, and to adapt our methods for calling structural variations which are particularly sensitive to artefactual variation in sequencing coverage.

7We have performed a number of analyses to make the Pf7 dataset as useful as possible for a broad range of users. These include descriptions of global population structure, geographic patterns of drug resistance, haplotypic analysis of drug resistance loci, hrp2 and hrp3 deletions that can cause failure of rapid diagnostic tests, and variation in the C-terminal of the csp antigen used in the most advanced malaria vaccines. These analyses are not intended to be comprehensive and technical users of the dataset can download the analysis-ready dataset for more specialised or detailed investigations.

8A high level view of the Pf7 dataset can be obtained from the data exploration tool at the MalariaGEN website. This shows the locations and years where the samples were collected, and the genotype-inferred drug resistance status of each sample. Most importantly it names the investigators who led the studies that contributed the samples at each location and thus made this global dataset possible.

9We used the Illumina platform to produce whole genome sequencing data on all samples and mapped sequence reads against the P. falciparum 3D7 v3 reference genome. The median depth of coverage was 107 sequence reads averaged across the whole genome and across all samples. We used an analysis pipeline for variant discovery and genotyping analogous to that used in Pf6, as outlined in the Methods section.

10In the first stage of analysis we discovered genomic variations in nearly half of the 23Mb P. falciparum positions (10,145,661 in total, Supplementary Table 3), including 4,397,801 single nucleotide polymorphisms (SNPs).

11For the analysis reported here, we excluded all variants in subtelomeric and internal hypervariable regions, mitochondrial and apicoplast genomes and applied stringent quality filters to the remaining variants as described in the Methods section. A total of 3,125,721 SNPs (of which 2,513,888 were biallelic) and 2,742,938 non-SNPs, i.e. short indels or SNP-indel combinations, passed all these filters. Some of the variant positions that were classified as SNPs in Pf6 are now classified as non-SNPs because they additionally include indel alleles.

12We performed quality control checks to remove samples with: (i) unverified or incomplete sample collection information; (ii) evidence of co-infection with other Plasmodium species; (iii) more than one technical replicate or time course sampling; (iv) low coverage; (v) a higher than expected number of singleton SNPs. In total, we retained 16,203 high-quality samples ( Table 1).

13This analysis-ready dataset with details of all participating partner studies and a python package providing convenience methods for accessing is available here.

14Unlike the previous version, nearly all samples that are new to this release (12,891/13,752, 94%) have been sequenced after undergoing selective whole genome amplification (sWGA). This process allows us to sequence samples collected as dried blood spots, which greatly simplifies many of the operational challenges in collecting venous blood.

15An artefact introduced by sWGA is high variability in coverage across the genome. This impacts on the use of local variation in genomic coverage as a way to identify large structural variations such as tandem duplications. We therefore developed a novel method based on GATK GermlineCNVCaller (gCNV) for typing duplications around mdr1 and plasmepsin 2–3 (associated with resistance to mefloquine and piperaquine, respectively) and deletions of hrp2 and hrp3 (associated with rapid diagnostic test failures). We started by compiling a list of observed breakpoints in and around the loci of interest. We then leveraged on the fact that the amplification bias introduced by sWGA, and the consequent variation in coverage, is relatively systematic and can be used for a cross-sample normalisation. Finally, we complemented the results with an analysis to detect presence of face-away reads around the known breakpoints and obtained a final set of calls. For plasmepsin 2–3 duplications, concordance between gCNV and the face-away reads methods was high, with 99% of samples called as duplication by gCNV also being called as duplication by the face-away method, and the remaining 1% all called as missing. For mdr1, concordance was significantly lower, with 19% of samples called as duplication by gCNV being called as no duplication by the face-away method. This could be explained by the fact that the set of breakpoints used is likely not exhaustive, and also by some duplications not being tandem duplications. For samples called as no duplication by gCNV, the vast majority were also called as no duplication by the face-away method (83%) or else missing (17%). For hrp2 and hrp3 deletions, we manually validated the results and identified evidence of breakpoints for all the deletion calls.

16To ensure that sWGA is not introducing biases in population structure, we analysed four sets of samples from the same location and time periods for which we had a substantial number of both sWGA and non-sWGA samples, and could detect no apparent stratification (Supplementary Figure 3).

17We grouped samples by location using the classification scheme known as first-level administrative division: we refer to these as sampling locations. Based on principal coordinate and neighbour-joining tree analyses of all samples, we identified ten major divisions of population structure: we refer to these as major sub-populations. We then determined the geographical range of each major sub-population by examining the sampling locations that it contained ( Figure 1 and Supplementary Figure 4). We identified four major sub-populations in Africa: AF-W (8,610 samples from western Africa), AF-C (573 samples from Kinshasa, DRC), AF-NE (316 samples from Sudan, Ethiopia, Uganda, and Kisumu county in western Kenya), AF-E (1,846 samples from east Africa). In Asia we identified four major sub-populations: two in South Asia (AS-S-E, 244 samples from the Indian states of Odisha and West Bengal, and AS-S-FE, 1,730 samples from Bangladesh and the far-eastern Indian state of Tripura), and two in Southeast Asia (AS-SE-W, 2,254 samples from the part of Southeast Asia west of Bangkok, and AS-SE-E, 4,620 samples from the eastern part). The two remaining major sub-populations were OC-NG (384 samples from Oceanian island of New Guinea) and SA (182 samples from South America).

18This geographical assignment of ten major sub-populations is a somewhat crude approximation of the underlying population structure, and it does not reflect international conventions for grouping countries or regions. However it provides a framework that allows a broad comparison of population genetic features between different parts of the world, such as the rate of decay of linkage disequilibrium (Supplementary Figure 5), nucleotide diversity (Supplementary Figure 6, lower panel) and complexity of infections (Supplementary Figure 6, upper panel). We also examined the fixation index between sub-populations (Supplementary Figure 7).

19We classified samples as resistant or sensitive to major antimalarials and combinations based on genotyping of known drug resistance alleles ( Table 2 - see here for details of the heuristics used). At a regional level, the frequency of samples classified as resistant to each drug is broadly consistent with known and previously reported geographical patterns, with the highest prevalence of multidrug resistance observed in Southeast Asia. Interestingly, in South Asia, we find that the frequency of resistance to chloroquine, sulfadoxine and pyrimethamine appears to be much higher in the far-eastern sub-population (Bangladesh and Tripura) than in the eastern sub-population (Odisha and West Bengal).

20Care is required when interpreting these findings, as most of the major sub-populations spanned a large geographic region, within which there could be considerable epidemiological diversity, and also because we aggregate samples that were collected over relatively long periods of time during which patterns of resistance may have changed (Supplementary Figure 8). To take West Africa as an example, if we consider samples collected between 2013 and 2016 ( Figure 2), we observed levels of chloroquine resistance varying from 0% in Volta, Ghana to 100% in Atlantique, Benin.

21Amplifications of the genes mdr1 and plasmepsin 2-3 are markers of resistance to mefloquine and piperaquine, respectively. Interestingly, two samples collected in 1993 in Cambodia have tandem duplications of both genes, an event which is relatively rare in more recent samples (only 21 samples in total out of 1,959 that have evidence of amplification of either gene). In addition to presence/absence of the amplification, we also provide details of the associated breakpoints for all samples which shows these two samples having two distinct breakpoints in plasmepsin, one of which is identical to that most commonly found in contemporary samples.

22Previous reports have shown that the current wave of multidrug-resistant P. falciparum in Southeast Asia is driven by the KEL1 lineage of the kelch13 artemisinin resistance locus – and is associated with multiple new mutations in the crt resistance locus. This dataset confirms the dramatic increase of KEL1 in Cambodia, Eastern Thailand, Laos and Vietnam that has occurred over the past ten years (Supplementary Figure 9). Analysis of the crt locus in samples with the K76T resistance variant reveals a major cluster of haplotypes on a common genetic background, the one observed in a widely used lab strain isolated in Asia in 1980 and commonly referred to as Dd2 ( Figure 3A, Supplementary Table 5). In addition to the original Dd2 haplotype, we observe 31 additional mutations. These are essentially restricted to eastern SE Asia with only four samples from outside this region, and they include mutations that have previously been associated with piperaquine resistance,. They are seen across all countries of eastern SE Asia and have risen rapidly in frequency leading us to consider them newly emerging Dd2-background haplotypes ( Figure 3B). Most have a single mutation on a Dd2 background, but we observe 13 haplotypes with two mutations and one haplotype (in a single sample) with three mutations (Supplementary Table 5). These findings highlight the value in retrospective analysis of drug resistance mutations, as most of these samples were collected and sequenced before the relevance of crt mutations to piperaquine resistance was appreciated.

23In addition to the selective pressures due to antimalarial drugs highlighted in the previous section, another area of interest is selective pressure due to vaccines. Having a baseline understanding of genetic variation in vaccine genes is likely to be valuable.

24The WHO has recently recommended the RTS,S vaccine for use in regions with moderate to high transmission which includes much of sub-Saharan Africa. The vaccine targets the gene csp and has a construct based on the 3D7 reference sequence of part of the central NANP repeat region where antibodies bind and the c-terminal region which contains T cell epitopes. Another vaccine based on the same region and sequence, R21, is also showing promise in early stage clinical trials. Vaccine efficacy is likely to depend on a number of factors, both host and parasite, and clinical trials show some variability between different locations in Africa. How similar the parasite is in the targeted region to the 3D7 sequence used in vaccine design could be a contributing factor to this variability. Genetic diversity in the construct region may or may not affect vaccine efficacy, and in order to understand this it will be important to monitor efficacy against diversity going forwards. Here we begin a systematic catalogue of population-level diversity in csp. While it is challenging if not impossible to genotype the central repeat region using short read data, we start here by looking at variation leading to amino acid changes in the c-terminal region of the protein.

25We identified all non-synonymous mutations in the csp c-terminal region and analysed the frequency of these in different populations. Interestingly, the vast majority of the samples across the globe carry non-reference alleles, i.e. different from the 3D7 sequence used in the vaccine design, at amino acids 301, 317, 318, 321 and 361 (Supplementary Figure 10). We found a total of 248 unique amino acid haplotypes of csp 277-397 out of 11,254 samples with no ambiguous calls. Amino acid haplotype sequences for the c-terminal region of csp for all samples can be found at the MalariaGEN website.

26Surprisingly, the most common haplotype in the dataset and the one with the second lowest number of differences from all other unique haplotypes is the one observed in lab strain Dd2, being found in 2,760/11,254 (25%) of the samples and having a mean of 4.7 differences from other haplotypes ( Figure 4a). In contrast, the 3D7 haplotype used for both csp-based vaccines is only found in 3% of samples and has on average 6.9 differences from all other haplotypes.

27Importantly, this striking difference also holds when examining each population separately, including in West Africa from where 3D7 is thought to have originated ( Figure 4b).

28Taken together, these results show perhaps surprising differences between the target haplotype used in the design of RTS,S and R21 and those circulating in natural parasite populations, and provide a systematic catalogue that can be used in future studies to elucidate any possible clinical significance of sequence diversity.

29Most widely used rapid diagnostic tests (RDTs) rely on detection of the products of the hrp2 and hrp3 genes, and deletions in these genes is known to lead to RDT failure. We used gCNV to call presence/absence of hrp2 and hrp3 deletions in 68% of QC passed samples. Frequencies of deletions vary greatly across countries, and deletions of hrp3 (1.9%) are more common than those of hrp2 (0.14%) (Supplementary Table 6). The only countries where we see deletions in both hrp2 and hrp3 that would cause HRP RDTs to fail are Peru (6/20 samples), Indonesia (2/115 samples) and Sudan (1/7 samples).

30There have been numerous reports of such deletions in recent years, but to date there has been little detail on the mechanisms causing such deletions. We manually inspected reads around the apparent breakpoints in order to classify the types of events driving these deletions. For hrp2, all deletion events can be explained by a process of telomere healing whereby the end of a chromosome is deleted and a telomere repeat sequence attached to the breakpoint, ( Figure 5). Telomere healing events can be determined with breakpoint precision and in almost all cases samples with the same breakpoints are from the same country (Supplementary Table 7). For hrp3, we also identified a number of telomere healing events, but also two other quite different types of event causing deletion of the gene ( Figure 5, Supplementary Table 7). In many cases a new hybrid chromosome appears to have been created by a recombination between chromosome 13 and 11 at a cluster of rRNA genes that have orthologous copies on both chromosomes. In other cases a recombination between chromosome 13 and an inverted section from within chromosome 5 containing the gene mdr1 can be identified. This remarkable event results in both deletion of hrp3 and duplication of mdr1.

31As with previous releases, in Pf7 we have created genotypes at SNPs genome-wide and CNVs in specific locations, but we intend to continue to expand the resource to consider other types of genetic variation. Various surface antigens, including vaccine candidate genes, have two distinctly different allelic forms. Often the two forms are so divergent that reads from a non-reference form will not map to the 3D7 reference genome, hence necessitating an alternative to the mapping-based genotype calling approach described above. As an example, the gene eba175 has two different allelic forms, known as the F- and C-types,. As a proof of concept for such a dimorphic gene, we used a novel kmer-based method to call these two types in each sample. We see that 7,380 (69%) samples have the F allele exclusively and 3,364 (31%) have the C allele. Although frequencies of each type vary between populations, we see >10% frequency of each type in all populations (Supplementary Figure 11). These results give weight to the argument that eba175 is under balancing selection, most likely negative frequency-dependent selection driven by interactions with the human immune system. Analyses of other such dimorphic genes will be left to future work, as will more detailed analysis of variation within these different allelic types.

32The Pf7 dataset increases the amount of curated open data for population genomic analysis of P. falciparum by almost threefold, and greatly increases the number of samples collected within the last five years. With denser geographical coverage it is possible to undertake higher resolution analysis of epidemiological variation within a region, e.g. we observe considerable heterogeneity of inferred chloroquine resistance in West Africa ( Figure 3), and it also allows us to identify new sub-populations with distinctive epidemiological features, e.g. we find two sub-populations in south Asia that have contrasting drug resistance profiles. There is useful historical information to be obtained from older samples that are included in this new data release, e.g. some samples collected in Cambodia in the early 1990s appear to be resistant to both piperaquine and mefloquine, which is highly relevant to the ongoing evolution of multidrug resistance to artemisinin combination therapy in Southeast Asia.

33An important technical advance is that Pf7 contains a large number of samples that were collected as dried blood spots in the field. We and others have previously described successful whole genome sequencing of P. falciparum from DBS after selective whole-genome amplification but it was unclear how well this methodology would perform at scale. Here we show that the SpotMalaria protocol for sWGA of DBS samples allows us to generate whole genome sequence data of sufficient quality to genotype the vast majority of SNPs with sufficient accuracy and reliability for large-scale population genomic analysis. We have introduced improvements to our pipelines for calling copy number variants, necessitated by the greatly increased heterogeneity of sequencing coverage following sWGA. There remain hypervariable gene families and other regions of the parasite genome that cannot be accurately genotyped in field samples using current methods, and these difficulties are compounded by sWGA, but by working on sequencing and analysis methods we aim to continually improve genome coverage in future releases.

34The knowledge that DBS samples can be used for whole genome analysis in large-scale studies is of practical importance, as it empowers field researchers and national malaria control programs to integrate population genomic information with other forms of epidemiological and public health data, and it paves the way to a global infrastructure for genomic surveillance of P. falciparum. Information about the processes and methods of the SpotMalaria Project can be obtained at the MalariaGEN website

35The Pf7 dataset includes a range of analyses and sample annotations that are intended to increase the utility of the data for researchers working on practical problems in malaria control. Compared to the Pf6 dataset, we have made improvements to methods for calling CNVs at the mdr1 and pm2 drug resistance loci and for calling hrp2 and hrp3 deletions that can affect rapid diagnostic tests. Other new analyses included in Pf7 include more detailed descriptions of: (a) hrp2 and hrp3 deletion breakpoints; (b) drug resistance locus haplotypes and in particular newly emerging crt haplotypes; (c) profiles of variation in the csp vaccine antigen and the vaccine candidate eba175. In future releases we aim to improve and expand analyses that are relevant to malaria control programmes and policymakers.

36The Pf7 dataset focuses entirely on genome sequencing data, but there is a growing body of data from amplicon sequencing and targeted genotyping approaches that is highly informative about multiple aspects of P. falciparum population genomics. For example, the GenRe-Mekong Project has used the SpotMalaria platform combined with amplicon sequencing to enable malaria control programmes in the Greater Mekong Region to conduct national genomic surveillance of multidrug resistance. In future data releases we aim to integrate data from these different sources to greatly increase sample size and geographical coverage, and thus improve the resolution of population genomic analysis.

37All samples in this study were derived from blood samples obtained from patients with P. falciparum malaria, collected with informed consent from the patient or a parent or guardian. At each location, sample collection was approved by the appropriate local and institutional ethics committees. The following local and institutional committees gave ethical approval for the partner studies: Walter and Eliza Hall Institute Human Research Ethics Committee, Australia; University of Antwerp, Belgium; Comite d'Ethique de la Recherche - Institut des Sciences Biomedicales Appliquees, Benin; Ministere de la Sante – Republique du Benin, Benin; Comité d'Éthique, Ministère de la Santé, Bobo-Dioulasso, Burkina Faso; Ministry of Health National Ethics Committee for Health Research, Cambodia; Institutional Review Board University of Buea, Cameroon; Comite Institucional de Etica de investigaciones en humanos de CIDEIM, Colombia; Research Ethics Committee of the Faculty of Medicine of the National University of Colombia; Comité National d'Ethique de la Recherche, Cote d’Ivoire; Comite d’Ethique Universite de Kinshasa, Democratic Republic of Congo; Armauer Hansen Research Institute Institutional Review Board, Ethiopia; Addis Ababa University, Aklilu Lemma Institute of Pathobiology Institutional Review Board, Ethiopia; Ghana Health Service Ethical Review Committee, Ghana; University of Ghana Noguchi Medical Research Institute, Ghana; Navrongo Health Research Centre Institutional Review Board, Ghana; Comite d’Ethique National Pour la Recherché en Santé, République de Guinée; KEMRI Scientific and Ethics Review Unit, Kenya; Ministry of Health National Ethics Committee For Health Research, Laos; College of Medicine, University of Lagos, Nigeria; Comité National d'Ethique auprès du Ministère de la Santé Publique, Madagascar; College of Medicine Regional Ethics Committee University of Malawi, Malawi; Faculté de Médecine, de Pharmacie et d'Odonto-Stomatologie, University of Bamako, Bamako, Mali; L'université des Sciences, des Techniques et des Technologies de Bamako, Mali; Ethics Committee of the Ministry of Health, Mali; Ethics committee of the Ministry of Health, Mauritania; National Bioethics Committees of Mozambique (CNBS); Institutional Review Board, Papua New Guinea Institute of Medical Research, Goroka, Papua New Guinea; PNG Medical Research Advisory Council (MRAC), Papua New Guinea; Institutional Review Board, Universidad Nacional de la Amazonia Peruana, Iquitos, Peru; Al Neelain University Institutional Review Board, Sudan; Federal Ministry of Health, Sudan; Ethics Committee of the Ministry of Health, Senegal; Medical Research Coordinating Committee of the National Institute for Medical Research, Tanzania; Ethics Committee, Faculty of Tropical Medicine, Mahidol University, Bangkok, Thailand; Gambia Government/MRC Joint Ethics Committee, Banjul, The Gambia; Liverpool School of Tropical Medicine, UK; London School of Hygiene and Tropical Medicine Ethics Committee, London, UK; Oxford Tropical Research Ethics Committee, Oxford, UK; University College London Hospitals Research Ethics Committee, UK; Walter Reed Army Institute of Research, USA; National Institute of Allergy and Infectious Diseases, Bethesda, MD, USA; University of Maryland School of Medicine IRB, USA; Ministry of Health Institute of Malariology-Parasitology-Entomology, Vietnam.

38Standard laboratory protocols were used to determine DNA quantity and proportion of human DNA in each sample as previously described,.

39Here we summarise the bioinformatics methods used to produce and analyse the data; full details are available here.

40Reads mapping to the human reference genome were discarded before all analyses, and the remaining reads were mapped to the P. falciparum 3D7 v3 reference genome using bwa mem. “Improved” BAMs were created using the Samtools FixMateInformation, Picard MarkDuplicates, and GATK base quality score recalibration. All lanes for each sample were merged to create sample-level BAM files.

41Putative variants were called in each sample independently using GATK (v4.1.4.0) HaplotypeCaller, then all samples were combined to jointly genotype the entire cohort using GATK GenotypeGVCFs.

42SNPs and indels were filtered using GATK’s Variant Quality Score Recalibration (VQSR). Variants with a VQSLOD score ≤ 2 were filtered out. Functional annotations were applied using snpEff version 4.3. Genome regions were annotated using BCFtools v1.10.2 ( http://www.htslib.org/doc/bcftools.html) and masked if they were outside the core genome. Unless otherwise specified, we used biallelic SNPs that pass all quality filters for all the analysis.

43VCF files were converted to zarr v2.4.0 format and subsequent analyses were mainly performed using scikit-allel v1.2.1 and the zarr files.

44We identified species using nucleotide sequence from reads mapping to six different loci in the mitochondrial genome, using custom java code. The loci were located within the cox3 gene (PF3D7_MIT01400), as described in a previously published species detection method. Alleles at various mitochondrial positions within the six loci were genotyped and used for classification.

45We created a final analysis set of 16,203 samples after removing samples with unverified identity, mixed species, replicate and low coverage samples, and samples with excessive numbers of singleton SNPs.

46We calculate genetic distance between samples using biallelic coding SNPs that pass filters using a method previously described.

47The matrix of genetic distances was used to generate neighbour-joining trees and principal coordinates. Based on these observations we grouped the samples into ten major sub-populations: South America, West Africa, Central Africa, Northeast Africa, East Africa, an eastern part of South Asia, a far-eastern part of South Asia, the western part of Southeast Asia, the eastern part of Southeast Asia and Oceania, with samples assigned to region based on the geographic location of the sampling site.

48F WS was calculated using custom python scripts using the method previously described. Nucleotide diversity (π) was calculated in non-overlapping 25kbp genomic windows, only considering coding SNPs to reduce the ascertainment bias caused by poor accessibility of non-coding regions. LD decay ( r 2) was calculated using the method of Rogers and Huff and biallelic SNPs with low missingness and regional allele frequency >10%. Mean F ST between populations was calculated using Hudson’s method.

49To call duplication genotypes around mdr1 and plasmepsin 2-3 from binned coverage, we adapted the GATK GermlineCNVCaller (gCNV) pipeline, which features a probabilistic Bayesian model that jointly infers both copy-number activity and a model for denoising sequencing systematics. After breakpoint genotypes were called, we performed an initial, permissive round of annotation-based call filtering, using hard cuts to identify failing samples and demarcate duplication and reference genotypes. This was followed by a final round of curation, based on manual inspection of the denoised copy ratios, to discard spurious duplication calls. The resulting filtered gCNV call set was integrated with an analogous call set based on consideration of face-away read-pair evidence, in which we set the breakpoint to be that with the highest proportion of face-away reads.

50Deletions in hrp2 and hrp3, genes which are located in subtelomeric regions of the genome with very high levels of natural variation, were identified using the same breakpoint-genotyping framework introduced above. As before, an initial round of permissive, annotation-based filtering was performed, followed by a final round of curation to discard spurious deletion calls. We identified deletion breakpoints by manual inspection of custom plots.

51The procedure used to map genetic markers to inferred resistance status classification is described in detail for each drug in the accompanying data release. In brief, we called amino acids at selected loci by first determining the reference amino acids and then, for each sample, applying all variations using the GT field of the VCF file. This same approach was used to identify haplotypes of csp 277-397. The amino acid and copy number calls generated in drug resistance genes were used to classify all samples into different types of drug resistance. Our methods of classification were heuristic and based on the available data and current knowledge of the molecular mechanisms. Each type of resistance was considered to be either present, absent or unknown for a given sample. eba175 F- and C-type calls are made by identifying samples that have 19bp kmers present that are unique to the C and F haplotypes.