Definition of Cis-eqtl

Definition of Cis-eqtl

We are searching data for your request:

Forums and discussions:
Manuals and reference books:
Data from registers:
Wait the end of the search in all databases.
Upon completion, a link will appear to access the found materials.

I have simple question, do cis-eQTLs use the SNPs that are within 1 Mb from the gene's TSS on the same chromosome ?

Kind of.

What you describe is a cis-eQTL, but the exact distance (1Mb in each direction, 1Mb total,… ) and the definition of the location of the gene (TSS, transcript borders,… ) depend on whoever calculates a cis-eQTLs and are not universal. The 'cis' part just means the SNP is somehow near the gene, there is no defined limit on the distance.

EQTL Analysis

Genetic variation in a population is commonly studied through the analysis of single nucleotide polymorphisms (SNPs), which are genetic variants occurring at specific sites in the genome. Expression quantitative trait loci (eQTL) analysis seeks to identify genetic variants that affect the expression of one or more genes: a gene-SNP pair for which the expression of the gene is associated with the allelic configuration of the SNP is referred to as an eQTL. Identification of eQTLs has proven to be a powerful tool in the study and understanding of diseases in human and other populations.

Using modern genotype and expression arrays, a typical eQTL analysis can involve millions of SNPs and tens of thousands of genes, making computation and multiple testing key challenges. Even local (cis) eQTL analyses that restrict attention to nearby genes and SNPs can involve tens of millions of gene-SNP pairs. Our initial work on eQTL analysis addressed fast computation of association statistics in homozygous populations, and subsequent testing. We are currently investigating the use of iterative testing methods to enhance the power of full (trans) eQTL analyses. Complementing eQTL testing, we have recently developed a simple log-of-linear model for assessing the effect size of an eQTL, an important problem that has not received much systematic attention in the literature.

To date, most eQTL studies have considered the effects of genetic variation on expression within a single tissue (typically blood). An important next step is the simultaneous analysis of eQTLs in multiple tissues. Multi-tissue analysis has the potential to improve the findings of single tissue eQTL studies by borrowing strength across tissues, and to address fundamental biological questions about the nature and source of the differences between tissues. An important feature of multiple tissue studies is that a SNP may be associated with the expression of a gene in some tissues, but not in others. Working with the NIH Genotype-Tissue Expression (GTEx) Consortium, we developed an empirical Bayes procedure, called MT-eQTL, for multi-tissue eQTL analysis. The procedure, which is able to test for complex patterns of association across multiple tissues, was one of two methods used for testing eQTLs in the Consortium’s recent Science paper. The MT-eQTL procedure is limited to nine or ten tissues, but we are currently working on extensions that will scale to as many as twenty or thirty tissues.


Strains used

The wild-types N2 and CB4856 and 54 RILs derived from a CB4856 x N2 cross were used (strains generated in [12]). For 49/54 of these strains low-coverage sequencing was applied to construct a more detailed genetic map (see also [22]). A matrix with the strain names and the genetic map can be found in Additional file 1.

Nematode culturing

The strains were kept on 6-cm Nematode Growth Medium (NGM) dishes containing Escherichia coli strain OP50 as food source [23]. Strains were kept in maintenance culture at 12 °C, the standard growing temperature for experiments was 20 °C. Fungal and bacterial infections were cleared by bleaching [23]. The strains were cleared of males prior to the experiments by selecting L2 larvae and placing them individually in a well in a 12-wells plate at 20 °C. Thereafter, the populations were screened for male offspring after 3 days and only the 100% hermaphrodite populations were transferred to fresh 9-cm NGM dishes containing E. coli OP50 and grown until starved.

Control, heat stress, and recovery from heat stress experiments for transcriptomics

The experiments were started by transferring a starved population to a fresh 9-cm NGM dish. This population was grown for 60 h at 20 °C to obtain egg-laying adults, which were bleached in order to synchronize the population. The eggs were transferred to a fresh 9-cm NGM dish. Three growing conditions were applied: (i) the control treatment was grown for 48 h at 20 °C, (ii) the heat-stress treatment was grown for 46 h at 20 °C followed by 2 h at 35 °C, and (iii) the recovery treatment was grown for 46 h at 20 °C, followed by 2 h at 35 °C and thereafter 2 h at 20 °C. Before the start of the treatment, the developmental stage of the population was determined by observing the developmental stage of the vulva in multiple individuals. Populations not consisting of L4 larvae were not isolated. Directly at the end of the treatment, the population was washed off the plate with M9 buffer and collected in an Eppendorf tube, which was flash frozen in liquid nitrogen. In this manner, 48 RILs per condition were assayed.

Genotypes and genetic map construction

Previously, 49 lines were sequenced and aligned. The single-nucleotide polymorphism (SNP) calls per strain were taken for constructing the genetic map [22]. The SNP density was determined per 10 kb bins and recombination events were recognized as transition of an area where there were no CB4856 SNPs in 10 consecutive bins into an area where there were CB4856 SNPs and the other way around. It was not allowed to have two recombination events within 10 consecutive bins (100 kb). The 10 kb bin where the first SNPs were detected was marked as the recombination event. Before use in mapping, the map was filtered for informative markers – that is - markers indicating a recombination event in at least one of the lines. This resulted in a map of 729 informative markers, each indicating the location of the recombination events within 10 kb (see the figure in Additional file 2).

The genetic map was investigated by correlation analysis to assess the linkage between markers. Markers on the centers of the chromosomes showed strong linkage (see also [24]). No strong in between chromosome correlations were found (see the figure in Additional file 3).

Transcript profiling

RNA isolation

The RNA of the RIL samples was isolated using the RNeasy Micro Kit from Qiagen (Hilden, Germany). The ‘Purification of Total RNA from Animal and Human Tissues’ protocol was followed, with a modified lysing procedure frozen pellets were lysed in 150 μl Rneasy Lysis Tissue buffer, 295 μl RNAse-free water, 800 μg/ml proteinase K and 1% ß-mercaptoethanol. The suspension was incubated at 55 °C at 1000 rpm in a Thermomixer (Eppendorf, Hamburg, Germany) for 30 min or until the sample was clear. After this step the manufacturer’s protocol was followed.

CDNA synthesis, labelling and hybridization

The ‘Two-Color Microarray-Based Gene Expression Analysis Low Input Quick Amp Labeling’ -protocol, version 6.0 from Agilent (Agilent Technologies, Santa Clara, CA, USA) was followed, starting from step five. The C. elegans (V2) Gene Expression Microarray 4X44K slides, manufactured by Agilent were used. Before starting cDNA synthesis, quality and quantity of the RNA were measured using the NanoDrop-1000 spectrophotometer (Thermo Scientific, Wilmington DE, USA) and RNA integrity was determined by agarose gel electrophoresis (3 μL of sample RNA on 1% agarose gel).

Data extraction and normalization

The microarrays were scanned by an Agilent High Resolution C Scanner with the recommended settings. The data was extracted with Agilent Feature Extraction Software (version 10.5), following manufacturers’ guidelines. For normalization, “R” (version 3.3.1 × 64) with the Limma package was used. The data was not background corrected before normalization (as recommended by [25]). Within-array normalization was done with the Loess method and between-array normalization was done with the Quantile method [26]. The obtained single channel normalized intensities were log2 transformed and used for further analysis.

Environmental responses

The transcriptional response to heat stress was determined by explaining the gene expression over the treatment with a linear model,

where y is the log2-normalized intensity as measured by microarray of spot i (i = 1, 2, . 45,220), and T is the treatment (either control, heat stress, or recovery from heat stress). This analysis ignored genotype.

The significances were corrected for multiple testing by applying the Benjamini Yekutieli method in p.adjust (R, version 3.3.1 Windows ×64) at FDR = 0.05 [27]. Thresholds of -log10(p) ≥ 2.87 for the control versus heat-stress treatment, −log10(p) ≥ 3.09 for the control versus recovery treatment, and -log10(p) ≥ 3.02 for the heat-stress versus recovery treatment were determined.

Developmental variation

Due to the setup of our experiment, potential variation in development could exist among the RILs and the treatments. The recovery animals were sampled two hours later than the control and heat-stress animals, furthermore, heat stress slows the developmental rate [19]. We estimated the relative age by using a set of

100 genes that show a strong, positive, linear response during development [9]. By setting the average age of the control RILs to 48 h we could estimate and compare the RILs in all treatments (Additional file 8).

Principal component analysis

A principal component analysis was conducted on the gene-expression data of the RILs over the three treatments. For this purpose, the data was transformed to a log2 ratio with the mean, using

where R is the log2 relative expression of spot i (i = 1, 2, . 45,220) in strain j (RIL) over all three conditions (n = 48 per condition), and y is the intensity (not the log2-transformed intensity) of spot i in strain j.

The transformed data was used in a principal component analysis, where the first six axes were further examined.

Expression quantitative trait locus analysis

EQTL mapping and threshold determination

The eQTL mapping was done in “R” (version 3.3.1 Windows ×64). The gene-expression data was fitted to the linear model,

where y is the log2-normalized intensity as measured by microarray of spot i (i = 1, 2, . 45,220) of RIL j. This is explained over the genotype (either CB4856 or N2) on marker location x (x = 1, 2, . 729) of RIL j.

The genome-wide significance threshold was determined via permutation, where the log2-normalized intensities were randomly distributed per gene over the genotypes. The randomized data was tested using the same model as for the eQTL mapping. This was repeated for ten randomized datasets. A false discovery rate was used to determine the threshold (as recommended for multiple testing under dependency) [27],

where FDS (false discoveries) is the outcome of the permutations and RDS (real discoveries) is the outcome of the eQTL mapping at a specific significance level. The value of m0, the number of true null hypotheses tested, was 45,220-RDS, and for the value of m, the number of hypotheses tested, the number of spots (45220) was taken. The q-value was set at 0.05. This yielded a threshold of –log10 (p) > 3.9 for the control, −log10 (p) > 3.5 for the heat stress, and–log10 (p) > 3.9 for the recovery treatment. For the analyses we used the most conservative thresholds measured, −log10 (p) > 3.9, for all the sets.

Statistical power calculations

In order to determine the statistical power at the set FDR threshold, QTL were simulated using the genetic map of the strains used per condition (n = 48 per condition). For each marker location, ten QTL were simulated that explained 20–80% of the variation (in increments of 5%). Random variation was introduced based on a normal distribution with sigma = 1 and mu = 0 and a peak of the corresponding size (e.g. a peak size of 1 corresponds to 20% explained variation) was simulated in this random variation. From the simulation, the number of correctly detected QTL, the number of false positives and the number of undetected QTL were counted. This was based on the thresholds determined in the permutations, −log10 (p) > 3.9. Furthermore, the precision of the effect-size estimation and the precision of the QTL location (based on a –log10(p) drop of 1.5 compared to the peak) were determined. A table summarizing the results can be found in Additional file 9.

EQTL analysis

The distinction between cis- and trans-eQTL was made on the distance between the physical location of the gene and the location of the eQTL-peak. For cis-eQTL the gene lies within 1 Mb of the peak or within the confidence interval of the eQTL. The confidence interval was based on a –log10(p) drop of 1.5 compared to the peak.

The amount of variation explained per microarray spot with an eQTL was calculated by ANOVA, by analysis of the gene expression explained over the peak-marker. For spots with multiple peaks, this analysis was conducted per peak, not using a full model, since a single-marker model was used in the analysis.

In order to identify trans-bands (an enrichment of trans-eQTL), a Poisson distribution of the mapped trans-eQTL was assumed (as in [28]). Therefore the number of trans-eQTL per 0.5 Mb bin were counted. Since trans-eQTL peaks were mapped to 107, 106, and 103 bins (respectively in control, heat stress, and recovery), it was expected that 9.16, 20.64, and 9.01 spots with a trans-eQTL were to be found at each of these markers. Based on a Poisson distribution, it was calculated how many trans-eQTL needed to be found to represent an overrepresentation. For example, for p < 0.001 there should be 20, 36, or 20 spots with a trans-eQTL at a specific marker (respectively in control, heat stress, and recovery).

To test for polymorphisms in genes with eQTL, we used the data from the CB4856 reference genome [22]. The genes with eQTL were matched to the polymorphisms. The frequencies of polymorphisms in each of the groups (genes with cis-eQTL, genes with trans-eQTL, and genes without eQTL) were counted and compared versus each other by a chi-squared test in “R” (version 3.3.1, ×64).

Detection of eQTL across treatments

Two criteria were used to detect the occurrence of eQTL over multiple treatments.

In the first criterion, it was tested whether or not an eQTL was mapped in treatment one versus treatment two, by simply comparing the tables listing the eQTL. This allowed for comparison of the actual mapped peaks and for comparison of eQTL effects of trans-eQTL regulated from different loci. In order to estimate the false-discovery rate associated with this comparison, the same analysis was applied to ten permutated datasets per condition, using the –log10(p) > 3.9 for eQTL discovery.

The second criterion compared the occurrence of eQTL at the exact same marker location. In this comparison, the eQTL mapped in one treatment were taken as lead for the occurrence of the same eQTL in the other two treatments. This comparison allowed for direct comparison of the eQTL effect at the locus. Based on observations on the effect distribution, this approach was used to estimate the number of trans-eQTL not detected due to statistical power or not detected due to absence of the eQTL in a treatment (see also text in Additional file 15).

Functional enrichment analysis

Gene group enrichment analysis was done using a hypergeometric test and several databases with annotations. The databases used were: the WS220 gene class annotations, the WS256 GO-annotation, anatomy terms, phenotypes, RNAi phenotypes, developmental stage expression, and disease related genes ( [29] the MODENCODE release 32 transcription factor binding sites ( [30, 31], which were mapped to transcription start sites (according to [32]) and the KEGG pathway release 65.0 (Kyoto Encyclopedia of Genes and Genomes, [33].

Enrichments were selected based on the following criteria: size of the category n > 3, size of the overlap n > 2. The overlap was tested using a hypergeometric test, of which the p-values were corrected for multiple testing using Bonferroni correction (as provided by p.adjust in R, 3.3.1, ×64). Enrichments were calculated based on gene names, not on spots.


A sequence variant that affects the expression of a gene on the same chromosome (a cis eQTL) should be detectable both as a local eQTL and by ASE. However, both local eQTL and ASE could be due to other causes besides a cis eQTL, for example local eQTL could be trans eQTL [5]. ASE could be due to parent of origin effects (i.e. imprinting) or nonsense mediated decay [5, 11]. Our statistical analysis attempts to distinguish cis eQTL from parent of origin as a cause of ASE by determining whether the direction of ASE is consistent with the parent of origin of the over-expressed allele or with the allele on the same chromosome at a nearby SNP (a dSNP). This analysis finds that PO-ASE is less common and has higher FDR than ASE due to the cis effect of a nearby SNP. The PO-ASE effects were also less often confirmed in multiple datasets. However, our analysis did find PO-ASE for several genes where imprinting has been reported. To avoid being misled by errors in the genome sequence of our cattle, we excluded complete mono allelic expression, and therefore would not find cases of complete imprinting where only one parental allele is expressed. Therefore, our PO-ASE tests are detecting “partial” imprinting. There are also some reports that imprinting can be tissue specific [13], which could explain why PO-ASE effects were not always confirmed in multiple datasets. However, the agreement of PO-ASE found in Holstein and Angus liver samples was also low perhaps because we had limited power to detect partial imprinting with only a small allelic imbalance.

Our analysis of ASE fits the effects of parent of origin and of a SNP in cis with the over-expressed allele jointly and therefore we can test for significance the cis effect of SNPs on gene expression, called here simply ASE. There are large numbers of SNPs associated with the expression of a nearby gene (local eQTL), or with ASE at that gene (Table 2). The SNPs associated with gene expression overlap significantly between ASE and local eQTL analyses and between datasets on different animals and different tissues and, among the SNPs that are significant in two different analyses, the same allele is associated with increasing gene expression most of the time (Table 4). These results support the conclusion that many of the SNPs associated with local eQTL and with ASE are cis eQTL.

Despite the significant overlap between local eQTL and ASE, only a small proportion of SNPs that are significant in one analysis are significant in the other. This could be due to lack of power in one or both analyses or systematic differences between ASE and local eQTL. When ASE in two different datasets is compared, there are on average 14.33 times more SNPs that are significant in both datasets than expected by chance (Table 4). Similarly, when local eQTL are compared in 2 datasets there are 12.87 times more significant SNPs than expected by chance (Table 4). However, when ASE in one dataset is compared to local eQTL in another, the fold enrichment is only 5.72 (Table 4). This indicates that ASE detects phenomenon that overlap with local eQTL, but have some systematic differences. One reason for these differences is that the local eQTL analysis uses all RNA reads from the gene whereas the ASE analysis uses only those containing the tSNP. When ASE in different datasets is compared using any tSNP within the gene, instead of the same tSNP in both datasets, the fold enrichment falls to 4.81 (Additional file 2: Table S4). Thus, some of the cis eQTL detected by ASE are tSNP specific probably because they are exon and splice variant specific [14]. Other reasons for systematic differences between ASE and local eQTL may include local eQTL being trans eQTL, nonsense mediated decay, and feedback mechanisms limiting gene expression. If one allele is more highly expressed than the other, feedback might limit expression of both alleles so that the total gene expression is much the same regardless of SNP genotype [8]. This would leave a significant case of ASE but no local eQTL.

Although there is some overlap between eQTL in different datasets, there are also systematic differences between datasets, breeds and tissues. The comparisons between ASE and local eQTL mapping showed within an RNA-Seq dataset (Additional file 2: Table S7), the SNPs found by ASE were more likely to be found in local eQTL mapping than expected by chance. The average fold enrichment for ASE and local eQTL in the same dataset was more than when comparing ASE in one dataset with local eQTL in another dataset, indicating that although some eQTL operate in both datasets, there are some dataset specific (tissue specific) eQTL. These results could be interpreted to mean that approximately half the cis eQTL in one tissue also operate in a second tissue which is the same conclusion as the GTex paper reached by different means [1]. Our results also show that when a cis eQTL operates in multiple tissues, it is nearly always the same allele that increases expression.

The systematic differences between datasets are due partly to differences in the tissue used, as shown by the similarities in local eQTL mapping and ASE between Holstein and Angus liver samples (Table 4). The ASE and local eQTL results were also more similar when comparing two Holstein or two Angus datasets than when comparing one Angus and one Holstein dataset (Table 4). The effect of breed might be due to differences between breeds in linkage disequilibrium, but could also be due to differences in gender, environment or physiological state between the Holsteins and Angus cattle we sampled.

These conclusions of the effect of tissue and breed are supported by comparing our results with ASE in 18 tissues of a Holstein cow. In these results, the percentage of SNPs showing the same direction of effect is higher when the tissues are the same and when the breed is Holstein in both cases. So, it seems that ASE and local eQTL are partially tissue specific which agrees with the findings of Chamberlain et al. (2015) [11]. However, there was overlap between all of our 4 datasets and many of the 18 tissues suggesting that some cis eQTL affect expression in many tissues. This is consistent with the GTex paper [1] that found the most cis eQTL either affected all 9 tissues or only one tissue.

There were significant (p < 0.05) overlaps between QTL influencing some complex traits and ASE or local eQTL (Table 5). Tenderness of the meat (MQLDPF) was associated with many SNPs that were also significantly associated with gene expression in multiple analyses and datasets. Fig 2 gives an example of a gene containing SNPs that significantly affect tenderness and expression of the gene calpastatin (CAST). The physiological role of CAST and the association between SNPs near CAST and tenderness are well known [15,16,17] but these results suggest that the known QTL is in fact an eQTL for calpastatin expression. Further examples include the calpain 1 gene (CAPN1), leptin receptor overlapping transcript-like 1 gene (LEPROTL1) and ligand dependent nuclear receptor corepressor-like (LCORL) (Fig. 3). CAPN1 influences meat tenderness [16]. LEPROTL1 is reported to influence body growth by negatively regulating leptin receptor (LEPR) cell surface expression, decreasing response to leptin and decreasing hepatic growth hormone action in mice [18]. LCORL has been reported to have an effect on feed intake, gain, meat and carcass traits [19] and its expression in muscle tissue has been associated with average daily feed intake in beef cattle [20]. In humans, there is a QTL for height near NCAPG and LCORL and it has not been possible to identify the causal gene. The result here suggests that LCORL, at least, is likely to affect growth.

There are many QTL for which we did not find a matching eQTL. This could be because the QTL acts through a different mechanism (e.g. a protein coding mutation) or due to the lack of power to find the relevant eQTL in our experiment. In our study gene expression was measured in just 3 tissues (muscle, liver and WBC) and only once. The eQTL in other tissues and their activity in other physiological or developmental states might underlie other QTL. The number of animals and the average transcript read depth in our study were also limiting factors in the power to detect eQTL. Therefore, if we want to find QTL influencing complex traits by controlling the expression of genes, it is important to sample the correct tissue at the appropriate time and to have sufficient sample size and transcript read depth for gene expression analyses.


Our genome-wide studies of expression levels of genes underlying cis- and trans-eQTLs provide strong support for the hypotheses that the categorisation of eQTLs has a genuine biological basis that can be detected in transcript expression levels, and that trans-eQTL clusters consist of functionally related and co-ordinately regulated transcripts.

We observed that patterns of correlation of expression profiles and of SDPs at the genetic location of the transcript whose expression was measured are strikingly different between the sets of cis- and trans-eQTLs in all four tissues ( Figures 1 – ​ 3). 3 ). Significant correlation between pairs of cis-eQTLs was found to be rare (ρ%, Table 2 ). A significant relationship has previously been shown between the absolute correlation coefficient of pairs of eQTLs derived from genes located on the same chromosome and the distance apart of those genes [10], [24]. Hence it was hypothesised that correlation of cis-eQTL genes can be explained in terms of linkage disequilibrium. We found that most of the observed correlations between cis-eQTL genes' expression profiles for co-localised genes can be explained by similarity of underlying genotypes. We also found that most significant correlations between cis-eQTLs located on different chromosomes can be explained by long-range allelic association. Such patterns of association have previously been observed in a SNP map of a mouse inbred strain panel [26], and hypothesised to have the potential to produce spurious associations between transcripts' expression profiles. Our findings suggest that correlation of expression profiles, in isolation, is not a suitable basis for the analysis of relationships between cis-eQTLs.

Much higher levels of significant correlation, encompassing 2.9�.9% of pairs, were observed between trans-eQTL genes. This was not found to be explained by correlation of the genotypes at the regions to which the transcripts map. However, when genotypes at the trans-eQTL are similar or the same, the gene pairs were disproportionately correlated (47.5 to 63.7% of significantly correlated pairs). This observation provides strong support for the hypothesis, formulated in studies of genetic regulation of gene expression in S. cerevisiae [19], that trans-eQTL clusters are indicative of co-regulation of the transcripts. In our dataset, the level of significant correlation within trans-eQTL clusters was found to be far higher than across the trans-eQTL dataset as a whole, averaging 83.5% across the 81 clusters ( Table 4 ). Additionally, we observed that much of the variability in the levels of trans-eQTL correlation between tissues can be explained by differences in the distribution and size of the trans-eQTL clusters (data not shown).

Functional investigation of the large trans-eQTL clusters showed significant over-representation of Gene Ontology (GO) terms in 80% of these clusters. While GO analysis may have limitations in scope and accuracy, owing to the annotation challenges posed by genes which often have complex and imprecisely defined roles and functions [29], the findings presented here are in agreement with those of Ghazalpour et al [25], who observed that genes in functionally related ‘pathway sets’ are typically highly correlated. KEGG pathway analysis of the trans-eQTL clusters was carried out through the DAVID interface, and pathways of potential interest were identified in a minority of the clusters (Table S6). It has previously been shown in a study of co-regulation in a mouse F2 intercross [20] that groups of highly correlated transcripts linked to the same genomic location can be identified in a functionally informed genome-wide correlation analysis. Here we show across multiple tissues that highly correlated expression profiles are a consistent motif of such trans-eQTL clusters. These results suggest a significant degree of functional relatedness of the genes making up the cluster ( Table 5 ).

Expression correlation between trans-eQTL cluster transcripts and cis-eQTL genes located in the linkage region has been previously described as a method of identifying candidate regulators and applied in yeast [30] and Arabidopsis [7]. Here we show a strong relationship between the distance of the cis-eQTL from the cluster peak of linkage and the strength of correlation ( Figure 5 ), which was not observed when unlinked transcripts located in the region were similarly tested (Figure S1). This, along with the observation that there is no association between distance from the peak of linkage and heritability of the cis-eQTL genes (Figure S2), suggests that genotype similarity at the linkage region, rather than the genetic influence on gene expression underlies the relationship. On this basis, correlation-based methods of prioritising candidates for genetic regulation of trans-eQTL clusters would be improved by taking into account the distance between the cluster linkage region and the map location of the candidate transcript. We therefore suggest that it may be of interest instead to identify outliers those cis-eQTL genes whose average correlation with the trans-eQTL cluster genes significantly deviates positively from the general negative regression trend. We were able to find 54 cis-eQTLs with significant positive Z-scores (ZϢ) (Table S4), and consider that these may be worthy of further investigation into the possibility that their correlation with the cluster transcripts has a biological basis.

In this study, we demonstrate the power of computational analysis of eQTL datasets across multiple tissues to provide new insights into the genome-wide correlation structure in gene expression data. Our findings consistently show that correlation of cis-eQTL genes' expression profiles is primarily indicative of similarity of genetic inheritance, as measured by the correlation of SDPs at the transcript location. Whereas, correlation between trans-eQTLs can frequently be explained in terms of co-regulation by a common linkage region even though correlated transcripts forming trans-eQTL clusters are located throughout the genome. The observation of functional enrichment within clusters is suggestive of a relationship between co-expression and function. Finally, we inform investigations of candidate regulators of trans-eQTL clusters by indicating that genetic linkage strongly influences co-expression of trans-eQTL cluster genes and candidate regulatory genes.


The precise spatial and temporal control of gene transcription is critical for biological processes, as evidenced by the causal role of gene expression perturbation in many human diseases [1]–[3]. Gene expression is controlled by regulatory proteins, RNAs, and the cell type specific cis-regulatory elements with which they interact. Genetic variation within cis-regulatory elements (CREs, e.g., transcription promoters, enhancers, or insulators) can affect gene expression in a cell type specific manner. An extensive body of work, performed in a variety of cell types in both humans and model organisms, has demonstrated that genetic variants that impact gene expression, or expression quantitative trait loci (eQTLs), are common and exist in both cis (local) and trans (over long genetic distances) [3]–[6]. Over of genotype-phenotype associations found in genome-wide association studies (GWAS) are with non-coding single nucleotide polymorphisms (SNPs), making their mechanistic interpretation more challenging. It is possible that these associated SNPs tag causal coding SNPs however, numerous compelling lines of evidence [2], [7]–[11] demonstrate that regulatory SNPs have causal roles in many complex human phenotypes. This is further supported by the finding that GWAS associations are enriched within DNase I hypersensitive (DHS) sites [12] and eQTL SNPs [13], [14], and by several elegant GWAS follow up studies that have mechanistically tied disease associations with SNPs that cause the misregulation of gene expression [15], [16].

Although eQTLs are increasingly used to provide mechanistic interpretations for human disease associations, the cell type specificity of eQTLs presents a problem. Because the cell type from which a given physiological phenotype arises may not be known, and because eQTL data exist for a limited number of cell types, it is critical to quantify and understand the mechanisms generating cell type specific eQTLs. For example, if a GWAS identifies a set of SNPs associated with risk of type II diabetes, the researcher must choose a target cell type to develop a mechanistic model of the molecular phenotype that causes the gross physiological change. One can imagine that the relevant cell type might be adipose tissue, liver, pancreas, or another hormone-regulating tissue. Furthermore, if the GWAS SNP produces a molecular phenotype (i.e., is an eQTL) in lymphoblastoid cell lines (LCLs), it is not necessarily the case that the SNP will generate a similar molecular phenotype in the cell type of interest. Furthermore, there are many examples of cell types with particular relevance to common diseases, for example dopaminergic neurons and Parkinson's disease, that lack comprehensive eQTL data or catalogs of CREs. The utility of eQTLs for complex trait interpretation will therefore be improved by a more thorough annotation of their cell type specificity.

While several studies have quantified the reproducibility of eQTLs within or between cell types derived from the same or different individuals [17]–[28] the determinants of eQTL cell specificity are still largely unknown. We address this need in this study by analyzing cell specific eQTLs collected from eleven studies performed in seven different cell types and by integrating these data with cell specific CRE data to mechanistically interpret cell specific eQTLs. We used Bayesian regression models to identify all cis-linked SNPs that are independently associated with each gene expression trait in each study. In an effort to identify the functional determinants of eQTL cell specificity, we quantified the associations between eQTL SNPs and CRE data sets, many of which were derived from the cell types used in eQTL discovery and are known to function in a cell type specific manner (e.g., transcription factor binding sites (TFBSs), DHS sites). We further considered the relationship between eQTL SNP-CRE overlap and the cell type specificity of eQTL replication. Lastly, we built a series of classifiers to predict the cell type specificity of eQTLs in the absence of additional gene expression data and to predict the function of GWAS SNPs with phenotypes relevant to cell types lacking eQTL data. We believe these predictive models will facilitate more substantial mechanistic analyses of individual SNPs by enabling the integration of disease genetics and regulatory SNPs with the thousands of genomic data sets that have been produced by projects like ENCODE [29], [30].


We have carried out a systematic analysis of eQTL networks constructed from cis- and trans-eQTL in 13 different tissues, using data available from the GTEx project. We have found that the structural properties of these networks provide functional insight into the regulatory roles of genetic variants across and within tissues. Using a community detection algorithm (30) to search for communities of densely connected SNPs and genes, we found that the eQTL network in each of the 13 tissues was organized into highly modular communities. When we examine the genes represented in each community, we find an enrichment for genes, located on multiple chromosomes, that share similar functions or are associated with coherent biological processes. While the FDR may not be well controlled for the large number of GO terms tested, our resampling analysis for GO terms shared across tissues suggests that the observed enrichment is unlikely to be due to the large number of tests. However, the possibility of post hoc plausibility explanations cannot be completely ruled out. Contrary to what one might expect, these communities are not driven by coexpression (excepting communities with very few genes), suggesting that it is the genetic influence of multiple cis- and trans-eQTL SNPs on functionally related groups of genes that drives the organization and structure of these communities.

When comparing communities between tissues, we find many communities with common patterns of functional enrichment across tissues, reinforcing the pleiotropic role of the SNPs in these communities. We do, however, also find TS communities that contain genes involved in TS functional processes such as cellular respiration in heart left ventricle or smooth muscle contraction in esophagus muscularis. There is a plausible mechanistic explanation for the tissue specificity of some of these communities: Using data in eight tissues from the Roadmap Epigenome Project we find that TS eQTL SNPs in these TS communities are enriched for active chromatin regions that are unique to that tissue. This suggests that the organization of these communities is driven by the epigenetic activation of chromatin regions surrounding specific SNPs and that these SNPs act in cis and trans to exert genetic effects on the expression of functionally related genes, genes with important roles in their respective tissue-level processes. In addition, these communities are not only enriched for specific tissue-relevant gene function they are also enriched for tissue-specific edges (eQTL), SNPs, and genes. This is relevant to the ongoing discussion of the tissue specificity of eQTL. Although most eQTL appear to be shared, TS eQTL emerge in concert with TS epigenetic changes and not only influence single genes, but also help coalesce TS gene expression into regulatory communities.

We find these 13 eQTL networks possess two informative types of hubs: community hubs or “cores,” which are SNPs highly connected to genes in their community, and global hubs, which are connected to many genes throughout the network. These two types of hubs have different biological properties across tissues: Community hubs are enriched for active chromatin regions close to the transcriptional start site, but not enhancers, while global hubs are enriched for distal elements such as nongenic enhancers. Moreover, community hubs are enriched for GWAS-associated variants, while global hubs are not. The degree distribution for trait-associated variants from the GWAS is also highly consistent across the 13 tissues: GWAS SNPs are enriched for intermediate network degree, depleted for low degree, and absent from global hubs. The significant overrepresentation of GWAS SNPs among the community cores provides another important insight. Across tissues, disease SNPs are those most likely to perturb groups of genes and, in doing so, may disturb important biological processes.

While the observed relationships between eQTL network properties and SNP/gene function are consistent across tissue types, we cannot rule out the possibility that the large number of statistical tests performed in the cis- and trans-eQTL analysis could lead to identification of some individual eQTL associations as significant when they are not. Although we cannot conclude, based on our analysis, that any individual SNP–gene association is correct, the consistency of our findings regarding the structure of the networks across multiple tissues and the consistent functional enrichment we observe for global and local hubs indicate that the higher-order structural organization of the networks likely provides a robust model of SNP regulatory effects. While one could imagine that the observed network patterns might be driven by unwanted systematic variations in the genotype and gene expression data, our identification of similar structural properties in an eQTL network derived using an independent chronic obstructive pulmonary disease (COPD) dataset (23) further supports the network structural associations we have described. Nevertheless, this possibility, along with more detailed analysis of specific network-prioritized SNPs, should be further investigated as additional TS gene expression and genotype data become available through the next release of GTEx and other large-cohort studies.

Our analysis of bipartite networks built from both cis- and trans-eQTL in 13 tissues provides important evidence about the collective role of eQTL in TS function and disease. The network communities reveal biological processes under the shared genetic influence of many variants, including both processes shared across tissues and those that are TS. The TS genetic regulation we observe is driven in part by SNPs that lie in TS active chromatin regions. This suggests that epigenetic profile analysis, applied to both genic and nongenic elements, will be essential for understanding the processes responsible for TS function. The eQTL networks also group together functionally related sets of variants, including GWAS SNPs, and the structure of the network provides a model of how multiple cis- and trans-acting variants can work together to influence function and phenotype. While the network models we describe do not fully resolve the question of how weak-effect variants determine complex traits and disease, this network approach provides a framework with distinct explanatory power that can serve as a basis for further exploration of the link between genotype and phenotype.

Author contributions

FXM, SM, and JB designed the experiments. FXM generated the yeast library. SC and SM generated the sequencing libraries. SC performed the ChIP-seq analysis. MRL and FXM designed and performed the RT–qPCR experiments. MAP and CTW designed, performed, and analyzed the growth experiments. AS designed, carried out, and analyzed the proteomics experiments. MC-Z and AB designed the analyses MC-Z and SR carried out the analysis. The manuscript was written by all authors. CTW, JB, and AB envisioned and supervised the project.

Cis-acting expression quantitative trait loci in mice

We previously reported the analysis of genome-wide expression profiles and various diabetes-related traits in a segregating cross between inbred mouse strains C57BL/6J (B6) and DBA/2J (DBA). By considering transcript levels as quantitative traits, we identified several thousand expression quantitative trait loci (eQTL) with LOD score >4.3. We now experimentally address the problem of multiple comparisons by estimating the fraction of false-positive eQTL that are under cis-acting regulation. For this, we have utilized a classic cis-trans test with (B6 x DBA)F(1) mice to determine the relative levels of transcripts from the B6 and DBA alleles. The results suggest that at least 64% of cis-acting eQTL with LOD >4.3 are true positives, while the remaining 36% could not be confirmed as truly cis-acting. Moreover, we find that >96% of apparent cis-acting eQTL occur in regions that do not share SNP haplotypes. Cis-acting eQTL serve as an important new resource for the identification of positional candidates in QTL studies in mice. Also, we use the analysis of the correlation structures between genotypes, gene expression traits, and phenotypic traits to further characterize genes expressed in liver that are under cis-acting control, and highlight the advantages and disadvantages of integrating genetics and gene expression data in segregating populations.


Informative SNP frequency across chromosome…

Informative SNP frequency across chromosome 16 between B6 and DBA. The horizontal line…

Results of cistrans…

Results of cistrans test on a subset of cis -acting eQTLs.…

Significant differences in the distribution…

Significant differences in the distribution of Pearson correlation coefficients involving genes with strong…

Scatter plot of the mlratios…

Scatter plot of the mlratios for the jumping translocation point ( JTP )…

Diagram explaining the overall negative…

Diagram explaining the overall negative correlation observed in Figure 3.

Highlighting what may be the…

Highlighting what may be the biologically relevant component of the correlation between the…

Utilization of genetics of gene…

Utilization of genetics of gene expression data in order to prioritize candidate genes…

Genes with eQTL that are coincident with the physical location of the 1810073K19Rik…

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

YH participated in the study design, produced RNA, analyzed allele specific expression and wrote the manuscript. FH aligned the RNA-seq data, created allele specific counts and manually curated the alignments when needed. LM bred the mice and provided the adipose tissue. AVN performed the F2 cross local-eQTL mapping. EE participated in study design and analysis. AJL and TD conceived the study, participated in study design and coordination and wrote the manuscript. All authors read and approved the final manuscript.

Watch the video: MIT CompBio Lecture 15 - eQTLs (October 2022).