Controlling for false positive findings of trans-hubs in expression quantitative trait loci mapping
- Jie Peng†^{1}Email author,
- Pei Wang†^{2} and
- Hua Tang^{3}
https://doi.org/10.1186/1753-6561-1-S1-S157
© Peng et al; licensee BioMed Central Ltd. 2007
Published: 18 December 2007
Abstract
In the fast-developing 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 (trans-hubs). In this paper, we develop statistical methods for eQTL mapping and propose a new procedure for investigating candidate trans-hubs. We use data from the Genetic Analysis Workshop 15 to illustrate our methods. After correlations among expressions were accounted for, the previously detected trans-hubs are no longer significant. Our results suggest that conclusions regarding regulation hot spots should be treated with great caution.
Background
The emerging microarray-based gene expression technology enables quantitative geneticists to conduct systematic, genome-wide linkage analysis to detect genomic loci that control gene-expression variations. One of the common features of expression quantitative trait loci (eQTL) studies is the detection of trans-hubs, "chromosomal regions that affect the expression of a much larger number of genes than expected by chance" [1]. However, a major concern in trans-hub detection is the high false-positive 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 co-expression is indeed co-regulation at this QTL.
In this paper, we first describe a robust score statistic method designed for three-generation pedigrees for linkage detection. We then introduce a new method for investigating candidate trans-hubs. To account for correlations among gene expressions, we treat the expression of additional genes as covariates in the variance-component 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 [2–4]. This data set consists of 14 three-generation Centre d'Etude du Polymorphisme Humain pedigrees. Genotypes of 2882 autosomal and X-linked SNPs are provided for each individual. Expression levels of ~8500 genes in lymphoblastoid cells were obtained using the Affymetrix Human Focus Arrays [4].
Methods
Variance-component model and score statistics
In this section, we derive robust score statistics under the variance-component model [5] to map QTLs. We assumed Hardy-Weinberg equilibrium and linkage equilibrium throughout. The single-locus additive model for a phenotype Y having a mean value μ is given by
Y = μ + α_{ m }+ α_{ f }+ e, (1)
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 identity-by-decent (IBD) sharing numbers under Eq. (1). For example, for two siblings i, j, $Cov({Y}_{i},{Y}_{j}|{\nu}_{ij}(\tau ))={\rho}_{s}{\sigma}_{Y}^{2}+{\alpha}_{0}({\nu}_{ij}(\tau )-1)$, where ν_{ ij }(τ) is the IBD sharing number between the two siblings at the trait locus τ, ρ_{ s }is the phenotypic correlation among siblings and ${\alpha}_{0}=E({\alpha}_{m}^{2})=E({\alpha}_{f}^{2})$ is the linkage parameter.
where ${\widehat{\nu}}_{j\kappa}$ 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 three-generation pedigree: sibship correlation (ρ_{ s }); grandparent-grandchild correlation (ρ_{ g }); and parent-offspring 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: ${I}_{\alpha \alpha}(t)={E}_{0}[{\ell}_{i}^{2}(t)|{Y}_{1},\cdots ,{Y}_{n}]={\displaystyle {\sum}_{i=1}^{n}{E}_{0}[}{\ell}_{i}^{2}(t)|{Y}_{i}]$. 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 $Z(t)=\ell (t)/\sqrt{{I}_{\alpha \alpha}(t)}$, 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 trans-hubs
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 trans-hubs, 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 variance-component model
Y_{ i }= μ + Y_{-i}β + α_{ m }+ α_{ f }+ e, (2)
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 over-fitting. 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 trans-hub 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 inter-marker distance is at least 0.1 cM to avoid linkage disequilibrium. The resulting map has an average inter-marker 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.
Number of expressions with at least one significant eQTL
Genome-wide significance level (point-wise significance level) | |||
---|---|---|---|
0.05 (4.7 × 10^{-6}) | 0.01 (3.2 × 10^{-7}) | 0.001 (6.7 × 10^{-9}) | |
Threshold by Gaussian approximation | 4.43 | 4.98 | 5.68 |
Number by chance (95% CI) | 177.7 (162.2, 193.2) | 35.6 (28.5, 42.6) | 3.6 (1.3, 5.8) |
Original data | |||
Threshold by permutation | 4.36 | 4.95 | 5.68 |
Number by permutation (95% CI) | 176.0 (139.0, 224.5) | 36.0 (24.0, 56.0) | 4.0 (0, 9.0) |
Observed number | 235 | 66 | 24 |
Sibship only | 173 | 40 | 13 |
Residual data (9p13.3) | |||
Threshold by permutation | 4.45 | 5.06 | 5.81 |
Number by permutation (95% CI) | 180.0 (154.0, 207.0) | 36.0 (25.0, 48.0) | 4.0 (0, 8.0) |
Observed number | 225 | 87 | 26 |
Residual data (14q32) | |||
Threshold by permutation | 4.46 | 5.07 | 5.85 |
Number by permutation (95% CI) | 178.0 (151.5, 205.0) | 36.0 (24.0, 47.0) | 3 (0, 8.0) |
Observed number | 212 | 70 | 23 |
Results
Expression phenotypes with the strongest evidence of linkage from genome scans
Point-wise p-value | Gene | Gene location | cis/trans | eQTL location |
---|---|---|---|---|
<10^{-15} | LRAP | 5q15 | cis | Chr 5 (99080578) |
<10^{-15} | HLA-DQB1 | 6p21.3 | trans ^{a} | Chr 6 (37592767) |
<10^{-15} | CHI3L2 | 1p13.3 | cis | Chr 1 (111704864) |
<10^{-15} | POMZP3 | 7q11.23 | cis | Chr 7 (75651464) |
<10^{-14} | CSTB | 21q22.3 | cis | Chr 21 (44061921) |
<10^{-14} | TBC1D8 | 2q11.2 | trans ^{a} | Chr 2 (108214542) |
<10^{-13} | DSCR2 | 21q22.3 | trans | Chr 9 (75300235) |
<10^{-13} | CRYZ | 1p31-p22 | trans ^{a} | Chr 1 (67949299) |
<10^{-11} | EGR2 | 10q21.1 | trans | Chr 20 (42643248) |
<10^{-11} | TM7SF3 | 12q11-q12 | trans ^{a} | Chr 12 (39239200) |
<10^{-11} | DDX17 | 22q13.1 | cis | Chr 22 (39410468) |
Number of cis-hits and trans-hits^{a}
Genome-wide significance level (point-wise significance level) | ||||
---|---|---|---|---|
0.05 (4.7 × 10^{-6}) | 0.01 (3.2 × 10^{-7}) | 0.001 (6.7 × 10^{-9}) | ||
cis-hits | Original (percentage) | 108 (0.72) | 74 (0.49) | 46 (0.30) |
Residual^{b} (percentage) | 182 (1.20) | 128 (0.85) | 68 (0.45) | |
trans-hits | Original (percentage) | 1166 (0.028) | 296 (0.0070) | 96 (0.0023) |
Residual^{b} (percentage) | 1057 (0.025) | 291 (0.0069) | 63 (0.0015) |
To further investigate our hypothesis, we adjusted for the expression correlations as described in the Methods on the trans-hub 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 genome-wide 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, HLA-DPB1, CSTB, BCKDHA, DSCR2, POMZP3, CHI3L2, HSD17B12, and TM7SF3). However, the hub phenomenon becomes much less obvious now: trans-eQTL hits were very evenly distributed along the genome (Figure 1, right panel). One explanation is that the overall pair-wise 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 trans-hits of one 5-Mb region at the genome-wide 0.05 significance level is only 4. The number of trans-hits at 9p13.3 drops to 3, while among the 500 permutation cycles performed on residual data, none has a maximum number of trans-hits smaller than 3. The same analysis is done for the trans-hub at 14q32 and the results are similar (Table 1). The number of trans-hits at 14q32 drops to 4.
Thus, we conclude that there is not enough evidence to claim either two candidate regions we examined as a trans-hub. However, we do find two statistically significant, trans-regulated, phenotypes: DSCR2 (Down Syndrome Critical Region gene 2) is the most significant (point-wise p-value 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 point-wise p-value of 3 × 10^{-7}. The correlation between expressions of DSCR2 and MAP3K6 is quite small: 0.074. These two trans-regulations 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 large-scale permutation studies. Using the original phenotypes, we identify two candidate trans-hubs, one at 9p13.3 and the other at 14q32. However, after accounting for the expression correlations in the linkage analysis, both trans-hubs 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 trans-regulation, 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.
Notes
Declarations
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/1753-6561/1?issue=S1.
Authors’ Affiliations
References
- Koning D, Haley C: Genetical genomics in humans and model organisms. Trends Genet. 2005, 21: 377-381. 10.1016/j.tig.2005.05.004.View ArticlePubMedGoogle Scholar
- 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: 422-425. 10.1038/ng1094.View ArticlePubMedGoogle Scholar
- 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: 1365-1369. 10.1038/nature04244.View ArticlePubMed CentralPubMedGoogle Scholar
- Morley M, Molony C, Weber T, Devlin J, Ewens K, Spielman R, Cheung V: Genetic analysis of genome-wide variation in human gene expression. Nature. 2004, 430: 743-747. 10.1038/nature02797.View ArticlePubMed CentralPubMedGoogle Scholar
- Almasy L, Blangero J: Multipoint quantitative-trait linkage analysis in general pedigrees. Am J Hum Genet. 1998, 62: 1198-1211. 10.1086/301844.View ArticlePubMed CentralPubMedGoogle Scholar
- Zou H, Hastie T: Regularization and variable selection via the elastic net. J Roy Stat Soc B. 2005, 67: 301-320. 10.1111/j.1467-9868.2005.00503.x.View ArticleGoogle Scholar
- Mallows L: Some comments on C_{ P }. Technometrics. 1973, 15: 661-676. 10.2307/1267380.Google Scholar
- Wang K, Huang J: A score-statistic approach for the mapping of quantitative-trait loci with sibships of arbitrary size. Am J Hum Genet. 2002, 70: 412-424. 10.1086/338659.View ArticlePubMed CentralPubMedGoogle Scholar
- Peng J, Siegmund D: Mapping quantitative traits under ascertainment. Ann Hum Genet. 2006, 70: 867-881. 10.1111/j.1469-1809.2006.00286.x.View ArticlePubMedGoogle Scholar
- Tang HK, Siegmund D: Mapping quantitative trait loci in oligogenic models. Biostatistic. 2001, 2: 147-162. 10.1093/biostatistics/2.2.147.View ArticleGoogle Scholar
- Abecasis G, Cherny S, Cookson W, Cardon L: Merlin-rapid analysis of dense genetic maps using sparse gene flow trees. Nat Genet. 2002, 30: 97-101. 10.1038/ng786.View ArticlePubMedGoogle Scholar
- 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-View ArticlePubMed CentralPubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.