Levels of human proteins in plasma associated with acute paediatric malaria

Sample cohort

Plasma samples from paediatric malaria patients, as well as community controls, were collected in the regional hospital as well as health centres in the catchment area of Nyagatare, Eastern Province, Rwanda. The patients were enrolled over a 3-year period, that begun in 2011 and ended in 2014. The attending paediatrician or research assistant recorded clinical and biometric parameters at the sample collection site. The clinical assessment of the disease stratifications was performed simultaneously as the samples were drawn, i.e. at the emergency room as the patients received initial triage and supportive care. The children enrolled in the study were aged from 3 months up to 6 years. Community controls were recruited from healthy siblings or other children in the villages where the patients lived. In total, 541 samples were analysed with associated clinical data, including 178 controls, 183 patients with mild malaria and 180 patients with severe malaria. Clinical patient characteristics and the different severe disease stratifications are present in the group with severe malaria as detailed in the supplementary information (Additional files 1, 2). A positive blood smear was used to confirm active infection with P. falciparum malaria for patients with both mild and severe disease. All healthy community control patients were screened and found negative for malaria. The subtypes of severe malaria included in the study were hyperparasitaemia, respiratory distress, hyperpyrexia, prostration, cerebral malaria, malaria with severe anaemia, febrile convulsions, or combinations thereof. This the classification of severity was based on the WHO criteria for the identification of mild and different manifestations of severe malaria [22]. HIV positivity was an exclusion criterion and all known HIV positive patients were excluded from the initial screening process. All patients were subsequently screened for HIV and the study subjects found positive were excluded from the study. The attending clinician duly denoted any suspected co-infections in the clinical data.

Blood samples were collected at admission in ethylenediaminetetraacetic acid (EDTA) tubes. The blood was processed immediately after sampling from the patient into separate blood constituents by gel flotation (Polymorphoprep, Alere Technologies AS). The plasma samples were subsequently snap frozen in liquid nitrogen and stored at below − 80 °C until proteomic analysis was performed [22]. Plasma samples from the patients were directly labelled using biotin for the bead array analysis (Additional file 3A).

Antibody selection

Antibodies for this screening were selected to validate previous findings by Bachmann et al. [21] in a new independent cohort, as well as for the discovery of new candidate protein markers to unravel the mechanism behind various stages of malaria infection. In addition to these previously suggested markers, additional proteins were selected from protein categories such as inflammatory proteins, erythrocyte surface proteins or proteins associated with microvesicles.

Overall, a total of 255 antibodies were included, targetting 115 unique human proteins, as well as assay specific controls. The majority of the antibodies used were polyclonal rabbit antibodies (Atlas Antibodies), originating from the Human Protein Atlas (HPA) (http://www.proteinatlas.org). In addition, 20 antibodies from R&D Systems, Abcam and Novus Biologicals were also included. Two assay specific controls were included in order to control the background; rabbit IgG (P120-301, Bethyl laboratories) for background binding to rabbit IgG molecules and a blocked bare bead (without coupled antibody) for background binding to the beads. Furthermore, two additional assay specific controls were used as sample loading controls; anti-human IgG (309-005-082, Jackson) and anti-human albumin (A0001, Dako). The selected antibodies were coupled onto carboxylated beads for the bead array assay (Additional file 3B).

Assay procedure and readout

The assay procedure was performed as described previously [23, 24]. Biotin-labelled samples were thawed at 4 °C, vortexed and centrifuged at 2000 rpm for 1 min (Beckman Coulter). The plasma samples were diluted by adding 1 μL labelled plasma to 50 μL assay buffer in 96-well microtitre plates (732-4828, Thermo Scientific) using a liquid handler (SELMA, CyBio). The assay buffer consisted of 1× PBS with 0.5% (w/v) polyvinyl alcohol (P8136, Sigma-Aldrich), 0.8% (w/v) polyvinylpyrrolidone (PVP360, Sigma Life Science), 0.1% casein (C5890, Sigma Life Science) and 0.5 mg/ml rabbit IgG (P120-301, Bethyl Laboratories). The diluted samples were subsequently heat-treated at 56 °C for 30 min in a water bath (TW8, Julabo) followed by 10 min cooling down at room temperature. The samples (45 μL) were transferred to 5 μL bead stock in a 384-well microtitre assay plate, (781 101, Greiner BioOne), using the liquid handler before overnight incubation at 650 rpm at room temperature. Afterward, the beads were washed three times with 60 μL PBS-T and incubated for 10 min at room temperature at 650 rpm with 50 μL 0.4% paraformaldehyde (PFA) (43368, Alfa Aesar) in PBS-T. Beads were washed three times with 60 μL PBS-T before incubation with 50 μL of 0.5 μg/mL RPE-labelled streptavidin (SA1004-4, Invitrogen) for 20 min at room temperature at 650 rpm. Finally, the beads were washed three times with 60 μL PBS-T, before they were re-suspended in 60 μL PBS-T and analysed on a Flexmap 3D instrument (Luminex corp.). Binding events were displayed as median fluorescence intensity (MFI) where at least 50 beads per bead ID were counted. To assess inter-assay reproducibility, the assay run was repeated using a re-labelled subset (n = 325) of samples.

Data analysis

Data analysis and visualization were performed with the open source software R (http://www.r-project.org) [25] and various R-packages, as described below.

As a part of quality control, inter- and intra-assay correlations were used to evaluate assay reproducibility using nonparametric Spearman’s correlation coefficients from the R-package ‘stats’. In order to explore sample heterogeneity, the distribution of the samples was illustrated with a nonlinear dimension reduction algorithm called t-distributed stochastic neighbour embedding (t-SNE) [26, 27]. The data was scaled and centred [28] before using the function Rtsne from the R-package ‘Rtsne’ with a perplexity of 30.

The univariate analysis included evaluation of individual protein’s ability to separate two sample groups. Student’s t-tests were used on log-transformed data for the analysis of the differences between controls and malaria cases (including both mild and severe malaria cases). Wilcoxon rank sum tests were used for the comparison between mild and severe malaria. False discovery rate (FDR) by Benjamini & Hochberg [29] was used to correct the p-values for multiple testing between controls versus malaria cases as well as mild versus severe malaria. Wilcoxon rank sum tests were used to calculate p-values for the pairwise comparisons of different severe disease stratifications (these p-values were not FDR adjusted). Boxplots of individual protein levels for different sample groups were created with R-package ‘graphics’ and ‘beeswarm’. Receiver operating characteristic (ROC) curves were used to demonstrate the classification power between two sample groups for the different proteins. Logistic regression was employed before ROC curves were created using log10-transformed data and the R-package ‘ROCR’.

Multivariate analysis was used to identify patterns in the data and to identify a potential combination of markers with higher separation power. Classification of protein level profiles was made using the unsupervised hierarchical clustering method Self-Organizing Tree Algorithm (SOTA) [30]. The data was scaled and centred [28] before using the R-package ‘clValid’ for SOTA clustering. In order to build a multivariate classification model the L1-penalized logistic regression, Least Absolute Shrinkage and Selection Operator (Lasso) [31] and logistic regression were combined. Lasso was applied using the R-package ‘glmnet’ to identify the driving variables in the dataset to evaluate potential multi-protein panels among all the antibodies cleared in the quality control. For the Lasso procedure, the MFI data were randomly divided into a training set and a test set. The training set consisted of 361 samples in total, including 119 control samples, 122 samples with mild malaria, and 120 samples with severe malaria. The test set consisted of 180 samples including; 59 control samples, 61 samples with mild malaria, and 60 samples with severe malaria. The tuning parameter lambda was determined by a five-fold cross-validation. The Lasso procedure determined the selection of contributing proteins using the value of lambda. Based on the output from Lasso, antibodies were selected to be included in a multi-protein panel. A model was built with this panel of antibodies using logistic regression on log10-transformed data. The classification power of the model was presented with ROC curves as described above and the area under the curve (AUC) was calculated, using the test set of samples. Heatmaps were used for exploratory analysis of samples and proteins using the R-package ‘pheatmap’ using hierarchical clustering with Euclidean distance on log10-transformed MFI data divided into two clusters.

Association networks between the proteins were explored using STRING v10.5 (https://string-db.org). The analysis was performed using all 61 presented proteins (Additional files 4, 5) as input and showed median confidence scores of minimum 0.4.

Study and experimental quality control

Sample heterogeneity was explored using t-SNE plots to visualize how clinical data correlated with the protein profiles. Six parameters were selected for visualization (Additional file 6), based on the presented proteins (Additional files 4, 5). As expected, a visual trend in separation could be observed between controls and malaria samples for the following parameters: disease subtype, body temperature and signs of dehydration. No separation was observed for sex, age or weight.

During the assay analysis, MFI was registered for each antibody-coupled bead ID and sample. Coefficients of variation (CV) were calculated for all replicated samples to assess the technical variability of the assay. All median CVs were below 10% (IQR = 5), both between replicates within the same labelling plate and between all replicates in the same assay plate. Furthermore, the median CV between assay plates in the first run was 8% (IQR = 2) and the median CV between the first and the second run was 10% (IQR = 3). Spearman’s correlations were calculated (Additional file 7) to further investigate intra- and inter-assay variability. The Spearman’s Rho for intra-assay correlation was 0.999 and for the inter-assay correlation 0.965. Signals from the two positive loading controls (anti-albumin and anti-human IgG) were used to identify wells to be excluded from further analysis. A stringent quality control of the antibodies was performed, including a Spearman’s correlation between assays above 0.85. Inter-assay correlations are shown for each presented antibody (Additional files 8, 9, 10). Signals correlating with negative controls (rabbit IgG and bare bead) were also excluded from further analysis. After the quality control, the dataset contained protein profiles of 70 proteins, targetted by 99 antibodies. Among the presented antibodies targetting the same protein, 13 proteins had antibody correlations above 0.7 (Additional file 11).

Protein profiles separate malaria cases from community controls

First, the relative protein levels were examined for all antibodies that passed the quality control to identify proteins with divergent amounts in community controls compared to malaria patients. Significant alterations in relative protein levels (p-values < 0.001) were observed for 76 antibodies corresponding to 57 proteins (Additional file 4).

The selected 27 proteins with the most significant differences were further characterized by SOTA clustering to reveal trends in protein level profiles with respect to disease progression (Fig. 2). From this clustering, 24 proteins showed a trend with higher levels in malaria patients compared to controls (Fig. 2a). In addition, an incremental trend in patients with severe disease compared to mild was observed. The remaining three proteins showed lower plasma levels in the malaria patient group compared to controls (Fig. 2b).

Furthermore, ROC curves were constructed for each respective protein, which together with significance levels further reveal the classification power of the selected proteins (Additional file 8). Six representative examples, showing both increased (ADSSL1, CEBPA, CRP, LBP) and decreased protein levels (CCL5, SPARC) in malaria patients compared to controls are shown in (Fig. 3).

Lasso regression analysis was used to identify whether combinations of proteins allowed an even higher discriminating power. Six proteins, covered by eight antibodies, were identified as the driving parameters that gave the optimal separation between the two sample groups in the training set. A multi-protein ROC curve based on logistic regression using all eight antibodies showed an AUC of 0.98 in the test set (Additional file 12). When including only one antibody per protein in the logistic regression model, this multi-protein 6-panel ROC curve still showed an AUC of 0.98 and the individual protein AUCs ranged between 0.68 and 0.95 (Fig. 4a). Out of the six proteins, three were identified as showing increased levels (CRP, CSF1, LBP) in malaria patients compared to controls while the remaining three displayed decreased levels (CCL5, CTSD, SPARC). Even when the classical inflammatory proteins CRP and LBP (also with the highest individual AUCs) were excluded, the joint AUC for the 4-panel was still 0.95 (Fig. 4b).

In addition, hierarchical clustering was used to identify potential sample clusters using all 541 samples and the 27 proteins (p-value < 10⁻¹⁴) (Fig. 5). Here, one cluster consisted of a majority of malaria cases (95% malaria cases and 5% controls) while the other consisted of a majority of control cases (78% controls and 22% malaria cases). Further scrutiny of the few malaria cases in the second cluster revealed the majority to be cases with mild disease.

Protein profiles in mild and severe malaria cases

Clinical manifestations of malaria infection are heterogeneous and cover a spectrum from mild disease to various degrees of severity. Less pronounced differences were found when comparing protein levels from patients with mild and severe malaria than between malaria cases and controls, which is to be expected. Hence, adjusted p-values < 0.05 were defined as significant for further analysis. Here, 37 proteins (number of UniProt IDs) were identified with the chosen significance level (Table 2 and Additional file 5). The identified proteins belonged to four major groups based on biological function: reactants of the acute phase response (e.g. CRP, LBP, TNF); erythrocyte specific (e.g. ANK1, BPGM, GYPC,); cell migration and adhesion (e.g. ADAMTS13, VWF); and tissue remodelling (e.g. AGT, CA2, MMP2) (Table 2). The relative differences for four representative proteins (Fig. 6) are displayed, out of the total set of 37 proteins (Additional file 9).

Primary biological process	Gene name	Uniprot ID	Gene description	p-value	Antibody
Cell migration and adhesion	ADAMTS13	Q76LX8	ADAM metallopeptidase with thrombospondin type 1 motif 13	2E−02	HPA042014
	CSF1	P09603	Colony stimulating factor 1	3E−02	HPA061864+
	IGFBP1	P08833	Insulin like growth factor binding protein 1	3E−02	MAB675 R&D
	ITGAV	P06756	Integrin subunit alpha V	2E−02	HPA004856
	VCAM1	P19320	Vascular cell adhesion molecule 1	2E−02	HPA069867+
	VWF	P04275	von Willebrand factor	4E−02	HPA002082
Erythrocyte associated	ANK1	P16157	Ankyrin 1	5E−02	HPA004842
	BPGM	P07738	Bisphosphoglycerate mutase	1E−02	HPA016493
	EPB41L2	O43491	Erythrocyte membrane protein band 4.1 like 2	3E−02	HPA005730
	GYPC	P04921	Glycophorin C (Gerbich blood group)	2E−02	HPA008965
	HBA1/2	P69905	Hemoglobin subunit alpha 1/2	1E−02	HPA043780
	MPP1	Q00013	Membrane palmitoylated protein 1	1E−02	HPA076675
Tissue remodelling	AGT	P01019	Angiotensinogen	3E−02	MAB3156 R&D
	CA2	P00918	Carbonic anhydrase 2	7E−03	HPA071085
	CALCA	P01258	Calcitonin related polypeptide alpha	4E−02	HPA064453
	ELANE	P08246	Neutrophil elastase	1E−02	MAB91671 R&D
	MMP2	P08253	Matrix metallopeptidase 2	2E−02	HPA001939
	MMP9	P14780	Matrix metallopeptidase 9	3E−02	HPA001238
Acute phase response	CD14	P08571	CD14 molecule	3E−02	HPA001887+
	CEBPA	P49715	CCAAT/enhancer binding protein alpha	3E−02	HPA052734
	CRP	P02741	C-reactive protein	3E−02	HPA027396
	LBP	P18428	Lipopolysaccharide binding protein	3E−02	HPA001508
	TNF	P01375	Tumour necrosis factor	5E−02	HPA077901+
Apoptosis	DAPK1	P53355	Death associated protein kinase 1	3E−02	HPA040472
	NGF	P01138	Nerve growth factor	5E−02	HPA063135
	RIPK2	O43353	Receptor interacting serine/threonine kinase 2	7E−03	HPA015764+
	TNFRSF1B	P20333	TNF receptor superfamily member 1B	2E−02	HPA004796
Axon associated	HAP1	P54257	Huntingtin associated protein 1	3E−02	HPA053019
	NEFM	P07197	Neurofilament medium	3E−02	HPA022845
Motor protein	MYL3	P08590	Myosin light chain 3	3E−02	HPA016564
	MYO15A	Q9UKN7	Myosin XVA	2E−02	HPA039770
T cell response	CD80	P33681	CD80 molecule	2E−02	HPA039851
	TNFSF13B	Q9Y275	TNF superfamily member 13b	4E−02	HPA030526
Cell cycle regulation	CDK14	O94921	Cyclin dependent kinase 14	3E−02	HPA015267
DNA replication	TIPIN	Q9BVW5	TIMELESS interacting protein	4E−02	HPA039704
Energy transduction	CKB/CKM	P12277, P06732	Creatine kinase B-type/M-type	4E−02	MAB5564 R&D

Information about the 37 proteins (number of UniProt IDs) showing differential levels between patients with severe and mild malaria. Proteins are grouped on primary biological function with displayed adjusted p-values based on the Wilcoxon rank-sum test. All these proteins showed higher levels in severe compared to mild malaria cases and include proteins, previously identified by us, as associated with severe malaria have been marked [21]. Proteins supported by multiple antibodies are marked + in the antibody column (see also Additional files 5, 11)

Dehydration is common among malaria patients and may present a potential experimental challenge as measuring relative protein levels may result in higher signals due to dehydration, regardless of the disease state. Therefore, clinically assessed hydration status of the patients were investigated to ensure that the increased protein levels in severe malaria cases compared to mild were not due to dehydration (which is more prominent among severe malaria patients). The samples were re-analysed, excluding dehydrated individuals in the two patient groups, to observe potential differences (Additional file 9). A total of 26 out of 37 proteins still showed a significant p-value (< 0.05) between the two sample groups.

Hierarchical clustering was used to visualize separation of mild and severe cases using these 37 proteins (Additional file 13). The two clusters revealed a minor separation between mild and severe disease, cluster one showing 44% of cases with mild malaria compared to 59% in cluster two. These findings highlight the complexity of the malaria disease pathogenesis.

In an effort to elucidate differences between the severe manifestations, exploration of proteins with divergent levels between the discrete groups was carried out. Five proteins present at different levels in severe cases presented with either cerebral malaria of febrile convulsions were identified (Table 3 and Fig. 7). This was evident despite the clinical overlap of these different patient groups [40]. The five identified proteins (ANK1, CD14, CDK14, ELANE, MPP1) were found to be elevated in the febrile convulsion group compared to cerebral malaria. In addition, all five proteins also showed significant differences when comparing malaria patients to controls as well as mild malaria patients to severe. The proteins were combined into a multi-protein signature ROC curve (AUC of 0.83) in order to further explore the discriminatory power of these proteins (Fig. 7f).

Primary biological process	Gene	Uniprot ID	Gene description	p-value	Antibody
Erythrocyte morphology	ANK1	P16157	Ankyrin 1	1e−02	HPA004842
	MPP1	Q00013	Membrane palmitoylated protein 1	1e−03	HPA076675
Acute phase response	CD14	P08571	CD14 molecule	4e−03	HPA001887+
Cell cycle regulation	CDK14	O94921	Cyclin dependent kinase 14	3e−03	HPA015267
ECM disassembly	ELANE	P08246	Neutrophil elastase	4e−03	MAB91671 R&D+

Information about the five proteins showing differential levels in patients with febrile convulsions compared to cerebral malaria, grouped on primary biological function, and unadjusted p-values based on Wilcoxon rank-sum test. All five proteins showed increased levels in patients with febrile convulsions compared to cerebral malaria. Proteins supported by multiple antibodies are marked + in the antibody column (see also Additional files 5, 12)

Overlap and interaction networks between human proteins of variable levels during different stages of malaria infection

Within this study, a total of 61 proteins were identified by significant p-values in any of the comparisons of different sample groups. For the 57 proteins showing differential levels in malaria cases compared to controls, 24 proteins were noted as only associated with malaria infection compared to controls. The remaining 33 proteins had also significant differences in the comparison between mild and severe malaria patients. For these 33 proteins, a rising-step trend of increment with disease severity could be noted (Fig. 8). Four proteins (CA2, CKB, CKM, DAPK1) were only found to be associated with severe compared to mild malaria.

STRING network analysis was used to identify possible connections between all the 61 proteins found as significant between any of the presented groups (Fig. 9). Protein interaction network associations were identified for the majority of the 61 proteins, with TNF and VWF as two central nodes connecting acute phase inflammatory proteins, cell migration/adhesion and tissue remodelling proteins. In addition, the erythrocyte membrane proteins EPB41L2, GYPC and MPP1, as well as ANK1 and ANK3, were shown to be associated.