Genetic effects on gene expression across human tissues.

Sample size and eGene discovery in the GTEx v6p studya, Illustration of the 44 tissues and cell lines included in the GTEx v6p project with the associated number of cis- (left) and trans-eGenes (right) and sample sizes. Each tissue has a unique colour code (defined in Supplementary Fig. 5). b, Fraction of genes that are eGenes across all tissues by transcript class. The three tissues highlighted are: testis, which has the highest proportion of trans-eGenes; skeletal muscle, which has the largest sample size; and fibroblasts, which have the highest proportion of cis-eGenes. Dark bars depict the fraction of all curated human genes in GENCODE v19. Light bars depict the fraction of genes expressed in one or more tissues. c, Proportion of expressed genes that are cis-eGenes (y-axis) as a function of tissue sample size (x-axis). Colours represent tissues, as in a. d, Number of trans-eQTLs (x-axis) per tissue (y-axis), with sample size indicated by point size.

Identification of cis-acting eQTLs using allele-specific expression at chromosome-wide distancesa, A logistic regression based model was developed to predict the probability of phasing error as a function of distance and variant minor allele frequencies. When applied to chromosome 2 of 1000 Genomes sample NA12878, this model had a receiver operating characteristic (ROC) area under the curve (AUC) of 0.87 using population phasing compared to transmission phasing. b, ROC when applying the beta-binomial mixture model to detect cis-acting regulation to the GTEx v6p subcutaneous adipose cis-eQTLs, with an AUC of 0.88. As the null, eGenes were shuffled with respect to eVariants. c, Power analysis using all nominally significant (P <1.0 ×10−5) linkage disequilibrium pruned associations within 100 kb of the TSS illustrating the number of eQTLs with nominally significant (P ≤ 0.01) evidence of cis-regulation as a function of phasing error and eQTL effect size. Expression QTL effect size was calculated using a companion method28, and uniform phasing error between 0 and 100% was introduced in silico. d, Proportion of nominally significant (P <1.0 ×10−5) linkage disequilibrium-pruned intrachromosomal eQTLs with nominally significant (P ≤ 0.01) ASE supported evidence of cis regulation in bins of increasing TSS distance. Observed indicates what is seen in the data, while Max Error indicates what would be expected in the worst-case scenario of phasing error (50%). e, Example of significant ASE supported cis-regulation at a distance of 52.7 Mb between eVariant rs17494053 and eGene ENSG00000108509 in whole blood. Each point represents allelic imbalance in a single eVariant homozygote (circle) or heterozygote (triangle).

Effects of sample size and assayed tissues on cis-eGene discoverya, Number of significant cis-eGenes at 5% FDR (y-axis) discovered in 44 GTEx tissues versus sample size (x-axis). b, Number of significant eGenes at FDR 5% (y-axis) discovered in nine GTEx tissues each subsampled to various sizes where possible (n =70, 100, 125, 150, 175, 200, 225, 250, 275, 300, 325, and 350; x-axis). We computed the number of cis-eGenes at each subsample size (circles connected by lines). We also plotted the number of cis-eGenes discovered with no subsampling of donors (diamonds). c, Number of significant cis-eGenes at FDR 5% (y-axis) as a function of sample size and number of tissues assayed (x-axis). Each tissue was subsampled to 70, 100, 150, 200, and 250 donors, and a forward search was used to assess sequential combinations of tissues that maximize the total number of unique cis-eGenes discovered.

Measuring cis-regulatory variation using ASEa, Proportion of protein-coding genes with ASE data in at least one tissue as a function of donors (top row) and with significant imbalance (binomial test versus 0.5, 5% FDR) stratified by ASE effect size (|log2(hapa count/ hapb count)|) deciles. Gene-level measurements of haplotype expression were calculated by aggregating counts per sample across all heterozygous variants with ASE data within the gene using population phasing. The following filters were applied on ASE data: total coverage ≥8 reads, no mapping bias in simulations76, UCSC mappability > 50, and no significant (FDR > 1%) evidence that variant is monoallelic in expression data75. b, log2 transformed cis-eQTL effect size (x-axis) versus log2 transformed ASE effect size (y-axis) for whole blood (Spearman’s ρ =0.82) and c, subcutaneous adipose (Spearman’s ρ =0.74).

Replication of cis-eQTLs between tissuesa, π1 statistics for cis-eQTLs are reported for all pairwise combinations of discovery (y-axis) and replication (x-axis) tissues. Higher π1 values indicate a stronger replication signal. Tissues are grouped using hierarchical clustering on rows and columns separately with a distance metric of 1 − ρ, where ρ is the Spearman’s correlation coefficient of π1 values. π1 is calculated only when the gene is expressed and testable in the replication tissue. b, c, π1 replication is reported between tissues subsampled down to 70 non-matched (b) and matched (c) donors. In (c), grey tiles indicate tissue pairs with fewer than 70 shared donors. d, Effect of sample size (x-axis) on average π1 replication (y-axis) across all replication tissues. e, Scatter plot of π1 scores among tissue pairs with matched (x-axis) and non-matched (y-axis) donors.

Tissue-specificity of cis- and trans-eQTLsa, Sharing of independently identified cis-eGenes across the 44 GTEx tissues (cis-eGenes are independently identified in each of the 44 tissues and then binned by the number of tissues in which they appear). b, Sharing of cis-eGenes across 44 GTEx tissues that were identified using the hierarchical multi-tissue analysis. c, The prior probabilities of having significant variant–gene association in different numbers of tissues, calculated using an empirical Bayes model. The prior probabilities are summed up on the basis of the Hamming weights of the corresponding cis-eQTL configurations. d–g, Meta-analysis performed using Meta-Tissue for trans-eGenes (50% FDR), randomly selected cis-eGenes (50 % FDR), and an equal number of the top cis-eGenes by P value. Distribution of the number of tissues that have Meta-Tissue m values greater than a given threshold (d, 0.5; e, 0.6; f, 0.9) across variant–gene pairs that have an effect (based on m value thresholding) in at least one tissue. g, The same distribution as d except that variant–gene pairs with predicted effect in zero tissues (based on the number of m values > 0.5) are included. Meta-Tissue predicts that many cis-eGenes (50% FDR) and trans-eGenes (50% FDR) will have an effect in zero tissues. The number of zero predictions is largely reduced for the top most significant cis-eGenes. h, Distribution of observed replication rate between pairs of tissues for trans-eQTLs (10% FDR) versus the predicted replication rate for trans-eQTLs (10% FDR) based on a negative binomial generalized linear model trained on cis-eQTLs (10% FDR0.1). This model directly accounts for effect size and minor allele frequency. Replication rates shown for a range of FDR thresholds in replication tissue. Box plots depict the IQR, whiskers depict 1.5× IQR.

Examples of trans-eQTLs shared across tissuesa, An example of a trans-eQTL (rs7683255–NUDT13) originally identified in sun-exposed skin (10%, P ≤ 1.1 ×10−10, indicated by asterisk) that has a global effect across tissues. The lines represent 95% confidence intervals of the effect sizes. A thicker line indicates that this variant–gene pair is called significant at P ≤ 0.05 in the corresponding tissue. b, An example of a trans-eQTL (rs60413914–RMDN3) that is genome-wide significant in putamen (basal ganglia) (10% FDR, P ≤1.2 ×10−13, indicated by asterisk) that has an effect in all five brain tissues tested but shows little effect in other tissues. c, Expression levels (RPKM) of RMDN3 in all donors across 44 tissues. Box plots depict the IQR, whiskers depict 1.5× IQR.

Sharing information across tissues for cis-eQTLsa, The proportion of expressed genes for which cis-eGenes are discovered in single tissues (5% FDR; origin) and the multi-tissue meta-analysis (m > 0.9; arrow), stratified by the sample size of individual tissues. In the meta-analysis, cis-eQTL discoveries are made using Meta-Tissue to identify tissues where the posterior probability a given cis-eQTL effect exists (that is, the tissue’s m value for the variant) is >0.9. b, The proportion of expressed genes that had a cis-eQTL in the subsampled data (n = 70) is shown on the y-axis, and the actual sample size of the tissue is shown on the x-axis. The proportion is shown for the tissue-by-tissue approach (5% FDR; origin) and using Meta-Tissue (m >0.9; arrow). c, For each of the subsampled tissue data sets (n = 70), we identified the additional K discoveries that were made using Meta-Tissue but were not significant at the 5% FDR threshold in the tissue-by-tissue analysis; we then identified the K most significant cis-eQTLs in the tissue-by-tissue analysis with a q value greater than 5% representing the additional discoveries we would make by simply relaxing the FDR. We then compared these two sets of K cis-eQTLs to the list of significant cis-eQTLs found in the full tissue-by-tissue analysis by calculating the proportion of the K cis-eQTLs that were significant in the full analysis (y-axis). The tissues are ordered along the x-axis by increasing sample size in the actual data set. d, The proportion of annotated and expressed genes that were found to be eGenes using the tissue-by-tissue approach and the hierarchical selection procedure implemented by TreeQTL. e, The number of cis-eGene discoveries per tissue (y-axis) against sample size (x-axis). The number of discoveries for the tissue-by-tissue approach are represented by the origin of each segment, while the number of discoveries from the hierarchical selection procedure are shown as arrows. As with Meta-Tissue, the hierarchical procedure improves cis-eGene discovery for tissues with low sample sizes, albeit fewer tissues overall have the benefits of this effect. Furthermore, an outcome of this procedure is that for tissues with high sample sizes, reported numbers of cis-eGenes are more conservative than those observed in the tissue-by-tissue analyses or Meta-Tissue. f, Improvement of cis-eGene discovery by incorporating genomic annotations. For the 26 tissues for which we can relate cell-type specific chromHMM annotations, we identify cis-eGenes accounting for the variant-level genomic annotations and corresponding enrichment estimates using the Bayesian FDR control procedure described previously80. For each tissue, the number of cis-eGenes identified by the Bayesian procedure (arrow) is plotted against the tissue-by-tissue results (origin).

cis effect size analysesa, For each autosomal protein-coding and lincRNA cis-eGene with eVariants discovered independently in at least five tissues, the median effect size was computed across these tissues. The empirical cumulative distribution function (CDF) of these median effect sizes is depicted. b, Normalized effect sizes of cis-eVariants located upstream of the gene, within the transcript, and downstream of the gene. c, cis-eQTL effect distributions stratified by discovery tissue. Tissues are sorted from largest sample size (muscle-skeletal, n = 361) to the smallest (uterus, n = 70). Box plots depict the IQR, whiskers depict 1.5× IQR.

trans-eQTL discovery restricted to informed subsetsa, Quantile–quantile (QQ) plot of trans-eQTL P values from all variants (x-axis) and the subset of trans-eQTL P values restricted to cis-eVariants (y-axis), illustrating enrichment of low trans-eQTL association P values for cis-eVariants. Data are plotted separately for skeletal muscle (pale blue) and adipose (red). trans-eQTLs with FDR ≤ 10% are shown as large circles, those with FDR > 10% are shown as small circles. b, QQ plot of trans-eQTL P values from all variants (x-axis) and the subset of trans-eQTL P values restricted to GWAS associated variants (y-axis), illustrating enrichment of low trans-eQTL association P values for cis-eVariants. Data are plotted separately for skeletal muscle (pale blue) and lung (green). trans-eQTLs with FDR ≤ 10% are shown as large circles, those with FDR > 10% are shown as small circles.

Replication of cis-eQTLs in TwinsUK dataa, All cis-eQTLs (5% FDR) from six tissues (adipose, subcutaneous; adipose, visceral omentum; cells, EBV-transformed lymphocytes; skin, not sun-exposed; skin, sun-exposed; and whole blood) were examined for replication in four closely matched tissues (LCLs, skin, whole blood, subcutaneous adipose) from the TwinsUK data. For each tissue pair (in facets), replication P value histograms illustrate strong enrichment of small P values. π1 statistics are provided for each tissue pair. b, All trans-eVariant–eGene pairs (10% FDR) from all tissues were examined for replication in four closely matched tissues (LCLs, skin, whole blood, subcutaneous adipose) from the TwinsUK data. Observed replication P values (y-axis) are plotted against the expected uniform P value distribution under the null hypothesis (x-axis). Replication P values are plotted separately for matched (grey) and unmatched (black) tissue pairs.

Patterns of tissue sharing of eQTL effectsa, Similarity (Spearman’s ρ) of Meta-Tissue effect sizes between tissues for cis- (upper triangle, 5% FDR) and trans- (lower triangle, 50% FDR) eQTLs. Tissues (by colours as in Fig. 1a) are ordered by agglomerative hierarchical clustering of the cis-eQTL results. b, The number of tissues in which a given eQTL is shared as a function of tissue sample size. For each tissue, we estimated the degree of sharing (number of tissues with m >0.9) for all eQTLs identified in that tissue at a 5% FDR. Tissues were then binned into quartiles on the basis of sample size. A higher proportion of eQTLs identified in tissues with small sample sizes have shared effects across multiple tissues compared with more deeply sampled tissues. This pattern inverts at higher sample sizes where more of the effects are tissue-specific. The median number of shared tissues is plotted for each quartile as a horizontal black line. c, Distribution of the number of tissues having Meta-Tissue m > 0.5 for the top variant for each trans-eGene at 50% FDR, and FDR-matched, randomly selected cis-eGenes (also 50% FDR). cis-eGenes were matched for discovery tissue to the trans-eGenes.

Disease gene enrichment for tissue-specific and shared cis-eQTLsEnrichment of shared and tissue-specific cis-eGenes in different disease gene data sets. Enrichments and 95% CI in each data set are calculated via Fisher’s exact test, and the odds ratio is plotted after log10 transformation.

General and replicated per-variant cis-eGene associations across tissuesa, Distribution of the number of unique cis-eGenes per variant within each tissue (points) or in the union of variant–eGene associations across all tissues. b, Proportion of variants with top-associated gene preserved between tissues for varying effect size thresholds, across all pairwise tissue comparisons. c, Proportion of variants with top-associated gene preserved between tissues for varying nominal P value thresholds, across all pairwise tissue comparisons. Red distributions include all variants with an associated cis-eGene in one of the two compared tissues, and blue distributions require the variant to have at least two associated cis-eGenes in each tissue. Distributions for which more than half of the pairwise comparisons (points) are empty are not shown.

Broad trans-regulatory locus 9q22 in thyroid tissuea, FOXE1 expression is thyroid-specific. b, Correlation between FOXE1 expression levels and thyroid PEER factors, compared to 100 random genes. For every gene, absolute correlation was sorted in decreasing order. The correlation of FOXE1 with the fifth, sixth, seventh and eighth PEER factors was significantly higher than the correlation of random genes at those rank ordered PEER factors (empirical P ≤ 0.05). c–e, Variants in the chr 9q22 locus were enriched for association with genes on other chromosomes in thyroid carcinomas compared to randomly selected variants nearby randomly selected genes. We used variants that were found within 35 kb upstream or downstream of the gene TSS. f, rs10759975 is associated with trans-eGene TMEM253. g, rs10759975 is associated with trans-eGene ARFGEF3. h, rs10759975 shows cis association with TRMO. i, rs10759975 is weakly associated in cis with FOXE1. Box plots depict the IQR, whiskers depict 1.5× IQR.

Trait-associated variants in skeletal muscle near interferon regulatory factor IRF1a, rs1012793 has broad regulatory impact in skeletal muscle. b, Gene set enrichment for potential trans-eGene targets (identified at P ≤ 0.001) of skeletal muscle 5q31 locus.

Functional characterization of cis-eQTLsa, Enrichment (x-axis) of eVariants in cis-regulatory elements (CREs) across 128 Roadmap Epigenomics project cell types, for each GTEx discovery tissue (y-axis). Enrichment estimated by comparing to random MAF- and distance-matched variants. Stronger enrichment was observed in matched tissues (coloured dots) than in unmatched tissues (box plots). b, Proportion of eQTLs shared between two tissues (m > 0.9) if the eVariant overlaps the same Roadmap annotation in both tissues (y-axis) or different annotations (x-axis). Points represent the mean across all tissues, coloured by the discovery tissue. c, Enrichment of eVariants (y-axis) in tissue-matched enhancers (black) and promoters (grey) for the first four conditionally independent eQTLs discovered for each eGene (x-axis). d, Proportion of eVariants overlapping tissue-matched DNase I hypersensitive sites (DHS; y-axis) as a function of the probability that a variant is causal (x-axis), coloured by the eQTL discovery tissue. e, Normalized absolute eQTL effect size (x-axis) for each eVariant annotation class (y-axis). f, Median (line) and interquartile range (shading) of normalized absolute eQTL effect size (y-axis), as a function of the number of tissues in which the eGene is expressed (x-axis). Box plots depict the interquartile range (IQR), whiskers depict 1.5× IQR. OR, odds ratio.

Functional characterization of GTEx trans-eVariantsa, Frequency distribution of Mendelian randomization t-statistic, for 296 cis–trans-eQTLs and matched background variants. b, CRE enrichment (y-axis) of trans-eVariants (10% FDR), cis-eVariants (10% FDR, to match trans-eVariants), and top most significant cis-eVariants. Box plots show promoter and enhancer enrichment (x-axis) in matched cell-type CRE annotations compared to MAF- and distance-matched background variants. c, Proportion (x-axis) of variants overlapping piRNA clusters, including randomly sampled background loci, trans-eVariants across all tissues, testis trans-eVariants, thyroid trans-eVariants, and trans-eVariants from all tissues other than testis and thyroid. Asterisks denote significant enrichment (permutation test, P ≤1.0 ×10−4). Box plots depict the IQR, whiskers depict 1.5× IQR.

Properties of cis-eQTL overlap with complex trait associated locia, Enrichment of tissue-specific and tissue-shared eGenes in disease and loss-of-function mutation intolerant genes. Tissue-specific and shared eGenes were defined as eGenes in the bottom and top 10% of the distribution of proportion of tissues with an eQTL effect. Bars represent 95% confidence intervals. b, Proportion of eQTLs (y-axis) discovered as a function of P cutoffs (x-axis). c, Proportion of variants (y-axis) with top associated protein-coding gene shared between tissues at varying P thresholds (x-axis). d, Number of GWAS loci (y-axis) and their co-localization results for each of 21 traits (x-axis), coloured by whether the eGene is the closest expressed gene to the lead GWAS variant. e, Proportion of GWAS loci (y-axis) with a significant co-localization for each of 21 traits (x-axis). Box plots depict the proportion explained in each of 44 tissues, red dots depict the proportion explained by the union of all tissues. Box plots depict the IQR, whiskers depict 1.5× IQR.

Characterization of complex trait-associated trans-eQTLsa, Association of rs1867277 with PEER-corrected TMEM253 expression (P ≤ 2.2 ×10−16). b, Quantile–quantile plot of associations between 19 variants in the 9q22 locus and all genes in GTEx thyroid gene expression levels, compared to 19 random variants from the same chromosome, and associations between 23 variants in the 9q22 locus and all genes in TCGA thyroid tumour expression data, compared to 23 random variants from the same chromosome. c, Network depicting cis and trans regulatory effects of rs1012793 mediated through interferon regulatory factor 1 (IRF1). Rs1012793 affects expression of IRF1 in cis and PSME1 and PARP10 in trans (box plots). IRF1 is significantly co-expressed with the trans-eGenes. Colours in scatter plots refer to genotype at rs1012793. d, cis and trans association significance of variants within 1 Mb of the IRF1 TSS in the chromosome 5 locus with cis-eGene IRF1 (blue) and trans-eGene PSME1 (brown), showing concordant signal across the locus. Box plots depict the IQR, whiskers depict 1.5× IQR.

Association of known covariates with expression components removed by PEEREach cell depicts the adjusted (R2) between a pair of variables. Scale bar specific to each panel is displayed at the bottom. Grey cells represent pairs of variables without sufficient data to estimate correlation. a, For each tissue, adjusted (R2) reflecting the proportion of variance explained by each covariate of the entire PEER component removed from the expression data. A selected set of the most relevant sample-specific covariates is shown here. b, For each tissue, adjusted (R2) reflecting the proportion of variance explained by each covariate of the entire PEER component removed from the expression data. A selected set of the most relevant donor-specific covariates is shown here. See Supplementary Information for complete set of covariates. c, Adjusted R2 between each PEER factor and known sample covariates in skeletal muscle. d, Adjusted R2 between each PEER factor and known donor covariates in skeletal muscle.

Controlling the discovery of cis- and trans-eGenesa, b, The number of unique cis- (a) and trans-eGenes (b) across all tissues at varying FDR thresholds. c, The per cent change in the number of cis-eGene discoveries comparing FDR (Storey’s q value) calculated across all tested genes in each tissue to FDR calculated only across autosomal lincRNA and protein-coding genes. Results are shown for each tissue and are stratified by gene type. The per cent change for each tissue and gene type is calculated as 100 ×(no. of cis-eGenes from q value on the restricted gene set – no. of cis-eGenes from q value on all tested genes)/(no. of cis-eGenes from q value on all tested genes). d, The per cent change in the number of cis-eGene discoveries comparing q value and the Benjamini–Hochberg (BH) procedure applied to all tested genes in each tissue. The percent change for each tissue and gene type is calculated as 100 ×(no. of cis-eGenes from BH – no. of cis-eGenes from q value)/(no. of cis-eGenes from q value). Box plots depict the IQR, whiskers depict 1.5× IQR.

Properties of cis-eGenes and non-eGenesIn these analyses, ‘non-eGenes’ refers to the set of genes with no significantly associated cis-eQTLs. a, Density of expression (mean across samples, median across tissues, in log10 scale) for cis-eGenes and non-eGenes. The difference in mean expression for cis-eGenes and non-eGenes is significant (Wald test; P <1 ×10−16). b, Box plots comparing the probability of being loss-of-function-intolerant (pLI) scores41 for cis-eGenes and non-eGenes, stratified by expression percentile across all genes (mean across samples, median across tissues, in log10 scale). On average, highly expressed non-eGenes have higher pLI scores than highly expressed cis-eGenes (t-test for the difference in mean pLI score between cis-eGenes and non-eGenes; BH adjusted P = 8.3 ×10−4 and 2.1 ×10−9 for (50,75]% and (75,100]%, respectively) and lowly expressed non-eGenes (t-test for the difference in mean pLI score between highly and lowly expressed non-eGenes; BH adjusted P ≤7.8 ×10−46, P ≤5.5 ×10−17 and P ≤2.1 ×10−3 for (75,100]% vs [0,25]%, (75,100 vs (25,50]% and (75,100]% vs (50,75]%, respectively). c, Gene Ontology (GO) analysis of tested protein-coding non-eGenes. We used the PANTHER overrepresentation test78 (release 20160715) against 20,972 human genes as background to test for enrichment in GO biological processes using the GO database release 2017-02-28. Significant GO IDs (Bonferroni adjusted P < 0.05) were selected for analysis with REVIGO79 to group similar ontological terms, which yielded 22 over-represented GO IDs. d, GO analysis of a more stringent set of protein-coding non-eGenes. Selected genes included those not tested in GTEx (532 genes) or those with a minimum nominal P value across tissues greater than 0.1 (692 genes). Of these stringent 1,224 non-eGenes, 808 were mapped in the GO analysis. Using a similar approach as in c, 20 over-represented GO IDs were identified. For both c and d, the x-axis represents the −log10 P value resulting from GO analysis. GO IDs are coloured by the broader enrichment category to which each corresponds. Box plots depict the IQR, whiskers depict 1.5× IQR.