Inferring the expression variability of human transposable element-derived exons by linear model analysis of deep RNA sequencing data

Recent years have witnessed the emergence of a few high throughput technologies, i.e. Affymetrix exon array and RNA-seq platforms, for generating large-scale exon-level gene expression profiling. While array-based data have been widely employed to study transcript splicing [1–3], there are inherent limitations of this technique, such as the potential cross-hybridization of the probes of one exon to the transcript of another exon, and the incomplete probe coverage [3, 4]. Especially, because the probe design for a gene is based on the exons included in one or multiple manually annotated or computationally predicted gene models [5], many infrequently used exons are never represented by any probe, and therefore their expression levels cannot be measured accordingly. This problem can be circumvented by RNA-seq technology, which provides hypothesis-free single nucleotide resolution of gene expression so that, theoretically, any expressed sequence can be detected and quantified, given appropriate computational/statistical methods and sufficient sequencing depth. RNA sequencing based gene expression data have become the major information source for detecting alternative splicing (AS), alternative cleavage, polyadenylation (APA) events, gene fusions, and splicing eQTLs (expression Quantitative Trait Loci) [4, 6–11].

Transposable elements (TEs) constitute approximately 44% of the human genome. In disease biology, the importance of TEs is highlighted by the potential association with genetic instability, one of the principal hallmarks and causative factors in cancer [12–15]. A recent study showed that Estrogen Receptor α (ERα), which is involved in human breast cancer, preferentially targets mammalian interspersed repeats (MIRs) transposons [16]. The exonization of mutated TE sequences has proven to be a significant mechanism for the creation of novel exons [17–20]. In the TranspoGene database, 1423 human exonized TEs (TE exons), involved in ~1700 RefSeq genes, have been collected [21]. Most TE exons are alternatively spliced during the post-transcriptional modification of RNAs [22]. Nevertheless, they are not infrequent in mature transcripts [4, 21, 23]. Past studies have suggested the diverse roles of TE alternative splicing in gene regulation [24–27]. Many human genetic diseases have been ascribed to TE exonization [28–30]. Analysis of single nucleotide polymorphisms (SNPs) also revealed that exonization of TEs can be population-specific, implying that exonizations may enhance divergence and lead to speciation [20].

The retention and splicing of a TE exon has been found related to multiple genomic factors, such as the TE category and exon length [21, 31]. These published results were mainly established by the analysis of EST and microarray data using descriptive statistical methods, where other factors cannot be considered simultaneously when inferring the effect of a specific one. In this work, we extended the previous studies by an integrative statistical analysis of multiple RNA sequencing data. In particular, we not only studied the digital expression levels of TE exons but also addressed the variability quantities of TE exon expressions across different biological samples and groups. Our findings provide advanced insights into the establishment level of TE exons and their roles in the evolution of human genome.

Computational identification of an exonized TE is usually based on the alignment with an exonic part of an EST (or cDNA) and the existence of canonical splice sites in the sequence. However, the functional importance of a TE exon is mainly determined by the expression level in a specific tissue and biological environment. In this study, we demonstrated that, with relatively high expression variability, most TE exons in mRNAs, especially those without exact counterparts in the UCSC RefSeq gene tables, demonstrate low but still detectable expression levels in most tissue samples. We also showed that the TE exons in CDSs are less likely to be constitutively expressed than those in 3′ (5′) UTRs. These results confirmed previous studies regarding the establishment of TE exons in the assembly of human mRNAs [4, 21, 22].

A recent study reported that MIR exons have highest inclusion levels among exonized TEs [23]. We found this result can be further generalized. That is, the exons derived from younger TE families (Alus, L1s) are less well established in the host genes than those derived from older TEs (MIRs, L2s) of the same super-classes, as indicated by the expression quantities including the measures of presence, level and variability. A deviation from the general association is that, although the activity burst of L1 families was 20–30 million years earlier than that of Alus [40, 41], the expression quantities of the exons derived from the elements of these two families are similar to each other (p & 0.05). This similarity may be due to their common consensus target sites [42, 43]. Another explanation is that the insertion of Alu and L1 elements frequently disturbs the stability of the target gene [14], and the human genome evolution process has likely formed a mechanism to remove the exons derived from such elements.

A previous study showed that the lengths of exons inversely influenced the TE inclusion levels in transcripts [44]. This relationship was confirmed by our analysis focused on exonized TEs. We are especially interested in this finding and wonder if it can be also derived by EST data analysis. Thus, we repeated the linear statistical analysis using EST-based inclusion level (retrieved from the TranspoGene/Human exonized Tes table) as the dependent variable in the place of rescaled RPKM. The result moderately supports the RNA-seq data based analysis in that the observed length effect is negative but not significant (p & 0.05). Such an observation even holds when the analysis is focused on the subset of Alu-derived exons. Herein, we speculated that, to which extent, the digital expression quantities of TE exons are consistent with the inclusion levels calculated using EST data [21]. By a simple regression analysis (Additional file 4), we found that 47.8% of variability across TE exons in the digital expression measure can be explained by the inclusion levels (p = 0). These results indicate that the EST data and digital expression data can be complementary in inferring the expression patterns of TE exons.

Our work also extended the previous studies by several novel findings. Among them, the finding that a high TE nucleotide proportion leads to lower retention rates of the expressed TE exons in mRNAs is intuitively correct. The other two as described below warrant further investigation.

First, the considered genomic factors impact the expression level and the expression variability of TE exons in inverse ways. For examples, compared to those derived from Alus, MIR-derived exons have high rescaled RPKMs but the variability measure (CV) across the samples is lower. The expression profile of MIR exons suggests that they may have become functionally essential in human.

Second, not only the exons derived from Alu elements but also the exons from the transposons (retrotransposons) of other families were preferentially established in C2H2 ZNF genes. In fact, such a preference is more remarkable in the exons derived from ancient TEs. For example, among the seven TE exons derived from L2 elements and hosted by ZNF genes, five are highly expressed. However, because most TE exons are derived from Alu elements, the roles of these elements in shaping ZNF genes and in the regulatory functions are most predominant.

Understanding the expression pattern of the TE exons not annotated in RefSeq is of special biological importance. For these exons, the distribution of rescaled RPKMs shown in Figure 1 (represents the 16 sample in cluster G1) can be regarded as a canonical profile. Herein, a challenging question is why the reads mapped to the genome regions are so rare to detect in the samples of cluster G2 (Additional file 1), which include three prostate tumors (prAd-1, 2, 3) and one normal prostate sample (prNorm-3). In the Results section, we explored this issue from biological aspects. Another potential explanation lies in the filtering process of the mapped reads. This can be scrutinized from two aspects. First, a relatively short read naturally has more chances to be aligned to multiple sites of the reference genome and thus to be filtered out as an “ambiguously mapped read”. Second, TE exon sourced reads are more likely to be “repetitively mapped” compared to others. As a result, the true expression level or un-expression ratio of TE exons may be underestimated or overestimated in G2 samples where the read length is short (33 bp). However, this elucidation is not supported by the information of two prostate (cancer) cell lines. That is, the samples LNCaP and PrEC have the read length (35 bp) close to G2 samples but the expression profile of TE exons similar to G1 samples. Apparently, this phenomenon warrants additional follow-up in independent laboratory studies.

As discussed above, the assembly of short sequencing reads derived from repetitive elements is potentially associated with a high unmapped fraction. The problem arises from exclusion of ambiguously mapped reads in data preprocessing, and is likely serious for a dataset generated by single-read sequencing methods. However, for a dataset generated by paired-end sequencing methods, in which the two ends of an mRNA fragment are measured, the uncertainty should be largely alleviated [45]. A plausible explanation is that, in mapping paired-end reads, the hit(s) of a read is determined not only by its alignment with the genome but also by the alignment of its mate read [46]. Usually, an mRNA fragment to be measured in paired-end RNA-seq experiments is long (at least 200 bp) [33, 35], sufficient to cover a sequence longer than a medium-sized exon. In this context, the probability that a paired read is mapped to multiple genomic positions is low. The observed repetitive hits can be largely attributed to the existence of paralogous genes. Accordingly, we can assume the chances for ambiguously mapped reads solely due to repetitive elements are rare in paired-end RNA-seq datasets.

It should be noted that in the linear model analysis, we included four independent variables whose quantities or categories can be directly inferred from the TranspoGene human exonized TEs table as the predictors of the retention level of TE exons. We believe that the quality of the splice sites of the exonized TEs is a promising predictor for their retention levels in mature mRNAs. However, the measurement of splice quality is yet an unsolved challenge, though many studies on the prediction of splice sites have been published [44, 47]. We have not devised an appropriate method to derive a measure that is related to the sequences of the splice sites and is in association with the retention levels. Solving this challenge is the goal of another study in our research agenda.

In this study, we employed a ratio quantity, i.e. the rescaled RPKM (see Methods Section), to infer the expression (retention) level of an exon. Theoretically, this estimate may be biased due to the potential existence of primary RNAs. However, due to following reasons, it does not impact the validity of our conclusion on the variability and influential factors for the TE exon retention. First, the reads mapped on the simulated exons, i.e. the intron sequences immediately downstream of the un-annotated TE exons, are very rare in the analyzed datasets except for the six ovarian cell line samples. Second, we used the rescaled RPKM of the simulated exons as the baseline to determine the “expressed” genes and to test the existence of the transcripts containing un-annotated TE exons. An alternative quantity (in the place of rescaled RPKM) is the inclusion level for an (TE) exon of interest. The value can be computed from the counts of the reads mapped to the flanking exon-exon junctions [6]. However, in this way, most reads in a dataset will be unused because the length of a junction must be shorter than the read length. As a result, given the limited sequencing depth of currently published RNA-seq datasets for individual samples, the inclusion level estimate of an alternative exon may be severely biased or difficult to explain because of the unbalanced read counts in flanking junctions as shown in the Figure 1 of [4]. Therefore, there is a dilemma on using either of the two methods to exactly model the splicing pattern of a specific exon, but this is beyond the scope of our study.

As mentioned in the Results section, we hypothesize that there is a tissue-specific mechanism underlying the association between the expression-based sample grouping (G1, G2 and G3) and the tissue/disease type-based classification. However, so far, we lack additional data to provide further evidence for this hypothesis. On the other hand, the differences among G1, G2 and G3 in the profiles of the retention rates of TE exons could be attributed to the variations of the used experimental protocols (especially those for cell culture) in generating the analyzed data. Previous studies have shown that pol III SINE transcripts increase when the cells are stressed [27, 38]. However, it seems that this finding cannot explain our observation. More specifically, many TE exons (including those derived from SINEs) demonstrate high inclusion levels in the six ovarian cell lines (G3), but the G3 cells were not stressed, i.e. they were cultured in normal medium and temperature [36]. Another related problem is that substantial pre-RNAs are observed present in the six ovarian cell lines (G3 class) but not in other samples (Figure 1, Additional files 1 and 2). One may speculate whether the variation of methods for extracting mRNAs may cause such a difference. We compared the experimental protocols of the 26 RNA-seq datasets listed in Table 1, and found that the method [36] used to extract mRNAs from G3 samples is the same as the method (using Invitrogen TRIzol) [32] employed for the liver tissue samples (of G1 class). This indicates that there could be other factors contributing to this issue.

A recent publication showed that Alu exonization is controlled by the competition of hnRNP C and U2AF65 genes [39]. We found that the RPKMs of hnRNP C gene in G2 and G3 samples are significantly (p < 0.01) lower than those in G1, while such a group specific characteristic is not observed for the expression level of U2AF65 gene (p & 0.05). This observation, together with the descending order of the digital expression levels of TE exons from G3, G1 to G2 modestly suggests that hnRNP C (but not U2AF65) regulate the exonization of TEs in a non-linear way. In other words, this implies that both too high (G3) and too low (G2) exonization levels of TEs are related to a lower expression level of hnRNP C gene.