iRNA-seq: computational method for genome-wide assessment of acute transcriptional regulation from total RNA-seq data

Abstract

RNA-seq is a sensitive and accurate technique to compare steady-state levels of RNA between different cellular states. However, as it does not provide an account of transcriptional activity per se, other technologies are needed to more precisely determine acute transcriptional responses. Here, we have developed an easy, sensitive and accurate novel computational method, iRNA-seq, for genome-wide assessment of transcriptional activity based on analysis of intron coverage from total RNA-seq data. Comparison of the results derived from iRNA-seq analyses with parallel results derived using current methods for genome-wide determination of transcriptional activity, i.e. global run-on (GRO)-seq and RNA polymerase II (RNAPII) ChIP-seq, demonstrate that iRNA-seq provides similar results in terms of number of regulated genes and their fold change. However, unlike the current methods that are all very labor-intensive and demanding in terms of sample material and technologies, iRNA-seq is cheap and easy and requires very little sample material. In conclusion, iRNA-seq offers an attractive novel alternative to current methods for determination of changes in transcriptional activity at a genome-wide level.

INTRODUCTION

Recently, RNA-seq has become the method of choice for assessing global changes of the transcriptome. This method provides several advantages over microarray analyses including higher reproducibility and greater dynamic range (1–7), reviewed in (8). The ease of performing RNA-seq has made it a standard technology to accompany almost any analysis of the genome. For example, it is used to correlate changes in the epigenome or binding of transcription factors to changes in the transcriptome. Here, changes in the transcriptome, as determined by RNA-seq, are used as a proxy for changes in transcriptional activity. Although this approximation works reasonably well, it is clearly far from perfect, since the amount of transcripts are also regulated at other levels than transcription. In addition, acute changes in transcriptional activity will be partially masked by existing mRNA transcripts. Thus, the transcriptome approximation is particularly problematic when addressing acute effects, and techniques that more directly assess transcriptional activity are required for such studies.

Over the past few years, multiple techniques have been developed to determine changes in transcription at a genome-wide level. These include RNA polymerase II (RNAPII) chromatin immunoprecipitation (ChIP)-seq (9,10), global run-on (GRO)-seq (11,12), 4-thiouridine (4sU)-RNA-seq (13) and bromouridine (Bru)-seq (14). These techniques, each have their strengths and weaknesses; however, all are significantly more laborious techniques and require much larger numbers of cells compared to RNA-seq. For example, RNAPII ChIP-seq involves DNA-protein cross-linking by formaldehyde and sonication of chromatin before immunoprecipitation and purification of RNAPII-associated DNA. GRO-seq requires isolation of nuclei, an in vitro transcription reaction, and subsequent pull down of labeled RNA, whereas techniques, such as 4sU-RNA-seq and Bru-seq, rely on metabolic labeling of cells prior to harvest and isolation of newly synthesized RNA. Thus, the applicability of these techniques for experiments with a restricted number of cells is limited, as is the potential for high-throughput experiments.

Here, we have developed a computational method, iRNA-seq (intron RNA-seq), for easy genome-wide determination of transcriptional activity based on total RNA-seq data sets. The method has been implemented in a Perl pipeline that quantifies differences in intron reads. We have used iRNA-seq to analyze our own unpublished data on the acute transcriptional response of human adipocytes to tumor necrosis factor (TNF) treatment, as well as data derived from the literature. We demonstrate that this new method is a sensitive, fast and easy way of simultaneously determining transcriptional activity and levels of mature transcripts at a genome-wide level from total RNA-seq data.

MATERIALS AND METHODS

Cell culture

Human SGBS cells were obtained from Dr. Martin Wabitsch, University of Ulm, Germany. Cells were passaged and differentiated to adipocytes as previously described (15).

RNA-seq

Following Isol® extraction and column purification of total RNA, ribosomal RNAs were removed using the Ribo-Zero® Human/Mouse/Rat kit (Epicentre). Library preparation was performed using TruSeq RNA Sample Preparation protocol according to the manufacturer's (Illumina) instructions.

cDNA synthesis and quantitative real-time polymerase chain reaction (qPCR)

cDNA synthesis and real-time qPCR were performed as previously described (16). Sequences of primers used for real-time PCR are available upon request.

ChIP-seq

ChIP experiments were performed according to standard protocol as described in (17). The RNAPII antibody used was from Diagenode (C15200004). Library preparation was performed as described in (18).

Additional data

Total RNA-seq data from TNF stimulation of human A549 cells (19) were downloaded from NCBI Sequence Read Archive (accession SRP020499). Total RNA-seq, GRO-seq and RNAPII ChIP-seq data from TNF stimulation of human IMR90 fibroblasts (20), 4sU-RNA-seq data from LPS stimulation of mouse dendritic cells (13), were downloaded from GEO data set browser (accession GSE43070 and GSE25432, respectively).

Data processing

All RNA-seq reads were mapped to their respective reference genomes with STAR (21) using default parameters. ChIP-seq and GRO-seq data were mapped to their respective reference genomes with STAR specifying –alignIntronMax 1 to avoid potentially aligning across exon–exon junctions.

Definition of unique intron, exon and gene regions

All RefSeq genes, exons and introns were extracted from the UCSC Genome Browser (22), and the gene lists were collapsed to the longest transcript for each gene. For each gene, regions overlapping another coding or non-coding gene were removed, so that only regions unique to a specific RefSeq gene were used for quantification. Lists of unique exon and intron regions were generated in a similar manner. Furthermore, for the intron list, all overlaps with genomic locations associated with mRNA sequences were subtracted. These regions were extracted from the UCSC Genome Browser (22), which uses all mRNA sequences submitted to the Genbank to create a list of genomic regions of origin of mRNA. For quantification of GRO-seq and RNAPII ChIP-seq, promoter proximal regions, i.e. regions from −1000 bp to +500 relative to transcription start sites were excluded to avoid quantification of stalled polymerase.

iRNA-seq pipeline

For read quantification and differential expression analysis, a Perl pipeline iRNA-seq was created that takes aligned RNA/GRO/ChIP-seq reads in either SAM or BAM format as input and uses featureCount (23) to quantify reads in all regions defined as unique introns, exons or genes. For each gene the sum of read counts in unique intron regions were used for quantification of primary transcripts (transcription), whereas unique read counts in exons were used for quantification of mature transcripts. iRNA-seq can then either analyze these summarized counts for differential expression by standard or blocked two-condition comparison using edgeR (24) or provide summarized non-normalized read counts for other purposes. iRNA-seq comes with gene, exon and intron lists for the human (hg19), mouse (mm9) and rat (rn5) genomes, as well as a script to generate custom list for other genome versions or organisms. The pipeline and instructions on how to use it is available at: http://www.sdu.dk/mandrupgroup.

Data access

The RNA-seq and RNAPII ChIP-seq data sets generated in this study have been submitted to the NCBI Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) under accession number GSE60462.

RESULTS

Building the iRNA-seq pipeline

Recent studies have demonstrated that information about nascent transcripts, co-transcriptional splicing and mRNA dynamics can be extracted from total RNA-seq data by analyzing intron reads (25,26). We therefore investigated if analysis of intron reads from total RNA-seq data, can be used for genome-wide assessment of acute transcriptional regulation of gene expression. For this purpose, we analyzed total RNA-seq and RNAPII ChIP-seq data obtained from in vitro differentiated human SGBS adipocytes treated with TNF or vehicle for 90 min. Inspection of the data in the UCSC genome browser indicated that at several loci, acute gene regulation could be detected at the level of intron reads as well as RNAPII occupancy, despite no apparent effect on transcript levels as determined by reads in exons (Figure 1A and B).

To investigate this at a genome-wide level, we designed a Perl pipeline, iRNA-seq, which for every gene quantifies and sums intron reads not overlapping with any other transcripts or non-coding RNA. Intron reads are counted and summed using featureCounts (23), and differential expression analysis is performed using edgeR (24) (Figure 2A). To minimize problems arising from inclusion of exon reads due to differential exon usage (Figure 2B, arrow A) and incomplete annotation (Figure 2B, arrow B), all regions associated with mRNA were subtracted from the regions subjected to counting (see Materials and Methods). In addition to the default quantification of intron reads, the iRNA-seq tool also allows quantification of reads in exons or reads in full-length genes. Thus, the iRNA-seq tool allows fast and easy parallel assessment of mature transcript levels as well as active transcription from one total RNA-seq experiment.

Evaluation of iRNA-seq

Application of the iRNA-seq pipeline for analysis of total RNA-seq data obtained from SGBS adipocytes stimulated with TNF or vehicle for 90 min, shows that significantly more genes are scored as differentially expressed by quantification of reads in introns as compared to reads in exons (Figure 3A–C). As would be expected, the difference is more pronounced for genes repressed by TNF than for genes activated by TNF, since reductions in transcription are partially masked by pre-existing mature mRNA (Figure 3C). This should be particularly true for transcripts with long half-lives, and to investigate this we obtained mRNA half-lives determined by an ActinomycinD-chase in human foreskin fibroblasts performed by the Wilusz laboratory (27). Indeed, genes that are detected as regulated by iRNA-seq have longer mRNA half-lives on average than genes that are differentially expressed at the exon level (Figure 3D). To confirm the differential regulation of genes that are detected by iRNA-seq, we performed qPCR of 3 induced and 4 repressed genes none of which could be detected as differentially expressed by quantification of exon reads. In all cases, qPCR using primers targeting introns revealed changes in expression that were not detected using primers targeting exons (Figure 3E). Importantly, fold change in transcription as determined by iRNA-seq correlates well with fold change as determined by RNAPII ChIP-seq (Figure 3F), whereas the correlation is weak between exon reads and RNAPII ChIP-seq data (Figure 3G).

To investigate the correlation between sequencing depth and iRNA-seq performance, we ran iRNA-seq on different subsamples of the RNA-seq data (Figure 3H). The dependency on sequencing depth is higher for iRNA-seq compared to quantification of exon reads due to the lower number of reads in introns (Figure 3I). Importantly, however, the number of significantly regulated genes detected by iRNA-seq is higher for all sequencing depths. Whereas biological duplicates markedly improved the performance of iRNA-seq, triplicates only modestly improved detection of significantly regulated genes in this system (Figure 3J). The same tendency was observed when analyzing exon coverage.

Time-course analysis with iRNA-seq

To investigate the relationship between treatment time and the benefit of iRNA-seq over regular RNA-seq, we downloaded data from a time-course experiment from the Brasier laboratory in which human A549 cells had been treated with TNF for different time points up to 6 h (19). These data were then processed through the iRNA-seq pipeline using either intron or exon reads. Similarly to what was observed for the SGBS adipocyte data, we identified markedly more differentially expressed genes using intron reads as compared to exon reads following 30 and 60 min TNF treatment; however, this difference levels off after 3 h (Figure 4A). Consistently, fold changes at the intron and exon levels correlate weakly at the 30- and 60-min time points, whereas there is a strong correlation at the 3- and 6-h time point (Figure 4B). Interestingly, acute transcriptional changes (30/60 min intron) show the strongest correlation with exon levels at the 3-h time point, whereas it weakens a little at the 6-h time point, possibly due to the effect of secondary transcriptional events at mRNA levels at later time points. In general, there is strong correlation between changes in primary transcript at a given time point with mature mRNA levels 2–3 h later, indicating that most transcriptional changes are translated into changes in mRNA levels within 2 h.

iRNA-seq performance is comparable to GRO-seq and RNAPII ChIP-seq

To test the iRNA-seq method in a different model system and to directly compare it with other approaches for analyzing acute transcriptional regulation, we downloaded total RNA-seq data from an experiment from the Ren laboratory, in which IMR90 fibroblasts had been treated with TNF or vehicle for 60 min (20). These data were processed through the iRNA-seq pipeline in intron and exon mode, and parallel GRO-seq and RNAPII ChIP-seq data were processed in the gene mode, with exclusion of promoter proximal regions in order to avoid quantifying stalled RNA polymerase. Consistent with the results obtained for SGBS adipocytes (Figure 3F and G), fold change determined by iRNA-seq correlates strongly with changes in RNAPII occupancy and change in GRO-seq signal (Figure 5A). By contrast, quantification of exons correlates weakly with RNAPII ChIP-seq and GRO-seq data. To directly compare the ability of these methods to detect changes in transcription, we determined the number of differentially transcribed genes using different subsamples of each data set (Figure 5B). At low sequencing depths GRO-seq and RNAPII ChIP-seq detect more differentially transcribed genes than iRNA-seq; however, the performance of iRNA-seq increases dramatically with sequencing depth. Above 30–40 mill reads, i.e. in the range recommended by ENCODE for RNA-seq (https://genome.ucsc.edu/ENCODE/protocols/dataStandards/RNA_standards_v1_2011_May.pdf), the performance of iRNA-seq is comparable to that of RNAPII ChIP-seq and GRO-seq, although GRO-seq remains the most sensitive method. In contrast, quantification by RNAPII ChIP-seq benefits less from increasing the sequencing depth above 10 mill reads, which is most likely due to lower signal-to-noise ratio of the RNAPII ChIP-seq method compared to the other methods.

We were unable to directly compare iRNA-seq with labeling-based techniques, such as 4sU-seq or Bru-seq, since we could not identify any parallel data sets. However, we downloaded 4sU-RNA-seq data from an experiment by the Regev laboratory, where mouse dendritic cells had been treated with LPS for 1 h (13). When these data were analyzed using the iRNA-seq pipeline, we identified significantly more regulated genes when restricting the analysis to only use the introns compared to using the full-length genes (Figure 5C), demonstrating that detection of regulated genes by 4sU-RNA-seq can be significantly improved by disregarding exon information in the iRNA-seq pipeline. This is particularly true for repressed genes, which may be due to the rather long labeling time used in this experiment (45 min (13)), resulting in inclusion of processed, stable transcripts, that were synthesized early in the labeling period. Thus, the levels of labeled transcripts do not reflect the transcriptional activities at the time of harvest but rather reflects the average transcriptional activity over the labeling period.

Taken together, these analyses demonstrate that with a sequencing depth in the range where total RNA-seq is usually performed there are sufficient reads in introns to quantify transcriptional activity with high sensitivity and accuracy using the iRNA-seq pipeline. Direct comparison with current more labor-intensive methods for assessment of transcription demonstrates that iRNA-seq can provide very similar output.

DISCUSSION

In this study we have developed a new pipeline, iRNA-seq, which counts and sums reads from total RNA-seq data mapping to uniquely annotated intron regions. We demonstrate that this is an accurate and sensitive method for genome-wide quantification of acute transcriptional changes from total RNA-seq data.

We used the iRNA-seq pipeline to analyze our own total RNA-seq data from TNF-treated adipocytes as well as to reanalyze data from the literature. These analyses clearly showed that while changes in mRNA levels may be an acceptable proxy for transcriptional changes over several hours, it is a poor approximation for detection of acute transcriptional changes (minutes to few hours). We demonstrate that such acute transcriptional changes are much more efficiently detected by iRNA-seq. Interestingly, however, we noted that there was a reasonable strong correlation between reads in introns at a given time point and reads in exons 2–3 h later. This is well in line with recent findings by the Lazar laboratory comparing transcriptional changes by GRO-seq with microarray data following exposure of 3T3-L1 adipocytes to rosiglitazone for different time periods (28).

Comparison of iRNA-seq with more labor-intensive methods, such as GRO-seq and RNAPII ChIP-seq, shows that iRNA-seq estimates similar fold changes in transcription and has similar sensitivity as these methods. Due to the low number of reads in introns, the sensitivity of iRNA-seq is more dependent on sequencing depth than the two other methods; however, this limitation is overcome at sequencing depths of 30–40 million reads, which is within the range of what is typically used and recommended by ENCODE for total RNA-seq experiments. Limitations of the iRNA-seq method include the risk of classifying exon or other types of non-intron reads as intron due to alternative splicing events, etc. not included in current genome annotations. The extent of this is hard to measure; however, given the comparable performance to RNAPII-ChIP-seq and GRO-seq, such cases must be rare. For correctly annotated genes, the iRNA-seq readout may even, in some cases, be more correct than that of the two other methods, the reason being that these techniques, in addition to elongating polymerase, measure engaged, paused polymerase at promoters and enhancers (29). Thus, changes in occupancy of paused polymerase at intron enhancers may be falsely quantified as changes in gene transcription. Another potential limitation of the iRNA-seq method is the fact that intron read density is influenced by processing rates in addition to pre-mRNA synthesis rate. However, recent investigations indicated that variation in processing rates is minor compared to synthesis, leaving intron read densities as an acceptable proxy for rate of transcription (26). Importantly, this is also reflected by the strong correlation between fold changes determined by iRNA-seq and those determined by RNAPII ChIP-seq and GRO-seq. Finally, iRNA-seq is not able to assess transcription for genes with no introns or few short introns. However, despite these limitations, we believe that for most purposes the benefits of this approach greatly outweigh the drawbacks.

The iRNA-seq method provides many technical advantages over other methods for the analysis of acute transcriptional changes (Figure 6). First, the amount of material required for total RNA-seq is much smaller than any of the current technologies for assessment of transcription (12,28,30,31,32,33). iRNA-seq therefore makes it possible to study acute transcriptional regulation in rare cell populations. In theory, since RNA-seq on single cells is becoming a standard protocol, it should also be possible to study transcriptional activity at a single-cell level using the iRNA-seq method (30). Second, the total RNA-seq procedure is much less laborious than the other methods and can be performed with standard technologies. In fact, the RNA-seq procedure is already fully automatized on multiple platforms (http://www.illumina.com/applications/sequencing/ngs-library-prep/automation.ilmn). While some automatized ChIP-seq solutions are also available, several steps in the ChIP-seq procedure (cross-linking, sonication, etc.) have to be optimized for each model system used. Third, the iRNA-seq methodology allows simultaneous assessment of primary transcript levels and mature mRNA levels from the same total RNA-seq experiment. This could be especially useful in time-course experiments, e.g. by allowing the dissection of the temporal relationship between active transcription and accumulations of mature mRNA levels in response to a given stimuli. Finally, iRNA-seq can be used to extract novel information about transcriptional activity from the large number of total RNA-seq data sets deposited at various databases.

In conclusion, we have developed a new method, iRNA-seq, for easy, sensitive and accurate genome-wide assessment of transcriptional activity. Despite minor limitations, this method offers several advantages over current methods in terms of time and material spent. The method will be of particular interest to researchers interested in determining acute changes in transcription, and researchers working with limited sample material, e.g. rare cell populations or even single-cell RNA-seq.

We thank members of the Mandrup laboratory, in particular, Rasmus Siersbæk, for helpful comments and discussions. Furthermore, we thank Per Sauerberg (Novo Nordisk A/S), for the kind gift of rosiglitazone and Martin Wabitsch (University of Ulm, Germany) for providing SGBS cells.

FUNDING

Danish Council for Independent Research|Health and Disease; Novo Nordisk Foundation; European Foundation for the Study of Diabetes; The VILLUM Foundation to the VILLUM Center for Bioanalytical Sciences at University of Southern Denmark. Funding for open access charge: Danish Council for Independent Research|Health and Disease.

Conflict of interest statement. None declared.

REFERENCES