Abstract
In the fastdeveloping field of expression quantitative traits loci (eQTL) studies, much interest has been concentrated on detecting genomic regions containing transcriptional regulators that influence multiple expression phenotypes (transhubs). In this paper, we develop statistical methods for eQTL mapping and propose a new procedure for investigating candidate transhubs. We use data from the Genetic Analysis Workshop 15 to illustrate our methods. After correlations among expressions were accounted for, the previously detected transhubs are no longer significant. Our results suggest that conclusions regarding regulation hot spots should be treated with great caution.
Background
The emerging microarraybased gene expression technology enables quantitative geneticists to conduct systematic, genomewide linkage analysis to detect genomic loci that control geneexpression variations. One of the common features of expression quantitative trait loci (eQTL) studies is the detection of transhubs, "chromosomal regions that affect the expression of a much larger number of genes than expected by chance" [1]. However, a major concern in transhub detection is the high falsepositive rate due to the complex correlation structure of gene expressions [1]. If a group of genes are highly correlated and a QTL is detected for one of them, then there is a large chance that other expression phenotypes in this group are also mapped to the same QTL, regardless of whether the reason for coexpression is indeed coregulation at this QTL.
In this paper, we first describe a robust score statistic method designed for threegeneration pedigrees for linkage detection. We then introduce a new method for investigating candidate transhubs. To account for correlations among gene expressions, we treat the expression of additional genes as covariates in the variancecomponent model of a target gene, and employ sparse regression techniques to remove the covariates' components before testing linkage. The effects of accounting for expression correlations in linkage analysis are illustrated in the Results.
We apply the proposed methods on Genetic Analysis Workshop 15 (GAW15) Problem 1 data [24]. This data set consists of 14 threegeneration Centre d'Etude du Polymorphisme Humain pedigrees. Genotypes of 2882 autosomal and Xlinked SNPs are provided for each individual. Expression levels of ~8500 genes in lymphoblastoid cells were obtained using the Affymetrix Human Focus Arrays [4].
Methods
Variancecomponent model and score statistics
In this section, we derive robust score statistics under the variancecomponent model [5] to map QTLs. We assumed HardyWeinberg equilibrium and linkage equilibrium throughout. The singlelocus additive model for a phenotype Y having a mean value μ is given by
where α_{i }= α_{i}(τ) denotes the additive genetic effect of allele x at locus τ (the subscript m or f denotes the parental origin of the allele). We assume that E(α_{m}) = E(α_{f}) = 0, that α_{m }and α_{f }are uncorrelated, and that the residual term e is uncorrelated with the explicitly modeled genetic effects. It is straightforward to compute the conditional covariances given the identitybydecent (IBD) sharing numbers under Eq. (1). For example, for two siblings i, j, , where ν_{ij}(τ) is the IBD sharing number between the two siblings at the trait locus τ, ρ_{s }is the phenotypic correlation among siblings and is the linkage parameter.
For each phenotype Y, the null hypothesis of no linkage can be written as H_{0}: α_{0 }= 0. The working assumption for the variancecomponent model is that the conditional distribution of the phenotypes in a pedigree given the pairwise IBD sharing information at a QTL is multivariate normal. At marker t (a putative trait locus), the score statistic for testing α_{0 }= 0 is then , where n is the number of pedigrees and
Here, Y_{i }denotes the vector of phenotypes in the i^{th }pedigree; A_{i}(t) is the IBD sharing matrix for the i^{th }pedigree at marker t:
where is the estimated IBD sharing number between the j^{th }and κ^{th }member at marker t; μ, σ_{Y}, Σ are the phenotypic mean, variance, and correlation matrix, respectively. We consider three different types of phenotypic correlations for the threegeneration pedigree: sibship correlation (ρ_{s}); grandparentgrandchild correlation (ρ_{g}); and parentoffspring correlation (ρ_{o}). All of these nuisance parameters are estimated by their corresponding sample estimators. We then standardize the above score by its conditional variance given the phenotypes: . Let a_{i,t}(j, κ) = A_{i}(t)(j, κ), then the calculation of the above quantity involves estimation of E[a_{i,t}(j, κ)a_{i,t}(j'k')]. We identify 11 different types of cross products for two pairs (j, κ) and (j'k') which have a nonzero expectation (e.g., (j, κ) and (j'k') are two sib pairs with one sibling in common). We pool the same type of cross products across all pedigrees and estimate the above expectation by sample mean.
We then define the robust score statistic at marker t as , which is asymptotically normally distributed with mean zero and variance one, no matter what the actual phenotypic distribution is. Because we do not know the location of the QTL τ, we scan the whole genome with the test statistic: Z_{max }= max_{t}Z(t), where the maximum is taken over all marker loci t throughout the genome.
Investigation of transhubs
When linkage exists between a genome region and an expression phenotype, the regulation can be "indirect" and act through one or more intermediate genes (that is, this region regulates some intermediate genes and their expression in turn regulate the phenotype). Such indirect regulations are usually less interesting. To detect biologically interesting transhubs, only direct linkage would be meaningful. In this section, we propose a method to distinguish direct and indirect regulations/linkages.
We will illustrate the idea through a simple example. Consider a system of three components: one candidate QTL (X) and two expression phenotypes: Y_{1}, Y_{2}. It can be shown that if both Y_{1 }and Y_{2 }are linked to X and if the linkage strength (defined as the proportion of variation explained by X) of Y_{2 }is no greater than that of Y_{1}, then the system will match to one of the two models: a) X regulates Y_{1 }and X regulates Y_{2 }(connection between Y_{1 }and Y_{2 }is allowed); or b) X regulates Y_{1 }and Y_{1 }regulates Y_{2}. (Due to limitation of space, detailed models and proofs are omitted.) We need to distinguish these two models in order to decide whether the linkage between X and Y_{2 }is direct (Model (a)) or indirect (Model (b)). This can be revealed through investigating the residual R_{21 }of the regression model Y_{2 }~ Y_{1}: under Model (a), R_{21 }links to X; while under Model (b), R_{21 }does not link to X. On the other hand, if we consider the regression model Y_{1 }~ Y_{2}, the residual R_{12 }will link to X under both models. However, the linkage might be weak. Therefore, in order to avoid performing unnecessary linkage tests on residuals, which decreases the power, we propose to first order the expression phenotypes with respect to their linkage strength at the candidate hub; and then for each expression trait, only those phenotypes with stronger linkage evidence are used as covariates to derive the corresponding residual in the model below. As a result, for any pair of expression traits, there is only one model having the two traits on the opposite sides of the equation.
According to the above discussion, we introduce the variancecomponent model
where a set of expression phenotypes other than Y_{i }are treated as covariates (Y_{i}). Define R_{i }= Y_{i } Y_{i}β. Model (2) becomes R_{i }= μ + α_{m }+ α_{f }+ e, for which the score statistics described previously can be applied to test linkage. Thus, the remaining task is to properly derive R_{i}: the residual of the regression model Y_{i }~ {Y_{i}}. Because of the high dimensionality of the expression phenotypes ({Y_{i}}), it is crucial to maintain sparsity in the regression models to avoid overfitting. For this purpose, we apply a sparse regression technique called elastic net [6] to derive R_{i}. Elastic net aims to minimize the loss function L(λ_{1}, λ_{2}, β) = Y  Xβ_{2}^{2 }+ λ_{2}β_{2}^{2 }+ λ_{1}β_{1}. The ridge penalty term encourages a grouping effect: strongly correlated predictors tend to be in or out of the model together; the lasso penalty term enables the algorithm to have a more sparse representation and thus serves as a model selection tool [6].
We propose the following procedure for investigating a candidate transhub region:
1. Order expression phenotypes according to the linkage strength to this region (based on the score statistics Z at the hub) from the largest to the smallest.
2. For the i^{th }ordered expression Y_{(i)}, perform Elastic net regression Y_{(i) }~ {Y_{(j)}}_{j<i }with λ_{2 }= 1 and maximum step k_{max}. Record the corresponding residue as R_{(i)}.
3. Perform linkage analysis on {R_{(i)}}_{i }using the robust score statistics.
4. An expression trait is considered to have a direct linkage to the candidate region if both the original expression Y_{i }and residual R_{i }show significant evidence of linkage.
For the GAW15 application, the maximum step for running elastic net is set to be k_{max }= 20, which is the mean optimal step chosen using Mallows' C_{P }criterion [7] among 100 randomly picked regression models (in each model, the expression of one randomly chosen gene is regressed on the expressions of all other genes).
Data analysis
We first performed an empirical normal quantile transformation for each gene's expression to make them marginally normal [8], for the purpose of improving power of linkage detection. We want to point out that the validity of our test statistic is robust to the distributional assumption of the phenotypes because it is standardized by the conditional variance of the score statistic [9,10].
For genotype data, 1197 SNPs were selected from 22 autosomes, such that the intermarker distance is at least 0.1 cM to avoid linkage disequilibrium. The resulting map has an average intermarker space of 3.1 cM, mean heterozygosity of 0.42, and mean missing rate of 3.89%. Merlin [11] was used for IBD inference based on the Rutgers sex average linkage map provided by Sung et al. [12]. Linkage tests were performed for the 3554 most variable genes selected by Morley et al. [4] and gender was used as a covariate.
The threshold at each genomewide significance level was estimated based on 500 permutations. In each permutation cycle, pedigrees of the same size were permuted; within a pedigree, phenotypes were permuted among individuals of the same type (i.e., generation and gender). We also calculated the thresholds according to a Gaussian approximation with skewness correction [10], which results in similar thresholds (Table 1). We observed that the confidence intervals derived from permutation are wider than what would be expected if all tests were independent. This also suggests that the correlations among gene expressions should not be ignored when examining false positives.
Table 1. Number of expressions with at least one significant eQTL
Results
The result of eQTL analysis are summarized in Tables 1 and 2. At the genomewide 0.001 significance level (pointwise pvalue < 7 × 10^{9}), we identify 24 expression phenotypes with at least one significant eQTL, among which five overlap with the genes reported in Table 1 of Morley et al. [4] (TM7SF3, HSD17B12, CHI3L2, DSCR2, and DDX17). We also performed genome scans using only sibship data, which, not surprisingly, resulted in a lower power compared to the analysis using whole pedigrees (Table 1).
Table 2. Expression phenotypes with the strongest evidence of linkage from genome scans
Following Morley et al. [4], we define cisregulators as the eQTL that map within 5 megabases (Mb) of the target gene and all other eQTL are defined as transregulators. To illustrate the proposed transhub investigation method, we examine the transeQTL events (based on the original expression data) at the genomewide 0.05 significant level to harvest enough eQTL hits for deriving candidate transhubs. The numbers of transhits and cishits at different significance levels are summarized in Table 3. The number of transhits dropped dramatically as the significance levels became more stringent, indicating an overall lower confidence of the translinkage detection. The genomic locations of the eQTLs detected at the 0.05 significant level suggest two possible hot spots, one at 9p13.3 (30.8 Mb to 35.1 Mb) with 13 transhits, and another at 14q32 (94.7 Mb to 98.2 Mb) with 11 transhits (Figure 1, left panel). The latter region, 14q32, has also been recognized as a candidate transhub by Morley et al. [4]. If regulators for expression phenotypes were independently and randomly distributed along the genome, the probability of the maximum number of hits being at least 13 or 11 are both less than 1 × 10^{6}. However, based on permutation results (in the original data), the chance of observing a transhub with at least 13 hits is about 11.2%, and the chance of at least 11 hits is as high as 22.4%. These numbers imply that the false detections of transeQTL are clustered instead of uniformly distributed along the genome. We hypothesize that this is partly due to the expression correlations.
Figure 1. Distribution of transhit along the genome. The xaxis represents the genome order of the 1197 markers. The yaxis represents the number of transhits in a 5Mb neighborhood region of each marker. Markers on different chromosomes are separated by vertical gray lines. The left panel is for the original expression data. The right panel is for the residual analysis with respect to 9p13.3. The positions of DSCR2 (21q22.3) and MAP3k6 (1p36.11), which show strong evidence of translinkage to 9p13.3 region, are indicated in the left panel.
Table 3. Number of cishits and transhits^{a}
To further investigate our hypothesis, we adjusted for the expression correlations as described in the Methods on the transhub at 9p13.3. The linkage results of residuals are also summarized in Table 1, which show comparable number of eQTL detections as before. At the genomewide 0.001 significance level, 26 expression phenotypes were identified, of which 9 overlap with the original 24 expression phenotypes mapping to the similar chromosomal regions (TBC1D8, HLADPB1, CSTB, BCKDHA, DSCR2, POMZP3, CHI3L2, HSD17B12, and TM7SF3). However, the hub phenomenon becomes much less obvious now: transeQTL hits were very evenly distributed along the genome (Figure 1, right panel). One explanation is that the overall pairwise expression correlations in the residual data are much smaller than those in the original data: the median absolute correlation is 0.052 and 0.139, respectively. The maximum number of transhits of one 5Mb region at the genomewide 0.05 significance level is only 4. The number of transhits at 9p13.3 drops to 3, while among the 500 permutation cycles performed on residual data, none has a maximum number of transhits smaller than 3. The same analysis is done for the transhub at 14q32 and the results are similar (Table 1). The number of transhits at 14q32 drops to 4.
Thus, we conclude that there is not enough evidence to claim either two candidate regions we examined as a transhub. However, we do find two statistically significant, transregulated, phenotypes: DSCR2 (Down Syndrome Critical Region gene 2) is the most significant (pointwise pvalue of 10^{12}) gene linked to the 9p13.3 region; MAP3K6 also shows strong evidence of linkage to the same region according to both the original phenotype and the residual phenotype. MAP3K6 encodes a member of the serine/threonine protein kinase family, and has a pointwise pvalue of 3 × 10^{7}. The correlation between expressions of DSCR2 and MAP3K6 is quite small: 0.074. These two transregulations may deserve further investigation.
Discussion
In this paper, we perform linkage analysis for GAW15 data using robust score statistics, which enjoy excellent computational efficiency (20 seconds for computing the score statistics for 3554 expressions on 1197 markers in R on a Thinkpad X40 laptop), and enable us to carry out largescale permutation studies. Using the original phenotypes, we identify two candidate transhubs, one at 9p13.3 and the other at 14q32. However, after accounting for the expression correlations in the linkage analysis, both transhubs disappear. This suggests that conclusions with regard to regulation hot spots should be interpreted with great caution.
Controlling false positives is one of the most important concerns in processing large high dimensional data sets. Without the controlling of false positives, power is not a meaningful quantity. In this paper, we focused on hubs of direct transregulation, which is conceptually different from the situation where both direct and indirect linkages are sought after. For this purpose, there are two types of false positives: i) the locus and the gene are not linked at all, while a linkage is claimed; ii) the locus and the phenotype are indirectly linked, while a linkage is counted as a direct regulation. The proposed method helps to prevent both types of false positives. Protection against the second type of false positives is discussed in the Methods. As to the first type of false positives, due to correlations among expressions, they do not randomly distribute along the genome. The proposed method also acts as a safeguard against detecting false hubs resulting from this source, since the residuals are usually much less correlated.
Competing interests
The author(s) declare that they have no competing interests.
Acknowledgements
The authors want to thank two anonymous reviewers whose comments led to significant improvements of the paper.
This article has been published as part of BMC Proceedings Volume 1 Supplement 1, 2007: Genetic Analysis Workshop 15: Gene Expression Analysis and Approaches to Detecting Multiple Functional Loci. The full contents of the supplement are available online at http://www.biomedcentral.com/17536561/1?issue=S1.
References

Koning D, Haley C: Genetical genomics in humans and model organisms.
Trends Genet 2005, 21:377381. PubMed Abstract  Publisher Full Text

Cheung V, Conlin L, Weber T, Arcaro M, Jen K, Morley M, Spielman R: Natural variation in human gene expression assessed in lymphoblastoid cells.
Nat Genet 2003, 33:422425. PubMed Abstract  Publisher Full Text

Cheung V, Spielman R, Ewens K, Weber T, Morley M, Burdick J: Mapping determinants of human gene expression by regional and whole genome association.
Nature 2005, 437:13651369. PubMed Abstract  Publisher Full Text

Morley M, Molony C, Weber T, Devlin J, Ewens K, Spielman R, Cheung V: Genetic analysis of genomewide variation in human gene expression.
Nature 2004, 430:743747. PubMed Abstract  Publisher Full Text

Almasy L, Blangero J: Multipoint quantitativetrait linkage analysis in general pedigrees.
Am J Hum Genet 1998, 62:11981211. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zou H, Hastie T: Regularization and variable selection via the elastic net.
J Roy Stat Soc B 2005, 67:301320. Publisher Full Text

Mallows L: Some comments on C_{P}.
Technometrics 1973, 15:661676. Publisher Full Text

Wang K, Huang J: A scorestatistic approach for the mapping of quantitativetrait loci with sibships of arbitrary size.
Am J Hum Genet 2002, 70:412424. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Peng J, Siegmund D: Mapping quantitative traits under ascertainment.
Ann Hum Genet 2006, 70:867881. PubMed Abstract  Publisher Full Text

Tang HK, Siegmund D: Mapping quantitative trait loci in oligogenic models.
Biostatistic 2001, 2:147162. Publisher Full Text

Abecasis G, Cherny S, Cookson W, Cardon L: Merlinrapid analysis of dense genetic maps using sparse gene flow trees.
Nat Genet 2002, 30:97101. PubMed Abstract  Publisher Full Text

Sung YJ, Di Y, Fu AQ, Rothstein JH, Sieh W, Tong L, Thomson EA, Wijsman EM: Comparison of multipoint linkage analyses for quantitative traits in the CEPH data: parametric LOD scores, variance components LOD scores, and Bayes factors.
BMC Proc 2007, 1(Suppl 1):S93. PubMed Abstract  Publisher Full Text  PubMed Central Full Text