BACKGROUND
the eventual rna transcript level of a gene is determined by regulation at multiple levels, including transcriptional initiation, elongation, splicing, export and degradation.
transcription initiation is regulated by complex interactions of sequence features, many of which involve chromatin modifications  <cit> .
although it is still unclear whether chromatin modifications are the cause or consequence of transcription, these chromatin modifications are often used to infer transcriptional regulation.
the chromatin modifications include several types of histone modifications, such as h3k <dig> tri-methylation , which is often observed around the transcriptional start sites of actively transcribed transcripts  <cit> .
in several large-scale projects they often used h3k4me <dig> sites as markers for active transcription, which allowed the characterization of transcriptionally active regions and estimation of transcript levels in a given cell at a given state  <cit> .
this partly reflects the fact that advances in next generation sequencing have enabled easy characterization of the sites bound by h3k4me <dig> sites using chromatin immunoprecipitation sequencing   <cit> .
indeed, a recent encode study conducted a large number of chip-seq experiments in difference cell types.
there have been several papers that modeled gene expression levels from chromatin features  <cit> .
it is evident that chip-seq data is not sufficient enough to model the steady-state rna expression levels for a number of genes, and regulatory mechanisms other than transcription initiation needs to be considered to understand the rna expression.
rna degradation is regulated by degradation factor, such as upf <dig>  exosc <dig> and stau <dig>  through rna-protein and protein-protein interactions.
upf <dig> is an essential mediator in nonsense-mediated mrna decay , in which aberrant rna containing a premature stop-codon  is recognized and degraded  <cit> .
in addition, recent genome-wide analyses by microarrays and rna-seq have suggested a regulatory role for upf <dig> in targeting 3-20% of bona fide mrna with full coding potential  <cit> .
upf <dig> is involved in other degradation pathways such as the staufen1-mediated mrna decay  and replication-dependent histone mrna decay  <cit> .
it was proposed that approximately 1% of human mrnas are regulated by staufen <dig> , suggesting that smd constitutes a significant post-transcriptional regulatory pathway  <cit> .
exosc <dig> is the essential component of the exosome complex, which functions in 3′ – 5′ rna degradation  <cit> .
however, even for these well-known factors, it is still unclear as to what extent they effect the eventual transcript levels.
in this study, we generated and integrated a dataset of bric-seq  <cit> , rna-seq and chip-seq, in order to uncover the contributions of rna decay in determining eventual genome-wide transcript levels  <cit> .
in bric-seq, the half-lives of transcripts are measured using 5′-bromouridine  based in situ labeling of rna.
bru added in culture medium is incorporated into cells, which convert it to brutp.
it is incorporated into nascent rna during transcription, and consequently, endogenous rnas are labeled with bru.
bru-labeled total rnas are isolated from cells at sequential time points after removal of surplus bru from the culture medium.
bru-labeled rnas are recovered by immunopurification followed by analysis by massive sequencing.
by this method, we can avoid artificial effects of the traditionally used transcriptional inhibitor, such as actinomyicin d, method, in which the physiology of the cell is known to be greatly affected  <cit> .
although 5′-ethynyl uridine labeling and 4′-thiouridine labeling methods have been used for measuring the transcriptome stability, these nucleotide analogues are more toxic than bru.
bru therefore has an advantage to determine rna stabilities in physiologically non-disturbed conditions.
it has been known that the rna abundance does not necessarily correlate with their transcription rates; however the reasoning behind the lack of correlation, have not been well characterized.
here we identified genes that have low rna abundance that could be explained by a particular rna half-life.
in addition, with the aid of computational simulation, we identified genes with rna abundance that was mediated by changes in rna stability by rna decay factors: upf <dig>  exosc <dig>  and stau <dig>
RESULTS
correlation between chip-seq and rna-seq data
first we analyzed the relationship between levels of the transcripts and the strength of active chromatin marks by performing chip-seq  analysis of h3k4me <dig> and pol ii on the illumina hiseq <dig> platform.
chip-seq peaks were called using a representative analytical program, macs  <cit> , using the default parameters .
for the transcript levels, we used rna-seq to determine the genome-wide gene expression in hela cells and we analyzed rna-seq by modeling the gene to the refseq gene model .
we identified a total of  <dig>   and  <dig>  genes that possessed h3k4me <dig> and pol ii “peak”, respectively, within the 3 kb regions around transcriptional start sites , out of  <dig>  refseq genes analyzed .
in the cases where the peaks were not identified, the chip-seq tags were mostly at the noise level; if peaks were identified, we associated the signal intensities of chip-seq of h3k4me <dig> and rna-seq gene expression levels  to the gene.
we observed genes possessing more than 1 ×  <dig> h3k4me <dig> sequence tags and pol ii peaks  and  <dig>   cases, among  <dig>  and  <dig>  genes with pol ii peaks respectively .table  <dig>
statistics of h3k4me <dig> chip-seq peaks against gene expression
total
no peak
low peak
high peak
number of genes
<dig>
<dig>
<dig>
<dig>
genes having rna level of > <dig> rpkm
<dig>
957
<dig>
<dig>
genes having rna level of > <dig> rpkm
<dig>
175
600
<dig>
peaks were called using macs and “peak present” represents genes with h3k4me <dig> and pol ii peaks within  <dig>  kb of the tss .
“no peak” represents genes without any h3k4me <dig> or pol ii peaks within  <dig>  kb of the tss.
we quantitatively analyzed the correlation between chip-seq and rna-seq data for genes with chip-seq “peaks”.
as shown in figure 1a, we observed a positive correlation; mrnas with higher expression levels were associated with higher chip-seq signal intensities, and we observed pearson’s correlation of r =  <dig>   with a log-transformed scatterplot .
when we examined individual genes, it was often possible to observe active transcription with large h3k4me <dig> chip-seq peaks, which is shown in figure 1c.
in contrast, we observed a number of cases in which rna expression was insignificant, despite significant chromatin marks, as shown in figure 1d.
a significant population deviated from the straightforward expected distributions in figure 1b.
when we set the threshold of more than 1 ×  <dig> for h3k4me <dig> chip-seq intensity, presence of polii chip-seq and less than  <dig> rpkm for rna expression, we identified  <dig> , genes /rna : upper left corner in figure 1b) in which significant levels of chip-seq peaks and low levels of rna-seq were detected, as shown in figure 1d.
in addition, for  <dig>  genes /rna : bottom right corner in figure 1b), although chip-seq intensities were less than 1 × 10− <dig> for h3k4me <dig> and no polii chip-seq peaks were detected, rna-seq indicated significant rna levels of more than  <dig> rpkm, as shown in figure 1e.
thus, we identified discrepancies between the chip-seq data and the rna-seq data for a significant population of genes.figure  <dig>
relationship between chromatin marks and eventual transcript levels.
boxplot of the intensities for the chip-seq peaks of h3k4me <dig>  intensities are plotted for the populations with different expression levels measured by rna-seq and are indicated on the x-axis .
asterisks indicate statistically significant differences, as evaluated by wilcoxon’s singed rank test .
scatterplot representing the chip-seq peak signal intensities of h3k4me <dig> on the y-axis and gene expression values on the x-axis .
pearson’s correlation co-efficiencies of the plots  are also shown in the graph area.
dotted lines on the x-axis show  <dig> rpkm and y-axis show 1 ×  <dig> h3k4me <dig> intensities.
labels on each quadrant of the graph  / rna ) are the names given to these set of genes, and are used continually throughout this manuscript.
graphical representation of the patterns of chip-seq  and rna-seq data for the rsp <dig> gene , ptgs <dig> gene  and nes gene .
arrows indicate the direction of transcription and n/a indicates the lack of recognized peaks by macs.
note that while the chip-seq of h3k4me <dig> and pol ii consistently indicate active gene expression of rsp <dig> and ptgs <dig> genes, essentially no gene expression was observed for ptgs <dig> .
for the nes gene, rna-seq tags were observed despite the lack of chip-seq signals of h3k4me <dig> and pol ii.
correlation among half-lives of the transcripts, chromatin marks and transcript levels
to examine the cause of the discrepancy, we focused on mrna stabilities.
we used our unique method, bric-seq, in which the nascent rnas are labeled with 5′-bromouridine  and subjected to massive sequencing analysis in a time-lapse manner.
by calculating the number of bru-labeled rna tags that remain in the population after a particular time duration, bric-seq can be used to measure each rna half-life at a genome-wide level  <cit> .
detailed sequencing statistics for representative cases are shown in additional file 1: figure s <dig>
we examined the relationship between the eventual mrna levels and the half-lives for genes with chip-seq “peaks”, which reflect active transcriptional initiation.
in these cases, we observed positive correlation , in which the half-lives of the transcripts were shorter in proportion to the decreasing expression levels.
however, we detected no correlations between the half-lives of the mrnas with the chip-seq intensities .
these results indicate that rna stability may be a contributing factor for the determination of eventual transcript levels.
furthermore, the mrna stability control is independent of transcriptional initiation, which is inferred by chromatin states.figure  <dig>
relationships between rna half-lives, transcript levels and chromatin marks.
boxplot chart of the rna half-lives of mrnas of varying expression levels.
boxplot charts of the signal intensities for the chip-seq peaks plotted for gene populations with different rna half-lives, as indicated on the x-axis.
boxplot chart of the half-life distributions in different chip-seq and rna-seq fractions.
labels on the x-axis refer to the quadrant from figure 1b.
for chip  / rna , ‘total’ indicates the total transcripts, and ‘×2’ and ‘× <dig> ’ indicates all the genes within 2-fold  and  <dig> -fold  of the least squared regression line, respectively.
asterisks in  indicate statistical significance by wilcoxon’s signed rank test .
scatterplot to show the distributions of chip  / rna  genes.
line is the least squared regression line.
red and green dots indicate genes within  <dig> -fold and 2-fold of the least squared regression line, respectively.
scatterplot representing the chip-seq peak signal intensities of h3k4me <dig> on the y-axis and gene expression values on the x-axis for transcripts with chip-seq peaks and measured half-lives .
genes having “short”  half-lives  are indicated in red dots and were statistically significant .
dotted lines show the  <dig> rpkm and 1 ×  <dig> peak intensities.
based on these observations, we speculated that control of the stability of mrnas might play a pivotal role in determining the eventual rna levels, particularly in case in which the chip-seq and rna-seq data were inconsistent /rna  population in figure 1b); where, the transcript levels may be suppressed at a low level, despite their active transcription, owing to fast rna turnover rates.
to examine this possibility, we compared the half-lives of mrnas between gene groups having h3k4me <dig> chip-seq intensities of more than 1 ×  <dig>  a presence of polii chip-seq peak, with gene expression levels below  <dig> rpkm /rna ) and those with h3k4me <dig> chip-seq intensities higher than 1 ×  <dig>  with a presence of polii chip-seq peak and gene expression values above  <dig> rpkm /rna ) or those with h3k4me <dig> chip-seq intensities lower than 1 ×  <dig>  with a lack of polii chip-seq peak and gene expression values above  <dig> rpkm /rna ), genes that half-lives could be measured.
as shown in figure 2c, we found that half-lives of the transcripts for chip /rna  genes were significantly shorter than those of chip /rna  genes  and chip /rna  genes .
from the correlations between rna-seq and chip-seq as shown in figure 1b, it was possible to infer the gene expression levels from the chip-seq intensities for those genes in which the rna-seq and chip-seq were consistent.
when we examined genes for which gene expression values were within a  <dig> or  <dig>  fold difference from those expected from the chip-seq intensities , we observed narrower distribution of the mrna half-lives with a median value of  <dig>  and  <dig>  hours, respectively .
these half-lives of chip /rna  genes may serve as the default half-lives of those genes, and if genes do exhibit this particular mrna half-life, the transcriptional initiation levels should be the major determinants of the eventual transcript levels.table  <dig>
statistics of the half-life associated with chip  and chip  genes
total
chip rna
chip rna
chip rna
chip rna
chip rna
chip rna
number of genes with a rna half-life
<dig>
<dig>
603
<dig>
<dig>
<dig>
187
median half-life
nd
nd
<dig>
<dig>
<dig>
<dig>
<dig>
in table  <dig>  we only considered genes with compatible half-life measurements from bric-seq.
chip: h3k4me <dig> intensities larger than 1 ×  <dig> and a presence of polii peak , chip: h3k4me <dig> intensities smaller than 1 ×  <dig> and an absence of polii peak.
rna: rpkm value larger than  <dig>  rna: rpkm value smaller than  <dig>  for chip/rna region; total: all genes, ×2: genes within two-fold of the regression line, × <dig> : genes within the  <dig> -fold of the regression line.
nd: not determined.
for further analysis into the correlation, we also selected mrna with ‘short’ half-lives, specifically a total of  <dig>  genes that had half-lives shorter than 4 hours.
from the default half-life of  <dig>  hours, we observed the standard deviation of the half-lives to be 3 hours, indicating that 4 hours was approximately the 95th percentile confidence level.
these particular mrnas are highlighted as red dots in figure 2e.
we observed the enrichment of mrnas with short half-lives in the chip+/rna- fraction of the scatterplot , in which gene expression values were lower than estimated by chip-seq data.
these results suggest that control of mrna stability is an important factor in determining the eventual mrna levels for this population and indicate that the expression levels of these  <dig> genes may be controlled by mrna stability in hela cells .
for validation, we used actinomycin d , a transcriptional inhibitor, and chased the rna decay by rt-qpcr .
bambi and med <dig> rna half-lives /rna/short rna half-life) determined by bric analysis were similar to those determined by actd analysis.
the transcripts of mmp <dig> and slc25a <dig> /rna/long rna half-life) were determined as stable rnas in both bric and actd chase analyses.
in contrast, the rna half-lives of znf <dig> and znf <dig> slightly varied between bric and actd chase analyses.
thus, most rna stabilities determined by bric analysis were confirmed by actd chase analysis.
to analyze the categories of genes that receive regulation at either transcriptional initiation or rna half-life levels, we ran go enrichment analysis.
among the chip /rna  genes, go terms associated with basic translation or transcriptional machineries were enriched .
for the chip /rna  genes, genes associated cytoplasm as a location were enriched , and for the chip /rna  genes, go terms associated with transcription factors were enriched, particularly among genes with a short half-life  .
these results suggest that different functional categories of genes are subjected to different modes of gene expression regulation.table  <dig>
list of go enrichments for transcripts in different chip-seq and rna-seq fractions
a) chip /rna  go enrichment
go:id
go: term
number of genes
false-discovery rate
go: 0044822
poly  rna binding
618
<dig> e-261
go: 0006412
translation
293
<dig> e-196
go: 0010467
gene expression
387
<dig> e-162
b) chip /rna  go enrichment
go: id
go: term
number of genes
false-discovery rate
go: 0005737
cytoplasm
138
<dig> e-2
go: 0070062
extracellular vesicular exosome
74
<dig> e-2
go: 0005635
nuclear envelope
15
<dig> e-2
c) chip /rna  go enrichment
go: id
go: term
number of genes
false-discovery rate
go: 0003677
dna binding
378
<dig> e-28
go: 0006355
regulation of transcription, dna-templated
266
<dig> e-17
go: 0006351
transcription, dna-templated
361
<dig> e-17
d) chip /rna /half-life < 4 h go enrichment
go: id
go: term
number of genes
false-discovery rate
go: 0003677
dna binding
216
<dig> e-58
go: 0006351
transcription, dna-templated
192
<dig> e-37
go: 0006355
regulation of transcription, dna-templated
147
<dig> e-34
identification of upf <dig>  exosc <dig> and stau <dig> as controlling factors for rna stabilities
to evaluate the potential contribution of known rna degradation factors to the control of global rna stability, we chose three representative factors for analysis: upf <dig>  exosc <dig>  and stau <dig>  it has been reported that upf <dig> regulates 3–20% of transcripts  <cit> , highlighting the potential importance of upf <dig> in regulating rna degradation and abundance.
exosc <dig> is an essential component of the exosome complex that is the major mrna degradation machinery in mammalian cells.
to analyze the alteration of global mrna turnover by perturbation of representative factors, we examined exosc <dig>  stau <dig>  which regulates around 1% of bona fide mrnas  <cit> , is a typical rna-binding protein involved in rna degradation.
we used data from bric assay in the cells depleted in upf <dig> by sirna .
as observed in a previous study  <cit> , mrna levels of the gadd45a gene, which is a known target of upf <dig>  were increased following upf <dig> knockdown , with an increase in half-life .
we then conducted similar experiments using exosc <dig> and stau <dig> knockdown cells and prepared a similar rna-seq and bric-seq dataset.
we observed an increase of fam120c mrna levels  and an increase in half-life  in exosc <dig> knockdown cells.
in stau <dig> knockdown cells, the mrna levels of cdkn2b were increased , with increased half-lives .figure  <dig>
identification of upf <dig>  exosc <dig> and stau <dig> as factors controlling rna half-lives.
chip-seq and rna-seq for the gadd45a , fam120c  and cdkn2b  genes, which showed significant changes in rna half-lives in upf <dig>  exosc <dig> and stau <dig> knocked-down cells, respectively.
upper two panels in : chip-seq data for h3k4me <dig> and pol ii.
lower two panels in : rna-seq data for control cells and the respective knockdown cells.
bric-seq normalized graph for gadd45a , fam120c  and cdkn2b  genes in control cells or cells knocked-down for upf <dig> , exosc <dig> , and stau <dig> .
scatterplot representing the chip-seq peak signal intensities of h3k4me <dig> on the y-axis and gene expression values on the x-axis with putative upf1-controlled genes, exosc5-controlled and stau1-controlled genes in blue, green and red dots respectively.
we next looked for genes that may be regulated by these factors, and identified  <dig>   <dig> and  <dig> genes where the mrna half-lives were extended by more than two-fold  or  <dig> -fold  and showed mrna expression increase by two-fold  or  <dig> -fold  in upf <dig>  exosc <dig> and stau <dig> knockdown cells, respectively.
as shown in figures 3a, 3b and 3c, we noticed that the transcripts that were not observed in the control knockdown cells appeared in the knocked-down cells in many cases.
we examined the distribution of the dots of these genes, whose transcripts were stabilized and increase in corresponding knockdown cells, in figure 1b and found that they are enriched in the upper-left corner of the plot.
in total, we identified  <dig>   <dig> and  <dig> genes  whose mrna half-lives are potentially controlled by upf <dig>  exosc <dig> and stau <dig> respectively, consisting of  <dig>   <dig>  and  <dig> % of the total of  <dig>  genes , rna , half-life < 4 h) in this area.
we also examined the overlap between the genes controlled by upf <dig> and exosc <dig> and found little overlap .
although we could identify upf <dig>  exosc <dig> and stau <dig> as control factors only for a limited population  by this approach.
these observations are only the first step in identifying the role of rna decay factors in determination of rna abundance through rna degradation, and further systematic analyses may facilitate identification of the complex regulatory mechanisms of mrna stabilities.table  <dig>
summary statistics used for the analysis of refseq transcripts
conditions
number of genes
active transcription  and half-life measured
<dig>
with short half-life
<dig>
with low expression  and high h3k4me <dig>
866
and upf <dig> target
<dig>
or exosc <dig> target
<dig>
or stau <dig> target
<dig>
computational modeling of the effect of mrna half-lives on eventual mrna levels
we predicted the rna abundance by normalizing the half-life to be  <dig>  hours, the estimated default half-life , and obtained the least squares regression line between the predicted rna levels and h3k4me <dig> intensities.
we made a threshold of ×  <dig>  and ×  <dig>  above and below the least-square regression line to define genes where the h3k4me <dig> intensities and rna abundance correlate.
we found that out of  <dig>  genes that were available from bric-seq dataset, we found  <dig>  and  <dig> genes that resided within ×  <dig> and ×  <dig>  of the regression line, respectively .
we checked the original gene expression of those genes prior to the simulation, and checked whether their measured rna abundance correlates with h3k4me <dig> intensities.
we found out of  <dig>  and  <dig> genes that resided within ×  <dig> and ×  <dig>  of the regression line from predicted gene expression, we found  <dig>  and  <dig> genes where measured to be outside of the threshold, respectively.
it means that the rna stability of those genes contributed to the rna abundance.
taken together, these results collectively support our claim that rna degradation significantly contributes in determining the eventual expression levels .figure  <dig>
computational simulations of the rna half-lives as a determinant for transcript level.
computational simulation into the effect of rna half-life on rna abundance.
each gene was simulated to have a half-life of  <dig>  hours, potential default half-life, and the predicted gene expression was plotted against the h3k4me <dig> intensities, along with the measured gene expression.
the dotted line indicates the least-squares regression line between the predicted gene expression and h3k4me <dig> intensities.
result of the computational simulation of the predicted rna expression when each rna abundance was predicted from rna half-life when upf <dig>   exosc <dig> and  stau <dig> were knocked-down, in comparison to each control knock-down conditions.
the predicted gene-expression of the decay factor knock-down and the measured gene-expression of control knock-down were plotted against h3k4me <dig> intensities.
the dotted line indicates the least-squares regression line between the predicted gene expression and h3k4me <dig> intensities.
summary statistics of computational simulation
condition
t1/2= <dig>  hours
siupf1
siexosc5
sistau1
additionally, we conducted chip-seq on h3k27ac, h3k27me <dig>  h3k36me <dig> to build a linear model as described by wang, c., et.
al.
<cit> ,.
we built one linear model incorporating h3k4me <dig>  h3k27ac, h3k27me <dig> and h3k36me <dig> intensities with half-life as an extra variable, and one without the half-life, to explain the rna abundance.
we found that fitting increased from  <dig>  to  <dig>  in r-value, which confirmed the previous finding by wang.
et.
al., <cit> .
we examined whether the changes in the rna half-lives from knockdown of upf <dig>  exosc <dig> or stau <dig> could explain the changes in the eventual transcript levels.
for this analysis, we conducted a computational simulation.
as shown in additional file 1: figure s <dig>  we found that the described theoretical model can predict the changes of eventual rna levels with pearson’s correlation co-efficiency of  <dig> ,  <dig>  and  <dig> , respectively.
overall, we demonstrate that the simple computational model could reasonably explain the changes of eventual mrna levels, thus supporting our idea that the major determinant of the eventual rna levels in these cases is at the level of rna stability.
we further simulated the rna abundance, from the changes in rna half-life, in relation to the chip-seq signal levels and we found  <dig>   <dig>  and  <dig> genes that were within ×  <dig>  of the regression line that lied outside of the threshold prior to simulation, for upf <dig>  exosc <dig>  and stau <dig>  respectively .
possible feedback between mrna turnover and transcriptional initiation
to analyze the possible feedback mechanisms between mrna turnover and transcription initiation, we used the genes in which both mrna half-lives and eventual transcript levels were increased more than two fold in knockdown cells.
in  <dig> and  <dig>  genes in upf <dig> and exosc <dig> knockdown cells, respectively, there were no significant changes in eventual transcript levels  observed despite remarkable changes to their rna half-lives .
we speculated that there might be a possible feedback between mrna turnover and transcriptional rate.
if mrnas for a particular group of transcriptional repression factors are included in the upf1/exosc <dig> targets and their stabilized mrnas result in increased protein levels of such transcriptional repression factors, thereby enhancing the repression activities on their target genes, it would explain unchanged balance of eventual transcript levels for these genes.
we examined whether any transcription factor binding sites were enriched in the upstream regions of the  <dig> and  <dig>  genes.
in the case of upf <dig>  we detected significant enrichment of the consensus binding site for hic <dig>  which is a transcription factor belonging to the zinc finger family.
we validated the changes of its expression level and half-lives in control and upf <dig> knockdown cells, and confirmed hic <dig> as a upf <dig> target .
we were unable to further validate direct binding of hif <dig> to target genes, since no effective antibodies are available.
also, it is possible that hic <dig> may not be the only candidate, which may contribute to the feedback regulation.
many zinc finger family transcription factors share consensus binding sequences.
the list of putative zinc finger family transcription factors that have significant homology to hic <dig> in their dna binding domains with extended half-lives  and increased eventual transcript levels upon upf <dig> knockdown are shown in additional file 2: table s <dig>  these factors may collectively enable elaborate regulation of gene expression.
identification of candidate genes controlled through rna stability in other cell types
to further extend our idea that controls at the level of mrna decay contribute to determining eventual mrna expression levels in other cell types, we analyzed the published encode data  <cit>  and dbtss  <cit> , which included chip-seq data of h3k4me <dig> and pol ii and rna-seq in a wide variety of cell types.
in addition, the rna-seq data of subcellular fractionated mrnas were included in the dataset.
we selected eight cell types for which all these datasets were available .
we retrieved and analyzed the chip-seq data and rna-seq data as performed with hela cells.
first, as shown in additional file 1: figure s <dig>  we observed weak correlations between chip-seq intensities of h3k4me <dig> and gene expression values as reported in hela cells.
since there were no bric-seq data in the encode dataset, we could not directly analyze the genes with short mrna half-lives.
nevertheless, we could select genes for which chip-seq tags of h3k4me <dig> and pol ii were associated, thus indicated as actively transcribed in the cell line, although their gene expression levels were not at the expected levels.
in the encode dataset, we also considered a positive signal of h3k36me <dig>  which is a chromatin mark for transcriptional elongation, to further assure active transcription  <cit> .
we identified an average of  <dig> candidate genes with active transcription and low rna abundance in each cell type , which should be regulated at the level of mrna decay.
we examined and identified go terms that were significantly enriched, depending on cell types, and the top enriched terms were predominantly associated with dna binding .
we also analyzed whether there are any cases for a particular gene to be selected as such a candidate in a cell type-preferred manner.
we found that  <dig>  potential controls at the level of rna half-life were observed in a single cell type .
we also examined if there was possible feedback between controls of rna half-lives and transcription initiations.
we found that several transcription factor binding consensus sequences of the ets family and arebp family genes are enriched in the promoters of the genes inferred to be regulated at the level of rna stability in gm <dig> cells and hepg <dig> cells, respectively.
in addition, we identified hic <dig> binding consensus sequences in the promoters of the genes inferred to be regulated at the level of rna stability in human h <dig> embryonic stem cells.
the feedback regulations may be common in various cell types with distinct responsive transcription factors depending on cell type .
we further analyzed the subcellular localizations of the mrnas of these genes using the corresponding encode data.
we found significantly enriched mrnas in the cytoplasm only in the embryonic stem cells   , suggesting there may be a characteristic regulatory mechanism for controlling rna stability in the cytoplasm of es cells.
on the whole, these data should provide an important complement to the encode annotations, which aim to generate a complete catalogue of genetic elements explaining gene expression regulation.figure  <dig>
screening of candidate genes that may be controlled at the level of rna degradation from the encode dataset.
number of candidate genes screened from the indicated cell types based on the encode data.
gene ontology  enrichment of genes in the chip  / rna  regions in the indicated cell types based on the encode data.
chip-seq signal intensities were comparable between different cell types, although rna-seq-based gene expression values were remarkably different.
the go term enrichment with the lowest p-value for each cell line is shown.
enriched consensus transcription factor binding sites for genes in the chip  / rna  region for encode dataset.
graphical representation of the enriched consensus binding sites in the promoter regions of the genes in chip  / rna  region for encode dataset.
the list of cell lines with consensus binding site enrichment from the transfac database.
the nuclear/cytoplasmic gene expression values for the genes where rna half-life may be the contributor to the rna levels and all other genes in h <dig> human embryonic stem cells .
statistical significance of the difference is indicated under the x-axis.
distinct controls of the rna stabilities of mrnas and non-coding rnas
to examine whether the regulations at the level of rna half-lives are observed in lincrnas, we conducted a similar analysis for lincrnas in hela cells, as shown in figures  <dig> and  <dig>  we tentatively defined the dataset of lincrna as that of lincdb  <cit> .
as shown in figures 6a and 6b, we associated the chip-seq peak intensities and gene expression values.
we unexpectedly observed distinct patterns from those of mrnas.
namely, among  <dig> lincrnas in hela cells,  <dig> lincrnas had “short ” rna half-lives.
of these,  <dig>  resided in the upper-left corner of the plot, suggesting that most of the lincrnas are controlled at the post-transcriptional level.
because there are reports on the possible involvement of nmd in regulating non-coding rna  <cit> , we examined the possible involvement of upf <dig> in the regulation of lincrnas.
interestingly, none of the  <dig> lincrnas were detected as potential upf1-controlled transcripts.
in this data set, upf <dig> may have only limited contribution.
by contrast, we found  <dig>  of the lincrnas  were regulated by exosc <dig> figure  <dig>
contribution of the rna half-lives to transcriptional regulations of lincrnas.
scatterplot showing the relationships between h3k4me <dig> intensities and gene expression for all lincrna  and those that were bric-seq compatible  .
examples of lincrnas regulated by upf <dig>  and exosc <dig> .
upper panels are h3k4me <dig> and pol ii peaks, and lower panels are rna-seq measurements for the control and indicated knockdown.
upf <dig>  and exosc <dig>  show the normalized decay curve from the bric-seq measurements for these examples.
red lines are siupf <dig>  or siexosc <dig>  data, respectively and blue lines are si_control data.
enhancer rnas  regulated by upf <dig>  or exosc <dig> .
the label shows the genomic coordinates of these ernas .
these show normalized decay curves for the bric-seq measurements for these examples as above.
summary statistics used for the analysis of lincrnas
conditions
number of genes
active transcription and half-life measured
141
with short half-life
103
with low expression  and high h3k4me <dig>
84
and upf <dig> target
0
or exosc <dig> target
<dig>
we next examined whether enhancer rnas , which facilitate the functions of the enhancers, are also regulated at the level of rna degradation.
we again used the encode data from hela cells.
we retrieved the chip-seq data of h3k4me and h3k27ac, which are representative chromatin marks of active enhancers.
among  <dig>  genomic regions having “peaks” of both h3k4me <dig> and h3k27ac, we identified  <dig> cases in which there were rnas in the overlapping regions ,  <dig>  kb away from any refseq gene body and their rna half-lives were extended by upf <dig> knockdown by more than two fold in bric-seq assay.
similarly, we identified  <dig> cases in which half-lives of rnas in overlapping regions were examined  and in which transcripts showed extended rna half-lives of more than two-fold in exosc <dig> knockdown cells.
although further detailed experimental validations are necessary, these results may indicate that controls mediated by rna stability are used in determining the transcript levels of non-coding rnas.
discussion
here, we have described a genome-wide correlation among the signal intensities of chip-seq, gene expression values measured by rna half-lives measured by bric-seq.
we identified that regulation at the level of rna degradation plays an important role in determining eventual rna levels.
we demonstrated that this control may exhibit a particularly large effect in cases in which the chip-seq data and rna-seq data are inconsistent.
indeed, in this study, we estimated that the abundance of  <dig> mrnas is regulated by rna degradation in hela cells.
furthermore, we applied similar approaches to analyze the public encode data and identified a total of  <dig>  candidate genes whose gene expression levels are likely to be controlled at the level of rna stability.
we also found that these controls appeared to vary among cell types.
to our knowledge, this is the first report that describes the integration of the chip-seq, rna-seq and rna half-life data in identifying genes that may receive post-transcriptional gene expression regulation.
in go analyses, there are some limitations because go terms may be loosely defined for a particular gene.
however, we first selected go terms with statistical significance , and we only used go term enrichment with more than  <dig> genes in a set to ensure minimum false positives, which will ensure that the go terms described are representative of the biological phenomenon.
we then found enrichment of transcription factors in transcripts with discrepancies between the chip-seq and rna-seq with short half-lives, and in particular, negative regulatory factors.
in a recent study, haimovich et al.
<cit>  indicated that some rna degradation factors play a role in transcription, implicating a feedback loop for gene expression.
our data suggest another mechanism by which this occurs, through a faster degradation rate of mrnas encoding transcription factors, particularly for those that negatively regulate gene expression, thus affecting the eventual rna levels.
previous studies conducted by wang et.
al., demonstrated that rna stability could in inferred from the residual errors in modeling the rna abundance from chip-seq data and they validated their claims from a half-life data; however, they used rna stability data from a different cell-line  <cit> , and they could not accurately predict the eventual rna abundance from the rna stability.
herein, by analyzing rna stability data from the same cell-line as rna-seq data, we were able to estimate the contribution of the rna stability to the rna abundance.
in addition we were able to estimate the rna stability contributions on the rna abundance upon upf <dig>  stau <dig> and exosc <dig> knockdown.
although we did not demonstrate how the rna half-lives are controlled in the current study and we identified upf <dig>  exosc <dig> and stau <dig> as control factors in some cases, they could explain at most 8% of the total mrna population.
even for the cases of candidate upf <dig>  exosc <dig> or stau1-controlled genes, it is possible that they may not be direct targets of these factors and that we may have picked up secondary or later effects as a consequence of upf <dig>  exosc <dig> or stau <dig> knockdowns.
another obvious drawback of our approach is that mrna half-lives were not directly measured by bric-seq for the encode dataset.
therefore, it is possible that they may be mediated by other regulatory mechanisms, rather than at the level of rna half-lives, such as rna halting and abortive transcriptions.
to minimize these possibilities, we selected the cases in which h3k36me <dig>  a marker of transcriptional elongation, should be significant in the transcript regions.
in spite of several drawbacks, we believe that genome-wide features of correlation among chip-seq, rna-seq and bric-seq should give an important starting point to further explore posttranscriptional regulatory mechanisms, for which only limited knowledge has been accumulated.
indeed, recent papers have begun to reveal many human diseases that are caused by malfunctions of rna decay pathway.
in particular, it has been made gradually clear that most immune-response mrnas are destabilized when they are not required via their cis-regulatory elements in the 3′ utr.
it is proposed that such rna-decay mechanisms collectively enable rapid up-/downregulation of gene expression in response to environmental changes.
the au-rich element  is one of such elements widely found in the 3′ utr of mrna of immune-related genes.
mice lacking are in the tnf-alpha mrna showed joint and gut-associated immunopathologies  <cit> .
the trans-acting factor regulatory rnase  <dig> , which is induced by toll-like receptor  ligands, interleukin -1β and mcp- <dig>  is involved in the destabilization of mrnas including il <dig> mrna.
regnase-1-deficient mice develop severe autoimmune disease because of excess production of cytokines  <cit> , highlighting the importance of rna degradation-mediated gene regulation.
in addition to immunological disorders, there are a growing number of cases with impaired rna decay regulation that cause disease, and they sometimes reveal unexpected connections between otherwise completely unrelated diseases.
perlman syndrome, an autosomal recessively inherited congenital overgrowth syndrome associated with high neonatal mortality, is an obvious example.
the survivors of this disease have a high risk of wilms tumor.
recently, it was reported that the responsible gene of this disease is the exoribonuclease dis3l <dig>  a homologue of exosome component dis <dig>  <cit> .
moreover, dis3l <dig> is mutated in approximately 3–6% of carcinomas  <cit> .
in this study, we have also analyzed the stabilities of non-coding rnas  and mrnas.
we and another study reported that stability of non-coding rnas is also tightly regulated, suggesting that the instability contributes to the dynamic nature of lincrnas  <cit> .
indeed, the stability of noncoding rnas has an impact on their biological function  <cit> , although the exploration into its relevance in human disease has just begun.
further enrichment of our knowledge on the control mechanisms on rna stability both for mrnas and non-coding rnas will shed new light on putative disease-associated genetic or somatic mutations.
CONCLUSIONS
by integrative analysis of chip-seq, rna-seq and our bric-seq, we showed that rna half-life may serve as an important post-transcriptional determinant of gene expression.
we suggest that upf <dig>  exosc <dig> and stau <dig> may play active roles in such controls.
in addition, we propose the linkage between transcription and rna decay through regulated degradation of mrnas encoding transcription factors to maintain the steady state level of rna abundance.

 the eventual rna transcript level of a gene is determined by regulation at multiple levels, including transcriptional initiation, elongation, splicing, export and degradation. transcription initiation is regulated by complex interactions of sequence features, many of which involve chromatin modifications. although it is still unclear whether chromatin modifications are the cause or consequence of transcription, these chromatin modifications are often used to infer transcriptional regulation. the chromatin modifications include several types of histone modifications, such as h3k tri-methylation, which is often observed around the transcriptional start sites of actively transcribed transcripts. in several large-scale projects they often used h3k4me sites as markers for active transcription, which allowed the characterization of transcriptionally active regions and estimation of transcript levels in a given cell at a given state. this partly reflects the fact that advances in next generation sequencing have enabled easy characterization of the sites bound by h3k4me sites using chromatin immunoprecipitation sequencing. indeed, a recent encode study conducted a large number of chip-seq experiments in difference cell types. there have been several papers that modeled gene expression levels from chromatin features. it is evident that chip-seq data is not sufficient enough to model the steady-state rna expression levels for a number of genes, and regulatory mechanisms other than transcription initiation needs to be considered to understand the rna expression. rna degradation is regulated by degradation factor, such as upf exosc and stau through rna-protein and protein-protein interactions. upf is an essential mediator in nonsense-mediated mrna decay, in which aberrant rna containing a premature stop-codon is recognized and degraded. in addition, recent genome-wide analyses by microarrays and rna-seq have suggested a regulatory role for upf in targeting 3-20% of bona fide mrna with full coding potential. upf is involved in other degradation pathways such as the staufen1-mediated mrna decay and replication-dependent histone mrna decay. it was proposed that approximately 1% of human mrnas are regulated by staufen, suggesting that smd constitutes a significant post-transcriptional regulatory pathway. exosc is the essential component of the exosome complex, which functions in 3′ – 5′ rna degradation. however, even for these well-known factors, it is still unclear as to what extent they effect the eventual transcript levels. in this study, we generated and integrated a dataset of bric-seq, rna-seq and chip-seq, in order to uncover the contributions of rna decay in determining eventual genome-wide transcript levels. in bric-seq, the half-lives of transcripts are measured using 5′-bromouridine based in situ labeling of rna. bru added in culture medium is incorporated into cells, which convert it to brutp. it is incorporated into nascent rna during transcription, and consequently, endogenous rnas are labeled with bru. bru-labeled total rnas are isolated from cells at sequential time points after removal of surplus bru from the culture medium. bru-labeled rnas are recovered by immunopurification followed by analysis by massive sequencing. by this method, we can avoid artificial effects of the traditionally used transcriptional inhibitor, such as actinomyicin d, method, in which the physiology of the cell is known to be greatly affected. although 5′-ethynyl uridine labeling and 4′-thiouridine labeling methods have been used for measuring the transcriptome stability, these nucleotide analogues are more toxic than bru. bru therefore has an advantage to determine rna stabilities in physiologically non-disturbed conditions. it has been known that the rna abundance does not necessarily correlate with their transcription rates; however the reasoning behind the lack of correlation, have not been well characterized. here we identified genes that have low rna abundance that could be explained by a particular rna half-life. in addition, with the aid of computational simulation, we identified genes with rna abundance that was mediated by changes in rna stability by rna decay factors: upf exosc and stau  correlation between chip-seq and rna-seq data first we analyzed the relationship between levels of the transcripts and the strength of active chromatin marks by performing chip-seq analysis of h3k4me and pol ii on the illumina hiseq platform. chip-seq peaks were called using a representative analytical program, macs, using the default parameters. for the transcript levels, we used rna-seq to determine the genome-wide gene expression in hela cells and we analyzed rna-seq by modeling the gene to the refseq gene model. we identified a total of and genes that possessed h3k4me and pol ii “peak”, respectively, within the 3 kb regions around transcriptional start sites, out of refseq genes analyzed. in the cases where the peaks were not identified, the chip-seq tags were mostly at the noise level; if peaks were identified, we associated the signal intensities of chip-seq of h3k4me and rna-seq gene expression levels to the gene. we observed genes possessing more than 1 ×  h3k4me sequence tags and pol ii peaks and cases, among and genes with pol ii peaks respectively.table statistics of h3k4me chip-seq peaks against gene expression total no peak low peak high peak number of genes genes having rna level of > rpkm 957 genes having rna level of > rpkm 175 600 peaks were called using macs and “peak present” represents genes with h3k4me and pol ii peaks within  kb of the tss. “no peak” represents genes without any h3k4me or pol ii peaks within  kb of the tss. we quantitatively analyzed the correlation between chip-seq and rna-seq data for genes with chip-seq “peaks”. as shown in figure 1a, we observed a positive correlation; mrnas with higher expression levels were associated with higher chip-seq signal intensities, and we observed pearson’s correlation of r =  with a log-transformed scatterplot. when we examined individual genes, it was often possible to observe active transcription with large h3k4me chip-seq peaks, which is shown in figure 1c. in contrast, we observed a number of cases in which rna expression was insignificant, despite significant chromatin marks, as shown in figure 1d. a significant population deviated from the straightforward expected distributions in figure 1b. when we set the threshold of more than 1 ×  for h3k4me chip-seq intensity, presence of polii chip-seq and less than rpkm for rna expression, we identified, genes /rna: upper left corner in figure 1b) in which significant levels of chip-seq peaks and low levels of rna-seq were detected, as shown in figure 1d. in addition, for genes /rna: bottom right corner in figure 1b), although chip-seq intensities were less than 1 × 10− for h3k4me and no polii chip-seq peaks were detected, rna-seq indicated significant rna levels of more than rpkm, as shown in figure 1e. thus, we identified discrepancies between the chip-seq data and the rna-seq data for a significant population of genes.figure relationship between chromatin marks and eventual transcript levels. boxplot of the intensities for the chip-seq peaks of h3k4me intensities are plotted for the populations with different expression levels measured by rna-seq and are indicated on the x-axis. asterisks indicate statistically significant differences, as evaluated by wilcoxon’s singed rank test. scatterplot representing the chip-seq peak signal intensities of h3k4me on the y-axis and gene expression values on the x-axis. pearson’s correlation co-efficiencies of the plots are also shown in the graph area. dotted lines on the x-axis show rpkm and y-axis show 1 ×  h3k4me intensities. labels on each quadrant of the graph / rna ) are the names given to these set of genes, and are used continually throughout this manuscript. graphical representation of the patterns of chip-seq and rna-seq data for the rsp gene, ptgs gene and nes gene. arrows indicate the direction of transcription and n/a indicates the lack of recognized peaks by macs. note that while the chip-seq of h3k4me and pol ii consistently indicate active gene expression of rsp and ptgs genes, essentially no gene expression was observed for ptgs. for the nes gene, rna-seq tags were observed despite the lack of chip-seq signals of h3k4me and pol ii. correlation among half-lives of the transcripts, chromatin marks and transcript levels to examine the cause of the discrepancy, we focused on mrna stabilities. we used our unique method, bric-seq, in which the nascent rnas are labeled with 5′-bromouridine and subjected to massive sequencing analysis in a time-lapse manner. by calculating the number of bru-labeled rna tags that remain in the population after a particular time duration, bric-seq can be used to measure each rna half-life at a genome-wide level. detailed sequencing statistics for representative cases are shown in additional file 1: figure s we examined the relationship between the eventual mrna levels and the half-lives for genes with chip-seq “peaks”, which reflect active transcriptional initiation. in these cases, we observed positive correlation, in which the half-lives of the transcripts were shorter in proportion to the decreasing expression levels. however, we detected no correlations between the half-lives of the mrnas with the chip-seq intensities. these results indicate that rna stability may be a contributing factor for the determination of eventual transcript levels. furthermore, the mrna stability control is independent of transcriptional initiation, which is inferred by chromatin states.figure relationships between rna half-lives, transcript levels and chromatin marks. boxplot chart of the rna half-lives of mrnas of varying expression levels. boxplot charts of the signal intensities for the chip-seq peaks plotted for gene populations with different rna half-lives, as indicated on the x-axis. boxplot chart of the half-life distributions in different chip-seq and rna-seq fractions. labels on the x-axis refer to the quadrant from figure 1b. for chip / rna, ‘total’ indicates the total transcripts, and ‘×2’ and ‘× ’ indicates all the genes within 2-fold and -fold of the least squared regression line, respectively. asterisks in indicate statistical significance by wilcoxon’s signed rank test. scatterplot to show the distributions of chip / rna genes. line is the least squared regression line. red and green dots indicate genes within -fold and 2-fold of the least squared regression line, respectively. scatterplot representing the chip-seq peak signal intensities of h3k4me on the y-axis and gene expression values on the x-axis for transcripts with chip-seq peaks and measured half-lives. genes having “short” half-lives are indicated in red dots and were statistically significant. dotted lines show the rpkm and 1 ×  peak intensities. based on these observations, we speculated that control of the stability of mrnas might play a pivotal role in determining the eventual rna levels, particularly in case in which the chip-seq and rna-seq data were inconsistent /rna population in figure 1b); where, the transcript levels may be suppressed at a low level, despite their active transcription, owing to fast rna turnover rates. to examine this possibility, we compared the half-lives of mrnas between gene groups having h3k4me chip-seq intensities of more than 1 ×  a presence of polii chip-seq peak, with gene expression levels below rpkm /rna ) and those with h3k4me chip-seq intensities higher than 1 ×  with a presence of polii chip-seq peak and gene expression values above rpkm /rna ) or those with h3k4me chip-seq intensities lower than 1 ×  with a lack of polii chip-seq peak and gene expression values above rpkm /rna ), genes that half-lives could be measured. as shown in figure 2c, we found that half-lives of the transcripts for chip /rna genes were significantly shorter than those of chip /rna genes and chip /rna genes. from the correlations between rna-seq and chip-seq as shown in figure 1b, it was possible to infer the gene expression levels from the chip-seq intensities for those genes in which the rna-seq and chip-seq were consistent. when we examined genes for which gene expression values were within a or fold difference from those expected from the chip-seq intensities, we observed narrower distribution of the mrna half-lives with a median value of and  hours, respectively. these half-lives of chip /rna genes may serve as the default half-lives of those genes, and if genes do exhibit this particular mrna half-life, the transcriptional initiation levels should be the major determinants of the eventual transcript levels.table statistics of the half-life associated with chip and chip genes total chip rna chip rna chip rna chip rna chip rna chip rna number of genes with a rna half-life 603 187 median half-life nd nd in table  we only considered genes with compatible half-life measurements from bric-seq. chip: h3k4me intensities larger than 1 ×  and a presence of polii peak, chip: h3k4me intensities smaller than 1 ×  and an absence of polii peak. rna: rpkm value larger than rna: rpkm value smaller than for chip/rna region; total: all genes, ×2: genes within two-fold of the regression line, ×: genes within the -fold of the regression line. nd: not determined. for further analysis into the correlation, we also selected mrna with ‘short’ half-lives, specifically a total of genes that had half-lives shorter than 4 hours. from the default half-life of  hours, we observed the standard deviation of the half-lives to be 3 hours, indicating that 4 hours was approximately the 95th percentile confidence level. these particular mrnas are highlighted as red dots in figure 2e. we observed the enrichment of mrnas with short half-lives in the chip+/rna- fraction of the scatterplot, in which gene expression values were lower than estimated by chip-seq data. these results suggest that control of mrna stability is an important factor in determining the eventual mrna levels for this population and indicate that the expression levels of these genes may be controlled by mrna stability in hela cells. for validation, we used actinomycin d, a transcriptional inhibitor, and chased the rna decay by rt-qpcr. bambi and med rna half-lives /rna/short rna half-life) determined by bric analysis were similar to those determined by actd analysis. the transcripts of mmp and slc25a /rna/long rna half-life) were determined as stable rnas in both bric and actd chase analyses. in contrast, the rna half-lives of znf and znf slightly varied between bric and actd chase analyses. thus, most rna stabilities determined by bric analysis were confirmed by actd chase analysis. to analyze the categories of genes that receive regulation at either transcriptional initiation or rna half-life levels, we ran go enrichment analysis. among the chip /rna genes, go terms associated with basic translation or transcriptional machineries were enriched. for the chip /rna genes, genes associated cytoplasm as a location were enriched, and for the chip /rna genes, go terms associated with transcription factors were enriched, particularly among genes with a short half-life. these results suggest that different functional categories of genes are subjected to different modes of gene expression regulation.table list of go enrichments for transcripts in different chip-seq and rna-seq fractions a) chip /rna go enrichment go:id go: term number of genes false-discovery rate go: 0044822 poly rna binding 618 e-261 go: 0006412 translation 293 e-196 go: 0010467 gene expression 387 e-162 b) chip /rna go enrichment go: id go: term number of genes false-discovery rate go: 0005737 cytoplasm 138 e-2 go: 0070062 extracellular vesicular exosome 74 e-2 go: 0005635 nuclear envelope 15 e-2 c) chip /rna go enrichment go: id go: term number of genes false-discovery rate go: 0003677 dna binding 378 e-28 go: 0006355 regulation of transcription, dna-templated 266 e-17 go: 0006351 transcription, dna-templated 361 e-17 d) chip /rna /half-life < 4 h go enrichment go: id go: term number of genes false-discovery rate go: 0003677 dna binding 216 e-58 go: 0006351 transcription, dna-templated 192 e-37 go: 0006355 regulation of transcription, dna-templated 147 e-34 identification of upf exosc and stau as controlling factors for rna stabilities to evaluate the potential contribution of known rna degradation factors to the control of global rna stability, we chose three representative factors for analysis: upf exosc and stau it has been reported that upf regulates 3–20% of transcripts, highlighting the potential importance of upf in regulating rna degradation and abundance. exosc is an essential component of the exosome complex that is the major mrna degradation machinery in mammalian cells. to analyze the alteration of global mrna turnover by perturbation of representative factors, we examined exosc stau which regulates around 1% of bona fide mrnas, is a typical rna-binding protein involved in rna degradation. we used data from bric assay in the cells depleted in upf by sirna. as observed in a previous study, mrna levels of the gadd45a gene, which is a known target of upf were increased following upf knockdown, with an increase in half-life. we then conducted similar experiments using exosc and stau knockdown cells and prepared a similar rna-seq and bric-seq dataset. we observed an increase of fam120c mrna levels and an increase in half-life in exosc knockdown cells. in stau knockdown cells, the mrna levels of cdkn2b were increased, with increased half-lives.figure identification of upf exosc and stau as factors controlling rna half-lives. chip-seq and rna-seq for the gadd45a, fam120c and cdkn2b genes, which showed significant changes in rna half-lives in upf exosc and stau knocked-down cells, respectively. upper two panels in: chip-seq data for h3k4me and pol ii. lower two panels in: rna-seq data for control cells and the respective knockdown cells. bric-seq normalized graph for gadd45a, fam120c and cdkn2b genes in control cells or cells knocked-down for upf, exosc, and stau. scatterplot representing the chip-seq peak signal intensities of h3k4me on the y-axis and gene expression values on the x-axis with putative upf1-controlled genes, exosc5-controlled and stau1-controlled genes in blue, green and red dots respectively. we next looked for genes that may be regulated by these factors, and identified and genes where the mrna half-lives were extended by more than two-fold or -fold and showed mrna expression increase by two-fold or -fold in upf exosc and stau knockdown cells, respectively. as shown in figures 3a, 3b and 3c, we noticed that the transcripts that were not observed in the control knockdown cells appeared in the knocked-down cells in many cases. we examined the distribution of the dots of these genes, whose transcripts were stabilized and increase in corresponding knockdown cells, in figure 1b and found that they are enriched in the upper-left corner of the plot. in total, we identified and genes whose mrna half-lives are potentially controlled by upf exosc and stau respectively, consisting of and % of the total of genes, rna, half-life < 4 h) in this area. we also examined the overlap between the genes controlled by upf and exosc and found little overlap. although we could identify upf exosc and stau as control factors only for a limited population by this approach. these observations are only the first step in identifying the role of rna decay factors in determination of rna abundance through rna degradation, and further systematic analyses may facilitate identification of the complex regulatory mechanisms of mrna stabilities.table summary statistics used for the analysis of refseq transcripts conditions number of genes active transcription and half-life measured with short half-life with low expression and high h3k4me 866 and upf target or exosc target or stau target computational modeling of the effect of mrna half-lives on eventual mrna levels we predicted the rna abundance by normalizing the half-life to be  hours, the estimated default half-life, and obtained the least squares regression line between the predicted rna levels and h3k4me intensities. we made a threshold of ×  and ×  above and below the least-square regression line to define genes where the h3k4me intensities and rna abundance correlate. we found that out of genes that were available from bric-seq dataset, we found and genes that resided within ×  and ×  of the regression line, respectively. we checked the original gene expression of those genes prior to the simulation, and checked whether their measured rna abundance correlates with h3k4me intensities. we found out of and genes that resided within ×  and ×  of the regression line from predicted gene expression, we found and genes where measured to be outside of the threshold, respectively. it means that the rna stability of those genes contributed to the rna abundance. taken together, these results collectively support our claim that rna degradation significantly contributes in determining the eventual expression levels.figure computational simulations of the rna half-lives as a determinant for transcript level. computational simulation into the effect of rna half-life on rna abundance. each gene was simulated to have a half-life of  hours, potential default half-life, and the predicted gene expression was plotted against the h3k4me intensities, along with the measured gene expression. the dotted line indicates the least-squares regression line between the predicted gene expression and h3k4me intensities. result of the computational simulation of the predicted rna expression when each rna abundance was predicted from rna half-life when upf exosc and stau were knocked-down, in comparison to each control knock-down conditions. the predicted gene-expression of the decay factor knock-down and the measured gene-expression of control knock-down were plotted against h3k4me intensities. the dotted line indicates the least-squares regression line between the predicted gene expression and h3k4me intensities. summary statistics of computational simulation condition t1/2= hours siupf1 siexosc5 sistau1 additionally, we conducted chip-seq on h3k27ac, h3k27me h3k36me to build a linear model as described by wang, c., et. al.,. we built one linear model incorporating h3k4me h3k27ac, h3k27me and h3k36me intensities with half-life as an extra variable, and one without the half-life, to explain the rna abundance. we found that fitting increased from to in r-value, which confirmed the previous finding by wang. et. al.,. we examined whether the changes in the rna half-lives from knockdown of upf exosc or stau could explain the changes in the eventual transcript levels. for this analysis, we conducted a computational simulation. as shown in additional file 1: figure s we found that the described theoretical model can predict the changes of eventual rna levels with pearson’s correlation co-efficiency of, and, respectively. overall, we demonstrate that the simple computational model could reasonably explain the changes of eventual mrna levels, thus supporting our idea that the major determinant of the eventual rna levels in these cases is at the level of rna stability. we further simulated the rna abundance, from the changes in rna half-life, in relation to the chip-seq signal levels and we found and genes that were within ×  of the regression line that lied outside of the threshold prior to simulation, for upf exosc and stau respectively. possible feedback between mrna turnover and transcriptional initiation to analyze the possible feedback mechanisms between mrna turnover and transcription initiation, we used the genes in which both mrna half-lives and eventual transcript levels were increased more than two fold in knockdown cells. in and genes in upf and exosc knockdown cells, respectively, there were no significant changes in eventual transcript levels observed despite remarkable changes to their rna half-lives. we speculated that there might be a possible feedback between mrna turnover and transcriptional rate. if mrnas for a particular group of transcriptional repression factors are included in the upf1/exosc targets and their stabilized mrnas result in increased protein levels of such transcriptional repression factors, thereby enhancing the repression activities on their target genes, it would explain unchanged balance of eventual transcript levels for these genes. we examined whether any transcription factor binding sites were enriched in the upstream regions of the and genes. in the case of upf we detected significant enrichment of the consensus binding site for hic which is a transcription factor belonging to the zinc finger family. we validated the changes of its expression level and half-lives in control and upf knockdown cells, and confirmed hic as a upf target. we were unable to further validate direct binding of hif to target genes, since no effective antibodies are available. also, it is possible that hic may not be the only candidate, which may contribute to the feedback regulation. many zinc finger family transcription factors share consensus binding sequences. the list of putative zinc finger family transcription factors that have significant homology to hic in their dna binding domains with extended half-lives and increased eventual transcript levels upon upf knockdown are shown in additional file 2: table s these factors may collectively enable elaborate regulation of gene expression. identification of candidate genes controlled through rna stability in other cell types to further extend our idea that controls at the level of mrna decay contribute to determining eventual mrna expression levels in other cell types, we analyzed the published encode data and dbtss, which included chip-seq data of h3k4me and pol ii and rna-seq in a wide variety of cell types. in addition, the rna-seq data of subcellular fractionated mrnas were included in the dataset. we selected eight cell types for which all these datasets were available. we retrieved and analyzed the chip-seq data and rna-seq data as performed with hela cells. first, as shown in additional file 1: figure s we observed weak correlations between chip-seq intensities of h3k4me and gene expression values as reported in hela cells. since there were no bric-seq data in the encode dataset, we could not directly analyze the genes with short mrna half-lives. nevertheless, we could select genes for which chip-seq tags of h3k4me and pol ii were associated, thus indicated as actively transcribed in the cell line, although their gene expression levels were not at the expected levels. in the encode dataset, we also considered a positive signal of h3k36me which is a chromatin mark for transcriptional elongation, to further assure active transcription. we identified an average of candidate genes with active transcription and low rna abundance in each cell type, which should be regulated at the level of mrna decay. we examined and identified go terms that were significantly enriched, depending on cell types, and the top enriched terms were predominantly associated with dna binding. we also analyzed whether there are any cases for a particular gene to be selected as such a candidate in a cell type-preferred manner. we found that potential controls at the level of rna half-life were observed in a single cell type. we also examined if there was possible feedback between controls of rna half-lives and transcription initiations. we found that several transcription factor binding consensus sequences of the ets family and arebp family genes are enriched in the promoters of the genes inferred to be regulated at the level of rna stability in gm cells and hepg cells, respectively. in addition, we identified hic binding consensus sequences in the promoters of the genes inferred to be regulated at the level of rna stability in human h embryonic stem cells. the feedback regulations may be common in various cell types with distinct responsive transcription factors depending on cell type. we further analyzed the subcellular localizations of the mrnas of these genes using the corresponding encode data. we found significantly enriched mrnas in the cytoplasm only in the embryonic stem cells, suggesting there may be a characteristic regulatory mechanism for controlling rna stability in the cytoplasm of es cells. on the whole, these data should provide an important complement to the encode annotations, which aim to generate a complete catalogue of genetic elements explaining gene expression regulation.figure screening of candidate genes that may be controlled at the level of rna degradation from the encode dataset. number of candidate genes screened from the indicated cell types based on the encode data. gene ontology enrichment of genes in the chip / rna regions in the indicated cell types based on the encode data. chip-seq signal intensities were comparable between different cell types, although rna-seq-based gene expression values were remarkably different. the go term enrichment with the lowest p-value for each cell line is shown. enriched consensus transcription factor binding sites for genes in the chip / rna region for encode dataset. graphical representation of the enriched consensus binding sites in the promoter regions of the genes in chip / rna region for encode dataset. the list of cell lines with consensus binding site enrichment from the transfac database. the nuclear/cytoplasmic gene expression values for the genes where rna half-life may be the contributor to the rna levels and all other genes in h human embryonic stem cells. statistical significance of the difference is indicated under the x-axis. distinct controls of the rna stabilities of mrnas and non-coding rnas to examine whether the regulations at the level of rna half-lives are observed in lincrnas, we conducted a similar analysis for lincrnas in hela cells, as shown in figures  and we tentatively defined the dataset of lincrna as that of lincdb. as shown in figures 6a and 6b, we associated the chip-seq peak intensities and gene expression values. we unexpectedly observed distinct patterns from those of mrnas. namely, among lincrnas in hela cells, lincrnas had “short ” rna half-lives. of these, resided in the upper-left corner of the plot, suggesting that most of the lincrnas are controlled at the post-transcriptional level. because there are reports on the possible involvement of nmd in regulating non-coding rna, we examined the possible involvement of upf in the regulation of lincrnas. interestingly, none of the lincrnas were detected as potential upf1-controlled transcripts. in this data set, upf may have only limited contribution. by contrast, we found of the lincrnas were regulated by exosc figure contribution of the rna half-lives to transcriptional regulations of lincrnas. scatterplot showing the relationships between h3k4me intensities and gene expression for all lincrna and those that were bric-seq compatible. examples of lincrnas regulated by upf and exosc. upper panels are h3k4me and pol ii peaks, and lower panels are rna-seq measurements for the control and indicated knockdown. upf and exosc show the normalized decay curve from the bric-seq measurements for these examples. red lines are siupf or siexosc data, respectively and blue lines are si_control data. enhancer rnas regulated by upf or exosc. the label shows the genomic coordinates of these ernas. these show normalized decay curves for the bric-seq measurements for these examples as above. summary statistics used for the analysis of lincrnas conditions number of genes active transcription and half-life measured 141 with short half-life 103 with low expression and high h3k4me 84 and upf target 0 or exosc target we next examined whether enhancer rnas, which facilitate the functions of the enhancers, are also regulated at the level of rna degradation. we again used the encode data from hela cells. we retrieved the chip-seq data of h3k4me and h3k27ac, which are representative chromatin marks of active enhancers. among genomic regions having “peaks” of both h3k4me and h3k27ac, we identified cases in which there were rnas in the overlapping regions,  kb away from any refseq gene body and their rna half-lives were extended by upf knockdown by more than two fold in bric-seq assay. similarly, we identified cases in which half-lives of rnas in overlapping regions were examined and in which transcripts showed extended rna half-lives of more than two-fold in exosc knockdown cells. although further detailed experimental validations are necessary, these results may indicate that controls mediated by rna stability are used in determining the transcript levels of non-coding rnas. discussion here, we have described a genome-wide correlation among the signal intensities of chip-seq, gene expression values measured by rna half-lives measured by bric-seq. we identified that regulation at the level of rna degradation plays an important role in determining eventual rna levels. we demonstrated that this control may exhibit a particularly large effect in cases in which the chip-seq data and rna-seq data are inconsistent. indeed, in this study, we estimated that the abundance of mrnas is regulated by rna degradation in hela cells. furthermore, we applied similar approaches to analyze the public encode data and identified a total of candidate genes whose gene expression levels are likely to be controlled at the level of rna stability. we also found that these controls appeared to vary among cell types. to our knowledge, this is the first report that describes the integration of the chip-seq, rna-seq and rna half-life data in identifying genes that may receive post-transcriptional gene expression regulation. in go analyses, there are some limitations because go terms may be loosely defined for a particular gene. however, we first selected go terms with statistical significance, and we only used go term enrichment with more than genes in a set to ensure minimum false positives, which will ensure that the go terms described are representative of the biological phenomenon. we then found enrichment of transcription factors in transcripts with discrepancies between the chip-seq and rna-seq with short half-lives, and in particular, negative regulatory factors. in a recent study, haimovich et al. indicated that some rna degradation factors play a role in transcription, implicating a feedback loop for gene expression. our data suggest another mechanism by which this occurs, through a faster degradation rate of mrnas encoding transcription factors, particularly for those that negatively regulate gene expression, thus affecting the eventual rna levels. previous studies conducted by wang et. al., demonstrated that rna stability could in inferred from the residual errors in modeling the rna abundance from chip-seq data and they validated their claims from a half-life data; however, they used rna stability data from a different cell-line, and they could not accurately predict the eventual rna abundance from the rna stability. herein, by analyzing rna stability data from the same cell-line as rna-seq data, we were able to estimate the contribution of the rna stability to the rna abundance. in addition we were able to estimate the rna stability contributions on the rna abundance upon upf stau and exosc knockdown. although we did not demonstrate how the rna half-lives are controlled in the current study and we identified upf exosc and stau as control factors in some cases, they could explain at most 8% of the total mrna population. even for the cases of candidate upf exosc or stau1-controlled genes, it is possible that they may not be direct targets of these factors and that we may have picked up secondary or later effects as a consequence of upf exosc or stau knockdowns. another obvious drawback of our approach is that mrna half-lives were not directly measured by bric-seq for the encode dataset. therefore, it is possible that they may be mediated by other regulatory mechanisms, rather than at the level of rna half-lives, such as rna halting and abortive transcriptions. to minimize these possibilities, we selected the cases in which h3k36me a marker of transcriptional elongation, should be significant in the transcript regions. in spite of several drawbacks, we believe that genome-wide features of correlation among chip-seq, rna-seq and bric-seq should give an important starting point to further explore posttranscriptional regulatory mechanisms, for which only limited knowledge has been accumulated. indeed, recent papers have begun to reveal many human diseases that are caused by malfunctions of rna decay pathway. in particular, it has been made gradually clear that most immune-response mrnas are destabilized when they are not required via their cis-regulatory elements in the 3′ utr. it is proposed that such rna-decay mechanisms collectively enable rapid up-/downregulation of gene expression in response to environmental changes. the au-rich element is one of such elements widely found in the 3′ utr of mrna of immune-related genes. mice lacking are in the tnf-alpha mrna showed joint and gut-associated immunopathologies. the trans-acting factor regulatory rnase, which is induced by toll-like receptor ligands, interleukin -1β and mcp- is involved in the destabilization of mrnas including il mrna. regnase-1-deficient mice develop severe autoimmune disease because of excess production of cytokines, highlighting the importance of rna degradation-mediated gene regulation. in addition to immunological disorders, there are a growing number of cases with impaired rna decay regulation that cause disease, and they sometimes reveal unexpected connections between otherwise completely unrelated diseases. perlman syndrome, an autosomal recessively inherited congenital overgrowth syndrome associated with high neonatal mortality, is an obvious example. the survivors of this disease have a high risk of wilms tumor. recently, it was reported that the responsible gene of this disease is the exoribonuclease dis3l a homologue of exosome component dis. moreover, dis3l is mutated in approximately 3–6% of carcinomas. in this study, we have also analyzed the stabilities of non-coding rnas and mrnas. we and another study reported that stability of non-coding rnas is also tightly regulated, suggesting that the instability contributes to the dynamic nature of lincrnas. indeed, the stability of noncoding rnas has an impact on their biological function, although the exploration into its relevance in human disease has just begun. further enrichment of our knowledge on the control mechanisms on rna stability both for mrnas and non-coding rnas will shed new light on putative disease-associated genetic or somatic mutations.  by integrative analysis of chip-seq, rna-seq and our bric-seq, we showed that rna half-life may serve as an important post-transcriptional determinant of gene expression. we suggest that upf exosc and stau may play active roles in such controls. in addition, we propose the linkage between transcription and rna decay through regulated degradation of mrnas encoding transcription factors to maintain the steady state level of rna abundance. 
rna binding proteins, such as upf stau and exosc may play active roles in such controls. here, we used stability data of whole transcriptome measured by 5′-bromouridine immunoprecipitation chase sequencing, which enabled us to determine the half-lives of whole transcripts including lincrnas, and we integrated bric-seq with chip-seq to achieve better estimation of the eventual transcript levels and to understand the importance of post-transcriptional regulation that determine the eventual transcript levels. encode data analysis supported the idea that rna stability control aids to determine transcript levels in multiple cell types. in addition, rna stability is an important determinant of eventual transcript levels. we identified upf exosc and stau well-studied rna degradation factors, as controlling factors for 8% of cases. we identified discrepancies between the rna abundance estimated by chip-seq and measured rna expression from rna-seq; for number of genes and estimated that the expression level of genes was controlled at the level of rna stability in hela cells. in addition, we propose a feedback circuit that includes the regulated degradation of mrnas encoding transcription factors to maintain the steady state level of rna abundance. integrative analysis of chip-seq, rna-seq and our bric-seq showed that transcriptional regulation and rna degradation are independently regulated. intriguingly, these regulatory mechanisms were distinct between mrnas and lincrnas. histone epigenome data determined by chromatin immunoprecipitation sequencing is used in identifying transcript regions and estimating expression levels. however, this estimation does not always correlate with eventual rna expression levels measured by rna sequencing. part of the inconsistency may arise from the variance in rna stability, where the transcripts that are more or less abundant than predicted rna expression from histone epigenome data are inferred to be more or less stable. bric-seqchip-seqintegrative analysisnext-generation sequencingrna stabilityestimation of transcriptional levelissue-copyright-statement© the author 2015 computational simulations reasonably explained the changes of eventual mrna levels attributable to the changes in the rates of mrna half-lives.
 histone epigenome data determined by chromatin immunoprecipitation sequencing is used in identifying transcript regions and estimating expression levels. however, this estimation does not always correlate with eventual rna expression levels measured by rna sequencing. part of the inconsistency may arise from the variance in rna stability, where the transcripts that are more or less abundant than predicted rna expression from histone epigenome data are inferred to be more or less stable. however, there is little systematic analysis to validate this assumption. here, we used stability data of whole transcriptome measured by 5′-bromouridine immunoprecipitation chase sequencing, which enabled us to determine the half-lives of whole transcripts including lincrnas, and we integrated bric-seq with chip-seq to achieve better estimation of the eventual transcript levels and to understand the importance of post-transcriptional regulation that determine the eventual transcript levels.  we identified discrepancies between the rna abundance estimated by chip-seq and measured rna expression from rna-seq; for number of genes and estimated that the expression level of genes was controlled at the level of rna stability in hela cells. encode data analysis supported the idea that rna stability control aids to determine transcript levels in multiple cell types. we identified upf exosc and stau well-studied rna degradation factors, as controlling factors for 8% of cases. computational simulations reasonably explained the changes of eventual mrna levels attributable to the changes in the rates of mrna half-lives. in addition, we propose a feedback circuit that includes the regulated degradation of mrnas encoding transcription factors to maintain the steady state level of rna abundance. intriguingly, these regulatory mechanisms were distinct between mrnas and lincrnas.  integrative analysis of chip-seq, rna-seq and our bric-seq showed that transcriptional regulation and rna degradation are independently regulated. in addition, rna stability is an important determinant of eventual transcript levels. rna binding proteins, such as upf stau and exosc may play active roles in such controls. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords bric-seqchip-seqintegrative analysisnext-generation sequencingrna stabilityestimation of transcriptional levelissue-copyright-statement© the author 2015 