Linkage and association analyses of principal components in expression data

Principal components theory

In microarray data, PCA provides orthonormal bases for the expression array profiles and the gene transcriptional responses. The details are well described elsewhere [4].

Linkage analyses

Our goal was to perform PCA on several different subgroups of expression phenotypes and then perform linkage analysis on the component scores. We performed the PCA in the R language [5]. In order to more closely match the sibling analysis in [1], we only used phenotypes from the sets of siblings in the third generation. In particular, we selected one sibling at random from each sibship and performed the PCA using those 14 independent individuals (assuming that the families are unrelated). Then the loadings were used to compute the components on the remaining siblings. This allowed us to perform linkage on a different set of individuals than those used for the PCA. To determine whether 14 individuals were sufficient to generate reliable components, we performed repeated samplings of 14 individuals at random from the third generation and re-generated component scores. We then examined how correlated the scores were on the overlapping individuals. We also examined the correlation with component scores generated using all 110 individuals from the third generation. Identity-by-descent matrices for linkage analysis were generated using Loki [6] on single-nucleotide polymorphisms (SNPs) thinned in order to avoid elevation from linkage disequilibrium [7]. Linkage analysis was performed using the variance-component (VC) methodology implemented in SOLAR [8].

We investigated the set of expression phenotypes linked to chromosome 14 by Morley et al. [1]. This phase consisted of 32 analyses: the expression levels of the 31 genes identified by Morley et al. [1], analyzed separately, and an analysis of the trait formed by the first component of the combined expression levels. We also examined the first component from the multivariate trait of all 3554 expression phenotypes. Finally, we used the provided gene expression data to examine PCA on three separate categories: 145 genes from the cytoskeleton category, 89 genes from the protein modification category, and 255 genes from the cell cycle category. The cytoskeleton and protein modification categories had been identified by the data providers as the most heavily represented categories among genes with highly variable expression [9]. The cell cycle category has been studied by other researchers using PCA techniques [2, 3].

Association analyses

For the strongest linkage signals, we downloaded data from the HapMap project for the CEPH families for SNPs within the two-LOD support [10]. We performed association analyses using an additive model to code the SNP genotypes using SOLAR. We also included age in the model when it was associated with the principal component. In this VC setting, age and SNP genotypes are treated as covariates predicting trait phenotype while a kinship matrix scaled by the trait heritability controls for the covariance from related individuals. We performed a two-stage analysis: first, we tested association with all HapMap non-synonymous SNPs; then, we tested association with all HapMap SNPs in the regions. We used a Bonferroni correction as well as a false-discovery rate (FDR) method (QVALUE [11]) to evaluate the significance of association in these two stages. The strategy of first focusing on non-synonymous SNPs has been applied to complex diseases [12]. Although mRNA expression levels of a single gene can be modified by polymorphisms across the gene footprint [13], in this context we are attempting to identify master regulators. We posit that non-synonymous changes may result in functional differences in proteins that regulate the mRNA expression of many other genes.

High heritability genes

For a number of phenotypes, SOLAR reported an estimated heritability of 1.0. Further investigation revealed that the individual expression levels as well as the component scores often showed an intraclass correlation (ICC) of more that 50% in sibships. Heritability is often estimated as twice the intraclass correlation in sibships; this would imply a heritability greater than one. This is clearly nonsensical, so we suspect that this indicates a shared environment with an effect on multiple expression levels. We performed a PCA on the 471 genes that had ICC in sibships greater than 0.5 (the "high heritability" genes). We also performed further tests to determine whether these genes contained an over-representation of ontological categories. In particular, we examined the 10 most represented ontological categories in the high heritability genes, and then performed 10,000 re-samplings of 471 genes from all 3554 to determine whether any of these 10 categories were over-represented.

Linkage analyses

We first examined the reliability of the PCA. Using the entire set of phenotypes, we found that the components generated from 14 unrelated individuals had a 98.7% correlation with the components generated from all 110 third-generation individuals. Furthermore, 10,000 repeated samplings of 14 unrelated individuals produced an average correlation of 96.1% with the 14 individuals initially chosen. We took this as evidence that the components were reliable, even using a relatively small number of individuals; we then used the same 14 individuals throughout. The PCA analysis of each of the six categories of gene expression then produced 14 principal components. Recently, Raiche et al. [14] developed a new technique to find the number of components to retain using an "acceleration factor," a measure of decrease in proportion of variance. For the cell cycle case, the analysis indicated that two components should be retained; in the other six cases, only the first component was retained. We then performed VC linkage analysis on these seven traits. We also performed linkage analysis for the separate expression traits reported as linked to chromosome 14 by Morley et al. [1] in order to replicate those findings.

Association analyses

For the association analysis, we focused on the five regions with linkage signals over 3.5 for the traits linked to chromosome 14, including the linked region on chromosome 14. We found the first principal component of these traits was strongly associated with age (p = 9.8 × 10^-6), so we included age in all SNP analyses. We identified 761 nonsynonymous SNPs in these regions, so we set our initial significance threshold to a conservative level of 0.05/761 = 6.6 × 10^-5. One SNP, rs10458896, is significant at this threshold (p-value of 4.5 × 10^-5). The corresponding q-value is 0.034; this is the only q-value less than 0.05. This is a non-synonymous SNP in KIF18A, a kinesin family member on chromosome 11. In the analysis of all the SNPs, we found 143,798 SNPs in these regions. We set our significance threshold to 0.05/143798 = 3.5 × 10^-7. One SNP in an intergenic region on chromosome 11, rs10768321, is significant at this threshold with a p-value of 5.8 × 10^-8 and a q-value of 0.008; again, this is the only q-value less than 0.05.

Ontological categories of "high heritability" genes

For the analysis of the top ten ontological categories represented by the 471 "high heritability" genes, we computed empirical p-values based on how frequently a random selection of 471 genes contained more genes of the category in question than in the high-heritability group. We observed an over representation of the categories "nucleus" (p = 0.0006), "nucleotide binding" (p < 0.0001), "ATP binding" (p = 0.0001), and "RNA binding" (p = 0.0002).