{'fraction_non_alphanumeric': 0.030770958504127126, 'fraction_numerical': 0.023770752615718645, 'mean_word_length': 4.907563025210084, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 7, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A novel locus of resistance to severe malaria in a region of ancient balancing selection

The high prevalence of sickle haemoglobin in Africa shows that malaria has been a major force for human evolutionary selection, but surprisingly few other polymorphisms have been proven to confer resistance to malaria in large epidemiological studies. To address this problem, we conducted a multi-centre genome-wide association study (GWAS) of life-threatening Plasmodium falciparum infection (severe malaria) in over 11,000 African children, with replication data in a further 14,000 individuals. Here we report a novel malaria resistance locus close to a cluster of genes encoding glycophorins that are receptors for erythrocyte invasion by P. falciparum. We identify a haplotype at this locus that provides 33% protection against severe malaria (odds ratio = 0.67, 95% confidence interval = 0.60-0.76, P value = 9.5 × 10(-11)) and is linked to polymorphisms that have previously been shown to have features of ancient balancing selection, on the basis of haplotype sharing between humans and chimpanzees. Taken together with previous observations on the malaria-protective role of blood group O, these data reveal that two of the strongest GWAS signals for severe malaria lie in or close to genes encoding the glycosylated surface coat of the erythrocyte cell membrane, both within regions of the genome where it appears that evolution has maintained diversity for millions of years. These findings provide new insights into the host-parasite interactions that are critical in determining the outcome of malaria infection.


A novel locus of resistance to severe malaria in a region of ancient balancing selection Malaria Genomic Epidemiology Network 1
The high prevalence of sickle haemoglobin in Africa shows that malaria has been a major force for human evolutionary selection, but surprisingly few other polymorphisms have been proven to confer resistance to malaria in large epidemiological studies [1][2][3] . To address this problem we conducted a multi-centre genome-wide association study (GWAS) of lifethreatening Plasmodium falciparum infection (severe malaria) in over 11,000 African children with replication data in a further 14,000 individuals. Here we report a novel malaria resistance locus close to a cluster of genes encoding glycophorins that are receptors for erythrocyte invasion by P. falciparum. We identify a haplotype at this locus which provides 33% protection against severe malaria (OR=0.67, 95%CI=0.60-0.76, P=9.5×10 −11 ) and is linked to polymorphisms previously shown to have features of ancient balancing selection, based on haplotype sharing between humans and chimpanzees 4 . Taken together with previous observations on the malaria-protective role of blood group O [1][2][3]5 , these data reveal that two of the strongest GWAS signals for severe malaria lie in or close to genes encoding the glycosylated surface coat of the erythrocyte cell membrane, both within regions of the genome where it appears that evolution has maintained diversity for millions of years. These findings provide new insights into the host-parasite interactions that are critical in determining the outcome of malaria infection.
In the discovery phase of this study we analysed GWAS data on 5,633 children with severe malaria and 5,919 population controls from The Gambia, Kenya and Malawi, and in the replication phase we analysed candidate SNPs in a further 13,946 case control samples from Burkina Faso, Cameroon, The Gambia, Ghana, Malawi, Mali and Tanzania (Extended Data Figure 1). The majority of samples used in the discovery phase have been analysed previously by lower resolution GWAS methods 3 . For this analysis we improved resolution by directly typing all samples at approximately 2.5 million SNPs using the Illumina Omni2.5M platform, followed by quality control (Extended Data Figures 1 and 2) and imputation of genotypes at over 10 million SNPs using haplotype data from the 1000 Genomes Project 6 . Imputation performance varied across populations, with accuracy decreasing as a function of the similarity between study and reference individuals (Extended Data Figure 3). When testing for genetic association, principal components analysis was used to correct for population structure (Extended Data Figure 4a-e) which reflected both geography and self-reported ethnicity. Similar results were obtained using a mixed-model approach (Extended Data Figure 4f).
To assess the evidence for association in the discovery phase we used an approach which allows for heterogeneity in the protective effect of an allele across different study sites. This could be particularly important in our data as high levels of genetic and ethnic diversity in Africa can result in variable patterns of linkage disequilibrium between study sites that can complicate GWAS analysis 7 . Other potential sources of heterogeneity include allelic heterogeneity and multiple independent origins of malaria resistance loci, as has been well documented at the HBB locus 1,3 , as well as the high levels of genetic diversity in the parasite 8 . Specifically, we used a Bayesian approach that combines evidence across multiple models of association by specifying a prior probability on the size and similarity of genetic effect across populations, as well as the mode of inheritance 1 . A single statistical summary of the signal of association was obtained by averaging the evidence across models, weighting each by its prior probability, and comparing the evidence to the null model of no association (model averaged Bayes factor; BF avg ). Having observed the data, a posterior probability was assigned to each model (Pr model ) conditional on it being a true association and the model assumptions, which are described in the Methods and Extended Data Figure  5. We replicated previously reported GWAS signals 2,3,9 at the HBB (BF avg =5.8×10 24 ), ABO (BF avg =6.7×10 9 ) and ATP2B4 (BF avg =4.4×10 5 ) loci, and a detailed analysis of key variants at these and other previously reported loci is presented elsewhere 1 . A previously reported association near the gene MARVELD3 2 did not replicate in this dataset (Supplementary Note 1). Genome-wide patterns of association with severe malaria at the 34 regions of the genome containing a variant with either a Bayes factor for the most probable model (BF max ) > 2.5×10 4 or with a model averaged Bayes factor (BF avg ) > 2.5×10 3 are summarised in Extended Data Figure 6 and Supplementary Table 1. Details of the evidence for association in these regions can be viewed online at www.malariagen.net/resource/14.
These data provide a rich resource of new candidate loci for further investigation. Here we focus on a region of chromosome 4 shown in Figure 1, where the strongest signal of association (at SNP rs184895969) is located between the gene FREM3 and a cluster of three glycophorin genes (GYPE, GYPB and GYPA). Glycophorins are sialoglycoproteins that are abundantly expressed in the erythrocyte membrane, providing a hydrophilic surface coat that is necessary for erythrocytes to flow freely in the circulation. A complex system of singlenucleotide and structural variants in this region determine the MNS blood group system 10 . These genes play a functional role in invasion of erythrocytes by P. falciparum. Glycophorin A is the receptor for the P. falciparum erythrocyte-binding ligand EBA-175 11 , and glycophorin B is a receptor for the parasite ligand EBL-1 12 . To follow up this observation, selected SNPs at this locus were genotyped by Sequenom iPLEX MassArray in the discovery and replication sample sets outlined above (Figure 1 and Extended Data Figure  7a). The combined dataset of 25,498 samples provided convincing evidence of association at rs186873296 by standard fixed-effect meta-analysis (P = 9.5×10 −11 ) as well as by the above Bayesian approach (BF overall = 1.3×10 8 ; Figure 2 and Methods). The minor allele frequency of rs186873296 was higher in East Africa than West Africa, and the greatest evidence of association was seen in Kenya where it was most common, with an allele frequency of approximately 10%. Using only replication data to avoid winner's curse, and assuming an additive model, we estimate that carrying one copy of the derived (non-ancestral) allele reduces the risk of severe malaria by about 40% in Kenya (OR = 0.60, 95% CI=0.46-0.79), with a slightly smaller effect across all populations (OR = 0.67, 95%CI 0.56-0.80 in frequentist fixed-effect meta-analysis). Further details are given in Supplementary Note 2.
The glycophorin gene cluster has a complex pattern of gene conversion and structural variation that has been previously noted; indeed it has been proposed that selective pressure due to pathogens, including malaria, has contributed to shaping diversity in this region 10,[13][14][15][16] . Using the human reference sequence and mapped sequence read data from the 1000 Genomes Project, we identified the boundaries of a 350kb region of sequence homology surrounding these genes as well as a set of segregating gene deletions (Extended Data Figure 8). The lead imputed marker (rs184895969) is located within 10kb of this complex region. Imputation accuracy within the region is low using current reference data, so it is possible that the causal variant lies within the glycophorin genes themselves but that this is obscured in the current imputed dataset. With this caveat, we computed a credible set of putatively causal variants in the region; this set includes both the lead imputed marker and a linked missense mutation (rs181620317) in FREM3 (Extended Data Figure 7b). We note that the protective allele at rs184895969 was associated with increased GYPA transcription in published gene expression data for HapMap lymphoblastoid cell lines 17 (P = 0.016; Extended Data Figure 9); other regional analyses are described in Supplementary Notes 1-2. Improved African genome variation reference panels 7,18 are needed to understand the complex patterns of variation in this region so that the causal variant of malaria resistance can be fine-mapped with greater confidence.
A striking feature of these data is that all of the loci which reach conventional criteria for genome-wide significance (P < 5×10 −8 using a fixed effects model) are in or near genes that play a key role in erythrocytes (HBB, ABO, ATP2B4, FREM3/GYPE), the primary host cell of P. falciparum. Other erythrocyte-related genes are identified in the discovery phase analysis but do not reach genome-wide significance, including EPB41, which encodes erythrocyte membrane protein band 4.1 and has been implicated as a possible receptor for P. falciparum invasion 19 (rs2985337: BF avg =3443; fixed-effect meta-analysis OR=1.16, 95%CI=1.09-1.23, P=1.2×10 −6 ; see Extended Data Figure 6).
An intriguing feature of the ABO locus is that it contains a number of polymorphisms that are shared between humans and other primates, and recent analyses of sequence variation across species indicated that some of these are ancient polymorphisms that have been maintained by balancing selection over millions of years 20 . The current findings are of particular interest since the FREM3/GYPE is one of the most prominent examples of putative ancient balancing selection in a genome-wide analysis of haplotype sharing between humans and chimpanzees 4 . The peak GWAS signal at this locus is less than 45kb away from the shared human-chimp haplotype, which falls within the region covered by the credible set (Extended Data Figure 7b), although it does not exhibit a strong association with severe malaria (Supplementary Note 3). To explore the genealogical relationship between the putative ancient balanced polymorphisms (ABPs) and SNPs associated with severe malaria in the FREM3/GYPE region, we inferred an ancestral tree 21 from the African (YRI + LWK) part of the 1000 Genomes data and used it to order haplotypes in the region, labelled with the positions of ABPs, malaria-associated SNPs, and other variants of interest ( Figure 3). All three haplotypes carrying the protective allele at the directly-typed marker with most evidence of association (rs186873296) carry the minor allele at the ABP markers (D'=1, P=0.017). By inferring the positions of putative causal mutations on the estimated genealogical tree 21 at the lead imputed marker (rs184895969) we found evidence for a single protective mutation in Kenya (log10 BF=3.09; OR=0.6) estimated to lie on the branch ancestral to the protective allele at the lead marker. Although the most likely position for an additional mutation was on the branch ancestral to the ABPs, a single haplotype explains most of the signal of association in this region.
These observations raise the question of whether malaria resistance loci are more likely to be found in regions of the genome that show evidence of ancient balancing selection. We therefore analysed the relationship between the regions of association in our GWAS and 125 regions of the genome found by Leffler et al 4 to contain haplotypes shared between humans and chimpanzees. The SNPs defining these haplotypes (ABPs) were not themselves enriched for association with severe malaria (P > 0.1, Methods). We used a simulation approach to assess the physical proximity of ABPs to the peak of association within the 34 strongest regions of association (tier 1) and 73 weaker signals (tier 2), and observed a significant relationship with tier 1 over a range of length scales (Extended Data Figure  10a,d). We also identified the nearest gene to the lead marker within each association region as well as the gene nearest to each ABP haplotype, and performed a gene-based test for genome-wide enrichment. Strong evidence for enrichment was seen at tier 1 loci (OR 41.0; P=4×10 −7 by Fisher's exact test, P sim =5×10 −4 using a simulation approach described in Methods, and P=1×10 −4 using the INRICH algorithm 22 ; Extended Data Figures 10b-d) and a weaker trend at tier 2 loci (OR 7.7, P sim =0.15). Apart from ABO and FREM3/GYPE, there were six other GWAS loci (4 in tier 1, 2 in tier 2) where the nearest gene to the lead marker was also the nearest gene to an ABP haplotype (DSCAM, NRG1, CNTNAP5, TMEM132C, CACNA2D1, RYR2). Although the current association evidence at these loci do not satisfy conventional criteria for genome-wide significance and they should be regarded as putative until convincingly replicated, it is noteworthy that they are all involved in key aspects of membrane biology (Supplementary Note 3).
In the largest genetic association study of malaria to date, we have identified a new locus of resistance to severe malaria that lies next to a cluster of glycophorin genes involved in erythrocyte invasion by P. falciparum, and that also overlaps a locus of putative ancient balancing selection identified by analysis of haplotype sharing between humans and chimpanzees. It is possible that malaria is not the cause of the ancient balancing selection, or that it is just one of a number of opposing evolutionary driving forces, as at ABO, where blood group O reduces the risk of severe malaria but increases the risk of severe cholera 23 . Nonetheless, these new findings raise the intriguing question of whether natural selection on malaria susceptibility has been shaping genetic diversity in humans and their ancestors for millions of years. P. falciparum is closely related to the chimpanzee parasite P. reichenowi and other parasites of African Great Apes [24][25][26] . It has been proposed that P. falciparum was introduced into human populations from chimpanzees or gorillas in the recent past, but this remains a matter of intense debate [25][26][27] . Population genetic data are consistent with an ancient origin followed by a marked expansion of the parasite population approximately ten thousand years ago, coincident with the introduction of agriculture 28 . The P. falciparum genome possesses a huge repertoire of polymorphism 8 and it is possible that the host and parasite genomes are engaged in a longstanding evolutionary arms race, each maintaining diversity to try to outflank the other 29 . Intriguingly, the parasite surface receptor EBA-175, which directly binds glycophorin A during red cell invasion, also contains structural polymorphisms with features of ancient dimorphism 30 . The present findings provide new leads both to investigate these evolutionary mechanisms and to discover further genetic determinants of resistance to severe malaria in African children.

Samples, ethics and clinical information
Samples were collected from nine partner projects from across sub-Saharan Africa (Extended Data Figure 1a) as described previously 1,3 . The studies and sample sets described in this manuscript form part of a larger ongoing project within the Malaria Genomic Epidemiology Network (www.malariagen.net). We used the World Health Organisation definition of severe malaria which comprises a broad spectrum of life-threatening clinical complications of Plasmodium falciparum infection 29,31 . Investigators from study sites worked together to agree on principles for sharing data and standardising clinical definitions, and to define best ethical practices across different local settings including the development of guidelines for informed consent. Relevant ethics committees are listed in Extended Data Figure 1a. Further information on policies, research and the consent process may be found on the MalariaGEN website (http://www.malariagen.net/community/ethics-governance).

DNA extraction and Sequenom typing
As described previously 1,3 , all samples submitted to the MalariaGEN Resource Centre underwent a standard set of procedures that included quantification using picogreen, genotyping of 65 polymorphisms (including HbS -rs334, and 3 gender-typing SNPs) on the Sequenom iPLEX MassArray platform and matching to baseline clinical data (e.g. gender, ethnic group and case-control status).

High-density genotyping
Three cohorts (Gambia, Malawi and Kenya) were genotyped on the Illumina HumanOmni2.5-4 (Kenya) and Illumina HumanOmni2.5-8 (Gambia, Malawi) platforms. As described previously 3 we used three different calling algorithms (Illuminus 32 , Illumina's Gencall algorithm as provided in BeadStudio, and GenoSNP 33 ), each of which uses slightly different information in the data. We formed final genotype calls by taking consensus between the three algorithms. Genotypes where any two of the three calling algorithms were discordant, and genotypes where fewer than two algorithms were confident enough to make Page 5 Nature. Author manuscript; available in PMC 2016 April 08. a call were treated as missing. This process showed improved calling, evaluated using Mendelian error counts in a subset of Kenyan samples, relative to each of the three algorithms separately (data not shown).
Following genotype calling, we aligned genotypes to the forward strand of the human reference sequence (GRCh37), using both the Illumina-supplied manifest and publically available strand files (www.well.ox.ac.uk/~wrayner/strand) obtained by mapping allele probes to the reference by BLAT 34 . We removed SNPs whose position or strand mismatched between the Illumina manifest and the strand file. To simplify analyses, we restricted attention to the set of SNPs having the same name, chromosome, position, strand, and probe sequences across the two genotyping platforms. Because the Omni2.5M contains multiple probes for some variants, we further removed SNPs to ensure positions were unique. In total we were left with 2322985 SNPs in each cohort across the autosomes and the X and Y chromosomes. We note that SNPs annotated as lying on the pseudo-autosomal region (PAR) of the X chromosome were not included, as these had position equal to 0 in the manifest for the HumanOmni2.5-4 array. Finally, we flipped alleles where necessary so that in downstream analyses the first allele always corresponds to the reference allele of the human genome sequence.

Sample Quality Control
We performed sample QC separately on each cohort by computing autosome-wide averages of normalised X channel intensity, normalised Y channel intensity, and heterozygosity and missingness based on the consensus call. To identify outlying samples we applied ABERRANT 35 , adjusting the lambda parameter per cohort to account for differences in genotyping quality (Extended Data Figure 1d-f). In Gambia, a tail of samples showing low heterozygosity but otherwise appearing to be well typed was apparent. We explicitly included these samples in downstream analyses.
To estimate genome-wide relatedness between samples, we selected a list of 178775 high quality SNPs satisfying the criteria missingness < 1%, MAF > 1% and thinned to be at least 0.005cM apart and to have pairwise r 2 < 0.3. Treating each cohort separately, we used SHELLFISH (www.stats.ox.ac.uk/~davison/software/shellfish/shellfish.php) to compute a matrix of pairwise relatedness values R = (r ij ) between samples, where r ij denotes the genome-wide average covariance of frequency-normalised genotypes in individuals i and j.
In samples with few close relationships, the value of r ij can be thought of as an estimate of kinship 36 with, for example, values close to 1 representing identity between samples, and values close to zero reflecting a lack of close relatedness (relative to the rest of the sample). To remove samples with duplicate typing and close relationships, we excluded one of each pair of samples with r ij > 0.2, taking all remaining samples through to phasing and imputation. Extended Data Figure 1c lists the number of samples before and after QC, and the number removed by each QC step; (*) denotes the number of samples removed after explicitly including samples with low heterozygosity in Gambia.
We used SHELLFISH to compute principal components (PCs) on the post-QC sample set. Consistent with our removal of poor quality, duplicate and closely related samples, principal components plots show no substantially outlying samples (see Extended Data Figure 4a-c).
The top few principal components in each cohort reflect substantial population structure, as evidenced by colouring by reported ethnicity. As found previously 3,9 , the top principal components are also significantly correlated with case-control status (Extended Data Figure  4d), indicating that population structure may act as a confounding factor in association analyses if not controlled.

Genotyping Quality Control
Treating each cohort separately, we used SNPTEST (mathgen.stats.ox.ac.uk/ genetics_software/snptest/snptest.html) to test for association at each autosomal SNP using a genotypic model of association which allows for different effects at heterozygote and homozygote genotypes, including the leading 5 principal components as covariates. Genotypic model association tests are particularly sensitive to confounding by genotyping error 37 , and we used the resulting P-values as a guide to finding appropriate QC criteria. To detect potential spurious genotypes due to batch effects, we modelled genotypes as predicted by the leading 5 principal components and case/control status in a linear regression framework in R 38 , and tested whether including an indicator of the plate on which each sample was genotyped contributed significantly to model fit. We refer to this as the "Plate test". For downstream analyses we excluded SNPs with minor allele frequency < 1%, missing data proportion > 5% (in Gambia and Malawi) or > 2.5% (in Kenya, which had fewer missing genotypes overall), Hardy-Weinberg P < 1×10 −20 in controls, or Plate-test P < 0.01. We inspected cluster plots of all remaining SNPs showing association test P < 1×10 −5 (Gambia) or P < 1×10 −4 (Kenya, Malawi) and excluded those with clear genotyping problems. Extended Data Figure 2a shows the number of SNPs excluded by each criterion. The post-exclusion genome-wide association analysis is shown in Extended Data Figure 2c.
SNP QC on the X chromosome was performed as on the autosomes, with a few differences as follows. We treated males and females separately, using a genotypic model of association in females and an allelic model in males (who have only one copy of the X chromosome). Because genotype calling was performed blind to the gender of samples, some males appear to be heterozygous at some X chromosome SNPs. We treated all such heterozygous calls as missing. We computed missing data proportion and plate test P-values in males and females separately, and tested for departure from Hardy-Weinberg in female controls and for differences in frequency between males and females. We excluded SNPs with MAF < 1%, missingness > 5% (Gambia, Malawi) or 2.5% (Kenya) in males or females, or plate test P < 0.01 in males or females. We further excluded SNPs with Hardy-Weinberg P < 1×10 −20 in female controls or showing significant difference in allele frequencies between males and females (P < 1×10 −20 ). Extended Data Figure 2b shows the number of SNPs excluded by each criterion.
Some regions of the X chromosome showed an elevated number of male heterozygote calls, contributing to the high number of SNPs excluded due to missingness in males. These included SNPs in the pseudo-autosomal regions at either end of the chromosome 39 (indicating that XY designation in the chip manifests does not adequately cover these regions) as well as SNPs within the X transposed region near the centromere.

Phasing and imputation
We phased genotype data in each cohort separately using SHAPEIT v2.r644 40 , specifying 200 hidden states and an effective population size of 17469, as recommended for African populations by the SHAPEIT documentation, and phasing each chromosome separately. We used IMPUTE v2.3.0 41,42 to impute phased genotypes into the 1000 Genomes 6 Phase I integrated variant set (version of 24th August 2012, as downloaded from the IMPUTE website, which we refer to here as the 1000 Genomes reference panel) in 5Mb chunks with a buffer region of 500kb. For phasing and imputation we used the combined HapMap recombination map in build 37 coordinates included with the 1000 Genomes reference panel. Unless otherwise stated, downstream analyses included only SNPs with minor allele frequency > 0.5% and impute info metric > 0.75.

Assessment of imputation
At each genotyped SNP, IMPUTE computes squared correlation (referred to here as accuracy) and concordance between typed genotypes and genotypes obtained by masking the SNP and re-imputing. To assess imputation performance, we plotted the cumulative distribution of concordance (Extended Data Figure 3a) and accuracy (Extended Data Figure  3b), accuracy by frequency (Extended Data Figure 3d) as well as average per-sample accuracy (Extended Data Figure 3c). We also assessed accuracy relative to direct typing on the Sequenom platform at variants typed in the FREM3/GYPE, ARL14 and INPP4B regions of association as described below.

Association testing and meta-analysis
Association testing-We used SNPTEST to test for association at approximately 38 million variants obtained through imputation, including 5 PCs as covariates to control for population structure separately in each cohort. SNPTEST uses a missing data likelihood to account for the uncertainty in genotypes at imputed SNPs. We fit additive, dominant, recessive, and heterozygote models of association and ran SNPTEST separately in each imputation chunk. Below, we refer to the estimated effect size for population i and mode of inheritance m as β i,m and its estimated standard error as SE i,m .
Frequentist meta-analysis-For each SNP and each mode of inheritance (additive, dominant, recessive, heterozygote) we computed the fixed-effect inverse variance-weighted meta-analysis effect size, standard error, and P-value. In this context, fixed-effect metaanalysis assumes a single true effect size which is identical between the three cohorts, and can be thought of as finding the maximum likelihood estimate of the effect size under the assumption that the likelihood in cohort i is proportional to the density of the normal distribution with mean β i,m and standard error SE i,m .
Bayesian meta-analysis-We have previously 1,3 used Bayesian meta-analysis techniques to allow for between-population heterogeneity of effect sizes in these three cohorts. Here, we applied this method to compute Bayes factors for association under four modes of inheritance and six different models of correlation between cohorts. In this framework, true effect sizes are modelled by a multivariate normal distribution centred on zero and with a given prior covariance matrix which can be written as σ 2 P, where σ 2 is a prior variance controlling the magnitude of plausible effects and P is a prior correlation matrix.
The correlation models we used were: Fixed effects (all elements of P equal to 1): as with frequentist fixed-effect meta-analysis, this assumes effect sizes are equal across the three cohorts.
Correlated effects (all off-diagonal elements of P equal to 0.96): this assumes effect sizes are similar but allows for some variability.
Independent effects (all off-diagonal elements of P equal to 0.1): this assumes there is little similarity between effect sizes in different cohorts.

Structured effects
We also considered models where effects in the two East African populations (Kenya and Malawi) were more similar to each other than to that in Gambia. We assumed either effects were fixed between Kenya and Malawi and correlated (ρ=0.96) between East Africa and Gambia (we refer to this as the fixed-structured-effect model), or that effects were correlated (ρ=0.96) between Kenya and Malawi and largely independent with Gambia (ρ=0.1) (referred to as the correlated-structured-effect model).
To form a single summary measure of evidence for association we formed a model-averaged Bayes factor (BF), referred to as the mean BF and denoted BF avg , as a weighted average of model-specific Bayes factors using the following weighting scheme. We assigned weights of 0.4 for additive mode of inheritance and 0.2 for dominant, recessive, or heterozygote modes of inheritance, reflecting a belief that variants which tag causal variants by LD are more likely to display additive effects. We assigned a weight of 0.4, 0.2, and 0.1 to fixed-, correlated-, and independent-effect models, and 0.2 and 0.1 to fixed-structured and correlated-structured effects models. Finally we assigned weight of 0.5 to small-effect models (σ 2 =0.2 2 ) and large-effect models (σ 2 =0.75 2 ). Overall prior weights were assigned by multiplying across these categories; for example, the model representing small effect size distribution, fixed-effect across cohorts with additive mode of inheritance (denoted smallfix-add) was assigned prior weight equal to 0.5×0.4×0.4=0.08, while the model representing small, correlated-structured, dominant effects (denoted small-cor-str-dom) was assigned a prior weight of 0.5×0.1×0.2=0.01. For each SNP we also recorded the model having the highest posterior weight, and refer to its Bayes factor as the max BF (denoted BF max ). Extended Data Figure 5a depicts slices through the combined prior on effect sizes across three cohorts for additive effect models. The mean BF behaves similarly to a minimum over all four fixed-effect meta-analysis P-values, but additionally captures effects that vary between cohorts (Extended Data Figure 5b).
As described previously 1,3 , to computed model-specific Bayes factors efficiently we used an approximation of the likelihood by the density of a normal distribution with the estimated Page 9 Nature. Author manuscript; available in PMC 2016 April 08. mean (β i,m ) and standard error (SE i,m ) in each cohort. Thus, overall, observed effect sizes are modelled as normally distributed around zero with covariance that depends on the prior covariance in true effect sizes and on model standard errors, where V is a diagonal matrix with ith diagonal entry equal to the squared standard error, SE 2 i,m . The approximate or asymptotic Bayes factor can then be computed by evaluating a ratio of two normal densities. In the univariate case this method is the same as described previously by Wakefield 43 . To facilitate working with genome-wide meta-analysis results, we wrote custom software to compute frequentist and Bayesian meta-analyses, and stored details of genotype counts, association model fit, and meta-analysis results directly in a SQLite database file.
Further discussion of the Bayesian approach can be found in Supplementary Note 4.

X chromosome association testing and meta-analysis
Association testing on the X chromosome was performed as for the autosomes with a few differences as follows. We ran SNPTEST separately in males and females, estimating effect sizes under additive, dominant, recessive and heterozygote modes of inheritance in females. In this usage, SNPTEST assumes a model of complete inactivation so encodes male genotypes as 0/1 and females as 0/0.5/1 for an additive mode of inheritance. We than metaanalysed the six gender-specific association analyses to produce frequentist fixed-effect and model-averaged bayesian meta-analyses.
For Bayesian analysis, in addition to summing over models of between-population heterogeneity, we adopted the view that differences in sex might lead to heterogeneity in effect. We therefore included models of heterogeneity between males and females as follows. Let ρ sex denote the correlation between effects in male and female samples within a single population. We included models where males and females have the same effect (ρ sex =1, termed fixed-sex model and given prior weight 0.45), correlated effects (ρ sex =0.96, termed correlated-sex model and given prior weight 0.225) or essentially independent effects (ρ sex =0.1, termed independent-sex model and given prior weight 0.225). Because some parts of the X chromosome escape inactivation 39,44 , we also included a model where the effect in females is twice that in males (given prior weight 0.1).
To fully specify prior correlation for each model, for each pair of populations A,B we also need to specify the correlation (denoted ρ cross ) between males in A and females in B. We set ρ cross = ρ pop × ρ sex where ρ pop is the chosen prior correlation between same-sex samples in A and B (i.e. ρ pop =1, 0.96 or 0.1 as defined above).
As above, we formed overall weights by multiplying across categories, so that for example the model of small, additive effects that are fixed across populations and across sexes (denoted small-fixsex-fix-add) had prior weight 0.5×0.4×0.4×0.45=0.036.

Page 10
Nature. Author manuscript; available in PMC 2016 April 08.

Linear mixed model analysis
To compare association test results for logistic regression as implemented in SNPTEST with the use of a linear mixed model, we reran association test scans in each discovery cohort using the program MMM 45 . We used the same genome-wide relatedness matrix as used to compute principal components (see above) and assumed an additive model of association. We plotted the -log10( P-value ) for MMM against the corresponding -log10( P-value ) based on the SNPTEST scan (which used 5 PCs as described above), for each discovery population and for fixed-effect meta-analysis (Extended Data Figure 4f). P-values under both methods at tier 1 loci are listed in Supplementary Table 1.

Lead SNPs, region and tier definitions
We formed a list of lead SNPs within approximately independent regions of interest as follows. We restricted to variants with IMPUTE info measure at least 0.75 across all three populations and ranked variants by the model-averaged Bayes factor (highest to lowest). We iteratively picked lead SNPs from the top of the list and excluded other variants within a recombination interval of 0.25cM±25kb centred at the lead SNP (referred to below as the association region), continuing until no more SNPs remained with BF avg >250 or BF max >2500.
We grouped lead SNPs into two tiers as follows: • tier 1: all lead SNPs with BF avg >2500 or BF max >25000 • tier 2: all lead SNPs not in tier 1 with BF avg >1000 or BF max >10000 In total, across the autosomes and the X chromosome there were 34 regions in tier 1 and 73 in tier 2, with association regions covering approximately 13Mb and 26Mb of the genome respectively.

Regional association analysis
For each tier 1 and 2 region, we examined the pattern of association in the region, generating a regional association plot for the region annotated with details of the meta-analysis for the lead SNP as follows.
LD computation-In each region we computed pairwise LD statistics between the lead SNP and surrounding SNPs using best-guess imputed haplotypes for control samples across three populations. To facilitate this computation, we stored imputed haplotypes in a SQLiteformat database allowing us effective random access to haplotype data. We used R to compute both Pearson correlation coefficient (r 2 ) and Lewontin's |D'| between the lead SNP and all other SNPs in the region.
Regional association plot-For each SNP in tier 1 and 2 we plotted log10(BF avg ) in the association region around the lead SNP plus 1Mb on either side, colouring points according to LD, with outer circles representing r 2 and inner circles representing |D'|. We further annotated SNPs that were typed in at least one cohort (using black plusses) and SNPs that had Sequenom genotype data available (with black triangles). Tables 4 and 5 of Leffler et al 4 ), eQTLs from the Genotype-Tissue Expression project 46 (GTEx), and previously reported 47 GWAS loci. Where SNPs matched an imputed or typed variant, we computed LD statistics and coloured these points as above. We plotted all RefSeq genes (as downloaded from the UCSC Genome Browser MySQL database 48 on 2013-03-18) in the region, annotating the direction of transcription. We further plotted local recombination rate estimates from the HapMap combined recombination map included with the IMPUTE haplotypes described above.

Regional annotations-We plotted all ABP variants (Supplementary
Meta-analysis-We annotated plots with effect sizes and confidence intervals for the lead SNP for each mode of inheritance (additive, dominant, recessive and heterozygote) that informed the model averaging. We also produced barplots showing the posterior distribution on models of association using the prior weights described above.
Cluster plots-Finally, for each region, we produced and inspected cluster plots for all typed SNPs with BF avg >= 10 that were not excluded by QC. As phasing fills in missing genotypes and potentially improves genotype calling based on LD with surrounding SNPs, we coloured plots based on genotype calls taken from the output of phasing. Website-Regional association plots and cluster plots can be viewed online at http:// www.malariagen.net/resource/14.

Validation and replication typing
Based on a preliminary version of the data presented here we selected SNPs for typing on the Sequenom platform across the whole MalariaGEN Consortial Project 1 sample set, which includes the discovery samples, further cases and controls in the same populations that were not included in the GWAS, and further large sample sets from five other populations from sub-Saharan Africa (see Extended Data Figure 1a). Data were available for SNPs tagging the lead markers in the FREM3/GYPE, INPP4B and ARL14 regions (r 2 >0.5 in controls, as estimated using the EM algorithm) as well as other regions not represented in tier 1.

Replication analysis
Replication analysis using Sequenom data was restricted to the set of samples with less than 10% missingness as measured across the set of 70 SNPs chosen for replication typing. In each population we conducted logistic regression in R, including five principal components (discovery samples) or reported ethnicity (replication samples) as covariates to control for population structure. For each GWAS lead SNP, we examined each Sequenom SNP in the region and computed r 2 and |D'| with the lead SNP. To allow for the effects of incomplete LD on the Bayesian model fitting, we recomputed the Bayesian analysis in the discovery samples based on Sequenom genotypes at each SNP, to obtain a Sequenom-based mean Bayes factor and a 'best model' Bayes factor at the model with highest posterior weight. Where LD is incomplete or where imputation is imperfect, this model may differ from the best model for imputed data.
For each SNP we computed fixed-effect meta-analysis across discovery samples, across replication samples, and across all samples. We also computed the Bayes factor for replication samples (BF replication ) for the model with highest posterior weight in the discovery samples. To compute an overall Bayes factor for association, we combined the discovery BF avg computed at the imputed lead marker with the replication bayes factor at the Sequenom SNP, as BF overall = BF avg x BF replication . We use the imputed lead marker here because the number of discovery samples directly typed was smaller than the number of imputed samples. Conditional on the lead imputed and replication markers reflecting the same signal of association, this BF overall represents an overall measure of the evidence for association at the locus that reflects all the samples in our study.
A discussion of the replication evidence in the FREM3/GYPE, INPP4B and ARL14 regions can be found in Supplementary Note 1.
MalariaGEN encourages individual study sites to perform more detailed analyses of local patterns of disease association, and a Tanzanian-focused analysis of FREM3 and other candidate SNPs that were genotyped as part of this study is reported elsewhere 49 .

Credible interval analysis
In a given region of the genome, under the assumption that exactly one variant (that is accessible to our typing or imputation) is causal, the posterior probability that each variant is the causal variant can be computed by a simple reweighting of Bayes factors 50 . Using this approach we computed 95% and 99% credible intervals (i.e. the smallest set of SNPs accounting for 95% or 99% of the posterior mass) for variants in the association region around rs184895969, plus a margin of 50kb at either end (Extended Data Figure 7b). We noted that rs181620317, which is annotated as a missense mutation for the gene FREM3, is within the 95% confidence interval in our data. We note that while this analysis is simple and appealing, its interpretation depends on assumptions about the true disease model, and on the behaviour of imputation in the region 50 ; in particular the difficulty of imputing variants around the three glycophorin genes (Extended Data Figure 8) might make this analysis fail to capture putatively important variation within or around GYPA, GYPB, or GYPE.

Sequence homology and alignability in the glycophorin region
To investigate the location of our GWAS signal with respect to the pattern of sequence homology around the three glycophorin genes, we generated a dot plot (Extended Data Figure 8a) showing co-occurrence of k-mers in the human reference sequence 51 in the region. Considerable sequence homology was observed over a region of about 350kb covering the three genes. Our lead GWAS marker and the ABPs lie just outside this region.
The high level of homology should affect our ability to align probes or sequences to this part of the genome. To confirm this, we also plotted the UCSC alignability track ('CRG alignability 100', Extended Data Figure 8c), which shows the degree to which the 100-mer starting at each position in the reference sequence is alignable (with a value of 1/n indicating that the 100-mer aligns to n positions across the genome, allowing up to two mismatches).

Page 13
Nature. Author manuscript; available in PMC 2016 April 08.
As expected, even short regions of shared kmers affect alignability considerably, with alignability dropping to an average of about 0.7 within the large region of homology. Similarly, we plotted imputation performance (as measured by the IMPUTE info score) and observed a marked drop within the region of sequence homology.

Structural variation in the glycophorin region
We attempted to identify structural variation in the glycophorin region by examining sequence read data generated by the 1000 Genomes' project, using the set of BAM files available from the 1000G project in October 2014, downloaded from ftp://ftp. 1000genomes.ebi.ac.uk/vol1/ftp/data. These data contain a mixture of read lengths, with most reads of at least 90bp and some reads of 76bp. We considered only reads with mappinq quality at least 20. For each sample we computed coverage across the glycophorin region, normalised by average coverage computed across chromosome 1, and refer to these values here as 'coverage'.
Coverage values across the region were correlated with genome alignability (as defined above), with coverage dropping substantially in regions of low alignability. We therefore restricted attention to the set of perfectly alignable positions, here defined as positions such that every 100-mer overlapping the position in the reference sequence aligns uniquely, allowing up to two mismatches (and computed using the CRG alignability 100 track described above). As expected, coverage at perfectly alignable positions showed little or none of the variation present in genome alignability (Extended Data Figure 8e-g).
Two large structural variants, present at frequencies of at least 1% in LWK+YRI, were evident in coverage data. Genotypes for both these variants were called by the 1000 Genomes phase I and are referred to as esv2668125 (also called MERGED_DEL_2_26708) and esv2662558 (also called MERGED_DEL_2_26722). Both deletions putatively represent deletions of all or part of GYPB. We noted some uncertainty as to the location of the breakpoints of both deletions -with 1000 Genomes' breakpoints differing from the breakpoints as they appeared in coverage data by at least 10kb. We further noted that three samples (NA18519, NA19185, NA19222) that were called as heterozygote for one or both deletions by the 1000 Genomes project, appeared to have homozygous genotypes (red, blue and green lines, Extended Data Figure 8e-g).

Expression QTL analysis
To investigate the effect of associated SNPs in the FREM3 region on expression of nearby genes, we downloaded publically available data on RNA expression levels 17 in HapMap samples, and plotted expression levels of genes in the region (Figure 1) against genotypes at the lead SNP and other SNPs of interest in African samples (YRI+LWK, Extended Data Figure 9). Data was available for genes INPP4B, USP38, GAB1, SMARCA5, GYPE, GYPB, GYPA, HHIP, ANAPC10 and ABCE1, with three probes available for GYPE, two for GAB1 and one probe for each of the other genes. Five samples in this dataset (NA19318, NA19324, NA19377, NA19429, NA19190) carry the protective allele at the lead marker rs184895969, while only two carry the protective allele at the directly-typed rs186873296. For each gene and SNP we tested for a trend of genotype on expression using linear Page 14 Nature. Author manuscript; available in PMC 2016 April 08.
regression. Extended Data 9 shows all probes for the glycophorins as well as other regional genes for which a P-value less than 0.05 was observed.

GENECLUSTER analysis
We used the program TREESIM 21 to estimate an ancestral recombination graph for the LWK and YRI individuals in the 1000 Genomes Project haplotype data in a region from 144.6Mb to 144.8Mb on chromosome 4 centred on the lead SNP in the FREM3 / GYPE region. We then ran GENECLUSTER 21 in the region, allowing it to assign either one or two causal mutations in each marginal genealogy. GENECLUSTER works by probabilistically assigning study individuals to the tips of the tree estimated by TREESIM, and attempts to explain case/control status by assigning causal mutations to the branches of the tree in a Bayesian framework. This analysis is somewhat similar in spirit to our marginal SNP analysis, but has important differences. In particular, GENECLUSTER may detect associated mutations anywhere on the tree (which may not correspond directly to any typed or imputed variant), and can assess models of association involving more than a single causal variant. However, GENECLUSTER does not take into account principal components or otherwise control for population structure, and for computational reasons we only included the Luhyan (LWK) and Yoruban (YRI) populations from the 1000 Genomes reference panel in the analysis.
To investigate the relationship between ancestry and variants of interest, we plotted IMPUTE reference panel haplotypes in the region ordered by the marginal tree at the position of the imputed lead marker (rs184895969), annotating the lead imputed and Sequenom-typed markers, ABPs, previously reported GYPE eQTLs 52 , common deletions, and variants determining the M/N and S/s blood groups as described in Supplementary Note 2 ( Figure 3).

Analysis of enrichment of malaria-associated loci in functional categories
Full details of enrichment analyses are provided in Supplementary Note 3.
Further software for generating key results for this paper will be made available at http:// www.malariagen.net/resource/14. Extended Data Figure 9.
Correlation between the genotypes at SNPs of interest within the GYPE/A/B locus and reported gene transcription levels in samples from the YRI and LWK HapMap cohorts 17 . P values are for a trend test of association where more than one genotype class is present. Only assays targeting the glycophorins, and those with a P-value below 0.05 are shown.

Page 27
Nature. Author manuscript; available in PMC 2016 April 08.

Author Manuscript
Author Manuscript Author Manuscript

Author Manuscript
Extended Data Figure 10.
Detail of enrichment analysis. a) Red histogram: the empirical distribution of the log 10 distance of observed tier 1 loci to the nearest ABP haplotype. Grey histogram: distribution of distances for 10,000 simulated tier 1 sets. b) The log 10 distance of tier 1 (filled red circles) and tier 2 (empty circles) loci to the nearest ABP, plotted against their rank in BF avg order (stronger signals have lower rank). Loci are annotated with the nearest gene where a gene exists within the association region. Asterisks denote nearest genes that are also the nearest gene to an ABP shared haplotype. c) Empirical null distribution of the odds ratio for the enrichment of tier 1 loci in the set of genes closest to an ABP shared haplotype, based on Page 28 Nature. Author manuscript; available in PMC 2016 April 08.
10,000 simulated SNP sets. The red asterisk and text indicate the odds ratio for the observed tier 1 loci. d) Distribution of the proportion of the genome which identifies a given gene as nearest, for genes in or not in the set annotated as nearest an ABP haplotype. Left: distribution of the length of the genome for which the given gene is unambiguously the closest gene. Middle: distribution of the number of SNPs in our study for which the given gene is the closest gene. Right: distribution of the number of SNPs in our study for which the given gene is the nearest gene within a recombination interval of 2.5cM±25kb around the SNP, as used to determine nearest genes to GWAS lead SNPs. e) Empirical P-values for enrichment of ABP haplotypes and coding SNPs in tier 1 and tier 2 GWAS regions. Second column: P-values for enrichment by gene overlap. Third to tenth column: P-values for enrichment by proximity at different length scales. †Results for simulations using SNPs frequency-matched to GWAS tier 1 loci in 1% frequency bins. † †Results for simulations excluding the regions of ABO, HBB, ATP2B4, FREM3, INPP4B, and HHIP-AS1. Signal of association with severe malaria across the FREM3/GYPE region. a) evidence for association (log 10 BF avg ) in the discovery data. Black plusses denote SNPs that were directly typed, and black triangles denote SNPs selected for typing on the Sequenom platform. Dotted red vertical lines indicate a region of 0.25cM±25kb centred at the lead SNP (rs184895969). Coloured circles denote the correlation (outer circles) and |D'| (inner circles) with rs184895969 in controls, computed from imputed haplotypes. b) Polymorphisms shared between humans and chimpanzees, eQTLs, and previously reported associations with other phenotypes. c,d) Genes in the region and the HapMap combined recombination rate.

Page 34
Nature. Author manuscript; available in PMC 2016 April 08.

Figure 2.
Evidence for association at SNPs in the FREM3/GYPE region assuming an additive model of association. a) Forest plot showing sample size, estimated odds ratio and 95% confidence interval for the lead imputed SNP in each population and under fixed-effect meta-analysis. The frequency of the protective allele in controls in each population is shown to the right. b) The posterior weight on different models of heterogeneity at rs184895969 under the prior used in the GWAS. Model names are described in Methods. c) Forest plot for the Sequenom-typed SNP rs186873296 in discovery and replication samples, with fixed-effect Page 35 Nature. Author manuscript; available in PMC 2016 April 08.