We investigated the transcriptomes of a diverse collection of human and mouse tissues and five breast and breast cancer cell lines that were recently sequenced at a depth of roughly 20 million short reads per sample using RNA-Seq protocols (Table S1). Gene expression was initially estimated by calculating read density as ‘reads per kilobase of exon model per million mapped reads’ (RPKM) [13]. These estimates are typically performed using common gene annotations (e.g., RefSeq) with the entire annotated transcript representing the ‘exon model’. These expression level estimates may however be confounded by the expression of shorter isoforms due to alternative cleavage and polyadenylation (Figure S1A and S1B). We found that excluding annotated 3′UTRs – which will sometimes vary between mRNA isoforms as a result of alternative cleavage and polyadenylation – enabled estimation of expression levels that correspond more closely with quantitative RT-PCR measurements (Figure S1C). We noted that removing the 3′UTR from calculation of gene expression yields a >2-fold change for over one thousand genes (Figure S1D), and that the effect of 3′UTRs on expression estimates does not seem to be a technical issue caused by secondary structure in the 3′UTR (Figure S2). We therefore advocate excluding UTRs from such estimates, and all subsequent gene expression estimates described here excluded 3′UTR regions.

We next sought to answer how many genes are expressed in a tissue or cell type. A comparison between the expression levels of exons and intergenic regions was used to first find a threshold for detectable expression above background (Figure 1A, algorithm in Figure S3), yielding a threshold RPKM value of 0.3 which balances the numbers of false positives and false negatives. For individual samples, we obtained threshold values between 0.2 and 0.8. As it is difficult to identify untranscribed DNA regions with confidence [20],[21], it is very possible that the background was overestimated. Applying the threshold 0.3 RPKM, the number of genes expressed in most human and mouse tissues varied from 11,000 to 13,000, corresponding to roughly 60–70% of RefSeq protein-coding genes (Table 1). These gene number estimates were stable across different sequencing depths (Figure 1B) and therefore represent bona fide tissue differences. Testis was a clear outlier, expressing more than 15,000 different genes (84% of RefSeq genes). As many as 7,897 genes (42%) were observed to be expressed in all tissues and cell lines (Dataset S1). The corresponding number for Ensembl annotation was 8,214, or 38% of protein-coding genes (Ensembl is an automated gene annotation system, whereas RefSeq is manually curated). Each ubiquitous gene was typically expressed at roughly the same order of magnitude in all tissues, suggesting that there were few problems with genes being considered ubiquitous when they were really specific to one or a few tissues but had a leaky, non-functional expression elsewhere (Figure S4). While we observed small numbers of reads for 8 genes known to have leaky transcription [22],[23] in several tissues, these genes were all too weakly or narrowly transcribed outside their main tissue to be detected as ubiquitous. The estimated number of ubiquitously expressed genes appeared to plateau as the number of samples used was increased to the full set of 24 (Figure 1C). The detection threshold used affects the number of genes detected (Table 1), and the number of detected ubiquitous genes can vary by up to ∼2,000 genes depending on threshold used. The number of samples is large enough that background is unlikely to cause relatively tissue-specific genes to be detected in every sample. These differences between thresholds therefore most likely reflect the presence of low-abundance RNA species. The number of ubiquitous genes we detected is much greater than the ∼1,000 shared genes identified by SAGE [8] and the 1–6% of genes from microarrays [9]–[11], but is in relatively good agreement with the ∼10,000 shared genes estimated by reassociation kinetics [6] and the 3,140 to 6,909 estimated from ESTs [24] (the higher number came from a cutoff of presence in 16 out of 18 tissues, used to remedy uneven EST sequencing across tissues). The increased number of ubiquitously expressed genes compared to SAGE and microarrays most likely results from the increased depth of mRNA-Seq data and improved detection of lowly expressed genes [22]. The number of genes expressed in a tissue ranged from 11,199 to 15,518 genes (Table 2), so a majority of the genes expressed in a specific tissue or cell type are ubiquitously expressed genes. These genes contribute ∼75% of the polyA⁺ RNA molecules in most tissues (Table 3), although this fraction was higher in the cancer cell lines, perhaps as a result of their elevated metabolic rate.

To characterize the set of ubiquitously expressed genes we had identified, we looked for functional enrichment compared to genes expressed only in a subset of the tissues analyzed (hereafter called non-ubiquitous). The protein products of human ubiquitously expressed genes were more likely to have intracellular localization and to be involved in metabolism and other core cellular functions such as macromolecule synthesis, general transcription and vesicles (Figure 1E). Genes that were expressed in only one or a few tissues were more often secreted or membrane-bound (Figure 1E; Dataset S2 and S3), suggesting that cellular contacts and communication are mediated more often by specialized tissue-specific components. Interestingly, an exception to this inside-outside rule was sequence-specific DNA binding proteins, which are nuclear yet seldom ubiquitously expressed. Among these transcription factors we found that POU, homeobox and forkhead genes had the fewest ubiquitously expressed members, consistent with roles in specifying cell and tissue identity [25], whereas e.g. basic-leucine zipper factors were more often ubiquitous (Table 4). Functional characterization of housekeeping genes has been done in the past [26],[27] (and indirectly by [28]), with comparable results, although transporters were found to be relatively tissue-specific in one study [26]. Rather than looking at ubiquitous expression, that study compared the mean number of tissues where the genes were expressed, which could explain the difference. Ubiquitous genes often had CpG islands near their promoters (Figure 1D), as has been observed previously for ubiquitous and developmental genes [29]. The set of ubiquitous genes with CpG-poor promoters were not enriched for any GO category compared to all ubiquitous genes, nor were those with CpG-rich promoters. These observations suggest that ubiquitous expression is a better indicator of housekeeping functions than promoter CpG content. Together, these analyses suggest that much of the internal cytoplasmic machinery and most nuclear functions are common to most or all tissues, and that a large portion of the differences between tissues lie primarily in expression of receptors and ligands that mediate communication, and in a subset of sequence-specific DNA binding transcription factors.

As RNA-Seq expression measurements are highly quantitative, we also explored tissue transcriptome composition in terms of mRNA abundance classes [1] and the extent to which mRNA populations are dominated by a few highly expressed genes. Genes were sorted according to their expression and the fraction of the total cellular polyA⁺ RNA pool devoted to the most highly expressed genes was determined. This analysis showed that mRNA expression in both tissues (Figure 2A) and cell lines (Figure 2B) followed a continuous distribution rather than separating into distinct abundance classes as reported in previous studies (e.g. [1],[6]).

In muscle and liver transcriptomes, a small number of genes contributed a large fraction of the total mRNA pool, e.g. the ten most highly expressed genes in liver and muscle made up roughly 20–40% of the mRNA population. Other tissue transcriptomes were more complex, with the ten most highly expressed genes contributing only 5–10% of the mRNAs in brain, kidney and testis. The remaining tissues had intermediate levels of complexity (Figure 2A). The breast cancer cell lines had similar or greater complexity than normal breast tissue (Figure 2B). Biological replicates in both human and mouse tended to have highly similar complexity distributions (Figure 2C, 2D). Mouse tissues had somewhat similar profiles to corresponding human tissues (Figure 2D), although a much higher expression of several acute-phase genes in both human liver samples shifted their curves toward lower complexity compared to mouse liver. We conclude that kidney, testes and brain tissues have more complex transcriptomes due to the expression of more genes and with less dominance of a few highly expressed genes, whereas liver and muscle tissues are the least complex and express fewer genes, with more dramatic contributions of highly expressed genes.

We next asked what fractions of total cellular mRNA are allocated to genes involved in different biological processes across the different tissues and cell lines. For this purpose, we developed a tool called FRACT (Functional Relative Allocation of Transcripts) that assesses relative gene expression from RNA-Seq read density for arbitrary sets of genes or broad gene ontology (GO) categories (results for a subset of tissues are shown in Figure 3A). This analysis provided a perspective on the functional priorities of cells in each tissue, since allocating a large fraction of the polyA⁺ RNA content in a cell (and likely of translational capacity) to one functional category represents a major investment of cellular resources. For some categories, including ‘metabolic process’, ‘transport’, and also ‘regulation of cell proliferation’, FRACT allocation varied relatively little across the tissues and cell lines (as measured by the coefficient of variation, CV, of the transcriptome fraction), consistent with the expected ‘housekeeping’ functions of these gene categories. Other categories had a far higher fraction of transcripts allocated to them in one tissue than in others, e.g. immune response (high in lymph node), muscle contraction, heart development and electron transport (all high in heart), and signal transduction and G protein-coupled receptor signaling (both high in brain). These examples, representing more specialized activities expected to be of increased importance in the corresponding tissues, provided a molecular-level validation of the integrity of the tissue samples and protocol used. In some cases, differences not readily apparent from the broad GO categorization shown in Figure 3A, could be detected by finer sub-classification of categories – an example is shown in Figure 3B.

We also investigated the expression of thousands of large non-coding RNAs (ncRNAs). These genes were found to contribute a small fraction of transcripts to polyA⁺ transcriptomes compared to mRNAs (Figure 4A) as a result of their considerably lower expression levels (Figure 4B). These levels are lower than for mRNAs for all degrees of tissue-specificity (Figure 4C).

Muscle and brain tissues from human and mouse were observed to have similar expression and FRACT distributions (Figure 2D and data not shown), raising the question of the extent of conservation of tissue-specific expression patterns. We compared global gene expression levels between human and mouse tissues and observed high correlations between expression of orthologous genes between human and mouse (Pearson correlation 0.76 for muscle, 0.77 for liver and brain). When different tissues were compared (e.g. human brain vs. mouse muscle) substantially weaker correlations were observed (Pearson correlations in the range 0.47 to 0.61). These observations indicate a fairly strong overall conservation of gene expression levels between mouse and man, consistent with previous studies based on microarrays [30].

The lengths of mRNAs were studied by mapping the reads to coding and untranslated regions. Using RefSeq annotations, the density of reads in untranslated regions was lower than in coding regions (Figure 5A), suggesting that expression of mRNAs with UTRs shorter than or distinct from those annotated in RefSeq is common. We therefore estimated the lengths of the UTRs as their relative number of reads to coding regions using the annotated coding region length. Mouse data from [13] was chosen for this analysis as this dataset had little 3′ bias (Figure S5). In all three mouse tissues studied, significant negative correlations were observed between expression level and transcript length (−0.31 in liver and muscle, −0.16 in brain; all tissues p<10⁻⁸⁷), showing that shorter mRNAs tend to be expressed at higher levels (Figure 5B). This result agrees with that from reassociation kinetics data [31]. Weighting each gene by the expression level to obtain length estimates for the bulk mRNA population in tissues to obtain the average mRNA length in each tissue, we found that brain mRNAs have longer 3′UTRs on average than liver and muscle mRNAs, by 300–400 nucleotides (Figure 5C).

To assess the protein functions encoded by transcripts with long or short UTRs, we calculated the median length of 5′ and 3′UTRs of genes associated with each GO biological process category (Figure 6B and data not shown). Transcripts coding for proteins involved in metabolism and RNA processing had the shortest UTRs (medians below 500 bp), while the longest median UTR lengths were observed in transcripts encoding proteins involved in development, morphogenesis and signal transduction (Figure 6A). The median lengths in the longest categories ranged between 1000 and 1500 nt, i.e. two- to three-fold longer than for typical metabolism- or RNA processing-associated transcripts. Some of these differences might reflect an increased role for 3′UTR sequences in localization of proteins to specific membrane locations, likely to be more common for proteins involved in signal transduction and morphogenesis than for metabolic or RNA processing-associated proteins, which are typically cytoplasmic or nuclear, respectively. These differences could also reflect differences in the complexity of translational regulation among these classes of genes.

We used short read data from human tissues from [14] (SRA002355.1) and [18], mouse tissues from [13] (downloaded from http://woldlab.caltech.edu/html/rnaseq), mouse embryonic cell and body data from [16] (http://grimmond.imb.uq.edu.au/mESEB.html) and cerebellum data from non-schizophrenic humans from [48]. See respective papers for details on library preparation, sequencing and general read mapping statistics. The data from [18] were mapped to build hg18 with bowtie [49] with setting –best and ambiguous reads were removed. Two human brain samples were used. The sample with lower sequencing depth from a mix of individuals was used in the comparison with RT-PCR data, while the deeper sample was used everywhere else.

We mapped read positions onto gene models and estimated gene densities as the number of reads divided by the number of read start positions. We used only reads that mapped uniquely to the genome, and only positions where a read could potentially map uniquely counted toward exon length. For testing different ways of measuring gene expression (by removing different parts of the gene structure), we selected a set of genes with >2 exons and only one annotated isoform in RefSeq whose expressions had been measured by the MicroArray Quality Control project [50] in the same two samples, UHR (universal human reference) RNA and brain. Cleavage and polyadenylation sites are from [51],[52]. RefSeq or Ensembl gene annotations without 3′UTRs were then used for all gene expression estimates. For genes with multiple splice variants, we fitted an RPKM value to each variant by least square regression and used the sum of the expression of all isoforms (Figure S6). Isoforms that did not overlap directly but were grouped only through overlap with a third isoform were not considered to represent the same gene. All Pearson correlations were calculated based on log-transformed expression. False discovery and false negative rates were estimated using the algorithm presented in Figure S3, which seeks to correct for the presence of spurious reads mapping to non-expressed genes. The extent of leaky ubiquitous transcription by comparison of the ubiquitous set of genes to shuffled controls (Figure S4).

For three mouse tissues, we calculated RPKM values in the same way as had been done for the human ones. Mouse genes were matched to human orthologs using Entrez Gene. A list of acute-phase genes was taken from http://www.informatics.jax.org. DAVID [53] was used for finding enriched gene ontology categories. Categorization of promoters by CpG content was performed as described in [29]. Transcription factor annotations are from [54].

RefSeq gene annotation was used for protein-coding RNA (i.e. accessions starting with NM_) and curated non-coding RNA (NR_). We used the liftOver tool from the UCSC genome browser to obtain human positions for lincRNA regions from [21].

GO annotations for Ensembl transcripts were downloaded from Ensembl (BioMart). The read density for each transcript in each tissue was distributed among its annotated GO categories (total transcript density/no. GO categories for the transcript). GO categories were sorted by the total transcriptome density across tissues and cell lines, and the 400 categories with greatest density (accounting for 94% of total density) were aggregated into 17 broad classes; the remaining categories (6% of total transcriptome density) were aggregated into an “other” class (see Dataset S4 for mappings). The total density of transcripts devoted to each class in each tissue was tabulated. The coefficient of variation in the fraction of each transcriptome devoted to different classes was computed, and a Z-score for each class was computed to identify particular tissues which devote a significantly different fraction of the transcriptome to particular classes (|Z-score|>2).

The UTR lengths were calculated as the number of reads in a UTR divided by the number of reads in CDS multiplied by the CDS length. For the expression weighted average gene lengths, we used the CDS length from Refseq gene annotation, but weighted according to the expression of each gene. To see the correlation between mRNA length and abundance, we took the CDS length from RefSeq annotation for gene isoforms and added UTR length according to the distribution of reads in the three regions. Only those expressed above 0.3 RPKM were included, in order to exclude genes with few reads that could drive an artificial correlation. To compare 3′ bias between samples, i.e. to what extent genes get more reads as you go in the 3′ direction, we plotted the average read density for all genes (weighted so that each gene contributed equally) across the coding region and fit a line y = kx+m where y = read density, x = location along coding region, and k/m is a measure of 3′ bias.