HLC and tissue collection.
The HLC was assembled from a total of 780 liver samples (1–2 g) that were acquired from Caucasian individuals from three independent liver collections at tissue resource centers at Vanderbilt University, the University of Pittsburgh, and Merck Research Laboratories (Table S1
). The Vanderbilt samples (n
= 504) included both postmortem tissue and surgical resections from organ donors and were obtained from the Nashville Regional Organ Procurement Agency (Nashville, Tennessee), the National Disease Research Interchange (Philadelphia, Pennsylvania), and the Cooperative Human Tissue Network (University of Pennsylvania, Ohio State University, and University of Alabama at Birmingham). The Pittsburgh samples were normal postmortem human liver and were obtained through the Liver Tissue Procurement and Distribution System (Dr. Stephen Strom, University of Pittsburgh, Pittsburgh, Pennsylvania). The University of Pittsburgh samples (n
= 211) were all postmortem, as were the Merck samples (n
= 65), which collected by the Drug Metabolism Department and reported previously [40
All samples were stored frozen at −80 °C from collection until processing for RNA and DNA; some samples had been stored for over a decade before being processed for this study. Demographic data varied across centers for these samples and were missing in many cases. In cases where age, sex, or ethnicity data were not available in the patient records, we imputed it from the gene expression and/or genotype data (see below). Of the 780 samples collected, high-quality DNA was isolated on 548 samples, and 517 of these were successfully genotyped on the Affymetrix genotyping platform (see Methods below). Of the 517 successfully genotyped samples, high-quality RNA was isolated and successfully profiled on 427 samples. This set of 427 genotyped and expression profiled samples comprised the HLC. Table S1
gives a summary of the demographics and other annotations on the 427 individuals that were successfully genotyped and expression profiled. All counts and descriptive statistics include the imputed data. All samples and patient data were handled in accordance with the policies and procedures of the participating organizations.
Mouse crosses and tissue collection.
C57BL/6J (B6) mice were intercrossed with C3H/HeJ (C3H) mice to generate 321 F2 progeny (161 females, 160 males) for the BXH wild type (BXH/wt). C57BL/6J (B6) mice were intercrossed with Castaneus (CAST) mice to generate 442 F2 progeny (276 females, 166 males) for the BXC cross. All mice were maintained on a 12 h light–12 h dark cycle and fed ad libitum. BXH mice were fed Purina Chow (Ralston-Purina) containing 4% fat until 8 wk of age. From that time until the mice were killed at 20 wk, mice were fed a western diet (Teklad 88137, Harlan Teklad) containing 42% fat and 0.15% cholesterol. BXC mice were fed Purina Chow until 10 wk of age, and then fed western diet (Teklad 88137, Harlan Teklad) for the subsequent 8 wk. Mice were fasted overnight before they were killed. Their livers were collected, flash frozen in liquid nitrogen, and stored in −80 °C prior to RNA isolation.
The BXH cross on an ApoE null background (BXH/apoE) was previously described [41
]. Briefly, C57BL/6J ApoE null (B6.ApoE–/–
) were purchased from Jackson Laboratory. C3H/HeJ ApoE null (C3H.Apo E–/–
) were generated by backcrossing B6.ApoE–/–
to C3H for ten generations. F1 mice were generated from reciprocal intercrossing between B6.ApoE–/–
, and F2 mice were subsequently bred by intercrossing F1 mice. A total of 334 (169 female, 165 male) were bred, and all were fed Purina Chow containing 4% fat until 8 wk of age, and then transferred to western diet containing 42% fat and 0.15% cholesterol for 16 wk. Mice were killed at 24 wk, and liver, white adipose tissue, and whole brains were immediately collected and flash-frozen in liquid nitrogen.
All procedures of housing and treatment of animals were performed in accordance with Institutional Animal Care and Use Committee regulations.
Microarray design, RNA sample preparation, hybridization, and expression analysis.
Array design and preparation of labeled cDNA and hybridizations to microarrays for the human liver cohort.
RNA preparation and array hybridizations were performed at Rosetta Inpharmatics. The custom ink-jet microarrays used in this study were manufactured by Agilent Technologies and consisted of 4,720 control probes and 39,280 noncontrol oligonucleotides extracted from mouse Unigene clusters and combined with RefSeq sequences and RIKEN full-length cDNA clones (Table S4
Liver samples extracted from the 427 Caucasian individuals were homogenized, and total RNA extracted using TRIzol reagent (Invitrogen) according to manufacturer's protocol. Three micrograms of total RNA was reverse transcribed and labeled with either Cy3 or Cy5 fluorochrome. Purified Cy3 or Cy5 complementary RNA was hybridized to at least two single microarrays with fluor reversal for 24 h in a hybridization chamber, washed, and scanned using a laser confocal scanner. Arrays were quantified on the basis of spot intensity relative to background, adjusted for experimental variation between arrays using average intensity over multiple channels, and fitted to an error model to determine significance (type I error), as previously described [42
]. Gene expression is reported as the mean-log ratio relative to the pool derived from 192 liver samples selected for sex balance from the Vanderbilt and Pittsburgh samples, because the RNA from the Merck samples had been amplified at an earlier date. The error model used to assess whether a given gene is significantly differentially expressed in a single sample relative to a pool composed of a randomly selected subset of samples has been extensively described and tested in a number of publications [42
The age, sex, race, center, alcohol use, drug use, and steatosis variables presented in Table S1
were tested for association to the gene expression traits. Only age, sex, race, and center were significantly associated with the expression traits beyond what would be expected by chance. As a result, all gene expression traits were adjusted for these covariates. The lack of association between the expression traits and alcohol use, drug use, and steatosis was somewhat surprising, but may be due to the sparseness of these data, resulting in a lack of power to detect significant associations.
Array design and preparation of labeled cDNA and hybridizations to microarrays for the mouse liver and adipose tissue samples.
RNA preparation and array hybridizations were again performed at Rosetta Inpharmatics. The custom ink-jet microarrays used in the BXH/wt, BXH/apoE, and BXC crosses were manufactured by Agilent Technologies. The array used for the BXH/apoE and BXH/wt samples consisted of 2,186 control probes and 23,574 noncontrol oligonucleotides extracted from mouse Unigene clusters and combined with RefSeq sequences and RIKEN full-length cDNA clones (Table S5
). The array used for the BXC cross consisted of 39,280 noncontrol oligonuceotides again extracted from the mouse Unigene clusters and combined with RefSeq sequences and RIKEN full-length cDNA clones (Table S6
Mouse adipose and liver tissues from all of the crosses were homogenized, and total RNA extracted using Trizol reagent (Invitrogen) according to manufacturer's protocol. Three micrograms of total RNA was reverse transcribed and labeled with either Cy3 or Cy5 fluorochrome. Labeled complementary RNA (cRNA) from each F2 animal was hybridized against a cross-specific pool of labeled cRNAs constructed from equal aliquots of RNA from 150 F2 animals and parental mouse strains for each of the three tissues for each cross. The hybridizations for the BXH/apoE cross were performed in fluor reversal for 24 h in a hybridization chamber, washed, and scanned using a confocal laser scanner. The hybridizations for the BXH/wt and BXC crosses were performed to single arrays (individuals F2 samples labeled with Cy5 and reference pools labeled with Cy3 fluorochromes) for 24 h in a hybridization chamber, washed, and again scanned using a confocal laser scanner. Arrays were quantified on the basis of spot intensity relative to background, adjusted for experimental variation between arrays using average intensity over multiple channels, and fitted to a previously described error model to determine significance (type I error) [42
]. Gene expression measures are reported as the ratio of the mean log10
DNA isolation. DNA isolation was performed at Rosetta Inpharmatics. DNeasy tissue kits from QIAGEN were used to carry out all DNA extractions. For each liver sample, 20–30 mg of liver was placed in a 1.5-ml microcentrifuge tube along with 80 μl buffer ATL and 20 μl proteinase K. The contents of each tube were then mixed thoroughly by vortexing, followed by incubation at 55 °C until the tissue was completely lysed. Transcriptionally active tissues such as liver and kidney contain high levels of RNA, which will co-purify with genomic DNA. Because RNA-free genomic DNA was required for processing, 4 μl RNase A (100 mg/ml) was added and mixed by vortexing, followed by incubation for 2 min at room temperature before continuing. Samples were then vortexed and 200 μl buffer AL was added to the sample and mixed thoroughly. After 10 min incubation at 70 °C, 200 μl ethanol (96%–100%) was then added and mixed again. The mixture was placed into the DNeasy Mini column and centrifuged at 6,000g (8,000 rpm) for 1 min. The DNeasy Mini spin column was then placed in a new 2-ml collection tube, and 500 μl buffer AW1 was added, followed by placement in a centrifuge for 1 min at 6,000g (8,000 rpm). The DNeasy Mini spin column was then placed in a new 2-ml collection tube again, and 500 μl buffer AW2 was added and centrifuged for 3 min at 20,000g (14,000 rpm) to dry the DNeasy membrane. Then the DNeasy Mini spin column was placed in a clean 1.5-ml or 2-ml microcentrifuge tube and 200 μl buffer AE was pipetted directly onto the DNeasy membrane. This was incubated at room temperature for 1 min and then centrifuged for 1 min at 6,000g (8,000 rpm) to elute. Two 200-μl elutions were performed followed by ethanol/sodium acetate precipitation and resuspension of the resultant pellet with TE buffer.
Genotyping data from the Affymetrix 500K panel. SNP genotyping was performed with the commercial release of the Affymetrix 500K genotyping array. The genotyping was carried out at the Perlegen genotyping facility in Mountain View, California. Genotyping was attempted on 548 samples. 18 samples were unable to be genotyped because of poor DNA quality, and an additional 13 samples were removed after genotyping because their overall call rate did not exceed the 90% cutoff we required. We then applied SNP-wise quality checks on the 517 samples that were successfully genotyped. The Affymetrix 500K array consisted of 500,568 SNPs in total, 429,545 SNPs provided quality data from the genotyping assay, and we rejected those SNPs with a call rate < 75%, resulting in a final panel of 393,494 SNPs. We further filtered out SNPs with minor allele frequencies < 4% (81,646 SNPs) or SNPs that deviated from Hardy-Weinberg equilibrium (p < 10−4; 1,104 SNPs). The resulting set of 310,744 SNPs were used to carry out tests for association to the liver gene expression traits in the HLC.
Genotyping data from the lllumina 650Y panel.
SNP genotyping was performed on the same set of samples that were genotyped on the Affymetrix 500K panel using the Sentrix humanHap650Y genotyping beadchip from Illumina. The genotyping was carried out at the Illumina genotyping facility in La Jolla, California. This chip consists of 655,352 tag SNP markers derived from the International HapMap Project (http://www.hapmap.org
) on a single BeadChip, with ~100,000 Yoruba-specific tag SNPs to provide more comprehensive coverage in African and African-American populations. Genotyping was attempted on 517 samples. A total of 497 samples were genotyped successfully, and 654,069 SNP assays genotyped successfully. The same genotype quality control measures applied to the Affymetrix 500K dataset were applied to Illumina HumanHap 650Y dataset to determine the analysis set. The sample set for analysis (n
= 397) was restricted to those identified or imputed as Caucasian. Of the 397 samples we attempted to genotype, 13 failed the Illumina genotyping assay (overall call rate < 75%), resulting in a set of 384 genotyped samples carried forward for the expression analysis. In total, 652,648 SNPs were called, with only two SNPs rejected because the call rate was <75%. We then sequentially removed 94,915 SNPs with MAF <4% and 491 SNPs that deviated from the Hardy-Weinberg equilibrium (p
). The resulting set of 557,240 SNPs was used to carry out tests for association to the liver gene expression traits in the HLC.
A total of 85,508 SNPs were represented in both the Illumina and Affymetrix SNP sets. Therefore, there were 782,476 unique SNPs successfully genotyped in the HLC such that the call rate was greater than 75%, the MAF >4%, and there was not significant deviation from Hardy-Weinberg equilibrium at the 0.0001 significance level. The sample set for analysis was restricted to the 427 HLC samples that had both genotype and gene expression data available, passed the criteria outlined above and those that were identified as Caucasian, or imputed to be Caucasian when data was missing (see below).
Sex confirmation. Sex identifiers were available for most of the liver samples obtained from the three study centers. We independently confirmed the sex of each individual providing a liver sample by two methods. First, we looked for expression of Y-specific genes in the liver gene expression based on three probes representing three distinct transcripts. Second, we scored heterozygosity of X-chromosome markers. We excluded any individual for which there was a discrepancy in any of the three measures of sex in order to ensure a coherent data set for analysis and that we had excluded as many potential cases of annotation or sample-handling errors as possible. For samples where sex was not noted in the records, we imputed the sex call if both the genotype and gene-expression data were concordant.
. Ethnicities were confirmed or imputed using STRUCTURE [45
]. A panel of 106 autosomal markers was randomly selected from around the genome to be unlinked and ancestry informative. Markers were selected from the HapMap data [46
] that were present on the Affy 500K panel such that the minor allele frequency was >0.05 and the absolute allele frequency difference in the Caucasians and African Americans ~0.5, with average minor allele frequency 0.5 (standard deviation = 12). Several K
were tested (K
= 1–6) with burn-in 100,000 and 100,000 reps of MCMC before any information was collected. In all cases, the greatest support was for K
= 2. Admixture was detected for some individuals in some runs and some individuals were reclassified. For those unknown and reclassified, population reassignment was made if the probability of group membership was >0.9 for that individual. This resulted in 469 individuals assigned to the Caucasian group, 28 individuals assigned to the African descent or African American group, and 18 individuals assigned as “unknown”. The data set for further analysis was restricted to Caucasian samples.
. Ages were imputed using the Elastic Net method [47
]. This method performs model selection and parameter estimation in a manner that is a combination of ridge-regression and the lasso. The prediction method is also explained in [47
]. For computational reasons, λ was set to zero, in which case the Elastic Net method reduces to the lasso method. For most applications, experience demonstrates that the optimal value for λ is zero or quite near zero.
Ages were imputed using separate models for each data source, due to evidence of a source effect, and each sex separately. In cases where the sex was missing or the reported sex was different from the sex implied by the expression data, the sex implied by the expression data was used. This was done so that in the case the annotation data and expression data were mismatched, the imputed age would correspond to the data used to predict it.
The 5,000 genes with the highest correlation to age were used as potential regressors. Cross-validation was used to select the number of steps in the model selection procedure. The number of predictors in the model was between 67 and 76 for the four different models. The percentage of variation explained in the training set is quite high (97%–99%) for three of the models. For the fourth, the model for Vanderbilt females, the percentage of variation explained was slightly lower, 0.92. This is a vast improvement over more naïve imputation methods that are used when adjusting for covariates with missing data, where mean values of the nonmissing data are used to fill in the missing values. Very few of the predictors we constructed were common between the different models. Given the number of predictors with high correlation to age, this is not surprising. Nonetheless, within a given data source (i.e., Pittsburgh or Vanderbilt samples), the male model is a reasonable predictor for the ages of the females and vice-versa. This same trend did not hold for predicting the ages of same-sex individuals across data sources.
Statistical and data visualization methods.
Expression trait processing. Expression traits were adjusted for age, sex, and medical center. Residuals were computed using rlm function from R statistical package (M-estimation with Tukey's bisquare weights). In examining the distributions of the mean log ratio measures for each expression trait in the HLC set, we noted a high rate of outliers. As a result, we used robust residuals and nonparametric tests to carry out the association analyses in the HLC. For each expression trait, residual values deviating from the median by more than three robust standard deviations were filtered out as outliers.
Genome-wide eQTL association analysis. The Kruskal-Wallis test was used to determine association between adjusted expression traits and genotypes. We chose this nonparametric method because of its robust nature to underlying genetic model and trait distribution. p
-Values were computed using nag_mann_whitney (for loci with two observed genotypes) and nag_kruskal_wallis_test (for loci with three observed genotypes) routines from NAG C library (http://www.nag.co.uk
). We used FDR for multiple-test correction. FDR was estimated as the ratio of the average number of eQTLs found in datasets with randomized sample labels to the number of eQTLs identified in the original data set. Since the number of tests was large (~1,010), we found the empirical null distribution was very stable and three permutation runs were sufficient for convergence to estimate FDR. FDR computation was performed separately for cis
(<1 Mb probe to SNP distance) and trans
associations resulting in nominal p
-value cutoffs of 5.0 × 10−5
and 1.0 × 10−8
Targeted set association analysis. The 3,346 SNPs identified in the first round of analysis as associating with expression traits in cis at an FDR < 0.1 were picked for a second round of analysis. To assess the significance of the resulting set of expression traits detected as associated with this set of SNPs, sets of randomly selected SNPs of size 3,346 with MAF distributions identical to the original set were generated. All sets of SNPS were then analyzed using the same method described above for genome-wide associations.
Identifying differentially expressed genes.
To assess whether a gene in a given sample was differentially expressed, we used a previously described and validated error model for testing whether the mean log ratio of the intensity measures between the experiment and reference channels was significantly different from zero [42
]. Based on this error model we obtained p
-values for each of the individual gene expression measures in each sample as previously described [33
]. We then computed the standard deviation of –log10
of the p
-value for each gene expression measure over all samples profiled for a given tissue, and then rank ordered all of the genes profiled in each tissue based on this standard deviation value (rank ordered in descending order). Genes that fall at the top of this rank ordered list can be considered as the most differentially expressed or variable genes in the study. We have previously shown that this type of ordering approach well captures the most active genes in a set of samples [33
]. For demonstrating the number of genome-wide significant eQTLs and eSNPs as a function of differential gene expression, we binned the expression traits into quartiles (Q1-Q4) based on the rank-ordered gene list, with each bin containing 10,025 genes and the bins increasing in significance with respect to differential expression, from Q1 to Q4.
Visualization of networks. Networks were visualized using the Target Gene Information (TGI) Network Analysis and Visualization (NAV) desktop application developed at Rosetta Inpharmatics. This tool enables rapid, real-time, graphical analysis of pathway network models built from a comprehensive and fully integrated set of public and proprietary interaction databases available through a back-end central database, described in detail in a separate report. Additionally, the TGI NAV tool supports experimentally generated systems biology data such as the statistical associations and causal relationships described here. TGI NAV enables integration and visualization of orthogonal data sets using network models as a framework and facilitates dissection of networks into smaller, functionally significant subnetworks amenable to biological interpretation.
To construct the local networks for H2-Eb1, Erbb3, and Rps26, the whole-gene probabilistic causal networks were loaded into the database and the TGI NAV tool was used to extract all edges from this network involving the central gene of interest. In the case of the Erbb3 network, the local network was expanded by extracting all additional edges involving any genes directly connected to Erbb3. Note that while the underlying networks describe causal relationships between transcripts, TGI NAV was used to translate this network into the space of genes using an integrated mapping database that clusters transcripts into gene models utilizing their genomic coordinates. As a result, multiple causal relationships between gene pairs can be observed in cases where multiple transcripts for a single gene were profiled. Visualization properties of nodes (e.g., color) are specified in TGI NAV either for individual nodes, or in a data-driven manner by associating attributes, such as KEGG pathway membership, with groups of nodes and mapping visualization properties to these attributes.