Using gene expression signatures to identify novel treatment strategies in gulf war illness
GWI cohort
As part of a larger ongoing study a subset of GWI male subjects (n?=?17) and healthy but sedentary Gulf War era veterans (n?=?22) were recruited from the Miami Veterans Administration Medical Centers, clinics
and the local veteran community between April 2006 and May 2008. All subjects were
comparable in age, body mass index (BMI), and ethnicity. Subjects were male and ranged
in age between 30 and 55. Inclusion criteria was derived from Fukuda et al.10], and consisted in identifying veterans deployed to the theater of operations between
August 8, 1990 and July 31, 1991, with one or more symptoms present after 6Â months
from at least 2 of the following: fatigue; mood and cognitive complaints; and musculoskeletal
complaints. Subjects were in good health prior to 1990, and had no current exclusionary
diagnoses 11]. Medications that could have impacted immune function were excluded. Use of the Fukuda
definition in GWI is supported by Collins et al.12]. Additional details may be found in Broderick et al.13].
All subjects signed an informed consent approved by the Institutional Review Board
of the University of Miami. Ethics review and approval for data analysis was also
obtained by the IRB of the University of Alberta.
Gene expression
Blood was drawn at rest at comparable times of day from each subject during the April
2006 to May 2008 period. At each blood draw three 8-mL tubes of blood were collected
in CPT vacutainers (B-D- Biosciences, San Jose, CA). The peripheral blood mononuclear
cells (PBMC) were isolated and stored in liquid nitrogen under conditions designed
to maintain viability. Specifically, whole blood was added to Ficoll-Paque, centrifuged
at 1000 g for 25 min. PBMC’s were isolated from the PBMC ring atop the Ficoll layer
into a separate tube, centrifuged at 300Â g for 10Â min, then re-suspended in PBS. Cells
were then counted using a Beckman Coulter viCell, and cryopreserved in freezing media
(temperature lowered 1
o
C per minute until -80
o
C).
Total RNA was extracted using TRI Reagent (Molecular Research Center, Cincinnati,
OH) following the manufacturer’s protocol. The quality and quantity of RNA was assessed
using the Agilent Bioanalyzer 2100 RNA 6000 Nano Kit (Agilent Technologies, CA). From
each sample, 300Â ng of total RNA was converted into cDNA by reverse transcription
using a T7-oligo(dT) primer and the Affymetrix 3? IVT Express Kit (Affymetrix, Santa
Clara, CA) according to standard manufacturer protocol. The generated cDNA was purified
using the GeneChip Sample Cleaning Module (Affymetrix) and labeled cRNA was generated
by in vitro transcription using the biotinylated nucleotide mix. This was then purified with
the Cleaning Module and quantified using the Nanodrop ND-1000 spectrophotometer (NanoDrop
Technologies, Inc., Wilmington, DE USA). In each preparation 11Â ?g cRNA was fragmented
in Fragmentation Buffer (Affymetrix) in a final reaction volume of 25Â ?l.
Hybridization, washing, staining and scanning were done using Affymetrix GeneChip
instruments (Hybridization Oven 640, Fluidics Station 450Dx, Scanner GCS3000Dx) and
Affymetrix Human U133 2.0 arrays (Affymetrix) as per manufacturer’s standards. Microarray
image files (.cel data) were generated using the Affymetrix GCOS software tool with
default microarray analysis parameters to provide overall within chip normalization
of the image intensity distribution. The quality parameters that were monitored besides
cRNA total yield and cRNA A260/A280 ratio included: (i) background noise (Q value),
(ii) percentage of present called probe sets, (iii) scaling factor, (iv) information
about exogenous Bacillus subtilis control transcripts from the Affymetrix Poly-A control
kit (lys, phe, thr, and dap), and (v) the ratio of intensities of 3? probes to 5?
probes for a housekeeping gene (GAPDH).
To generate a broad comparison group the Gene Expression Omnibus (GEO) DataSets 14], 15] and data from Suthram et al.16] were used to obtain a set of gene expression profiles describing a number of human
disease conditions. We restricted our selection to include only those sets in which
both disease and a corresponding healthy control group were measured in the same cell
type or tissue in the same experimental conditions. Sets that included different exposure
times, exposure concentrations or multiple cell/tissue types were each treated as
a separate disease condition. For consistency, and to avoid complications arising
due to cross platform comparisons, datasets were restricted to the Affymetrix Gene
Chip Human Genome U133A, U133 Plus 2.0, U133A 2.0 and U95 Version 2 arrays, to align
with our GWI gene expression data. All diseases affecting male subjects in the GEO
database meeting these criteria were included in this study. Overall, this resulted
in 101 human disease gene expression datasets (Additional file 1: Table S1). The disease profiles selected provide a comprehensive set of diseases
and include various cancers, neurodegenerative diseases, autoimmune illnesses, chemical
exposures, viral infections, neurological and neuromuscular disorders.
Data transformation and normalization
Gene expression data was Log2 transformed then normalized using a Z-score transformation
for each microarray sample to allow for the direct comparison of values across various
samples and diseases.
Gene functional modules
Here we use 4,620 functional modules defined by Suthram et al.16] from the human protein-protein interaction network. Individual genes i in module
j were compared between GWI subjects and healthy controls using an unpaired t-test to generate a T-Score, G
ij
. Differential module activity (MA
j
) was determined using the average of the absolute T-Score values for all N
j
genes in a given module, j, such that:
MAj=?i=1NjGiNj(1)
This is a modification of the module activity described in 16] were the signed T-score was used. When averaging over signed T-scores equal but opposite
scores will tend to cancel, resulting in no module activity, when in fact there are
genes significantly over or under expressed. Here the absolute value was used in the
differential module activity to capture both up and down-expressed genes. A threshold
of MA?=?1.5, corresponding to a p-value of ~0.1, was chosen as a liberal cutoff value
to identify differentially expressed gene modules. Only modules with MA values above
threshold were taken to be differentially expressed in GWI and considered for further
analysis, all others were discarded.
Individual gene expression analysis
The 202 genes in the 19 modules with MA values above threshold were individually compared
between GWI and controls. Individual genes i were compared between GWI subjects and
healthy controls using an unpaired t-test. To account for multiple comparisons false discovery rates (FDR) were then calculated
for each comparison from these resulting p-values using the procedure introduced by
Storey 17]. Genes with FDR of less than or equal to 0.05 were taken to be significantly different
in GWI compared to control. Fold change was calculated by calculating the ratio of
average gene expression between GWI and health controls for non-Log2 transformed data.
Positive values were taken to indicate the fold-increase while ratios that were less
than 1 were inverted and given a negative value to denote the fold-decrease. An absolute
fold change of 1.5 was taken as a cutoff 18].
Mining the pharmacogenomic knowledge base (PharmGKB) database
Genes from modules affected in GWI were screened against the PharmGKB database (8.1.2015)
19] to find gene-drug and gene-disease relationships supported with pharmacogenomics
research reported in the literature.
Pathway-based functional annotation of gene modules
Functional annotation of gene modules with gene-drug relations was performed using
the ConsensusPathDB 20]–22] to provide biological pathway information for each gene set. Over-representation
analysis 20] incorporating the Kyoto Encyclopedia of Genes and Genomes (KEGG) (73.0) 23], Netpath (1.1.2015) 24], the Small Molecule Pathway Database (SMPDB) (8.1.2015) 25], the Integrating Network Objects with Hierarchies (INOH) (1.1.2015) 26], Biocarta (2009_05_12) 27], Humancyc (18.5) 28], Signalink (8.1.2015) 29], Edinburgh human metabolic network (Ehmn) (1.1.2015) 30], Reactome (51) 31], PharmGKB (8.1.2015) 19], Wikipathways (9.1.2015) 32] and the Pathway Interaction Database (PID) (2014_02_14) 33] pathway sets was used to interpret the function of druggable gene modules. Here the
significance of the observed overlap between the gene module and the members of known
pathways, compared to random expectations, was calculated based on the hypergeometric
distribution. A minimum overlap of 2 genes between the gene module and the pathway
set at a p-value cutoff 0.01 was required. Specifically, the p-value is calculated
as the probability of randomly finding k or more successes from the population in
N total draws. Thus, small p-values indicate a greater overlap than expected by chance.
Pathway sets containing the majority of the druggable genes, the highest number of
module genes overall, and the lowest p-value were taken as the functional annotation
of the module. Pathway annotation was performed only to provide biological pathway
information for each gene module set. All subsequent analysis was performed on these
gene module sets constructed on the bases of human protein-protein interactions 16] and not on the bases of known pathway membership.
Gene module alignment across illnesses
To provide a quantitative measure of similarity in module expression between GWI and
another disease a partial Spearman correlation was calculated between the T-Score
values of genes compared between select diseases (Additional file 1: Table S1) and their controls within each of the modules preferentially expressed
in GWI 16]. Partial correlation measures the degree of association between two random variables,
while removing the effect of a set of controlling random variables. Partial Spearman
correlation was used to assess disease similarity while explicitly factoring out dependencies
between different gene expression experiments 16]. The Partial Spearman correlation coefficient between GWI (G) and a disease condition (D), r GD,C
, was conditioned on the T-Score values between the control samples (C), corresponding to each comparator disease, such that:
rGD,C=rGD?rGCrDC1?rGC21?rDC2(2)
with the generic Spearman correlation coefficient, r GD
, given as:
rGD=?i=1NjGi?G¯Di?D¯?i=1NjGi?G¯2?i=1NjDi?D¯2(3)
where G i
the ranked TScores between GWI and controls, D i
the ranked TScores between a given human disease condition and controls, and C i
the ranked TScores between between GWI and human disease condition control groups
for the i th
gene in module j containing N j
genes. Note that the terms in the numerator of Equation (2) are given by Equation (3) with substitutions for G, D, and C appropriately as labeled. Here in Equation (3) the barred terms (
G¯,D¯,C¯
) denote the averaging of all genes in a given module. Correlations with p-values???0.01 were considered statistically significant. Additionally, as we were
looking for similarity between modules, inverse correlations were not considered.
Partial correlations were calculated between GWI and the 101 disease conditions (Additional
file 1: Table S1) in the GWI affected modules.
Significance of correlation
A random background distribution of disease correlations across each module was created
to determine if module correlations between GWI and other diseases were significantly
different from random chance. Gene to module assignments were first randomized, while
preserving the total number of modules, the number of genes per module, and the number
of modules to which each gene belongs as in Suthram et al.16]. Then disease and control datasets for the disease conditions were randomly shuffled
before calculating partial correlations with intact GWI datasets. This process was
performed 100 times for each module to generate the random background distributions.
The cumulative distribution of the random background correlations per module was used
to calculate the p-value of the true correlations in each module. FDR was calculated
for each correlation from these resulting p-values using the procedure introduced
by Storey 17]. Correlation values with FDR of less than or equal to 0.05 were taken to be significantly
different from random background. All others were removed.