PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of elifeeLifeRecent contentAbout eLifeFor authorsSign up for alerts
 
eLife. 2013; 2: e00299.
Published online Jul 16, 2013. doi:  10.7554/eLife.00299
PMCID: PMC3713456
Integrative genomic analysis of the human immune response to influenza vaccination
Luis M Franco,1,2 Kristine L Bucasas,1 Janet M Wells,3 Diane Niño,3 Xueqing Wang,4 Gladys E Zapata,4 Nancy Arden,5 Alexander Renwick,6 Peng Yu,1 John M Quarles,5 Molly S Bray,4 Robert B Couch,2,3 John W Belmont,1,4* and Chad A Shaw1
1Department of Molecular and Human Genetics, Baylor College of Medicine, Houston, United States
2Department of Medicine, Baylor College of Medicine, Houston, United States
3Department of Molecular Virology and Microbiology, Baylor College of Medicine, Houston, United States
4Children’s Nutrition Research Center, Baylor College of Medicine, Houston, United States
5Department of Microbial and Molecular Pathogenesis, Texas A&M University System Health Science Center, College Station, United States
6Interdepartmental Program in Structural and Computational Biology and Molecular Biophysics, Baylor College of Medicine, Houston, United States
Emmanouil T Dermitzakis, Reviewing editor
Emmanouil T Dermitzakis, University of Geneva Medical School, Switzerland;
*For correspondence: jbelmont/at/bcm.edu
These authors contributed equally to this work.
School of Public Health, University of Alabama at Birmingham, Birmingham, United States.
Received October 2, 2012; Accepted April 25, 2013.
Identification of the host genetic factors that contribute to variation in vaccine responsiveness may uncover important mechanisms affecting vaccine efficacy. We carried out an integrative, longitudinal study combining genetic, transcriptional, and immunologic data in humans given seasonal influenza vaccine. We identified 20 genes exhibiting a transcriptional response to vaccination, significant genotype effects on gene expression, and correlation between the transcriptional and antibody responses. The results show that variation at the level of genes involved in membrane trafficking and antigen processing significantly influences the human response to influenza vaccination. More broadly, we demonstrate that an integrative study design is an efficient alternative to existing methods for the identification of genes involved in complex traits.
Research organism: Human
Vaccines increase resistance to disease by priming the immune system to respond to specific viruses or microorganisms. By presenting a weakened (or dead) form of a pathogen, or its toxins or surface proteins, to the immune system, vaccines trigger the production of antibodies against the virus or microorganism. If a vaccinated individual then encounters the pathogen, their immune system should be able to recognize and destroy it. Many vaccines also include a secondary agent, known as an adjuvant, to further stimulate the immune response.
Influenza, an RNA virus commonly referred to as the ‘flu’, is an infectious disease that affects both birds and mammals. Seasonal epidemics occur each year affecting 2–7% of the population. According to the World Health Organization, influenza leads to nearly 5 million hospitalizations each year and causes up to half a million deaths. Vaccination is a primary strategy for the prevention of seasonal influenza, but responses to the vaccine vary markedly, partly because of variation in the genetic makeup or genotype of individuals. However, the details of how genes influence response to vaccination, and indeed susceptibility to influenza, remain unclear.
To investigate the genetic basis of variation in the immune response of healthy adults to the seasonal influenza vaccine, Franco et al. combined information about the genotypes of individuals with measurements of their gene transcription and antibody response to vaccination. They identified 20 genes that contributed to differential immune responses to the vaccine. Almost half of these encode proteins that are not specifically associated with the immune system, but have more general roles in processes such as membrane trafficking and intracellular transport.
Focusing on these genes may enable researchers to spot those individuals who are less likely to respond to a vaccine. It could also open up new avenues of research for vaccine development: rather than designing adjuvants that target known immune mechanisms, researchers should develop adjuvants that target the proteins encoded by these 20 genes.
Research organism: Human
Influenza remains one of the major threats to human health worldwide and is responsible for an estimated 250,000–500,000 deaths each year (World Health Organization, 2009). Attempts at immunization pre-dated the isolation of the virus from humans in 1933 (Smith et al., 1933) and vaccination remains the cornerstone of prevention strategies. Since 1977, strains of influenza A (H3N2), influenza A (H1N1), and influenza B have been responsible for the majority of documented human infections and trivalent vaccines are updated annually to contain the circulating strains. Animal models have demonstrated that immune responses and susceptibility to influenza infection can be strongly influenced by host genetic factors (Trammell and Toth, 2008; Srivastava et al., 2009). As with viral infection, variability in the immune response to vaccination is likely to be influenced by genotype. Accordingly, twin and sibling studies have shown heritability estimates as high as 45% for a varicella vaccine (Klein et al., 2007) and 90% for a measles vaccine (Tan et al., 2001). Studies investigating influenza vaccine immunogenicity in humans have consistently shown large inter-individual variability, but the genetic contribution to this variability remains poorly understood.
Gene expression is strongly controlled by common genetic variants (Morley et al., 2004; Stranger et al., 2007) with both broad (Bullaughey et al., 2009) and tissue-specific effects (Innocenti et al., 2011; Rotival et al., 2011), referred to as expression quantitative trait loci (eQTL). Moreover, genome-wide association studies have identified hundreds of variants associated with human disease risk that are also eQTL, implying that the mechanism by which they influence risk involves variation in transcriptional responses (Emilsson et al., 2008; Cookson et al., 2009; Naukkarinen et al., 2010; Nicolae et al., 2010; Rotival et al., 2011; Barreiro et al., 2012) Finally, integrative genomic studies in model organisms (Schadt et al., 2005; Amit et al., 2009) have demonstrated that the combination of genetic and transcriptional information can allow direct tests of causal mechanisms in controlled experiments. We hypothesized that integrating genome-wide genotype data with serial measurements of the transcriptional and humoral responses to an influenza vaccine in a clinical study could be used to identify loci that influence vaccine responsiveness and subsequent immunity to influenza in humans.
We immunized an ethnically homogeneous group of 119 healthy adult male volunteers with licensed trivalent influenza vaccine. DNA was obtained from peripheral blood and genome-wide SNP genotyping was performed. We also measured global transcript abundance in peripheral blood RNA specimens before and at three time points (days 1, 3, and 14) after vaccination. Type-specific antibody measurements (H1N1, H3N2, and FluB) were made in serum samples before and at two time points (days 14 and 28) after vaccination. An identical study was then carried out with an independent validation cohort of 128 ethnically homogeneous healthy adult female volunteers. This experimental design allowed us to search for loci that show evidence of a transcriptional response to vaccination, genetic regulation of gene expression (cis-acting eQTL), and correlation between gene expression and the magnitude of the antibody response.
Multiple genes show evidence of a transcriptional response to the vaccine and genetic regulation of expression
We performed mixed model regression analysis with SNPs located in 1-Mb intervals around each expression reporter sequence. We began by identifying SNP-transcript pairs with both significant evidence of a cis-acting eQTL and significant changes in gene expression in response to vaccination. Thresholds for local significance were initially explored, since only SNPs flanking each reporter sequence were tested for cis association. In the discovery cohort, 3229 SNP-transcript pairs, corresponding to 408 unique genes, exhibited significant genotype-expression association (genotype effect p<1 × 10−4) and concomitant evidence of a transcriptional response to the vaccine (day effect p<0.01). Of these, 2606 SNP-transcript pairs, corresponding to 256 genes, were validated in the independent cohort of female volunteers (genotype effect p<0.05 and day effect p<0.01). When more stringent thresholds were applied, 756 SNP-transcript pairs, corresponding to 114 unique genes, exhibited significant genotype-expression association (genotype effect p<5 × 10−8) and concomitant evidence of a transcriptional response to the vaccine (day effect p<0.01) in the discovery cohort. Of these, 654 SNP-transcript pairs, corresponding to 93 genes, were validated in the second cohort (genotype effect p<0.05 and day effect p<0.01). A majority of these (467 SNP-transcript pairs, corresponding to 78 unique genes) would pass equally stringent thresholds in both cohorts (genotype effect p<5 × 10−8, day effect p<0.01). A Manhattan plot of these results is presented in Figure 1. Data for the individual SNP-transcript pairs that passed equally stringent thresholds in both cohorts, including results of significance testing and gene identifiers, are provided in Table 1 via the Interactive Results Tool (which is also available to download from Zenodo and shown within Supplementary file 1).
Figure 1.
Figure 1.
Multiple genes show both a transcriptional response to the vaccine and evidence of genetic regulation of gene expression (cis-acting eQTL) in both cohorts.
At some loci, the genetic effect is enhanced or only apparent after the experimental perturbation
We hypothesized that, at some loci, the magnitude of the genetic effect could be different before and at different time points after vaccination. This type of effect, which would not be observed in a cross-sectional study design, could be directly examined with our serial expression data. We analyzed the additive effect of genotype on expression at each day in the study. Using a cis-effect significance threshold of p<1 × 10−4 in the discovery cohort and p<0.05 in the validation cohort, this analysis identified 5155 validated eQTL SNP-transcript pairs (3011 at baseline and 3417, 2496, and 3043 at days 1, 3, and 14, respectively). These SNP-transcript pairs correspond to 543 unique genes. We then identified the SNP-transcript pairs in which the expression variance explained was most strongly increased after vaccination (highest change in genetic variance explained, which we termed delta-Rg2). This analysis revealed multiple loci at which the genetic effect was either enhanced or only apparent after the experimental perturbation. An example is presented in Figure 2A, which displays local Manhattan plots for the NECAB2 locus before and 3 days after vaccination in both cohorts.
Figure 2.
Figure 2.
At some loci, the magnitude of the genetic effect changes after the experimental perturbation.
Theoretically, the observed temporal changes in the estimated genotype effect after vaccination could be driven by an increase in the effect size, a relative decrease in the variability within genotype strata, or both. We analyzed all SNP-transcript pairs for loci at which we observed both a strong cis-acting eQTL and a transcriptional response to vaccination, calculating the relative magnitude of slope and within-genotype variance between the pre-vaccination and maximal Rg2 time points. Figure 2B shows that an increase in the strength of the genotype effect (slope of the additive association) was the main driver for the observed change in Rg2, and that this amplitude change was a general feature of the loci in which we observed both a strong cis-acting eQTL and a transcriptional response to the vaccine stimulus. The delta-Rg2 values were consistent between the cohorts when evaluated by Spearman’s rank correlation analysis using all SNP-transcript pairs (Cor = 0.25, p<2 × 10−16). To select a conservative set of candidate loci based on this property for further analysis, we identified the SNP-transcript pairs that were in the top 1% of the delta-Rg2 distribution and also showed evidence of a strong cis-acting eQTL (genotype effect p<5 × 10−8), in both cohorts. Data for the resulting set of 146 SNP-transcript pairs, including Rg2 values, are provided in Table 2 via the Interactive Results Tool (which is also available to download from Zenodo and shown within Supplementary file 1).
Content analysis shows enrichment for genes involved in membrane trafficking, antigen processing, and antigen presentation
Of the 78 genes that had the strongest validated evidence of a genotype effect and a transcriptional response to the vaccine, 14 were also in the list of 34 genes with the strongest evidence of an increase in the magnitude of the genetic effect after vaccination. Content analysis on the union of the two sets (98 genes) showed significant enrichment for genes involved in antigen processing and presentation, cytotoxic T-lymphocyte-mediated apoptosis of target cells, dendritic cell maturation and function, and membrane trafficking (Figure 3).
Figure 3.
Figure 3.
Content analysis shows enrichment for genes involved in membrane trafficking, antigen processing, and antigen presentation.
Integration of genotype, expression, and antibody titer data identifies 20 genes with the strongest evidence for genetic variation influencing the humoral immune response to influenza vaccination
We and others have shown that for some transcripts there is significant correlation between the magnitude of the transcriptional and antibody responses to the vaccine stimulus (Zhu et al., 2010; Bucasas et al., 2011; Nakaya et al., 2011) In a combined analysis of the two cohorts in the present study, 301 transcripts were found to correlate with the magnitude of the antibody response (Figure 4). Additional details of these 301 transcripts, including correlation coefficients and days of maximum correlation, are provided in Table 3 via the Interactive Results Tool (which is also available to download from Zenodo and shown within Supplementary file 1). We imposed an additional selection threshold based on this correlation, and identified 20 genes that show evidence of significant genotype-expression association (genotype effect p<5 × 10−8), a significant correlation between the transcriptional and antibody responses (expression-antibody effect p<0.05), and either a transcriptional response to the vaccine (day effect p<0.01) or evidence of a change in the magnitude of the genetic effect after vaccination (top 1% of the delta-Rg2 distribution) in the two independent cohorts. These loci have the strongest evidence of genetic variation influencing the immune response to the vaccine, and include TAP2, SNX29, FGD2, NAPSA, NAPSB, GM2A, C1orf85, JUP, FBLN5, CHST13, DIP2A, PAM, D4S234E, C3AR1, HERC2, LST1, LRRC37A4, OAS1, RPL14, and DYNLT1. Remarkably, seven of these encode proteins involved in intracellular antigen transport and processing (Figure 5).
Figure 4.
Figure 4.
Gene expression at specific loci correlates with the antibody response to vaccination.
Figure 5.
Figure 5.
Genetic variation in intracellular antigen transport and processing influences the human immune response to influenza vaccination.
We determined genetic associations to the antibody response using 137 eQTL SNPs from these 20 loci. The quantile-quantile plot from the association tests performed on these SNPs shows marked deviation from the empirical null distribution for QTL associations (Figure 6), supporting the idea that these loci are enriched for true genetic associations.
Figure 6.
Figure 6.
SNPs at the 20 loci identified show evidence of association with the antibody response to the vaccine.
The study design permits causal and reactive model analyses
We explored three types of associations in our work: genotype to gene expression (eQTL), gene expression to antibody titer, and genotype to antibody titer (QTL). We now considered alternative models for the relationships between these distinct types of association (Figure 7A), and we evaluated our data to determine which of these alternatives appears most consistent with our observations. The alternative models considered were: (i) genotype association with gene expression is independent of genotype association or trends of association with antibody response (independent model); (ii) genotype association or trends of association with antibody response are mediated by gene expression patterns that are strongly correlated with genotype (causal model); and (iii) genotype associations to antibody response are not mediated by expression, but instead gene expression patterns are a response to the antibody trait or its early correlates (reactive model). To perform a comparative analysis of these alternatives we extended the framework for causal modeling (Pearl, 2010) in eQTL data recently developed by others (Millstein et al., 2009) and applied the method to our time-course gene expression study. We used the 137 eQTL SNP-transcript pairs from the 20 loci with the strongest evidence of genetic variation influencing the immune response to the vaccine, as described above. We found that the patterns in the data trend toward the causal model compared to the reactive model (Figure 7B), but a power analysis based on the distribution of the empirical effect sizes of our observed associations also indicates that our sample size is too modest to support definitive conclusions (Figure 7C).
Figure 7.
Figure 7.
The study design permits causal and reactive model analyses.
The results provide an unbiased integrative survey of the genetic and transcriptional components of the humoral immune response to influenza vaccination in humans. They suggest that variation at the level of genes involved in antigen processing and intracellular trafficking is an important determinant of vaccine immunogenicity. Even in healthy, young individuals, there are a significant number of people who do not develop a protective antibody response after influenza vaccination. If these individuals could be identified prior to vaccination, modifications to the type or dose of vaccine could be attempted, with the goal of reducing the number of unprotected vaccinated individuals. The genes identified in this study as playing a role in variation in the humoral response to vaccination would be a logical starting point for the development of DNA- or RNA-based predictive biomarkers. Prospective evaluation of such biomarkers would be the next step towards clinical implementation.
Understanding the mechanisms that underlie variation in response to the vaccine may also direct modification of factors that enhance the response. Most of the efforts to date have focused on vaccine adjuvants that activate known immunologic mechanisms. Surprisingly, many of the genes identified in this study encode proteins that are not specifically immune but play a more general role in membrane trafficking and intracellular transport. Interventions aimed at increasing vaccine antigen affinity to these proteins or altering their intracellular concentrations could represent new avenues in vaccine development.
More broadly, the results demonstrate that a longitudinal, integrative genomic analysis study design, applied to a clinical intervention, is an efficient alternative to cross-sectional methods for the identification of genes involved in medically relevant complex traits. By making repeated measurements on the same individual over time after a controlled experimental perturbation, we were able to account for individual variation in a way that would not have been possible otherwise. The dynamic nature of the measurements also allowed us to uncover genetic effects that are either enhanced by or only evident after the experimental perturbation. The specificity of gene identification in this study emerges from the genome’s acute response to the perturbation, which cannot be assessed by a cross-sectional eQTL analysis or a genome-wide association study. This approach could be used for a broad variety of medically important problems whenever there is the opportunity to test a well-controlled intervention such as drug, dietary, or vaccine responses.
Several limitations of the study are worth noting. First, we studied two samples of healthy young adults, thereby excluding the segments of the population that are most likely to have a poor response to influenza vaccination: children, the elderly, and individuals with severe illnesses. Second, in order to minimize the risk of false associations related to population stratification, we studied an ethnically homogeneous group of individuals. Third, while an interesting aspect of our study design is that it could open the door for direct comparisons of causal and reactive models, the sample size in this study was not sufficient to establish whether or not there is a causal relationship between the loci for which an association was identified and the antibody response to the vaccine. Finally, while antibody titers have historically been used to evaluate vaccine responsiveness, it is clear that they do not capture the complexity of the human immune response to vaccination. Additional studies would be necessary to determine whether the genes identified are also related to variation in influenza vaccine responses in groups other than the one chosen for this study, whether there is a causal relationship between these genes and the antibody response, or whether they also influence the cell-mediated immune response to the vaccine.
Visual summary
A visual representation of the study design, the resulting data sets, and the integrative analysis scheme, is presented in Figure 8.
Figure 8.
Figure 8.
Study design and integrative analysis scheme.
Study subjects
Healthy volunteers ages 18 to 40 years were enrolled. Individuals who were known to have received an influenza vaccine in the previous 3 years or who had signs or symptoms of an active infection at the time of enrollment were excluded. To minimize false-positive results related to population stratification, enrollment was limited to individuals of self-reported Caucasian ancestry. Enrollment, vaccination and sample collections were conducted at a university campus. The initial (discovery) cohort was restricted to males. The validation cohort was enrolled approximately 18 months after the initial cohort and was restricted to females. The protocol was approved by the institutional review boards of all participating institutions. Informed consent was obtained from each subject prior to enrollment.
Vaccine
Study participants were immunized on day 0. Those enrolled in the initial cohort received the 2008–2009 inactivated trivalent influenza vaccine (A/Brisbane/59/2007[H1N1], A/Brisbane/10/2007[H3N2], B/Florida/4/2006; Sanofi-Pasteur, Lyon, France). The validation cohort received the 2009–2010 vaccine, which came from the same manufacturer and included (A/Brisbane/59/2007), (A/Brisbane/10/2007[H3N2]), and (B/Brisbane/60/2008) strains.
DNA samples
Whole blood samples (7 ml) for DNA purification were collected in Vacutainer acid citrate dextrose (ACD) tubes (Beckton-Dickinson, Franklin Lakes, NJ), on day 1 after vaccination. DNA was purified using Qiagen Gentra Puregene Blood Kits (Qiagen Sciences, Germantown, MD). Quantitation and quality control were performed with a NanoDrop-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA) and using the Quant-iT PicoGreen dsDNA Reagent (Life Technologies, Carlsbad, CA) in a Tecan GENios microplate reader (Tecan Group, Mannendorf, Switzerland).
RNA samples
Peripheral blood samples for RNA purification were obtained immediately before (day 0) and on days 1, 3, and 14 after vaccination. To minimize changes in gene expression induced by sample handling and processing, whole-blood samples (2.5 ml) were collected in PAXgene RNA stabilization tubes (Qiagen Inc., Valencia, CA) and frozen at −80°C. RNA purification was performed using the PAXgene Blood RNA system (Qiagen Inc.). All RNA samples from an individual were consistently purified in the same batch. Spectrophotometry (NanoDrop-1000 Spectrophotometer; Thermo Fisher Scientific) and microfluidic electrophoresis (Experion Automated Electrophoresis System; Bio-Rad Laboratories, Hercules, CA) were used for quality control.
Whole-genome genotyping
Genotyping was performed on Illumina HumanOmniExpress microarrays (Illumina, Inc., San Diego, CA) following the manufacturer’s instructions. The arrays include 730,525 SNP markers. Basic quality control of the genotyping data was performed on GenomeStudio software, version 2010 (Illumina, Inc.). All microarrays had call rates >0.99.
Measurement of global transcript abundance
In vitro transcription was performed using Ambion Illumina TotalPrep RNA Amplification Kits (Applied Biosystems/Ambion, Austin, TX). cRNA was hybridized onto Illumina HumanHT-12v3 or HumanHT-12v4 Expression BeadChips (Illumina, Inc.), following the manufacturer’s protocol. All samples for a given individual were processed on the same slide. The arrays have 48,742 (v3) and 47,301 (v4) reporters, representing approximately 25,000 genes and non-annotated gene candidates.
Serum samples and antibody titer measurements
Whole blood (10 ml) was collected in Vacutainer Serum Separator Tubes (Beckton–Dickinson). Serum was separated by centrifugation prior to storage at −20°C. Hemagglutination inhibition (HAI) tests were performed as previously described (Dowdle et al., 1979), except for a starting serum dilution of 1:4 and the use of turkey red blood cells. HAI test antigens were allantoic fluid harvests from infected embryonated hen’s eggs (whole-virus antigens). Neutralizing antibody tests were performed as previously described (Frank et al., 1980) except that hamster serum was not included. Test strains were the same as those used in the vaccine.
Antibody titer data analysis
HAI and neutralizing antibody titers were measured for each viral antigen included in the vaccine. For all antigens, the change in antibody titer post-vaccination is known to be negatively correlated with the pre-vaccination titer. The responses to the individual viral antigens were correlated within individuals. For each individual, we computed a Titer Response Index (TRI), as previously described (Bucasas et al., 2011). The TRI characterizes the magnitude of an individual’s antibody rise accounting for their pre-vaccination titer and integrating the responses across the three vaccine components to provide an overall measure of vaccine responsiveness.
Expression data processing and quality control
Raw signal intensity data from all individuals and all time points in the discovery cohort were first processed in a single batch. Background adjustment, variance stabilization transformation (Lin et al., 2008) and robust spline normalization were performed using the lumi package (Du et al., 2008) in R (R Development Core Team, 2009). Eight individuals had missing expression data from two or more time points and were excluded. We required a detection p-value of <0.01 in at least 80% of the samples for a transcript to be considered detected. We also aligned the entire set of expression reporter sequences to the human genome reference sequence (Build 36 [March 2006]/hg18) by applying the BLAT algorithm in BlatSuite34 software (Kent, 2002), and excluded any reporters that did not map or mapped to more than one region. Using these two thresholds for the data in the discovery cohort, the final data set included 9809 detected transcripts. This data set was then used for eQTL analysis in the discovery cohort.
Once data generation for the validation cohort was completed, the expression microarray signal intensity data from all individuals and all time points in that cohort were processed in a single batch. Background adjustment, variance stabilization transformation, robust spline normalization, and the application of detection thresholds were performed identically to the discovery cohort. Five individuals had missing expression data from two or more time points and were excluded. Because two different array versions were used (HT12-v3 and HT12-v4), unique reporter identifiers (ProbeID and nuID) for the 9809 reporters selected in the discovery cohort were used to subset the data from the validation cohort. This data set was then used for eQTL analysis in the validation cohort.
The analysis of correlation between gene expression and antibody titer in the discovery cohort was previously published (Bucasas et al., 2011). As part of the integrative genomic analysis described in the present study, we performed a similar analysis of this expression/titer correlation, but included expression data from both cohorts. For this purpose, the two data sets described above were combined, and an additional quantile normalization step was performed to account for batch effects between cohorts.
Genotyping data processing and quality control
Array quality was initially assessed using GenomeStudio software (Illumina, Inc.). Default algorithms were used to normalize, generate SNP clusters, and make genotype calls. SNPs with minor allele frequency (MAF) <0.05 and Hardy-Weinberg Equilibrium (HWE) χ2 < 1 × 10−7 were removed prior to analysis. For local eQTL mapping, we restricted the analyses to SNPs within 1 Mb (500 kb upstream or downstream) of each unambiguously mapped expression reporter sequence. We used the genotype information to reliably establish that the individuals in each cohort were indeed unrelated and ethnically matched. For this, we estimated pairwise identity-by-descent (IBD) metrics on the basis of identity-by-state (IBS) information from their genome-wide SNP genotyping results, followed by multidimensional scaling (MDS) analysis on the resulting matrix of pairwise distances, using PLINK (Purcell et al., 2007). The first three components of the MDS analysis on our study samples and the founders from six HapMap populations were plotted with R package Scatterplot3D (Ligges and Machler, 2003) and are shown in Figure 9. As expected, our study samples from both cohorts cluster with the HapMap CEU population. In the male cohort, we identified four pairs of individuals with cryptic or undisclosed familial relationships (pi-hat ≥ 0.125). 12 outliers were identified by component 2; two of the four pairs of cryptic or undisclosed relatives were among the outliers. One individual from each related pair and all remaining outliers were removed prior to analysis. In the female cohort, we identified one pair of individuals with cryptic or undisclosed familial relationship (pi-hat ≥ 0.125). 11 outliers were identified by component 2 and one by component 1. One individual from the related pair and all outliers were removed prior to analysis.
Figure 9.
Figure 9.
Study samples cluster with the HapMap CEU population.
Integrative genomic analysis
Analyses were carried out separately for each cohort and in the combined data set. Because we have repeated gene expression measurements on the same individuals over a series of time points, we performed a random effects linear model analysis. For each SNP-transcript pair, we fit a model with terms for day, additive effects for genotype and day–genotype interactions, and a random effect for each person:
equation m1
Where Y = observed expression value, i = day (0, 1, 3, 14), j = individual, D = day effect, g = genotype effect, P = person effect (random), and e = random error. We fit this model using the R package nlme and the estimation method of the same name. The day effect (D) is a shift term which considers time-ordered change in expression. The genotype effect (g) is a main effect for genotype that considers the expression differences between genotypes to be additively related to the number of copies of a reference allele. The day–genotype interaction allows the additive effect to vary by day (see Figure 10 for a visual representation of the analysis scheme). The person effect (P) is modeled as a random effect; this term helps to efficiently account for the fact that the same individual appears in each time point, without overfitting the expression model to the particular sample ascertained in our study. We explored different parametrizations, such as allowing the person effect to vary by day, and found that the simplest model with a single term for each person was effective. This person effect can account for inter-individual or batch differences which are stable within individuals for the duration of the study.
Figure 10.
Figure 10.
Genetic and transcriptional analysis on a prospective cohort.
In addition to the random effects linear model, we analyzed each expression trait with respect to genotype at each time point using
equation m2
where Yj denotes the matrix of expression traits of individuals j at a given time point (days 0, 1, 3, 14) and Gi,k is a matrix of genotypes for individuals j at SNP locus k such that each element is assigned 0, 1 or 2 according to the number of minor alleles at the kth locus of the jth individual. This allowed estimates of the proportion of expression variance explained by individual genotypes from the model’s coefficient of determination (R2g). This model permitted a more detailed examination of the changes in the strength of the genotypic association with expression at each time point. To detect SNP-transcript pairs where the magnitude of the genotype effect varied after immunization, we took the difference in R2g measures (ΔR2g) from pairwise linear models fitted after (days 1, 3 and 14) and before vaccination (day 0). We then took the top 1% ΔR2g values between day 0 and a later time point. This cis-acting eQTL subset represents the loci that show significant changes in the genetic effect as a result of the vaccine.
Correlation of transcripts with the humoral immune response
The titer response index was treated as a dependent variable and modeled as a linear combination of expression values for each transcript across the time course. A single F-statistic p-value was determined to evaluate the explanatory value of the expression data using all days. Separate evaluations were made in each cohort and in the combined data. The 301 genes displayed in Figure 4B had F-statistic p-values in the combined data <0.01 and absolute values of the maximum average correlation between expression and titer response in the two cohorts >0.15.
Content analyses
Enriched biological and functional pathways were analyzed using DAVID Bioinformatics Resources (Dennis et al., 2003), and Ingenuity Pathway Analysis (IPA) software.
Causal and reactive model analyses
To evaluate the relationship between the associations identified in our study, we extended a recently published analysis framework for causal modeling in eQTL data (Millstein et al., 2009). The relationship is modeled as
equation m3
where T is trait, E is gene expression, and G is genotype at a locus. In our experiment, the term E includes a separate value for each day. The term G has a separate effect on the trait under the independent null model, while under the causal alternative, where it acts through gene expression, its influence is entirely transmitted through E. The statistic used to assess the influence of G is the partial F reported by the ‘aov’ function in R (version 2.15.1). The analyses utilize an equivalence test approach, where the null hypothesis is that the second independent variable has an effect on the response conditional on the first variable, while the alternative is complete co-linearity in the associations. The null distribution of the F statistic–which under parametric assumptions would have a non-central F distribution—is derived using a permutation test procedure, as follows: The relationship between E and T is decoupled by regressing E on G, permuting the resulting residuals, then adding the permuted residuals to the predicted values to arrive at E*, which is independent of T but maintains the marginal variance and G-correlation of E. This procedure enforces the assumption of the independent model, while leaving other properties unchanged. The partial F statistic for G in
equation m4
represents a sample from the distribution under the independent null hypothesis. This process, repeated 1000 times, provided a null distribution. The p-value for the observed F is the proportion of the null distribution that is below the observed value.
The reactive model is evaluated using the same approach, but with the expression and trait terms swapped, that is
equation m5
This scheme tests the extent to which the trait value mediates the relationship between genotype and gene expression. In this analysis, the permutation step decouples the expression-to-trait relationship from the eQTL association and asks if the eQTL association disappears after accounting for the expression-to-trait association.
Power analyses for these models apply the same scheme to simulated data. The term G is simulated as binomial, while E and T are normal, with parameters chosen to match the minor allele frequency, variance, and correlations observed in our data set. The power to reject the false (independent) model is shown using a threshold of p<0.05. The power properties for the tests used to evaluate both the reactive and causal models are equivalent because the model structures and r2 relationships are symmetric between the two alternatives.
Acknowledgements
The clinical trials were conducted at the Beutel Health Center at Texas A&M University with the strong support of Dr Martha Dannenbaum, Scott Draper, and Richard Darnell. Mr Uri Kelman of Kelman Design contributed to the generation of Figure 5.
Funding Statement
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Funding Information
This paper was supported by the following grants:
  • US National Institutes of Health - National Institute of Allergy and Infectious Diseases NO1-AI-030039 to Robert B Couch, John W Belmont.
  • US National Institutes of Health - National Institute of Allergy and Infectious Diseases 1K23AI087821-01 to Luis M Franco.
  • US National Institutes of Health - Ruth L Kirschstein National Research Service Award T32 GM07526 to Luis M Franco.
  • US National Institutes of Health - Ruth L Kirschstein National Research Service Award F31AI071372-03 to Kristine L Bucasas.
Additional information
Competing interests
The authors declare that no competing interests exist.
Author contributions
LMF, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article.
KLB, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article.
JMW, Acquisition of data.
DN, Acquisition of data.
XW, Acquisition of data.
GEZ, Acquisition of data.
NA, Acquisition of data.
JMQ, Acquisition of data.
AR, Analysis and interpretation of data.
PY, Analysis and interpretation of data.
MSB, Acquisition of data, Drafting or revising the article.
RBC, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article.
JWB, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article.
CAS, Analysis and interpretation of data, Drafting or revising the article.
Ethics
Human subjects: The protocol was approved by the institutional review boards of all participating institutions. BCM IRB identifiers: H-28980, H-25352. Informed consent and consent to publish was obtained from each subject prior to enrollment.
Additional files
10.7554/eLife.00299.013
Supplementary file 1.: Compilation of the tables and figures (XLS).
This is a static version of the Interactive Results Tool, which is also available to download from Zenodo (see major datasets).
Major datasets
The following datasets were generated:
Franco LM, Bucasas KL, Wells JM, Niño D, Wang X, Zapata GE, , et al., 2013, Time series of global gene expression after trivalent influenza vaccination in humans, GEO SuperSeries GSE48024; http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE48024, Publicly available at GEO (http://www.ncbi.nlm.nih.gov/geo/).
Franco LM, Shaw CA, Belmont J, 2013, Time series of global gene expression after trivalent influenza vaccination in humans (male cohort), GSE48018 (SubSeries of GSE48024); http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE48018, Publicly available at GEO (http://www.ncbi.nlm.nih.gov/geo).
Franco LM, Shaw CA, Belmont J, 2013, Time series of global gene expression after trivalent influenza vaccination in humans (female cohort), GSE48023 (SubSeries of GSE48024); http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE48023, Publicly available at GEO (http://www.ncbi.nlm.nih.gov/geo).
Franco LM, Bucasas KL, Wells JM, Niño D, Wang X, Zapata GE, , et al., 2013, ”Interactive” version of data associated with the eLife paper “Integrative genomic analysis of the human immune response to influenza vaccination”, 10.5281/zenodo.6960; https://zenodo.org/record/6960?ln=en#.UdPjW1OItzY, Publicly available at Zenodo (https://zenodo.org/).
Franco LM, Bucasas KL, Wells JM, Niño D, Wang X, Zapata GE, , et al., 2013, Adult influenza vaccine genetics, phs000635.v1.p1; http://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000635.v1.p1, Publicly available at dbGAP (http://www.ncbi.nlm.nih.gov/gap).
  • Amit I, Garber M, Chevrier N, Leite AP, Donner Y, Eisenhaure T, et al. 2009. Unbiased reconstruction of a mammalian transcriptional network mediating pathogen responses. Science 326:257–63. doi: 10.1126/science.1179050. [PMC free article] [PubMed] [Cross Ref]
  • Barreiro LB, Tailleux L, Pai AA, Gicquel B, Marioni JC, Gilad Y. 2012. Deciphering the genetic architecture of variation in the immune response to Mycobacterium tuberculosis infection. Proc Natl Acad Sci USA 109:1204–9. doi: 10.1073/pnas.1115761109. [PubMed] [Cross Ref]
  • Bucasas KL, Franco LM, Shaw CA, Bray MS, Wells JM, Nino D, et al. 2011. Early patterns of gene expression correlate with the humoral immune response to influenza vaccination in humans. J Infect Dis 203:921–9. doi: 10.1093/infdis/jiq156. [PMC free article] [PubMed] [Cross Ref]
  • Bullaughey K, Chavarria CI, Coop G, Gilad Y. 2009. Expression quantitative trait loci detected in cell lines are often present in primary tissues. Hum Mol Genet 18:4296–303. doi: 10.1093/hmg/ddp382. [PMC free article] [PubMed] [Cross Ref]
  • Cookson W, Liang L, Abecasis G, Moffatt M, Lathrop M. 2009. Mapping complex disease traits with global gene expression. Nat Rev Genet 10:184–94. doi: 10.1038/nrg2537. [PubMed] [Cross Ref]
  • Dennis G, Jnr., Sherman BT, Hosack DA, Yang J, Gao W, Lane HC, et al. 2003. DAVID: database for annotation, visualization, and integrated discovery. Genome Biol 4:P3. doi: 10.1186/gb-2003-4-5-p3. [PMC free article] [PubMed] [Cross Ref]
  • Dowdle WR, Kendal AP, Noble GR. 1979. Influenza viruses. In: Lennette EH, Schmidt NJ, editors. editors. Diagnostic procedures for viral, rickettsial and Chlamydial infections. 5th ed Washington, D.C.: American Public Health Association; p. 603–5.
  • Du P, Kibbe WA, Lin SM. 2008. lumi: a pipeline for processing Illumina microarray. Bioinformatics 24:1547–8. doi: 10.1093/bioinformatics/btn224. [PubMed] [Cross Ref]
  • Emilsson V, Thorleifsson G, Zhang B, Leonardson AS, Zink F, Zhu J, et al. 2008. Genetics of gene expression and its effect on disease. Nature 452:423–8. doi: 10.1038/nature06758. [PubMed] [Cross Ref]
  • Frank AL, Puck J, Hughes BJ, Cate TR. 1980. Microneutralization test for influenza A and B and parainfluenza 1 and 2 viruses that uses continuous cell lines and fresh serum enhancement. J Clin Microbiol 12:426–32. [PMC free article] [PubMed]
  • Innocenti F, Cooper GM, Stanaway IB, Gamazon ER, Smith JD, Mirkov S, et al. 2011. Identification, replication, and functional fine-mapping of expression quantitative trait loci in primary human liver tissue. PLOS Genet 7:e1002078. doi: 10.1371/journal.pgen.1002078. [PMC free article] [PubMed] [Cross Ref]
  • Kent WJ. 2002. BLAT–the BLAST-like alignment tool. Genome Res 12:656–64. doi: 10.1101/gr.229202. [PubMed] [Cross Ref]
  • Klein NP, Fireman B, Enright A, Ray P, Black S, Dekker CL. 2007. A role for genetics in the immune response to the varicella vaccine. Pediatr Infect Dis J 26:300–5. doi: 10.1097/01.inf.0000257454.74513.07. [PubMed] [Cross Ref]
  • Ligges U, Machler M. 2003. Scatterplot3d - an R Package for Visualizing Multivariate data. J Stat Softw 8:1–20.
  • Lin SM, Du P, Huber W, Kibbe WA. 2008. Model-based variance-stabilizing transformation for Illumina microarray data. Nucleic Acids Res 36:e11. doi: 10.1093/nar/gkm1075. [PMC free article] [PubMed] [Cross Ref]
  • Millstein J, Zhang B, Zhu J, Schadt EE. 2009. Disentangling molecular relationships with a causal inference test. BMC Genet 10:23. doi: 10.1186/1471-2156-10-23. [PMC free article] [PubMed] [Cross Ref]
  • Morley M, Molony CM, Weber TM, Devlin JL, Ewens KG, Spielman RS, et al. 2004. Genetic analysis of genome-wide variation in human gene expression. Nature 430:743–7. doi: 10.1038/nature02797. [PMC free article] [PubMed] [Cross Ref]
  • Nakaya HI, Wrammert J, Lee EK, Racioppi L, Marie-Kunze S, Haining WN, et al. 2011. Systems biology of vaccination for seasonal influenza in humans. Nat Immunol 12:786–95. doi: 10.1038/ni.2067. [PMC free article] [PubMed] [Cross Ref]
  • Naukkarinen J, Surakka I, Pietilainen KH, Rissanen A, Salomaa V, Ripatti S, et al. 2010. Use of genome-wide expression data to mine the “Gray Zone” of GWA studies leads to novel candidate obesity genes. PLOS Genet 6:e1000976. doi: 10.1371/journal.pgen.1000976. [PMC free article] [PubMed] [Cross Ref]
  • Nicolae DL, Gamazon E, Zhang W, Duan S, Dolan ME, Cox NJ. 2010. Trait-associated SNPs are more likely to be eQTLs: annotation to enhance discovery from GWAS. PLOS Genet 6:e1000888. doi: 10.1371/journal.pgen.1000888. [PMC free article] [PubMed] [Cross Ref]
  • Pearl J. 2010. An introduction to causal inference. Int J Biostat 6:1–59. doi: 10.2202/1557-4679.1203. [PMC free article] [PubMed] [Cross Ref]
  • Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. 2007. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet 81:559–75. doi: 10.1086/519795. [PubMed] [Cross Ref]
  • R Development Core Team 2009. R: a language and environment for statistical computing. 2.9 ed.Vienna, Austria: R Foundation for statistical computing. doi: 10.1007/s10985-007-9065-x. [Cross Ref]
  • Rotival M, Zeller T, Wild PS, Maouche S, Szymczak S, Schillert A, et al. 2011. Integrating genome-wide genetic variations and monocyte expression data reveals trans-regulated gene modules in humans. PLOS Genet 7:e1002367. doi: 10.1371/journal.pgen.1002367. [PMC free article] [PubMed] [Cross Ref]
  • Schadt EE, Lamb J, Yang X, Zhu J, Edwards S, Guhathakurta D, et al. 2005. An integrative genomics approach to infer causal associations between gene expression and disease. Nat Genet 37:710–7. doi: 10.1038/ng1589. [PMC free article] [PubMed] [Cross Ref]
  • Smith W, Andrewes CH, Laidlaw PP. 1933. A virus obtained from influenza patients. Lancet 222:66–8. doi: 10.1016/S0140-6736(00)78541-2. [Cross Ref]
  • Srivastava B, Blazejewska P, Hessmann M, Bruder D, Geffers R, Mauel S, et al. 2009. Host genetic background strongly influences the response to influenza a virus infections. PLOS ONE 4:e4857. doi: 10.1371/journal.pone.0004857. [PMC free article] [PubMed] [Cross Ref]
  • Stranger BE, Forrest MS, Dunning M, Ingle CE, Beazley C, Thorne N, et al. 2007. Relative impact of nucleotide and copy number variation on gene expression phenotypes. Science 315:848–53. doi: 10.1126/science.1136678. [PMC free article] [PubMed] [Cross Ref]
  • Tan PL, Jacobson RM, Poland GA, Jacobsen SJ, Pankratz VS. 2001. Twin studies of immunogenicity–determining the genetic contribution to vaccine failure. Vaccine 19:2434–9. doi: 10.1016/S0264-410X(00)00468-0. [PubMed] [Cross Ref]
  • Trammell RA, Toth LA. 2008. Genetic susceptibility and resistance to influenza infection and disease in humans and mice. Expert Rev Mol Diagn 8:515–29. doi: 10.1586/14737159.8.4.515. [PubMed] [Cross Ref]
  • World Health Organization 2009. Seasonal Influenza. World Health Organization Factsheet No. 211. Available from: http://www.who.int/mediacentre/factsheets/fs211/en/index.html.
  • Zhu W, Higgs BW, Morehouse C, Streicher K, Ambrose CS, Woo J, et al. 2010. A whole genome transcriptional analysis of the early immune response induced by live attenuated and inactivated influenza vaccines in young children. Vaccine 28:2865–76. doi: 10.1016/j.vaccine.2010.01.060. [PubMed] [Cross Ref]

Decision letter
Emmanouil T Dermitzakis, Reviewing editor
Emmanouil T Dermitzakis, University of Geneva Medical School, Switzerland;
 
eLife posts the editorial decision letter and author response on a selection of the published articles (subject to the approval of the authors). An edited version of the letter sent to the authors after peer review is shown, indicating the substantive concerns or comments; minor concerns are not usually shown. Reviewers have the opportunity to discuss the decision before the letter is sent (see review process). Similarly, the author response typically shows only responses to the major concerns raised by the reviewers.
Thank you for choosing to send your work entitled “Integrative genomic analysis of the human immune response to influenza vaccination” for consideration at eLife. Your article has been evaluated by a Senior editor and 2 reviewers, one of whom, Manolis Dermitzakis, is a member of our Board of Reviewing Editors.
The Reviewing editor and the other reviewer discussed their comments before we reached this decision, and the Reviewing editor has assembled the following comments based on the reviewers' reports.
The manuscript describes a project that aims to integrate genetic, gene expression, and vaccination data, together with immunological traits. The design is innovative and promises interesting discoveries, some of which are described in the submission. Some key concerns were raised that need to be addressed before a final decision can be made:
1) The eQTL discovery was performed in a cohort comprised of males only and the replication was performed in a cohort comprised of females only. What was the purpose of this design? The authors do not address the issue that discovery and replication are confounded by sex.
2) We suspect that an interaction term in the models would actually be more efficient in detecting genetic effects modified by the time point after vaccination. Why was this not used (instead a simple comparison among time points was made)?
3) The fact that there is enrichment in the QQ plots on Figure 6A does not provide evidence about causality of SNPs to vaccine response but simply of correlation among the two effects. More complex models need to be implemented to show causality.
4) The study design used two different array versions: Illumina HT-12v3 and HT-12v4. The authors should describe how they combined data from these arrays, particularly for transcripts where the probe sequence differs between versions.
5) The expression data processing should be described in more detail to make it clear how normalization, and so on, was performed given the different time points, array versions, and discovery/replication cohorts.
6) The authors mention that the person effect within their model would account for differences in cell populations among individuals, given that these are stable within individuals during the sampling timeframe. No data or references are provided to support this claim. Immune cell relative proportions can change rapidly.
7) The authors use the delta r2 metric to identify eQTLs that differ between time points or in response to the vaccination. It seems that the authors could have performed the analysis using the delta of transcription(time point x)−transcription(baseline) as a trait. This normalizes to each individual's starting point. Furthermore, a difference of r2 between time points can have multiple causes (as the authors point out). It seems to us that what they are most concerned with is cases where the effect size (beta) of a given eQTL differs between time points, or where the actual response (transcription time point x−transcription baseline) differs between individuals. The authors reasoning behind these choices should be made clearer.
8) For the trans-eQTL analysis, the authors consider as genome-wide significant any SNP-transcript pair with p<5×10-8. Indeed for a GWAS of common variation, this may be appropriate for a single trait, but the authors should account for testing nearly 10K transcripts.
[After resubmission, the editors also asked for the following comments to be addressed prior to acceptance.]
A) The revision is not satisfactory for point 7, above, on the delta r2. The case is still not convincing because only 146 of the 541 cases were validated. It is not clear why the difference in r2 is necessarily interpreted as a significant change when only 30% of the cases are validated.
B) The modelling as previously discussed in point 3 is somewhat weak and causality is not as well demonstrated, as the authors want to claim. The fact that the association with the residuals is lower than with the original values is not direct but indirect evidence that there are casual relationships. Also, how were the “random” SNPs chosen? Establishing a relevant null is important here. Also, direct causal models would help in establishing the relationships rather than the indirect approach the authors use.

 
1) The eQTL discovery was performed in a cohort comprised of males only and the replication was performed in a cohort comprised of females only. What was the purpose of this design? The authors do not address the issue that discovery and replication are confounded by sex.
Many loci show sex-specific effects at the level of gene expression. The study design for the discovery cohort, therefore, sought to minimize the effects of suspected variables, such as previous vaccine exposure, ethnicity, age, and sex. The follow-up study was designed to fill the resulting gap in information about responses to the seasonal influenza vaccine in women. We agree that the consequence of this feature of the data is confounding of potential sex-specific effects with other technical variables in the trials. However, it should not result in confounding with other biologically important variation, particularly in the validation of eQTLs. The literature is inconsistent with respect to terminology, but we prefer to make the distinction between the terms “validation” and “replication” (Igl et. al Hum Hered 2009; 67: 66-68).
Validation involves a confirmation attempt on a sample from a different population, where systematic variability (in this case, the gender difference that was noted by the reviewers) would be expected in addition to random variability. If a genetic association can be validated, it shows greater generalizability than a replicated association, in which the confirmation sample is drawn from the same population as the discovery sample. We agree that future studies should include both men and women as a way to identify potential sex-specific differences in transcriptional response to vaccine. Additional trials will also be required to assess the effects of age, other risk-increasing clinical variables, such as compromised immunity, and ethnic/racial origins.
2) We suspect that an interaction term in the models would actually be more efficient in detecting genetic effects modified by the time point after vaccination. Why was this not used (instead a simple comparison among time points was made)?
We agree with the reviewers and did attempt to investigate this. The model we employed for the analysis is as follows:
equation m6
This model includes interaction terms (βigj) for the change in the coefficient for the additive effect of genotype by day. In the loci showing significant transcriptional response to the vaccine, we found that the main effect of genotype was far stronger than these interaction effects. The sign of the additive effect does not change across the days, but the magnitude tends to increase after vaccination. We therefore speculate that our sample size does not provide sufficient power to exploit the interaction terms. Future larger studies may allow more precise estimates of the time × gene interaction and allow identification of additional vaccine response loci beyond those detected here.
3) The fact that there is enrichment in the QQ plots on Figure 6A does not provide evidence about causality of SNPs to vaccine response but simply of correlation among the two effects. More complex models need to be implemented to show causality.
We agree with the reviewers that the enrichment shown in the QQ plots does not provide evidence of causality. As noted, the results do support the possibility that the selected markers are enriched for true associations to the humoral immune response to influenza vaccine. We have revised the manuscript to clarify this point and have separated Figure 6A from Figure 6B to help reinforce the difference in interpretation.
A potential causal relationship between the genes identified and the antibody response is supported by the conditional independence analysis shown in Figure 7 (formerly Figure 6B). For this analysis, we first computed the R2 between SNP (genotype at SNPs that were validated eQTL and where gene expression was associated with the magnitude of the antibody response) and antibody response using an additive model. We then performed a regression analysis predicting antibody response from gene expression and we used the residuals from that analysis to perform another association analysis. We observed that the association between SNP and antibody, as measured by non-zero R2 values, disappeared after accounting for gene expression. This reduction in R2 was not observed in randomly-sampled sets of SNPs. The preceding data showed that the SNPs used are eQTL and that gene expression changes (e.g., day 1 post-vaccine) precede the change in antibody titers (day 28). Taken together, these observations support the interpretation that the predictive information in the genotype for the change in antibody titer is mediated by gene expression. This analysis supports a causal role for these genes in the change in antibody titer.
4) The study design used two different array versions: Illumina HT-12v3 and HT-12v4. The authors should describe how they combined data from these arrays, particularly for transcripts where the probe sequence differs between versions.
The sub-heading “Expression data processing and quality control” in the Materials and methods section has been re-written to more clearly explain how the expression microarray data were processed. In summary, for the initial part of the study (the eQTL component), the two cohorts were analyzed independently. Because of the difference in array versions, unique reporter identifiers (Illumina ProbeIDs and sequence-based nuIDs) were used to ensure that the same set of reporters was being analyzed in each of the cohorts. The two resulting expression data sets, which included data from the same set of 9,809 reporters, were then combined for the integrative genomic analysis component of the study.
5) The expression data processing should be described in more detail to make it clear how normalization, and so on, was performed given the different time points, array versions, and discovery/replication cohorts.
The sub-heading “Expression data processing and quality control” in the Materials and methods section has been re-written to more clearly explain how the expression microarray data were processed. The data for all individuals and all time points for a particular cohort were processed together. In the eQTL components of the analysis, the two cohorts were analyzed independently. Background adjustment, variance stabilization transformation, and normalization (VST-RSN) were carried out using identical parameters for each cohort. When the two data sets were combined for the integrative genomic analysis component, an additional quantile normalization was performed to adjust for batch effects between cohorts.
6) The authors mention that the person effect within their model would account for differences in cell populations among individuals, given that these are stable within individuals during the sampling timeframe. No data or references are provided to support this claim. Immune cell relative proportions can change rapidly.
We agree with the reviewers. Our model has a random-effect term to account for the fact that individuals were serially sampled. This person effect could account for systematic variation that was stable within an individual over the duration of the study (days 0–14), but would not account for changes in cell composition following vaccine administration. The corresponding paragraph in the manuscript has been corrected to reflect this.
7) The authors use the delta r2 metric to identify eQTLs that differ between time points or in response to the vaccination. It seems that the authors could have performed the analysis using the delta of transcription(time point x)−transcription(baseline) as a trait. This normalizes to each individual’s starting point. Furthermore, a difference of r2 between time points can have multiple causes (as the authors point out). It seems to us that what they are most concerned with is cases where the effect size (beta) of a given eQTL differs between time points, or where the actual response (transcription time point x−transcription baseline) differs between individuals. The authors reasoning behind these choices should be made clearer.
We attempted analyses with an approach similar to that proposed by the reviewers in the course of our work, but found that the approach that was ultimately chosen to explore the day differences was the most powerful. Specifically, utilizing a pairwise comparison (difference) between time points as the substrate for the eQTL analysis would lead to an increase in the technical variance of the phenotype, as the sum of two independent (technical) errors has twice the variance of an individual measurement. We agree with the reviewer that our interest is in investigating the source of changes in R2 explained on each day, and we analyzed this change in R2 to show that (a) in principle, R2 for the explanatory power of the genotype increases in the days after vaccination, and (b) this change is due to an increase in the magnitude of the additive effect, and not in a change in the dispersion of the errors across days.
8) For the trans-eQTL analysis, the authors consider as genome-wide significant any SNP-transcript pair with p<5×10-8. Indeed for a GWAS of common variation, this may be appropriate for a single trait, but the authors should account for testing nearly 10K transcripts.
In this study we performed only cis-eQTL analysis within an interval of 1 Mb centered on the reporter sequence. Positive associations in this analysis are most likely to result from direct allele-specific effects on transcription rate or transcript stability. A p-value threshold of 5 x 10-8 would, therefore, be conservative based on either number of SNPs or LD structure around the transcript. Future studies with much larger sample sizes might allow a full trans-eQTL analysis in which one would seek to identify genetic variants at other loci that indirectly affect transcript abundance measures. We have removed the term “genome-wide significance”, in reference to the p-value, from the relevant section of the Results.
[After resubmission, the editors also asked for the following comments to be addressed prior to acceptance.]
A) The revision is not satisfactory for point 7, above, on the delta r2. The case is still not convincing because only 146 of the 541 cases were validated. It is not clear why the difference in r2 is necessarily interpreted as a significant change when only 30% of the cases are validated.
We chose to highlight the observed phenomenon of a change in the magnitude of the genetic effect after the experimental stimulus (delta-Rg2) because it is evident from our data (as exemplified by the NECAB2 locus in Figure 2A) and it is a general property of the observed eQTLs (Figure 2B). Importantly, this observation could not be made by cross-sectional study designs such as genome-wide association studies or single-time-point eQTL analyses.
The Reviewing editor has raised concern about the validation of the delta-Rg2 property in the two cohorts. We agree that the metric is not well defined in terms of the type I error, and that we should avoid the term ‘validation’ in this setting. The manuscript has been modified to reflect this change. However, our data shows that the delta-Rg2 is a stable property of SNP-transcript pairs in both cohorts. To more objectively illustrate this, we evaluated the Spearman's rank correlation of the delta-Rg2 between the two cohorts using all SNP-transcript pairs (Cor = 0.25, p < 2x10-16). This information has been added to the manuscript. To select a conservative set of candidate loci based on this phenomenon, we chose to impose a stringent intersection rule, requiring a threshold of the top 1% of the delta-Rg2 distribution in both cohorts, in addition to evidence of a strong cis- acting eQTL (genotype effect p < 5 × 10-8) in both cohorts. The manuscript was modified to clarify that the 146 SNP-transcript pairs that were selected for further analysis comprise the set that passed this strict intersection rule.
B) The modelling as previously discussed in point 3 is somewhat weak and causality is not as well demonstrated, as the authors want to claim. The fact that the association with the residuals is lower than with the original values is not direct but indirect evidence that there are casual relationships. Also, how were the “random” SNPs chosen? Establishing a relevant null is important here. Also, direct causal models would help in establishing the relationships rather than the indirect approach the authors use.
We agree with the Reviewing editor and, as stated above, spent considerable time reconsidering our approach to the evaluation of causal versus alternative association model comparisons.
Our work considered three distinct types of associations in the vaccine response data: genotype to expression (eQTL), expression to antibody response and genotype to antibody response. An important advantage of experimental designs such as ours is the opportunity to consider relationships between these distinct types of association. To respond to the reviewing editor's questions and to more formally evaluate alternative models, we have now adapted the approach recently published by Eric Schadt and colleagues (Millstein et al., 2009), which provides a statistical testing framework for so-called causal versus reactive model analyses. This work is related to ideas and methods developed by Judea Pearl and others, but more specifically applied to the context of eQTL data. Under this approach, a formal statistical hypothesis test is developed for evaluating three possible scenarios, referred to as independent, causal and reactive models.
This concept is now depicted in Figure 7A and described in the text of our paper. This formulation incorporates the three pairs of associations considered in our study: locus-to-expression, expression-to-trait, and locus-to-trait, together with a more explicit analysis of conditional independence between two of these variables given the third, such as the association (and change in association) between locus and trait conditional on gene expression, which is the focus of the causal model. Two separate p-values are generated, considering both the causal and reactive scenarios against the null model of independence. Details of our approach are described in the Methods section of the paper and the results of the analyses are presented in Figure 7B and described in the text. In keeping with the Reviewing editor's comments, and in an attempt to more globally evaluate the trend observed in our data, we performed a large-scale simulation under a null model that used permutation to sever the connection between eQTL and expression-to-trait associations. In each permutation, we selected a set of candidate SNP-transcript pairs that passed eQTL and expression-to-permuted-phenotype association. We examined the causal versus reactive p-value distributions for these simulation-based results. In no instance did we observe the trend towards the causal versus the reactive model that was observed in our data.
Finally, we performed a power analysis to examine the sample size that would be required to establish statistical significance under this methodology when effect sizes and errors are consistent with our data. These analyses suggest that a sample size of 200 – which corresponds to the combined sample size of our two cohorts – has 50–60% power to detect deviations from the independence model if they are present (Figure 7C). This result suggests that our sample size, while giving us some opportunity to discriminate these alternatives, is still too small, and that a sample size at least twice as large would be needed for more clarity.
Articles from eLife are provided here courtesy of
eLife Sciences Publications, Ltd