Methods and Protocols

This study was a randomised, open‐label, single‐centre, descriptive study ({"type":"clinical-trial","attrs":{"text":"NCT02080559","term_id":"NCT02080559"}}NCT02080559) conducted by the Oxford Vaccine Group, University of Oxford. One hundred and eighty‐seven healthy infants aged 8–12 weeks who had not yet received their routine infant immunisations were enrolled. Infants were randomly assigned (1:1) to a control group (to receive their routine vaccines according to the UK immunisation schedule, Fig 1, i.e. PCV13 [Prevenar®, Pfizer] and DTaP‐IPV‐Hib [Pediacel®, Sanofi Pasteur MSD] at 2 and 4 months of age) or a test group (4CMenB [GSK Vaccines] plus aforementioned routine/control vaccines at 2 and 4 months of age). Control vaccines were administered in the anterolateral left thigh, while 4CMenB was administered intramuscularly in the anterolateral right thigh. Both test and control groups also received oral rotavirus vaccine (Rotarix, GSK) at 2 and 3 months; and MenC‐TT at 3 months (NeisVac‐C, Baxter Vaccines; Fig 1 and Appendix Fig S1). PCV13 and DTaP‐IPV‐Hib vaccines contained aluminium phosphate adjuvant; MenC‐TT and 4CMenB vaccines contained aluminium hydroxide adjuvant. In this study, paracetamol was not given prophylactically but post‐vaccination paracetamol/ibuprofen was administrated at the parent's/guardian's discretion.

All participants were Caucasian infants (defined as having two Caucasian parents). Ethnicity is a factor that is known to influence baseline gene expression characteristics, and contemporary gene expression analysis methods frequently exclude data from different ethnic groups (to reduce data heterogeneity) from downstream analysis (Spielman et al, 2007). Therefore, it seemed unethical to recruit participants whose data were not likely to be included in differential gene expression analysis. Blood samples for transcriptomic and proteomic analyses were taken at 4 months of age, 7 days or fewer before a second dose of 4CMenB in the test group, then either 4 h, 24 h, 3 days or 7 days post‐vaccination (Fig 1). This sub‐allocation was to minimise the number of blood tests each infant received. Sub‐allocation was dependent on parental availability; however, re‐allocation was possible if an earlier visit was unsuccessful. Blood samples for immunogenicity were taken at 5 months of age. Written informed consent was obtained from the parent or legal guardian, after a detailed explanation of the study. A national research ethics committee (South Central ‐ Oxford A, 14/SC/0077) approved this study.

Plasma samples collected 1 month after a second dose of 4CMenB in the test vaccine group were assessed for capsular group B meningococci SBA activity (plasma or serum can be used in this assay) (Borrow et al, 2005). SBA using human complement was performed at Vaccine Evaluation Unit, Public Health England, Manchester using capsular group B meningococci strain H44/76‐SL (Borrow et al, 2005). The lower limit of quantification for the SBA assays was 2; samples without detectable SBA activity were assigned an arbitrary value of 1.

The clinical haematology laboratory (Oxford University Trust Hospital NHS foundation trust) performed the full blood counts (FBC).

The clinical biochemistry laboratory (Oxford University Trust Hospital NHS foundation trust) performed the C‐reactive protein (CRP) assessment using spectrophotometry.

The body temperature of participants at the time of their 4‐month immunisations was recorded by using a wireless continuous temperature monitoring system (iButton®). The iButton® is a non‐invasive device for human skin temperature measurements, that has the capacity to record a temperature reading each minute over a 24 h period. In this study, the iButton® was applied to the infant's abdomen and used to record temperatures in the first 24 h post‐vaccination. Fever was defined as any iButton® recording of ≥ 38°C within the first 24 h of vaccination. The chi‐square test was used to compare the iButton® fever rates between the control vaccine and 4CMenB + control vaccine groups. Time to first fever in both groups was illustrated using Kaplan–Meier failure curves, and the log‐rank test was used to compare the groups. In the transcriptomic analysis, only iButton® datasets with < 90 min of missing data were considered “complete” and included in the reactogenicity analysis (to minimise misclassification).

All procedures were performed in accordance with the terms of the UK Home Office Animals Act Project License. The University of Oxford Animal Care and Ethical Review Committee approved procedures. Mice were immunised intramuscularly under general anaesthesia. Cardiac bleeds were performed under general anaesthesia followed by cervical dislocation. A mouse immunisation model was designed to recapitulate the human infant study, previously described (Sheerin et al, 2019). Groups of six 6‐ to 8‐week old female C57BL/6 mice (Harlan, UK) were immunised intramuscularly with 1/15 of the human dose (to comply with the maximum volume of vaccine allowed per mouse outlined in the project licence) for vaccine combination comparisons—4CMenB only, control group (routine vaccines only), 4CMenB + control vaccines (test group) and phosphate‐buffered saline (PBS control)—or 1/5 of the human dose for single vaccine/antigen comparisons—E. coli LPS (Invivogen, France, tlrl‐eblps) in alum (Invivogen, France, vac‐alu‐250), ultrapure E. coli LPS (Invivogen, France, tlrl‐3pelps) in alum, E. coli peptidoglycan (Invivogen, tlrl‐pgneb) adsorbed on alum, combination of alum‐adsorbed peptidoglycan and LPS, alum only or no immunisations at all (naïve control; Dataset EV1). Blood samples were taken 24 h after the second dose (day 22) and stored in RNAprotect Animal Blood Tubes (QIAGEN) containing RNA‐stabilising reagent and incubated at room temperature for 2 h to lyse blood cells.

Peripheral blood (up to 1.5 ml) was collected into a reduced volume PAXgene™ RNA stabilisation reagent (ratio of blood to PAXgene equivalent to manufacturer's specifications). Total RNA was extracted using the Blood RNA Kit (PreAnalytiX, Switzerland), using the automated protocol (QIAcube instrument, QIAGEN, Germany). The ribodepleted and globin depleted fraction were selected from the total RNA using Ribo‐Zero™ Gold rRNA removal Kit (Illumina, USA). RNA was converted to cDNA, second‐strand cDNA synthesis incorporated dUTP. The cDNA was end‐repaired, A‐tailed and adapter‐ligated and prior to amplification, samples underwent uridine digestion (to ensure strand‐specific sequencing). The prepared libraries were size selected, multiplexed and quality‐controlled before 100bp paired‐end sequencing (HiSeq4000). The multiplexing blocking strategy is available with the raw sequencing data (Gene Expression Omnibus, {"type":"entrez-geo","attrs":{"text":"GSE131929","term_id":"131929"}}GSE131929). Sequencing was conducted at the Wellcome Trust Centre for Human Genetics (Oxford, UK).

The sequencing data (fastq) files were aligned against the whole human (Homo sapiens) genome build GRCh38 (https://ccb.jhu.edu/software/hisat2/index.shtml), using HISAT2 (version 2.0.5) (Pertea et al, 2016). Gene counting was conducted using the HTSeq (version 0.9.1), utilising human gene annotation gtf (General Transfer Format) version GRCh38.88 (www.ensembl.org). To remove genes with very low counts across most libraries, only genes with an abundance of more than three counts per million in nine or more samples were carried forward. Genes assigned a “gene biotype” of ribosomal RNA (rRNA), sex chromosome genes, mitochondrial RNA or haemoglobin were excluded from downstream analysis. Human leucocyte antigen typing of RNA‐sequencing data using HISAT‐genotype (version 1.0.1‐beta) was used to check correct pairing of pre‐ and post‐vaccination samples (preprint: Kim et al, 2018).

Plasma cytokines were measured at baseline and 4 and 24 h post‐vaccination, in all participants for whom plasma was available. Prior to proteome measurement, plasma samples were thawed at room temperature and then clarified by spinning at 10,000 g for 15 min at 4ºC to remove any residual platelets and debris.

Twenty‐six cytokines were measured in multiplex using Luminex® technology (MILLIPLEX® Multiplex Assays, Merck, USA). The Human cytokine/chemokine panel (cat # HCYTOMAG‐60K‐23C) and human CVD panel 3 (cat # HCVD3MAG‐67K) were used to measure the following plasma proteins: L‐selectin, epidermal growth factor (EGF), transforming growth factor alpha (TGF‐α), granulocyte colony‐stimulating factor (G‐CSF), granulocyte–macrophage colony‐stimulating factor (GM‐CSF), fractalkine (FKN), interferon‐γ (IFNγ), GRO, interleukin‐1α (IL‐1α), interleukin‐1β (IL‐1β), interleukin‐2 (IL‐2), interleukin‐3 (IL‐3), interleukin‐4 (IL‐4), interleukin‐5 (IL‐5), interleukin‐6 (IL‐6), interleukin‐8 (IL‐8), interleukin‐10 (IL‐10), interleukin‐13 (IL‐13), interleukin‐17A (IL‐17A), interleukin‐1 receptor antagonist (IL‐1RA), IP‐10/CXCL10, tumour necrosis factor alpha (TNF‐α) and soluble CD‐40 ligand (sCD40L).

Samples and standard curves were run in duplicate. Mean fluorescence intensity (MFI) was read on a MagPix® (Luminex Corporation, USA) instrument. MFI was converted to concentration based on standard curves via the xPonent version 4.2 software (Luminex Corporation, USA) using default settings (logistic 5P weighted curve). Average concentrations were calculated for samples run in duplicate. A two‐sample Wilcoxon rank sum test was applied to compare post‐vaccination plasma protein levels between the concomitant 4CMenB and control vaccine groups.

RNA was extracted from whole blood samples using a Mouse RiboPure™‐Blood RNA Isolation Kit (Ambion, USA). Samples were depleted of α‐ and β‐globin messenger RNA (mRNA) transcripts using a GLOBINclear™ Mouse/Rat Kit, (Ambion, USA). Polyadenylated mRNA transcripts were selected by oligo (dT) beads, uridine digested, converted to complementary DNA, amplified and labelled using a TotalPrepTM‐96 RNA Amplification Kit (Illumina, USA). The prepared libraries were size selected for 75 bp fragments and multiplexed before paired‐end sequencing using an Illumina HiSeq4000 (Illumina, USA), at the Wellcome Trust Centre for Human Genetics (Oxford, UK). Reads were aligned to the mouse (M. musculus) GRCm38 reference assembly and annotation (release 93, http://www.ensembl.org) using STAR v2.6 (Dobin et al, 2013). Reads per gene were determined simultaneously by using the STAR “‐quantMode GeneCounts” option. Counts were estimated at the gene level. Lowly expressed genes (fewer than 0.5 reads per million in 2 samples) were filtered from the data.

Differential gene expression was undertaken using the R Bioconductor packages “edgeR” and “limma” (Robinson et al, 2010; McCarthy et al, 2012; R core team, 2013; Ritchie et al, 2015). RNA‐sequencing data were normalised for RNA composition using trimmed mean of M‐value (TMM) method (Robinson & Oshlack, 2010). Data were transformed using the limma “voom” function. A linear model was fitted to the data using the limma “lmFit” function using the empirical Bayes method to borrow information between genes (Ritchie et al, 2015). For the human study, paired analysis was conducted to compare pre‐ and post‐vaccination samples at each of the study time points and the statistical cut‐off for significance was set at false discovery rate (FDR) < 0.01. For the murine study, unpaired analysis was conducted to compare each vaccine group with a group of PBS‐immunised mice, and the statistical cut‐off for significance was set at FDR < 0.01.

Cell composition of infant whole blood samples was evaluated using the CIBERSORTx method (Newman et al, 2019). A filtered (as described above), non‐log space reads per kilobase million (RKPMs) sample gene matrix and a “signature” gene files (LM22; 22 immune cell types and immunoStates; 20 immune cell types (Vallania et al, 2018)) were used to deconvolute cell abundances. CIBERSORTx was run in relative mode, with B‐mode batch correction with 1,000 permutations, using cell type‐specific gene expression profiles for the analysis (Newman et al, 2019). Cell composition of mouse samples was evaluated using the ImmQuant tool (Frishberg et al, 2016). A filtered (as described above), log₂ transformed sample gene matrix and an Immunological Genome Project mouse gene expression dataset were used to infer immune cells based on signature markers (Heng et al, 2008). Predicted relative abundances were calculated relative to the PBS‐immunised mouse samples.

Gene set enrichment analysis was performed on differentially expressed genes (FDR < 0.01) employing XGR, using the functional categories of Gene Ontologies (GO; biological process, cellular component and molecular function) (Fang et al, 2016; The Gene Ontology Consortium, 2017).

Blood transcriptional module analysis was undertaken using the “tmod” R package on genes ranked by their log‐ratio (LR) value; statistical testing for module expression was evaluated using the “tmodCERNOtest” function, which is a non‐parametric test working on gene ranks (Yamaguchi et al, 2008; Weiner, 2016).

To identify baseline transcriptome predictors that can discriminate between infants that develop fever following 4CMenB vaccination, we applied sequential, iterative modelling “overnight”, SIMON, as described previously (preprint: Tomic et al, 2019). Briefly, in the first step of the SIMON analysis, all data are centred and scaled; then, the dataset is divided into training (75%) and test sets (25%). The same training and test sets are used to evaluate the performance for each of the 215 algorithms tested and to select the best performing models. To evaluate model performance and test the validity of class predictors, we implemented the two‐step procedure in SIMON. The accuracy of the predictors was first tested by 10‐fold cross‐validation on 75% of the data from the initial dataset (training set only). The cross‐validation process is repeated five times, and cumulative error rate is calculated. To prevent optimistic accuracy estimates resulting from overfitting, in the last step, each model is evaluated on the withheld test set (Kohavi, 1995). The performance of classification models was determined by calculating the area under the receiver operating characteristic curve (AUROC) for training (train AUROC) and test set (test AUROC). Out of all 215 algorithms tested, only 33 algorithms (15%) successfully built models, and only 10 models had good performance (training AUROC value above 0.7; Dataset EV13). Three of these models were overfitted, having AUROC value < 0.7 as evaluated on the test set. In the final step, SIMON calculated the contribution of each feature to the model as variable importance score (scaled to maximum value of 100), as described in the R package “caret” (Kuhn, 2017). Since variable importance score for each feature is determined by using the model information, this ultimately incorporates the correlation structure between the predictors and the importance calculation.

Similarly, a supervised machine learning algorithm was used to determine a predictive model of quantitative post‐vaccination (1 month after 4CMenB 4 months immunisation) MenB SBA titres based on pre‐vaccination gene expression profiles. For this purpose, we utilised a radial basis function kernel support vector regression (SVR), a support vector machine (SVM) that is appropriate for regression analysis (Smola & Schölkopf, 2004). Again, a partition was created with 75 and 25% of pre‐vaccination (test group only) RNA‐sequencing data in a training dataset and test dataset, respectively. Highly correlated genes (R > 0.9) within the training dataset were removed from both datasets, which were then scaled and centred. SVR was used to tune model parameters, using 25 bootstrapped iterations of the training dataset with log₁₀ MenB SBA titres as the outcome measure, utilising the R package “caret” (Kuhn, 2017). The summary metric used to select the optimal model was the root‐mean‐square error (RMSE). The final model from the training set was then used to assess its ability to predict log₁₀ MenB SBA titres in the test dataset.