Gene coexpression network analysis for family studies based on a meta-analytic approach

Abstract

For a better understanding of the biological mechanisms involved in complex traits or diseases, networks are often useful tools in genetic studies: coexpression networks based on pairwise correlations between genes are commonly used. In case of a family-based design, it can be problematic when there is a large between-family variation in expression levels. We propose here a gene coexpression network analysis for family studies. We build a coexpression network for each family and then combine the results. We applied our approach to data provided for analysis in the Genetic Analysis Workshop 19 and compared it to 2 naïve approaches—ignoring correlations among the expressions and decorrelating the gene expression by using the residuals of a mixed model—and a single-probe analysis. Our approach seemed to better deal with heterogeneity with regard to the naïve approaches. The naïve approaches did not provide any significant results, while our approach detected genes via indirect effects. It also detected more genes than the single-probe analysis.

Background

Weighted gene coexpression network is a widely used method for studying biological networks based on pairwise correlations. This method provides more insight in the underlying biological mechanisms and offers a tool for dimension reduction by summarizing identified modules (clusters) of genes [1, 2]. How to perform such an analysis for family data is an open question. For family data Kraft et al. [3] noted that testing association between expression levels and traits without taking into account the family structure can lead to spurious results, especially when the number of families is small and in the presence of large between-family variation. In this paper, we propose a novel strategy for network analyses in a small set of relatively large families.

For this family-based approach, we first construct family-specific coexpression networks and test for association between the modules and the traits of interest. A common set of genes for all families was obtained by using the intersection and the union of family-specific modules. We compare this family-based approach with 2 naïve approaches: namely, one using the gene expression of the families directly (ignoring correlation) and one that first decorrelates the gene expressions and then applies the standard approach. We also compare our results with single-probe analyses.

Methods

Study sample

The gene expression data set is composed of 647 individuals from 17 large families. These samples are from the data set described in Almasy et al. [4]. Here, we focus on the largest 5 families: namely families 2, 5, 6, 8, and 10 with 65, 55, 45, 62, and 49 family members, respectively. The total number of individuals is 276 and the total number of probes from which gene expression is available is 20,364. We used the simulated quantitative phenotypes systolic blood pressure (SBP) and the phenotype Q1 at time point 1 as outcome variables. The simulation model of SBP comprises 15 genes and that of Q1 does not contain any of these genes. SBP, Q1, and all probes were corrected for age and sex by regressing out covariates and using residuals.

To decorrelate the gene expressions, we fitted for each probe a linear mixed model: X _ij  = μ + u _ij  + v _i  + ε _ij , with X _ij the value of the probe for the individual j in family i, u _ij a normally distributed random genetic effect: u _ij  ~ N(0, S) where S = 2 * K * s _g with K kinship matrix and s _g genetic variance, v _i a normally distributed random effect representing shared environmental effects, and ε _ij a normally distributed residual. To obtain the residuals X ^* _ij of this model, we used the function lmekin, which fits linear mixed models with specific structure of the variance-covariance matrix from the package coxme [5] in R.

Single-probe analysis

For the single-probe analysis the following mixed model was used:

$$ {Y}_{ij} = \mu +{u}_{ij}+{v}_i+\beta {X}_{ij}+{\upvarepsilon}_{ij} $$

(1)

with Y _ij the value of SBP or Q1 and X _ij the value of the probe for individual j of family i. The random effects u _ij , v _i and ε _ij are the genetic effect, the shared environmental effect and residuals respectively. The parameter β represents the effect of the probe on the outcome variable.

Network constructions

Coexpression networks were built on the data set without correction for family structure based on X _ij (naïve approach), the data set adjusted for family structure based on X ^* _ij (naïve decorrelated approach), and on the data sets from the 5 families separately.

We used signed coexpression networks. The adjacency matrix A = [a _lk ] of each network was computed as follows: a _lk  = |0.5 + 0.5 × cor(x _l , x _k )|^γ, with cor(x _l , x _k ) the correlation between x _l the values vector of probe l and x _k the values vector of probe k. The parameter γ is acting as a soft threshold in the adjacency matrix, when we increase the value γ the coefficient of the adjacency matrix will tend toward zero except for values really close to 1. We used the biweight midcorrelation based on the median, which is more robust than the Pearson correlation. The co-expression networks were constructed with the R package WGCNA (weighted gene correlation network analysis) [6]. For each obtained module, the first principal component (eigengene) was computed.

Phenotype analysis

From all modules and all families, the following models were fitted:

$$ {Y}_j = \mu +{u}_j+\beta eigengen{e}_j^k+{\upvarepsilon}_j $$

(2)

where Y _j is the outcome, u _j the random genetic effect, and eigengene_j ^k the value of the eigengene of module k of family member j. Let E ^M _F2 to E ^M _F10 be the most significant eigenvalues of the family specific networks (N_F2 to N_F10) and let E ^M _F be the most significant eigenvalue of these 5 eigenvalues and M ^M _F be the corresponding module. Identify the modules of the family-specific networks, which have the highest overlap with M ^M _F (denoted as M ^O _F2 to M ^O _F10 ). Next, 2 common sets of genes for all families were obtained by taking the intersection (M_F = M ^O _F2  ∩ M ^O _F5  ∩ M ^O _F6  ∩ M ^O _F8  ∩ M ^O _F10 ) and the union (M_F = M ^O _F2 ∪ M ^O _F5 ∪ M ^O _F6 ∪ M ^O _F8 ∪ M ^O _F10 ) of the family specific modules. The first principal components of the 2 common sets were computed. The principal component that explained most of the variance of the corresponding set of genes was used as the eigengene E_F of the family-based approach.

The eigengenes of the naïve approach (E_N), the naïve approach after decorrelation (E_ND), and the family-based approach (E_F) are tested for association with the 2 phenotypes SBP and Q1. Here, the following mixed model was used:

$$ {Y}_{ij} = \mu +{u}_{ij}+{v}_i+\beta eigengen{e}_{ij}^k+{\upvarepsilon}_{ij} $$

(3)

with Y _ij the phenotype value for individual j of family i and eigengene ^k _ij the value of eigengene of module k of individual j of family i. And u _ij , v _i and ε _ij are the genetic effect, the shared environmental effect and residuals, respectively. The parameter β represents the effect of the eigengene k on the outcome variable.

Finally, because spurious associations are especially expected in the presence of large between-family heterogeneity [1], we also performed a network analysis using the subset of 25 % most heritable probes when performing the network analysis (n = 4911 probes with heritability between 0.33 and 0.88).

To test for significance we used a nominal alpha level of 0.05 and the Bonferroni correction was applied to take into account multiple testing.

Results

Results obtained with all probes

For per family analysis, the module that showed the highest correlation with the SBP was the magenta module obtained in family 8 (M ^M _F8 ) (β = 2.52, p = 0.0011). M ^M _F8 comprises 710 genes. Table 1 gives for each family the number of genes of the module with the highest overlap. The intersection and the union of these 5 family modules, comprises 62 and 1746 probes, respectively. The first principal component (eigengene) of the probes in the intersection set explained more than 50 % of the variance for each family, while for the union set the eigengene explained only between 23 % and 31 % of the variance of the expression levels. Therefore the eigengene of the intersection set was used as summary for the family approach (E_F). In Table 2, for each family the effect of E_Fi on SBP (β of model (2)) is given. For families 2 and 8, the eigengenes (E_F2 and E_F8) were significantly associated with SBP.

Table 1 Module size of M ^O _F2 to M ^O _F10 and overlap size with M ^M _F in the all-probes analysis

Table 2 Parameter estimates of the association between eigengenes and Q1 and SBP

When analyzing all families together, none of the approaches provided significant results. The joint analysis of the families using E_F as eigengene in model (3) did not provide a significant association SBP (β = −0.13, p = 0.49). For the naïve approach, the eigengene of the module magenta (E_N) had the smallest p value (β = −3.21, p = 0.01). For the naïve approach using the decorrelated data set, the eigengene of the module grey60 (E_ND) had the smallest p value (β = −3.03, p = 0.0061). After multiple-testing correction (between 43 and 50 modules in each network) none of the results were significant. Finally the single-probe analysis preformed in the 5 families by using model (1) provided 1 significantly associated probe with SBP (CRIP2; β = −13.68, p = 1.7e-06).

The intersection module of the family based approach did not contain any of the 15 genes used for the simulation. Also, the identified gene of the single-probe analysis is not among these 15 genes. We hypothesized that correlation might exist between E_F2, E_F8, and the gene expression of CRIP2 on one hand and the set of 15 genes on the other hand. Indeed, E_F2 showed significant correlation with PSMD5 (p = 0.004) and GTF2IRD1 (p = 0.007), and E_F8 showed significant correlation with ZNF443 (p = 5e-05), PSMD5 (p = 3e-05), and ABTB1 (p = 6e-05). When the presence of these 15 genes in the modules was investigated, it appeared that they were in different modules (Table 3). The gene CRIP2, which was significant in the single-probe analysis, showed significant correlation with the gene KRTAP11-1 (p = 3.1e-03).

Table 3 List of the top genes involved in the simulation model and their module color in each network

Analysis of Q1

The results of the analysis of Q1 are also given in Table 2. For the family approach, none of the modules obtained in family-specific network analysis was significantly associated with Q1 and no common set could be defined. In Table 2, the estimates of strongest associated modules E ^M _F for each family are given. For the naïve approach, the module red (E_N) (β = 2.75, p = 0.00056) was significant and for the naïve approach using the decorrelated data the module green (E_ND) (β = −1.41, p = 0.00094) was significantly associated with Q1.

Results obtained with the 25 % most heritable probes

For the naïve and the family approaches, the results of the network-based analyses using only the gene expressions of the 25 % most heritable probes (n = 4911 probes with heritability between 0.33 and 0.88) are also given in Table 2. None of the 15 genes used in the simulation model for SBP was among this set of most heritable probes. For Family 6 the E_F6 was significantly associated with SBP (p = 0.01). The association of E_F in the 5 families with SBP was also significant (p = 0.02). For Q1, none of the approaches provided significant results. With regard to the single-probe analysis, no probe other than CRIP2 was significantly associated.

Discussion

In this paper, we have proposed a novel strategy to perform a coexpression network analysis with family data: building a network for each of the large pedigrees and defining a common module by taking the intersection of family-specific modules. We compared our family-based approach with 2 naïve network approaches and a single-probe analysis. All analyses were performed in a small set of 5 relatively large families. None of the 15 genes in the simulation model was identified in this small data set. However, the family-based approach identified significant associations between the eigengene and SBP in 2 families. This eigengene was significantly correlated with 4 of the 15 genes. When analyzing all families jointly the family eigengene was not significant. Also the naïve network approaches did not provide any significant result. The single-probe analysis provided 1 significant gene, which was correlated with 1 of the 15 genes. To study the performance of the methods with regard to false-positive findings, we also analyzed the trait Q1 for which no gene expressions were included in the simulation model. The family approach did not provide a significant finding, whereas both naïve approaches identified a significant module for Q1. The result in the naïve approach based on gene expression (X _ij ) is in line with the findings of Kraft et al. [3]. We did not expect to have a false-positive finding when using the decorrelated data (X ^* _ij ) as input for our network analysis. Possible explanations for this finding are the fact that the correlation based on the kinship coefficient might not be appropriate for gene expressions, and randomness. In addition to the set of all probes, networks were also built using only the 25 % most heritable probes. Especially for these variables that show large between-family variation spurious associations might occur when the family structure is not taken into account. This was not confirmed in our analysis. More research is needed to study the sensitivity of the methods for between-family variation.

We did not know the answers when we developed the family-based approach and analyzed the data. The simulation model used to create the data sets may not be well suited to pick up the 15 genes directly by network analysis. The 15 genes present in the model were in different pathways: they were not correlated. Moreover our approach was able to identify indirect effects: that is, the eigengenes were correlated with the 15 genes. Thus the significant association of the family-based network approach represented the largest number of genes from the simulation model. We expect that especially in the presence of large between-family variation our approach would perform best. A thorough simulation study is required to investigate the performance of our method further.

Network analysis provides a tool to reduce the number of tests by first summarizing the data in sets of genes with correlated gene expressions and summarizing the gene set by the first principle component. Another obvious reduction step is to only consider the heritable probes for the analysis. It appeared that by using the heritable probes the results across the families were less heterogeneous. The family approach, as well as the naïve approach using decorrelated data, provided significant results for SBP.

Conclusions

In this paper we combined the family-specific modules by taking the intersection of the modules which showed most overlap. This approach worked well for the relatively small set of 5 families. When we applied our method to the 6 largest families, similar results were obtained (data not shown). However, intersection might not be the most appropriate approach to combine modules across families, because the intersection set becomes too small. Alternative approaches have to be developed. For example, LASSO (least absolute shrinkage and selection operator)-type methods can be used to select probes from the union sets. Development of a method for constructing a common set from the family specific modules is a topic for future research.

Finally more research is needed to evaluate the performance of our method with regard to false-positive and false-negative findings in relationship to heterogeneity, family size, the number of families, and the heritability of gene expressions.