New insights about host response to smallpox using microarray data

Active Modules

The goal of our analysis was to identify differentially expressed modules in monkeys with early death, i.e., died in the first three days, and monkeys with a later death (ED and LD, respectively). Our hypothesis is that these modules play crucial roles in the host response to smallpox infection. The gene groups were defined according to the KEGG database (updated on 04/17/2006) [8], as described in the Methods section.

Initially, this analysis was based on the ED and LD monkey groups using all KEGG Orthology (KO) categories (at level 4), indicating three highly active modules in LD monkeys: Cytokines and Cell Adhesion Molecules (CAMs), which are consistent with the classical model of immune and inflammatory responses, where an activation and proliferation of leukocytes take place. The Alzheimer's Disease module also appeared as differentially activated, most likely due to its hybrid composition (20/22 genes that compose the AD module are present in the microarray platform used by Rubins et al [1], and 14/20 are also members of other modules).

In addition to the ED and LD analysis, modules were also analyzed during the disease course. We have found 27 modules presenting differential activation in at least one point (Figure 1). Several of those modules are related to immune response, like CD molecules, CAM ligands, antigen processing and presentation and the previously considered activated in a late death Cytokines and Cell adhesion molecules.

A novel finding highlighted by our systems biology approach for all macaques (Figure 1) was an altered expression in the Natural Killer (NK) mediated cytotoxicity module. NK cells play a major role in the host-rejection of virally infected cells. Due to their strong cytolytic activity and the potential for auto-reactivity, they are tightly regulated. The activation of NK cells requires the action of proinflammatory cytokines in combination with cell surface receptors [9]. Under normal conditions, inhibitory NK cell receptors engage MHC class I (MHC-I) molecules. Such inhibitory activity normally predominates over that provided by the activating receptors, which bind to ligands expressed on virus-infected and other stressed cells [10]. Cells with normal levels of MHC-I are generally protected from NK cell-mediated lysis. Virally infected cells often express reduced levels of MHC-I and this inability to engage the inhibitory receptors makes the cells susceptible to NK cell attack [10]. The down regulation of Natural killer cell mediated cytotoxicity module reported here could indicate a susceptibility to a NK cell mediated lysis (days 1–3, Figure 1).

The contributions of individual molecules for each significantly altered module are available in the additional material [see Additional files 1 and 2], represented by scores and the respective p-values.

Relevance Networks (RNs)

To understand the profiles of interaction between the genes of the groups studied, we used a procedure known as relevance networks (RNs), initially proposed by Butte et al [11]. Interaction here represents significant correlation between two variables (i.e. molecules), which may (but not necessarily will) physically interact with each other.

Cytokines and cell adhesion molecules (CAMs) modules were highly expressed in LD monkeys in the early phase of the infection (Figure 1) and were further explored in more detail by using the RN approach. A metamodule (CAMs and ECM receptors combined) was created to gain a better understanding of leukocyte migration to the infection locus. The connections between ICAM2 and ITGAV and ITGA6 shows a significant change (p value < 10^-3) when comparing day 0 versus the earliest infection groups (days 1–3), see Figure 2A.

To establish a successful infection, pox viruses affect the release of immune modulators by the infected cells, so it is crucial to understand the dynamics of the host immune system [6]. The population of the immune cells varies during the infection, as do their exposed membrane proteins. Therefore, the Relevance Networks approach can reveal a more detailed perspective on the dynamics of the molecules involved in the disease. A better assessment of the variation of those molecules during the onset of the infection could give a better understanding of the host response to smallpox, and ultimately of what may be the key factors that set the 30% mortality rate.

To reach the infection focus and kill the pathogen, the immune cells have to adhere and overpass the endothelium. These events are possible due to the interaction between the integrins and the extra cellular matrix receptors and are modulated by cytokines. More specifically, ESAM is an endothelial cell adhesion molecule and ITGAV, ITGB1 and ITGB5 are integrins, which are heterodimeric integral membrane proteins composed of an alpha chain and a beta chain. ITGAV encodes an alpha chain; ITGB1 and ITGB5 encode beta chain integrins, known to participate in the cell-surface mediated signaling. ICAM2 is an intercellular adhesion molecule that promotes an inflammatory response and the activation of AKT (a kinase). AKT phosphorylates components that block the apoptosis pathway [12, 13]. To explore further the graphs produced by RNs for the ECM receptors + CAMs the relations concerning ICAM2 and ESAM are detailed. The complete network is available in the additional material [see Additional files 3 and 4]. RNs of day 0 and days 1–3 show statistically significant alterations in the interactions between ICAM2 and the integrins ITGAV and ITGA6 (Figure 2A), which are important extracellular matrix receptors, involved in the adhesion, migration, differentiation, survival and proliferation of cells and cytoskeleton organization. Day 0 versus days 4–6 RNs revealed a significant relation between ICAM2 and other molecules that play crucial role in immune response against pathogens, PLXNC1 [14] and CDH5 [15] (Figure 2B). The relevance networks produced for cytokine group is available in the additional material [see Additional file 5].

ESAM shows a strong relation with ITGAV, ITGB1 and ITGB5 with variable patterns along the time course of infection. The quartet of ESAM, ITGAV, ITGB1 and ITGB5 reveals the relationship between alpha and beta chains in the composition of integrins and the relation of such heterodimeric proteins with ESAM. It is worth mentioning that ITGB3, ITGB6 and ITGB8, elements that interact with ITGAV, are present in the microarray platform design, but were excluded from analysis in the filtering steps (see Methods section).

Active modules and the RNs results pointed out drastic changes in the expression levels and interactions of important genes involved in essential cellular processes of the host response, such as cellular adhesion, apoptosis and immune response. These results constitute interesting new biological hypotheses about smallpox virus infection that could be further explored by expert virology groups.

Overlapping with Rubins et al and Jahrling et al

The data analyzed in our manuscript was previously analyzed by Rubins et al [1] and Jahrling et al [16]. The former focused on the microarray data and the later on the clinical findings. Rubins et al [1] reported the expression profile in IFN response, cell cycle/proliferation response, TNF-α/NF-κ B regulated genes, dose response and lymphocyte and Ig responses. Meanwhile Jahrling et al [16] discuss the hemorrhagic diathesis, cytokines activation (IL6, IL8, INFG, CCL2 and CCL4) and lymphocyte depletion.

Rubins et al [1] performed a detailed analysis of IFN response genes. The KEGG module of cytokines includes IL6, IL8, IFN -α, -β and -γ, CCL2 (MCP-1), CCL4 (MIP1β), TNF-α and NF – κ B. The influence of each molecule in the modules activation was also computed [see Additional files 1 and 2]. Once KEGG separates the molecules into functional classes, it is not possible to define cascading effects, like Rubins et al did [1]. Despite this fact, we manually compiled a module with the IFN-responsive genes detailed by Rubins et al [1] (in Figure 2A). For the ED and LD analysis, this module is significantly down regulated in monkeys that died early and up regulated in those which died later (see Methods for details). In addition, for the time course analysis, the IFN-responsive module is down regulated for day 0 and up regulated for days 1–3 and inactive for the other days (data not shown). Therefore, our results suggest that activation of IFN-responsive genes may play an important role in a successful immune response upon infection, as previously reported by Rubins et al [1].

The decrease in the relative abundance of B cells and MHC II reported by Rubins et al [1] are corroborated by the profile of module Antigen processing and presentation (Figure 1). The Cell cycle, Complement and coagulation cascades, Cytokines and Hematopoietic cell lineage modules identified here were also reported by Rubins et al [1] as differentially expressed using gene by gene approaches.

Altered D-dimer levels were interestingly investigated by Jahrling et al [16], which reflect the serious influence of smallpox virus inhibitors in the complement and coagulation cascades. Our active modules are also in agreement with such findings. Our results suggest an unbalanced cytokine levels, which is strongly associated with sepsis and organ failure death in infected animals.

In addition, we have also observed alterations in the hematopoietic cell lineage module which could be a consequence of the influx of immature leukocytes into the bloodstream, particularly neutrophils (aka "left shift", using medical jargon).

Rubins et al [1] observed a considerable depletion of T cells, and a large number of apoptotic cells in the splenic periarterial lymphatic sheath, possibly due to the viral replication in antigen presenting cells in the lymphatic tissues, which could lead to apoptosis in the lymphocytes. The lysis of vulnerable cells induced by NK cells (described earlier) and the down regulation of the NK cell mediated cytotoxicity module supports such hypothesis.

Data description

In this work we used cDNA microarray data containing 37,632 elements that represents approximately 11,000 unique genes and approximately 12,000 genes without identification [1, 17]. The slides were produced with human cDNA clones and the test samples were extracted from primate models (cynomolgus macaques). The samples of interest were always labeled with Cy5 and the reference samples were labeled with Cy3. The original dataset was obtained from the SMD [18, 19] and is associated with the Experimenter and ExptDate keys "KRUBINS" and "2001-08-05", respectively. In this work, only data associated to studies 1 (Harper or India 7124 variola strains through a combination of aerosol exposure [5 × 108 plaque-forming units (pfu)] and high-titer i.v. inoculation [109 pfu]) and 2 (India 7124 strain via i.v. route alone [109 pfu]) [1] were used. Our goal is to analyze the host response to a successful infection; therefore the third study (varied viral dose) data was not used here.

The raw data was submitted to exploratory analyses. The main problem found was a bias of log ratio values from left to right in the majority of the chips, see Figure 3A. To address this problem, we used a normalization procedure proposed by Futschik and Crompton [20] known as OLIN. The OLIN method has some clear advantages: (i) it relies on strong statistical formalism; (ii) has the capacity to remove spatial and intensity (systematic) errors in microarray data, which is the main problem we found in the original dataset; (iii) Futschick and Crompton [20] precisely addressed the importance of optimization of model parameters on microarray data normalization (instead of simply accepting the method default values); (iv) the method performance was conveniently compared to other proposed methods [20]. Basically, the effect of log ratio values in relation to mean expression values (M and A, respectively) was is estimated using LOWESS regression [21]with an optimization step for the selection of the parameter α _A that control the smoothing of the regression. Subsequently, a new LOWESS regression was performed with log ratio values based on the x and y directions on the slide, and again, the optimization step of the parameters α _x and α _y was performed. This model successfully corrected the spatial-dependent artifact observed in the dataset (Figure 3B). Also, in a filtering step, we used only spots with signal intensity greater than background intensity in any sample and all spots with signal intensity greater than 2.5 times background intensity in 80% of the samples.

Active modules

In order to find significantly altered gene groups we used the "active modules" approach, originally proposed by Segal et al [22]. The rationale of this method is to evaluate differential activation or repression of gene modules, which represent the biological processes.

The first step was to organize the data in a matrix E _ij , i = 1,...,n _s , j = 1,...,n _l , where e _ij represents the gene i in the sample j. In the present study, 199 groups were selected according to the level 4 of KEGG [8] (updated on 04/17/2006) [see Additional file 6], so each group has n _g genes, g = 1,...,199.

After the data normalization described earlier, expression values were standardized by subtracting the average expression calculated over all the arrays. This process adjusted the mean values of each gene to zero, which were then used to build a discrete version of E _ij , named E _ij *, by applying the indicator function e _ij * = I(e _ij ≥ 1) - I(e _ij ≤-1) for all i and j (an indicator function assumes value 0 or 1 for false or true, respectively). Thus, E _ij * is a matrix where each element can assume the discrete values 1 (induction), 0 (equal expression) or -1(repression).

For each of the g groups of interest, each array j was then splitted in two disjoint sets of genes: one with exactly the n _g elements of the group, and other constituted by the remaining n _s -n _g genes. One should also consider two random variables Y_1j^(g)and Y_2j^(g)representing the number of induced genes in the first and second group, respectively.

Under the null hypothesis of no biological activity for g in sample j, Segal et al [22] propose that Y_1j^(g)~hypergeometric(n _g ; n _s -n _g ; q_1j) and Y_2j^(g)~hypergeometric(n _g ; n _s -n _g ; q_2j). Intuitively, the absence of biological activity makes the probability of an induced (or repressed) gene to belong to the first group equal to 1/2 and p-values associated with the observed values y_1j^(g)and y_2j^(g)can be easily calculated. This procedure defined a new matrix P _gj where the element p _gj represents the fraction of induced genes of group g in the array j, or zero for the non-significant modules (we assumed a significance level α = 0.05). In the cases where the significance is for the repressed genes, we just made the fraction p _gj negative. The arrays were then categorized in biological types t = 1,, t and for each pair of gene set g and biological type t, we count the number of arrays of type t that showed gene set g induced (repressed). We then used the hypergeometric distribution to classify each gene group as significantly induced or repressed across the biological types (p value < 0.05).

The biological types are defined as following: n_t = 1 means that the macaque was not infected, n_t = 2 means death 1 to 3 days after infection, n_t = 3 means death 4 to 6 days after infection and n_t = 4 means death after the 6^th day. n_t = 2 was also labeled Early Death (ED) and n_t = 3 and n_t = 4 were labeled Late Death (LD). The contribution of molecules for activation of a significantly altered module was computed for time course and "ED and LD" analyses [see Additional files 1 and 2]. Using these values, it is possible to evaluate the influence of each gene in the module activation or repression (p-value).

Relevance Networks (RNs)

For a better understanding of the Relevance Networks (RN) method, a definition of the term interaction is required. Here we are interested in a more general definition of protein interaction, which goes beyond physical interaction and includes also functional or phylogenetic relationships. Interaction here is inferred by significant correlation between two variables (i.e. molecules), which may (but not necessarily do) physically interact with each other. Aiming to identify such correlations, the Relevance Networks (RN) method, initially proposed by Butte et al [11], was employed. The method computes the squared linear Pearson correlation coefficient, $r_{ij}^{\text{2}}$, between all gene pairs for each biological type (as described in the active modules section of methods), defining a fully connected graph. Using a re-sampling method, the algorithm estimates a cut-off value, c, to split the graph into small sub-graphs where $r_{ij}^{\text{2}}$ > c. These sub-graphs are the relevance networks (RNs).

In order to identify significant differences between two RNs, we adapted this method to contrast two distinct conditions. Therefore, the linear Pearson correlation coefficients between every pair of genes in each condition, $r_{ij}^{\text{1}}$ and $r_{ij}^{\text{2}}$ are computed and then a Fisher's Z-transformation [23] is used to translate the probability distribution of the random variable $r_{ij}^{\text{1}}$ - $r_{ij}^{\text{2}}$ into an approximated standard normally distributed random variable, which permits the identification of pairs of genes that displayed significant (positive or negative) differences between $r_{ij}^{\text{1}}$ and $r_{ij}^{\text{2}}$ (p < 0.001). This type of analysis was accomplished for the 199 gene groups used in the active modules procedure, linking the findings with biological knowledge.