High-Resolution Mapping of Expression-QTLs Yields Insight into Human Gene Regulation

The Distribution of cis-Acting eQTLs

For each of the 11,446 genes, we tested for putative cis-acting eQTLs by regressing measured mRNA levels against SNP genotypes, independently for each SNP in the cis-candidate region, using a standard linear regression model. Consistent with previous reports [20], we found a substantial number of genes with strong evidence for containing at least one eQTL. A total of 744 genes (6.5%) had at least one SNP with a p-value <7×10⁻⁶. If the smallest p-value in each gene is treated as a summary statistic, this threshold yields a gene-level false discovery rate of 5% [28].

We also observed that, in many cases, the SNPs most strongly associated with mRNA levels for a particular gene lie in a restricted region, allowing relatively precise localisation of eQTLs. Figure 1 plots examples of p-values in three genes, illustrating both the strong association signal that is often achieved, and the relatively localised nature of many of the signals (Figure S5).

Encouraged by the potential for these data to localise eQTLs, we next examined the distribution of the physical location of putative eQTLs within the cis-candidate region. For each gene with an eQTL (defined as having at least one p-value <7×10⁻⁶) we took the position of the most significant SNP as an estimate of the location of the functional site. In practice, we expect that the most significant SNP will sometimes be the actual functional site, but usually it will not since (1) HapMap contains only ≈1/3 of common SNPs [25]; (2) some eQTLs may be due to SNPs in LD with nearby copy number variants, though in practice few of the copy number variants known to be associated with expression are well-tagged by SNPs in these data (data not shown; [19]); (3) a non-functional SNP in strong LD with the functional site may have a smaller p-value by chance. Using simulations we estimate that the median distance between the functional SNP and the most significant SNP in our data is 7.5 kb (Figures S6 and S7). As expected, local recombination rates are strongly inversely correlated with the distance between the functional SNP and the most significant SNP (Figure S8).

Figure 2 shows histograms of the locations of the most significant eQTL SNPs, as a function of gene size. (The plots incorporate a correction factor for the possibility of spurious signals due to undetected SNPs in the expression probes; see Methods.) Several interesting features emerge. First, the distribution of the most significant eQTL SNPs is roughly centered on the transcribed region. Second, nearly all such eQTL SNPs lie close to genes: we find relatively few that are >50 kb from the corresponding gene. Third, as shown in Figure S9, there is a significant enrichment of eQTL SNPs in exons compared to introns. We will return to this observation later in the paper.

Finally, for all three gene sizes, the highest density of eQTLs is around the TSS and immediately upstream of the TES, as reported previously in yeast [22]. The TSS peak was reported in a previous plot of these data [20], but in that previous analysis the TES signal peak was concealed due to the variability of gene lengths (see Figure S10). The TES signal is quite asymmetric: among genes with an eQTL, 10% (75) have the most significant eQTL in the 4 kb upstream of the TES, compared with just 4% (29) in the 4 kb immediately downstream.

A Hierarchical Model of eQTLs in the cis-Candidate Region

While Figure 2 reveals the broad distribution of eQTLs and makes few modeling assumptions, it does not easily enable formal model testing about which aspects of gene structure (or other sequence features) are most important for generating eQTLs. Moreover, since the most significant SNP is not always close to the functional site, this approach can be expected to flatten out the true peaks of eQTLs and to increase the numbers of eQTLs that appear to lie far from the target genes.

Consequently, we next developed a Bayesian hierarchical modeling approach that solves many of these problems (see the Methods for further details). We considered a collection of models in which the parameters predict the prior probability that any given SNP in the cis-candidate region will be an “eQTN” (i.e., the functional nucleotide that creates an eQTL). Each model incorporates information about the physical locations of SNPs and, in some of our models, additional functional annotation of the SNPs. (Our calculations assume that the actual functional site is included in the HapMap genotype data; see below for further discussion.) The model parameters are estimated by maximizing the overall likelihood of the expression data, across all genes.

To implement our hierarchical approach, we switched to using Bayesian regression to test for association between SNPs and gene expression [29] (Methods). For each SNP in the cis-candidate region around a gene, we computed a Bayes factor that measures the relative support for the alternative hypothesis (the SNP is an eQTN) compared against the null (the SNP is independent of gene expression). For these data, the Bayes factors are highly correlated with p-values from standard linear regression. However, a key advantage of Bayes factors is that, combined with the prior probabilities specified by the model, they allow us to compute the posterior probability that each SNP is the actual eQTN.

The hierarchical model shares information across all genes about the distribution of signals and this in turn allows better weighting of which SNPs in individual genes are most likely to be eQTNs. For example, consider a hypothetical gene in which two SNPs that are associated with expression are in perfect LD (r ² = 1). Suppose that one SNP is very close to the TSS, and the other is 30 kb upstream. In the p-value analysis, we would assign each of these SNPs 50% weight. In contrast, the hierarchical model downweights the upstream SNP because it is apparent from the overall data that eQTNs are much more abundant near the TSS, suggesting that the SNP near the TSS is much more likely to be responsible for the signal. Simulations show that the hierarchical model provides a more accurate profile of the distribution of eQTNs (see Figures S5 and S11).

Of course, some degree of complication is added by the fact that current HapMap data do not yet contain all SNPs. Therefore, the sites that we infer to be “eQTNs” in this study surely include many SNPs that are tags of nearby functional SNPs that are not in HapMap. This effect will systematically reduce our estimates of the importance of any particular factor in predicting eQTNs. In the case of factors relating to physical location (such as distance from the TSS) simulations show that this has a modest impact on spreading out the signal peaks that we observe, and that the overall distribution of signals is still estimated very well (see Figure S5, S11, and S12). In contrast, in the case of factors relating to functional categories (e.g., whether a SNP lies in a conserved element) we would expect the impact to be much more serious since functional elements are usually small and tag SNPs are unlikely to fall within the same element as a functional site. A second complication is caused by the possibility that undetected SNPs in the expression array probes might create spurious signals [30]. Our hierarchical model includes an explicit correction for this, using the 634 genes with a known SNP in the probe as training data.

Distribution of eQTNs with Respect to the Transcribed Region

We first set out to get a more refined view of the distribution of eQTNs across the cis-candidate region. The basic versions of our hierarchical model described the positions of SNPs relative to a single “anchor” point such as the TSS. SNPs were grouped into discrete bins based on their distance upstream of the anchor, or downstream (treated separately). The bins nearest the anchor point were just 1 kb wide, to accommodate rapid changes in the rate of eQTNs, while more distant bins were wider (this improves the parameter estimates since the distant bins generally contain few eQTNs). Each bin was associated with a single parameter that relates to the proportion of SNPs in that bin that are eQTNs. The rationale for this model was that it would provide a good description of the data if, for example, the abundance of regulatory elements could be well predicted by distance from the TSS alone.

We also considered models with pairs of anchor points (e.g., the TSS and the TES). In those models, each SNP belonged to two bins, each corresponding to the distance from one anchor point. This model treats the probability that a SNP is an eQTN as the sum of an effect due to the first anchor plus an effect due to the second anchor. Recall that our gene set includes only genes with a single annotated transcript, so that this analysis does not incorporate alternative transcription start or end sites.

Table 1 compares eight different models using either a single anchor point (e.g., TSS or TES), or pairs of anchors (TSS and one other anchor). We used AIC (Akaike Information Criterion) to penalize the two-anchor models for the extra parameters that they use.

In summary, the results provide compelling support for a model including both the TSS and TES over all other models (Table 1). Two other two-anchor models (namely TSS+probe location, and TSS+coding sequence end) also performed well, presumably because the Illumina probes and the coding sequence end positions are usually near to the TES. However, given that the TSS+TES model had by far the strongest support, we use this model in our subsequent analyses.

We next replotted the locations of eQTNs, using the posterior probabilities estimated by the hierarchical model (Figure 3). Compared to the p-value-based analysis, the two strong peaks of signal near the TSS and TES are considerably strengthened. Also, in the hierarchical model, the level of background signal upstream and downstream of the gene is greatly reduced, presumably because most of the background signal in the p-value analysis can be explained as resulting from LD with SNPs near the TSS and TES. The hierarchical model estimates that the total number of eQTLs is considerably larger than the number that we detected by linear regression at the rather stringent false discovery rate of 5% (1586 vs. 744). This difference is partly because the hierarchical model adds fractional probabilities for eQTLs that have only partial support for being true eQTLs, and partly because the hierarchical model is more sensitive to signals in locations that are likely a priori.

Another view of the hierarchical model results is shown in the cumulative plots in Figure 3, which plot the cumulative distribution of eQTNs across the gene region. Most eQTNs lie close to the gene, with less than 7% of the detected cis-eQTNs located more than 20 kb outside the gene. Overall, there are about 3-fold more eQTNs in the upstream region of the gene (5′ of the TSS) than downstream (3′ of the TES) (30% vs. 9%).

We next investigated the peaks of signal near the TSS and TES in more detail, using a finer bin partition close to the TSS and TES (see Figure 4A and Methods). At this finer scale, the TES signal is extremely sharply peaked over a region of just ∼100 bp immediately upstream of the TES. The data strongly reject a model in which the signal is symmetric around the TES (p = 3×10⁻⁷). In contrast, the TSS signal is more spread out, and spans both sides of the TSS. There is no evidence of asymmetry in the TSS signal (p = 0.34).

We also observed that the TSS and TES peaks both correspond with two parts of the typical gene structure that, averaging across all 11,446 genes, tend to be highly conserved across the mammalian phylogeny (Figure 4B). The correspondence of the two eQTN peaks with the peaks of conservation suggests that there may be a causal link between these two types of signals. We propose that the sequence conservation reflects, at least in part, the roles of these two locations in regulating mRNA levels, though further work will be needed to verify the connection.

Similarly, the TSS peak also matches up closely with the peak binding densities of a collection of transcription factors that are involved in transcription initiation (reported previously by the ENCODE group, based on ChIP-chip data collected for a set of regions spanning ∼1% of the genome [4]). As might be expected, the ENCODE data identified almost no transcription factor binding near the TES. We return to these latter observations in the Discussion.

Distribution of eQTNs with Respect to Functional Annotation

We next used our hierarchical model to examine the impact of various types of functional annotation on the probability that a SNP is an eQTN. We first classified SNPs that lie inside genes into categories based on the exon/intron structure (e.g, first, coding and last exons; first, internal, and last introns; Figure S13). In order to make the model fully identifiable, we estimate the effect of each annotation relative to the abundance of eQTNs in internal introns (as this category has the greatest number of SNPs). Since gene position is highly predictive of eQTN abundance, we controlled for SNP position using the TSS+TES model. In effect, the hierarchical model now tests whether the annotation adds any predictive value beyond the basic position information. As noted above, incomplete SNP ascertainment in HapMap means that we will generally underestimate–perhaps substantially–the impact of relevant annotations.

The main result of this first analysis is that internal introns have a deficit of eQTNs, compared to coding exons, as well as first and last exons and introns (Figure 5, Table S1). For example, SNPs in coding exons are ∼2-fold more likely than SNPs in internal introns to be eQTNs. First introns are also relatively enriched for eQTNs compared to internal introns (controlling for position). However, since the total amount of sequence contained in introns vastly exceeds that in exons, 53% of genic eQTNs lie in internal introns compared to 10% in coding exons (see Table S1). SNP density differs slightly between exons and introns, but not nearly enough to generate a 2-fold difference in eQTN abundance (Table S2). Overall, the hierarchical model that includes the gene structure annotation as well as position effects relative to the TSS and TES is substantially better than the TSS+TES-only model (by 7 units of AIC).

We then considered the impact of a variety of other types of SNP annotation (see Methods and Figure S14). None of these annotations shows convincing signals of enrichment (Table S3). We estimate a 1.9-fold enrichment of eQTNs inside conserved noncoding elements, as might be expected if these identify functional elements, however the 95% confidence interval narrowly overlaps 1. We also tested for an enrichment of eQTNs at computationally predicted microRNA binding sites, reasoning that SNPs in these binding sites might affect mRNA degradation. We found a suggestive, but non-significant, enrichment of eQTNs in these sites (1.4-fold). It is unclear whether the absence of significant effects in these analyses indicates that these types of annotation are not strongly associated with eQTNs or instead reflects the incompleteness of HapMap and the limitations of current functional annotations.

Finally, based on ENCODE results showing that the promoter regions of genes with CpG islands tend to have more accessible chromatin and greater occupancy by transcription factors [4], we predicted that CpG status might also provide relevant annotation. Indeed, we find that genes with a CpG island spanning the TSS are expressed at higher average levels, and are ∼50% more likely than genes without a CpG island on the TSS to have a cis-eQTN (15% vs 11%). This effect is consistent with the observation that genes with CpG islands are more likely to be expressed in a wide range of tissues than are genes without CpG islands [31]. After adjusting for the different overall rates of eQTNs, the distribution of signal locations in the two classes of genes is very similar (Figure S15).

Genotype Data

We analyzed genotype and expression data from 210 unrelated individuals studied by the International HapMap project [24],[25]. These include 60 Yoruba (YRI) and 60 CEPH (CEU) parents, and 45 unrelated Chinese (CHB) and 45 unrelated Japanese (JPT) individuals. We used the HapMap Phase II genotype data, release #21 (phased and with missing data imputed). We used data from the 22 autosomal chromosomes only, giving a total of 3,304,587 SNPs. Since allele frequencies in CHB and JPT are extremely similar [24], these two samples were treated as a single analysis panel of 90 Asians (“ASN”).

Gene Expression Data

We used gene expression levels that were measured previously in lymphoblastoid cell lines from all 210 unrelated individuals, using Illumina's human whole-genome expression array (WG-6 version 1) [19]. We downloaded the data that were normalized first by quantile normalization within replicates and then median normalized across all HapMap individuals [19] [ ftp://ftp.sanger.ac.uk/pub/genevar/].

Since mean expression levels at many loci differ between the HapMap populations [26],[42],[20],[27], there is a potential for spurious eQTLs in the combined sample due to population structure. To control for this effect, we applied a normal quantile transformation to the data for each gene, within each HapMap population (ASN, CEU, YRI), prior to combining the samples. That is, for each gene, separately in each population, we transformed the rth largest gene expression value to the (r−0.5)/nth quantile of the standard normal distribution, where n is the number of individuals with gene expression data from that population [29]. By forcing each population to have the same distribution of expression values, we avoid concerns about spurious associations due to allele frequency differences between the HapMap populations. (Note that the overall results within populations are very similar; Figures S2, S3, and S4.) This normalization also reduces the effect of outlying expression values on the regression [29].

Selection of Genes and Probes

We used BLAT [43] to map the 47,294 Illumina array probes onto human RNA sequences from RefSeq (hg18) [44]. The accession numbers of the RNA sequences were mapped against the Entrez Gene database and all probes that mapped with greater than 90% identity to multiple genes were discarded. Of the remaining probes we retained only those with exact matches to a unique gene, leaving us with 19,536 valid probes. Of these, we kept the 13,244 probes for which the gene has a single RNA accession in the RefSeq database. This was done to simplify the analysis by avoiding genes with multiple splice forms or multiple annotated start sites, etc. These 13,244 probes map to 12,227 unique autosomal genes.

Of these 12,277 genes, 85% contained exactly one probe. For the genes with multiple probes, we analyzed only a single probe, selecting the probe nearest to the 5′ end of the gene. We selected this probe because overall the probes are strongly biased towards the 3′ end of the gene, and we wanted to reduce this bias as far as possible. Then, we removed 634 genes for which there was at least one HapMap SNP inside the probe since it is known that such SNPs can impact the measured expression level [30]. Finally, 147 very large genes (size greater than 500 kb) were discarded, leaving our core data set of 11,446 genes.

Gene Structure and SNP Annotation

Gene structure annotation was obtained from the RefSeq gene table [44] for human genome build 35 (hg17). For each gene the TSS and the TES genomic locations were obtained from the fields “Transcription start position” and “Transcription end position” of the RefSeq table, respectively. We checked the genomic positions of the TSSs against dbTSS, a database of experimentally-determined TSSs, [45] and found no differences among the 84% of gene transcripts in our data set that are also in dbTSS. We defined the CDS (coding sequence) to be everything between the translation start and stop positions defined by the fields “cdsStart” and “cdsEnd”, respectively, of the RefSeq Table. We then assigned every genic SNP to one of 8 mutually exclusive gene-related annotations (see Figure S16):

• First (non-coding) exon. If the gene has at least 2 exons, this is the part of the first exon that is not located inside the CDS. If the gene has only one exon, we do not consider it to have a first exon.

• First intron. If the gene has at least 2 exons, this the intron following the first exon, provided that it is not located inside the CDS. Otherwise there is no first intron.

• Noncoding exon. This is any part of an exon located outside the CDS region and excluding the first and last exons.

• External intron. This is an intron located outside the CDS region and excluding the first and the last introns.

• Coding exon. This is any part of an exon located inside the CDS region. Note that exons containing the translation start or stop generally contain both coding exon and noncoding (or first/last) exon. Coding SNPs were further subdivided into synonymous and nonsynonymous, according to their annotation in dbSNP.

• Internal intron. This is an intron located inside the CDS region.

• Last intron. If the gene has at least 2 exons, this is the intron preceding the last exon, provided that it is not located inside the CDS. Otherwise there is no last intron.

• Last (noncoding) exon. If the gene has at least 2 exons, this is the part of the last exon that is not located inside the CDS. Otherwise there is no last exon.

We also included annotations that indicate whether a SNP is in the following special categories: SNP is in a (1) CpG island; (2) conserved noncoding region; (3) predicted cis-regulatory module; (4) predicted micro-RNA binding site; or that (5) a predicted binding site of the CTCF insulator protein lies between the SNP and the TSS. See the Supplementary Methods (Text S1) for further details.

Finally, note that in our analysis design, each SNP is tested for association with every gene that is within 500 kb. This means that typical SNPs contribute data to multiple genes. Our analysis treats these multiple tests as independent, which is likely a good approximation since we identified only five SNPs that are eQTLs for > one gene in cis.

Statistical Analysis

P-Value Method

In the first part of the paper we used standard linear regression to test the gene expression data at each gene for association with SNPs in the cis-candidate region, as follows. The effect of individual i's genotype at SNP j (g_ijk) on his/her gene expression level (y_ik) is assumed to follow an additive linear model:

(1)

where μ is the mean expression level at that gene for individuals with g = 0, where a_jk is the additive effect of the minor allele at SNP j and ε_ijk is the residual. A standard p-value from a 1 df test can then be obtained for the hypothesis that SNP j is an eQTN for gene k (a_jk ≠ 0).

We used the following procedure to generate the results plotted in Figure 2. For each gene with expression data we assigned each SNP in the cis-candidate region to a single bin (see below). Let m be the total number of SNPs that fall into bin b, summing across all genes. (Note that most SNPs are in the cis-candidate regions of multiple genes and hence can contribute data to multiple bins.) Next, for each gene, we tested every SNP for association with gene expression. If the p-value of the most significant SNP was <7×10⁻⁶ then we considered this to be one “signal” in the bin that this SNP lies in. (Note that the results are robust to the choice of the p-value cutoff; Figures S17 and S18.) For genes in which the smallest p-value was shared by n>1 SNPs, we considered that the signal was divided equally among the n most significant SNPs (i.e., a fraction 1/n of the signal was assigned to each SNP). Suppose that, by this way of counting, there are s signals in bin b.

Prior to reporting the data, we also applied a correction for the possibility of spurious signals due to ungenotyped SNPs in the expression array probe. We used the 634 genes with a known HapMap SNP inside the probe to create a profile of the abundance of spurious signals as a function of distance from the probe. This profile was used to adjust the observed number of signals, s, to a corrected number s′, that removes the predicted number of spurious signals in each bin (see Figure S19 and Text S1 for details). In practice, we estimate that the contribution of spurious signals does not substantially change the overall uncorrected distribution of signals. Finally, we computed the fraction of most significant SNPs in bin b as s′/m.

Hierarchical Model

We present here an overview of the hierarchical model. Complete details on the models are provided in the Supplementary Methods section (Text S1).

Bayesian Regression Model

The hierarchical model applies the Bayesian regression framework of Servin and Stephens [29]. The effect of individual i's genotype at SNP j (g_ijk) on his/her gene expression level (y_ik) is assumed to follow a linear model:

(2)

where μ is the mean expression level at that gene for individuals with g = 0, and where a_jk and d_jk are the additive and dominance effects of the minor allele at SNP j. The residual, ε_ijk, is assumed to be N_(0,1/τ) and independent for each y_ik, where 1/τ is the variance of expression levels within each genotype class. The indicator function I(g_ijk = 1) is defined as 1 if the genotype is heterozygous (g_ijk = 1) and 0 otherwise.

Let denote the probability of the expression data Y_k under the null hypothesis that there are no cis-eQTNs in gene k (i.e., a_jk = d_jk = 0 for all j). Similarly, let denote the probability of the expression data Y_k assuming that SNP j is the eQTN. In this case, the effect sizes a_jk and d_jk are modeled as being drawn from mixtures of normal distributions centered on 0 (see Text S1 for details). The Bayes factor (BF) for SNP j in gene k is defined as

(3)

and measures the relative support for the hypothesis that SNP j is an eQTN for gene k, versus the null hypothesis. We use priors on effect sizes that allow the BF to be calculated analytically (see Text S1).

The Hierarchical Model

We describe first the basic version of our hierarchical model. All the results presented in this paper additionally include a correction for the possibility that genes might show signals due to undetected SNPs in the probe. We describe that extension later in the Methods, briefly, and in detail in the Supplementary Methods (Text S1).

Our basic model assumes that there are two mutually exclusive categories of genes. With probability Π₀ there is no eQTN in the cis-candidate region, and with probability Π₁ = 1−Π₀ there is a single eQTN. Then the likelihood of the expression data at gene k is

(4)

where denotes the probability of the expression data Y_k given that there is no eQTN in gene k and denotes the probability of the expression data given that there is exactly one eQTN. Note that our model allows for at most one eQTN per gene. If in fact there is more than eQTN, our model will usually assign the signal to the strongest of these. In practice, we see little variation in average effect size as a function of location, so this modeling simplification is unlikely to seriously distort the results.

Given that there is a single eQTN in gene k, the probability of the observed expression data, , can be written as

(5)

where is the probability of the expression data given that SNP j is an eQTN, and π_jk is the prior probability that SNP j is an eQTN, given that exactly one SNP in gene k is an eQTN.

A key feature of the hierarchical model is that the probability that SNP j is an eQTN, π_jk, is allowed to depend on the physical location of SNP j relative to one or more “anchor” points, and other relevant annotations (see Text S1). Suppose that we consider L different kinds of annotation, and let the indicator δ_jkl equal 1 if SNP j at gene k has the lth annotation, and equal 0 otherwise. Then define

(6)

where Λ = (λ₁,…,λ_L) is a vector of annotation effect parameters. We use a logistic model to relate π_jk to these annotation indicators, namely,

(7)

As detailed in the Supplementary Methods (Text S1), we parameterized the effect of distance from the anchor locations using a series of discrete bins that represent absolute physical distance from the relevant anchor. The bins nearest to the anchor are 1 kb wide, and increase in width to 10 kb and finally 100 kb with increasing distance from the anchor. For the two-anchor models, each SNP belongs to two position bins, each of which indicates distance from one anchor.

Likelihood for the Hierarchical Model

Substituting the above expressions for into (4), the likelihood for the hierarchical model is

(8)

(9)

where Θ denotes the model parameters and BF_jk is the BF from the Bayesian regression (3). To be explicit, the model parameters Θ include the annotation parameters Λ, the proportion Π₀ and other parameters related to the Bayes factor computation (see Text S1). The likelihood of the entire data set is the product of (9) across all K genes. We fit the hierarchical model by maximizing the log-likelihood

(10)

with respect to the model parameters Θ. (Note that the first term, involving does not depend on Θ, and so need not be evaluated.)

Accounting for the Effects of SNPs in Probes

Since undetected SNPs in the probe sequence sometimes generate eQTLs, the results that we report include a modification to account for this effect. We used the 634 genes that have a known SNP in the probe region as training data to help parameterize the model. We assume that these represent ∼1/3 of all probes with common SNPs [25].

Suppose that with probability there is a gene inside the probe sequence (this is set to 1 for the training data), and suppose that when there is a SNP in the probe, there is a probability Π_s that this generates a spurious signal. Then let be the probability of a spurious signal. We consider that we are only interested in real signals if there is no spurious signal, so we write the probability of the data as

(11)

where the first term is the likelihood when there is no spurious signal (as in Equation 4), and where the second term gives the likelihood () when there is a spurious signal.

Posterior Probabilities

Once the likelihood has been maximized, we can compute the posterior probability of a given SNP j to be an eQTN for gene k. In the case without spurious signals this is

(12)

and the general version is given in the Supplementary Methods (Text S1).

Sequence Conservation and Transcription Factor Binding

To compute the average sequence conservation as a function of position for Figure 4B, we estimated the average number of substitutions per site across the phylogeny of seven mammalian species (human, chimpanzee, macaque, mouse, rat, dog, and cow), using data and alignments from the UCSC browser. This was done for the main set of 11,446 genes analyzed in this paper. For each gene, 5 kb on each side of the TSS (and separately for the TES) was split into non-overlapping 50-bp bins. We then concatenated all the sites across all genes that lay in the same bin. After excluding sites in coding exons we estimated the average number of substitutions at each site using baseml, a program in the PAML package [47].

We obtained results on transcription factor binding density using ChIP-chip data collected by the ENCODE project (4). We used data for eight transcription factors that showed large numbers of binding fragments at a 1% false discovery rate in the ENCODE study. The left-hand panel of Figure 4C is essentially a replotting of data presented in Figure 5 of (4), while the right-hand panel shows analogous data plotted with respect to the TES.