CART Classification of Human 5′ UTR Sequences

Class I Transcripts Have Long 5′ UTRs Filled with Stable Secondary Structures, uAUGs, and uORFs

Kozak (1991) presented a comprehensive review on 5′ UTR features involved in translation control and predicted that many of the growth-related proteins would be poorly translated. Substantial experimental evidence has been accumulated in recent years that supports this prediction. Some of the well-studied transcripts that are poorly translated because of the presence of stable secondary structures or the presence of uAUGs in the 5′ UTR are ornithine decarboxylase (ODC), TGF-β3, β1,4-galactosyltransferase (β4GalT-I), cyclin D1, p53, AdoMetDC, RARβ2, and potassium channel ROM-K3. Our classification model classified all these transcripts (not included in the training set) in class I along with many other that have highly stable secondary structures and uAUGs. In most of these cases, translation occurs by the cap-dependent scanning model (Merrick and Hershey 1996). The cytoplasmic cap-binding protein, eIF-4E, participates in unwinding the secondary structures, and hence, its availability is crucial for the translation of these highly structured transcripts. When the availability of active eIF-4E is limiting, these transcripts are poorly translated. One way to overcome this problem is the overexpression of eIF-4E. Elevated levels of eIF-4E have been found in many tumor cell lines and almost all breast carcinomas. As a consequence, some of these poorly translated transcripts in class I might be efficiently translated in cells with eIF-4E overexpression. ODC is a good example for this as its levels were found to be drastically increased in eIF4-E transformed cells (Shantz et al. 1996).

One of the other ways these poorly translated transcripts can get rid of these inhibitory features is by a shift in the transcription start site and alternative splicing. TGF-β3 and β4GalT-I are good examples for this mechanism. Enhanced translational efficiency of TGF-β3 was observed in human breast cancer cells, and its 5′ UTR lacks the 5′ end of ∼870 nucleotides (Arrick et al. 1994) that contained inhibitory secondary structure. The β4GalT-I gene results in two transcripts with different 5′ UTRs. Charron et al. (1998) showed that mammary gland-specific β4GalT-I transcript, with truncated 5′ UTR that lacks extensive secondary structure, was efficiently translated both in vitro and in vivo. Both these transcripts that were not included in the training set have been successfully classified as class III transcripts by our classification model.

Class II and III Transcripts Have Small 5′ UTRs Free from Stable Secondary Structures, uAUGs, and uORFs

Class II mRNAs contain a 5′ TOP that regulates the translation of these transcripts in a growth-dependent manner. 5′ UTRs of this class were relatively short and almost completely free from the inhibitory features that were commonly observed in class I 5′ UTRs. However, we found a few transcripts from the other two classes that contain this regulatory element. Avni et al. (1997) showed that elongation factor 2 (EF2) and β1-tubulin, which contain 5′ TOP, are not regulated in a growth-dependent manner but regulated in a cell type-specific manner. They showed that the downstream sequences suppressed the regulatory features of the 5′ TOP and suggested that the mRNAs with longer 5′ UTRs might not be regulated in the same way as ribosomal proteins. Our classification model correctly classified all those transcripts of class I and III even though some of them contained 5′ TOP.

In the classic review, Merrick (1992) suggested the optimal characteristics for efficient mRNA translation, and most of the transcripts in class III have the favored characteristics. Hence, we suggest that most of these mRNAs are likely to be efficiently translated or, at least not repressed at the translational level. For a definite proof, we would have to wait for the experimental results. In a personal communication, Dr. David Morris provided the list of highly translated genes in human foreskin fibroblast cell lines by using a method called sucrose gradient analysis (Zong et al. 1999). Some of the genes in the list that were not in our training set are vimentin, desmin, CD59, caveolin-1, decorin, Ku80, and cytokeratin 8. Our classification model was able to correctly classify all these into class III.

Why CART Is Good for the Present Analysis

We analyzed large multivariate data that included both continuous and categorical variables. The CART technique is particularly applicable for studies like this, in which many of the variables considered do not seem to follow any particular distribution. In other words, we didn't make any parametric assumptions regarding the distributions of the variables under study. Moreover, our analysis was pattern driven rather than model driven; rather than building a coherent global model that includes all variables of interest, our classification algorithm produced a set of statements about local dependencies among predictor variables (in rule form with yes or no answers).

Also, CART uses predictor variables independently. That is, initially the entire data is partitioned into two subgroups according to the variable that produces the best split, for example, presence of TOP. Then, in each of the resulting strata, the process is repeated recursively until none of the selected variables shows significant influence on the split or the size of the subgroup is too small. In the final process, subgroups of cases that do not differ in any of the characteristics under study are joined together to form homogeneous classes.

CART also picks the best discriminating variables and ranks all the variables according to the relative discriminating power. We tried other classical methods such as linear discriminant analysis (LDA) and quadratic discriminant analysis (QDA) by considering the top three relevant variables picked by CART. The models of LDA and QDA didn't seem to give any better prediction than the CART model (data not shown). This might be due to the non-normality of the data and forcing parametric assumptions that didn't seem to exist.

Our experience shows CART as a useful data-mining tool for analyzing large data with many variables, where conventional statistical methods like LDA or QDA are not effective.

Limitations of CART

CART exhibits its greatest strengths in classification trees with a highly nonlinear structure (e.g., the 5′ UTR data in the present study). The closer the model is to linear, the less useful CART will be. When data exhibit a genuinely linear structure, CART is not a particularly useful analytical technique. Another important problem with CART is heteroscedasticity (within class variance). If the cases within a node genuinely belong together but have a high variance because of heteroscedasticity, CART may select a spurious split to partition the data. Although cross validation is designed to protect against the retention of such splits, some do survive the pruning process. In the present study class II is highly homogeneous, followed by class II and class III. The classification model could clearly segregate class II from the other two. However, 6.6% of class I were misclassified in class III, and 13.2% of class III were misclassified in class I (Tables 2 and 3). This might be due to the heteroscedasticity present in classes I and III. On the other hand, experimentalists are encouraged to look into these misclassified transcripts for proper reasons for their misclassification.

Conclusion

We made a comprehensive analysis of a large collection of 5′ UTRs and broadly classified the data into three functional groups. The class I transcripts seem to be very poorly translated under normal conditions, and those from class III might be candidates for efficient translation. Our classification model and the data we have generated may provide valuable information for experimentalists engaged in translational control and regulation studies. For example, the next natural step is to look for IRES or internal entry points within class I 5′ UTRs experimentally as well as computationally.

One of the main goals of our study is to develop a complete gene prediction system. As a first step toward this goal, we recently added a 3′-terminal exon-recognition module (J. Tabaska, R.V. Davuluri, and M.S. Zhang, in prep.) to the internal exon finder, MZEF (Zhang 1997). Our next step in achieving this task is the development of a 5′-terminal exon-detection module. We are presently working on the extension of the results of this work to the 5′-terminal exon-prediction program. The feature variables studied here would be valuable to identify the correct start site and separate 5′ UTR from the coding region. The drastic differences in 5′ UTR features observed between the three classes indicate that distinct models could be used to predict 5′-terminal exons. Hence, the classification of the 5′ UTRs into homogenous classes would facilitate the building of separate models for each class so that the overall quality of the 5′-terminal-exon prediction is expected to be higher than the prediction based on a single mixture model.

5′ UTR Database

A set of 954 human 5′ UTR sequences was obtained from the 5′-end-enriched cDNA library (Suzuki et al. 1997, 2000) with their mRNA start sites. The 5′-end-enriched cDNA library was constructed to isolate the mRNA start site of long mRNA, by using a method called oligo-capping (Maruyama and Sugano 1994) with some modifications. We collected 5′ UTRs from this library as follows: First, cDNA sequences were clustered with DYNACLUST (Dynacom) after removing the oligo-capped 5′-oligonucleotide sequence from each 5′ end.DYNACLUST is a database management software, which clusters the sequences using BLAST with the score of e − 40 for 400 bp. The position of the translation start site (ATG) was marked for each sequence according to the annotation in GenBank. Then, the sequence between the oligo-capped 5′-oligonucleotide sequence and the translation start site (ATG) was extracted from each cluster. If alternative mRNA start sites or translation start sites were observed, then the cDNAs containing the longest 5′ UTRs at both the 5′ and 3′ boundaries were selected as the representative. (for details, see Suzuki 2000).

The experimentally derived set of 954 5′ UTRs was augmented with a second set of 1613 full-length 5′ UTR sequences retrieved from the UTRdb (Pesole et al. 2000) database. Only those sequences with UT feature tag as complete 5′ UTR are considered. These sequences were extensively verified by going through their corresponding GenBank records, and only those records with evidence = experiment were considered. All the redundant and ambiguous sequences were eliminated, and finally, a nonredundant set of 2312 5′ UTR sequences was prepared for the analysis. A sequence was considered redundant if it has 90% similarity and 90% overlapping with a larger sequence in the database. However, there may be more than one 5′ UTR for some genes because of alternative splicing and usage of different transcription start site. From this database the following three classes of 5′ UTRs were considered for analysis: Class I, the first class, consists of 5′ UTRs of growth factors, their receptors, transcription factors, proto-oncogenes, cytokine receptors, and tumor suppressor genes. Most of these are understood to be translationally repressed mRNAs. Class II, the second class, consists of TOP mRNAs. TOP mRNAs are vertebrate transcripts with a C residue at the cap site, followed by an uninterrupted stretch of 4–13 pyrimidines, called 5′ TOP, encode for ribosomal proteins and elongation factors 1α and 2α. The translation of this class of mRNAs is regulated in a growth-dependent manner. Class III, the third class, consists of 5′ UTRs of highly expressed genes, tubulins, globins, globulins, myosins, caseins, glycolytic enzymes, β-actin, γ-actin, and histones. The expression of these genes is controlled mainly at the transcriptional level, and their transcripts are believed to be efficiently translated. In other words, these genes are either translationally efficient or (at least) not repressed at the translational level. In contrast, the first two classes of genes are tightly regulated at the translational level in stringent ways. There are 226 5′ UTRs in the first class, 70 in the second, and 76 in the third class. The complete data set is available at ftp://cshl.org/pub/science/mzhanglab/ramana.

Data Analysis

Procedure for Constructing a CART Tree

The key components of tree-structured data analyses are tree growing, tree pruning, and optimal tree selection. Tree growing depends on splitting rules and stopping criteria. CART begins with all the data points in the learning sample, L. The CART classification tree initially consists of one node—the parent node of the tree, which contains all the points in L. The CART program searches through all possible values of all the variables, looking for the split that best separates the classes. The first split creates two child nodes. CART takes each of the child nodes and recursively partitions each child node in the same way that it partitioned the parent node. CART evaluates the goodness of any candidate split using an impurity function. A node that contains members of only one class is perfectly pure, and the node that contains an equal proportion of every class is least pure. Given a node t with estimated class probabilities p(j/t), j = 1,…,J, and a measure of node impurity, CART searches for the split that most reduces node, or equivalently, tree impurity. CART provides different impurity functions, for example, Gini Measure, Twoing criterion, etc. (for more details, see Breiman et al. 1984). We tried all the impurity functions and finally selected Gini Measure, which best suited our data. For a node t with estimated class probabilities p(j/  t), j = 1, … , J, Gini Measure is defined as

The next important step in tree growing is stopping. However, stopping is not an essential component of CART. CART grows the tree until no further growth is possible, that is, terminal nodes have only one case, or terminal nodes with more than one case are identical on the classification variable. The resulting maximal tree is called T_max. Because T_max usually overfits the data (i.e., being overly sensitive to irregularities in data) especially when it is noisy, CART applies a pruning and evaluation procedure to find the optimum-sized classification tree,T_opt. The pruning procedure generates a nested sequence of smaller and smaller subtrees, {T_max,…,T₂,T₁}, from which CART selects the T_opt that has the lowest or near lowest cost of misclassification as determined by a cross validation procedure. In other words, CART performs postpruning, whereby as many child nodes as possible are eliminated as long as the overall estimated accuracy of the tree is not significantly reduced.

Feature Variables

The following variables were used as predictor variables in CART analysis:

1. |UTR length − 80|: Kozak (1991) has shown that the rate of mRNA translation increases proportionally as the length of the 5′ UTR is increased from 17 to 80. Hence, taking 80 as the optimal 5′ UTR length for efficient translation, we considered |UTR length − 80| as one of the variables.

2. Free energy estimate of secondary structure (ΔG): The latest version of the mfold (Mathews et al. 1999) program was obtained from Michael Zuker (mfold v.3.1). The mfold program minimizes a free energy function, which sums contributions from stacking, loop lengths, etc. It actually estimates the difference between the free energy of the unfolded state and folded state. For any given RNA sequence length, the lower the energy estimate the more stable the predicted secondary structure. Local free energy estimates were worked out for each 5′ UTR. This was done by cutting the longer 5′ UTRs into sequences 200 bp in length with an overlapping window of 100-bp size between successive sequence fragments. The most stable secondary structure predicted in each sequence fragment was considered, and its free energy estimate is recorded in the database.

3. Number of stable folds: The number of predicted folds with minimum folding energy of less than −50 kcal/mole were counted.

4. GC percentage: G + C percentage was calculated for each 5′ UTR.

5. G/C ratio: The absolute value of G/C − 1 was calculated. If the value of G/C is in the neighborhood of 1, then the chance of forming stable secondary structures is high, and it is less otherwise.

6. A/T ratio: Similar argument holds good for the absolute value of A/T − 1.

7. Number of uAUGs: The most important positions for efficient translation are a purine at the −3 position and a G at position +4, where A of the AUG codon is position +1 (Kozak 1997). If either of these main features is unfavorable, positions +5 and +6 can influence the start site efficiency, with a A or C being preferred at position +5 and a U at position +6 (Boeck and Kolakofsky 1994). An AUG is said to be in good context if it satisfies these criteria, and all those uAUGs that were in good context are counted for each 5′ UTR.

8. Number of uORFs: Some mRNAs with a long 5′ UTR contain one or more uORFs, and these uORFs are often inhibitory for the translation of the downstream coding region. All those uORFs with good initiation context, as defined above, were counted for each 5′ UTR.

9. TOP: This is a categorical variable with value 1 if 5′ TOP exits and 0 otherwise.

10. Start site context: This is another categorical variable taking value 1 if the start codon (AUG) is in good context and 0 otherwise. Along with these 5′ UTR features, we have also included the following two more variables that may also influence the mRNA translation:

11. CDS length: The length of the coding sequence in terms of the number of amino acids was calculated.

12. Codon bias: Codon bias was calculated according to the Karlin formula (Karlin and Mrazek 1996) with slight modifications.

where p(a) are the amino acid frequencies, N is the total amino acid count and B_g(a) are calculated as follows:

with g^(r)_xyz being the relative codon frequencies for the gene g_r and g_xyz the relative codon frequencies of a reference set. Here the reference set constitutes a set of 30 genes with poor translational efficiency.