A hidden Markov model approach for determining expression from genomic tiling micro arrays

Background Genomic tiling micro arrays have great potential for identifying previously undiscovered coding as well as non-coding transcription. To-date, however, analyses of these data have been performed in an ad hoc fashion. Results We present a probabilistic procedure, ExpressHMM, that adaptively models tiling data prior to predicting expression on genomic sequence. A hidden Markov model (HMM) is used to model the distributions of tiling array probe scores in expressed and non-expressed regions. The HMM is trained on sets of probes mapped to regions of annotated expression and non-expression. Subsequently, prediction of transcribed fragments is made on tiled genomic sequence. The prediction is accompanied by an expression probability curve for visual inspection of the supporting evidence. We test ExpressHMM on data from the Cheng et al. (2005) tiling array experiments on ten Human chromosomes [1]. Results can be downloaded and viewed from our web site [2]. Conclusion The value of adaptive modelling of fluorescence scores prior to categorisation into expressed and non-expressed probes is demonstrated. Our results indicate that our adaptive approach is superior to the previous analysis in terms of nucleotide sensitivity and transfrag specificity.


Background
Tiling micro arrays query genomic sequence in a manner not biased towards coding transcripts and has proven a valuable resource in the exploration of genomic expression. To date, the approach has been applied to Human, Arabidopsis, and Rice [3][4][5][6][7][8][9][10][11] and has been reviewed recently [12,13]. The challenge of tiling micro array analysis is to categorise genomically consecutive probes as either expressed or non-expressed. The Affymetrix tiling array analyses [1,3,4] invoke extensive experience from gene arrays. Here, chip noise between replicates is dealt with using a quantile-normalisation procedure [14]. Background and probe sequence-specific signal is normalised using the MAS 5.0 method [15]. In order to convert expression values into transcribed fragments (transfrags) Affymetrix uses a smoothing window approach, whereby expression values within a genomic neighbourhood are averaged. A region is labelled as expressed if all probes exceed a threshold value. The threshold is determined from expressed RNA spiked into the micro array experi-ment. Neighbouring regions closer than a distance threshold are then joined, forming transfrag predictions. We will use the collective acronym ANTM for for Affymetrix's normalisation and transfrag methods. For a review of existing algorithms for tiling micro array analysis see Royce et al. [16].
Previous related work includes analysis of ChlP-chip experiments using Affymetrix tiling arrays to predict transcription factor binding sites [17]. These ChlP-enriched regions on the tiling arrays are predicted using log-odds scores for each probe generated by a two-state hidden Markov model. Tiling array intensities have been used together with a rudimentary HMM gene finder to establish likely exon-intron transcripts in windows of significant expression intensity [18]. This work, however, does not treat the probe intensities as observables of the HMM but as separate evidence supporting prediction of expression.
The approach to tiling array analysis presented here employs a hidden Markov model (HMM), trained directly on the correspondence between tiling array fluorescence scores and annotation of expression. The model is then used to categorise consecutive probes as expressed or nonexpressed. Training the model using annotation of expression presents a problem since, in many cases, annotation does not represent actual expression. This may be due to a number of factors, such as non-expressed genes, alternative splicing, un-annotated exons, antisense transcription, and non-coding RNAs. This type of incomplete information presents a problem to machine learning approaches as missing information compromises the training of the model unless properly addressed. We describe a two-step modelling and training approach that allows the training algorithm to down-weight dubious annotation in the context of the relevant cell line. The adaptive nature of our method avoids the ad hoc thresholds characteristic of previous work. As a case study we analyse Human tiling array data of a single cell-line from Cheng et al. [1].

Data
The data we use is derived from the tiling array experiments on Human chromosomes 6,7,13,14,19,20,21,22, X, and Y [1]. The tiling arrays use 25-mer probes querying genomic positions at five nucleotide intervals on the genomic sequence. Only non-repetitive sequence is tiled resulting in a coverage of approximately 30%. From the Cheng et al. raw array data we analyse the subset querying poly-adenylated RNA extracted from Human neuroblastoma cells (SK-N-AS cell line).

Array and probe normalisation
The Affymetrix MAS 5.0 method addresses probe and array normalisation by, in addition to perfect match (PM) probes, also designing mismatch (MM) probes, where the complementary base at the 13th position of the 25 mer probe is used instead of the actual base. The mismatch probe is expected to capture much of the probe sequence and array-specific contributions to the signal as well as problematic probe eccentricities such as secondary-structure and cross-hybridisation effects while capturing little real signal. The difference in fluorescence score for the PM probe and an Idealised-MM probe score is then used to estimate expression signal from the queried genomic position.
We tested a number of strategies for array and probe normalisation. These included PM -MM with and without quantile normalisation and a novel approach inspired by Naef and Magnasco [19] that shares similarities with GC-RMA [20,21].
For the normalisation inspired by Naef and Magnasco sequence and array specific signal is modelled in a way that allows this to be subtracted from the raw signal. For each chip a large sample of probe scores are selected and used to fit position-specific estimates of nucleotide contribution to the signal. Three free parameters for each position times the 25 parameters for each probe sequence position gives rise to 75 parameters. These are modelled by the following linear system: where each nucleotide type (A, C, G or T) in position k in a 25-mer probe is modelled by a triple in the indicator function a i (a 3k-2 , a 3k-1 , a 3k ). I.e. an A at position k has corresponding coefficients (a 3k-2 , a 3k-1 , a 3k ) = (1, 0, 0), a C has corresponding coefficients (0, 1, 0), a G has corresponding coefficients (0, 0, 1), and a T has corresponding coefficients (-1, -1, -1). A least-squares solution to this system of equations for x i provides estimates of the expected contribution to the expression for any nucleotide at position k. That is, the estimated contribution to the expression level of an A at position k is given by x 3k-2 , a C at position k is given by x 3k-1 , a G at position k is given by x 3k , and a T at position k is given by -x 3k-2 -x 3k-1 -x 3k . To obtain a representative expression level over different replicates the mean log-likelihoods are computed: An illustration of this strategy and comparison with unnormalised distributions is shown in Figure 1 The performance of each normalisation strategy was evaluated (see discussion) and was chosen as the more appropriate.

Training data
From Ensembl (Ensembl Mart version 30) [22] we extract all annotated Human RefSeq transcripts [23] for the relevant ten chromosomes. Transcripts with no introns are discarded. To obtain a one-to-one correspondence between annotation and sequence, overlapping annotation of expression is collapsed, forming a single set representing all annotated expression. The training set consists of RefSeq genes each represented by the corresponding sequence of genomically consecutive probe scores. Using the 5' position of probes on the genome each probe is labelled as expressed or non-expressed in accordance with the collapsed RefSeq annotation (see Figure 2). The training set includes 6411 sequences of labelled probe scores that each represent the transcribed part of a RefSeq gene. The entire set is used for training the hidden Markov model (HMM). For evaluation purposes a ten-fold crossvalidation is used.

HMM Architecture and training
The correspondence between scores and annotated expression, presented by the training set, is captured by an HMM, see e.g. Durbin et al. [24]. An HMM is a probabilistic model that consists of a set of connected hidden states that each emit observables. In this case, states represent expression and non-expression and these states emit scores that constitute the sequence of probe scores over expressed or non-expressed regions. First order emission probabilities are used to model the dependency of scores for neighbouring-probes.

A B
Labelling of probes Figure 2 Labelling of probes. Each probe is labelled as either expressed (red) or non-expressed (blue) based on expression annotation of the genomic position corresponding to the 5' end of the probe.

Expressed Non−expressed Non−expressed
The HMM is shown in Figure 3. Two loop states (E 0 and N 0 ) model the bulk of probe scores in expressed and nonexpressed regions, capturing the score distributions of the two cases. The loop states are connected by states that capture score gradients induced by genomically consecutive probes that overlap borders of expressed regions. The loop probabilities are tied during training of the HMM and thus estimated as one parameter. This is done in order to avoid bias of predictions towards typical exon length. The HMM used is discrete and expression scores are binned. Considering the distributions of scores for probes anno-tated as expressed and non-expressed, we choose 13 bins uniformly covering the scores from -50 to 50. Another two bins capture extreme positive and negative values, (see Figure 4).
The type of HMM used is a Class HMM [25]. The collection of states in a Class HMM is divided into classes, each designated by a label. The observables (i.e. probe scores) in the training set are each assigned a label corresponding to a class. An observable with a given label can only be emitted by a state belonging to the class designated by that label. Hence, the labelling of states and observables in the training set determines which parts of the training data are used to estimate which parts of the model. In our case, the labels are expression and non-expression and correspond to the red and blue colours in Figure 3. The advantage of the Class HMM is that it allows all model parameters to be estimated simultaneously from sequences of labelled probe scores. Alternatively, the individual parts would have to be estimated separately from expressed and nonexpressed probe scores. The model is trained using the Baum-Welch algorithm [26].

Self-supervised re-training
As noted above, the correspondence between labelling of probe scores and actual expression in the relevant cell line is at best an approximation. This issue is addressed by retraining after adding a parallel shadow model that mirrors the core model shown in Figure 3. The combined model is shown in Figure 5. Each state in the shadow model, , shares emission probabilities with a corresponding state, X i , in the core model but does not contribute to the estimation of these parameters. The labelling of states in the shadow model differs from those of the core model in that the states model non-expressed probes labelled as expressed and the state models expressed probes labelled as non-expressed. Additional states impose a minimum on the number of consecutive probe scores that can be captured by the shadow model. Without this constraint the model is not able to distinguish natural variation and noise from dubious annotation. The shadow component models regions annotated as expressed using the score distribution initially learned from nonexpressed regions. Analogously regions annotated as nonexpressed are modelled using the score distribution initially learned from expressed regions. Re-training after adding the shadow model allows the Baum-Welch algorithm to weight the contribution of training information to parameter estimation in proportion to how characteristic this is: In the E-step of the Baum-Welch algorithm the posterior probabilities of each emission from state X i are This figure shows how normalised expression scores are dis-cretised prior to HMM training and decoding    Core and shadow model  training, however, produces a starting point that allows the model to reach a desirable local optimum.
Once the parameters of the HMM have been estimated the parallel shadow architecture is removed. The sequences of probe scores corresponding to the tiling of each chromosome are then decoded using an N-best algorithm [25]. In this step the most likely sequences of contiguous expressed and non-expressed probes are established. In addition, the forward-backward algorithm [27] is used to assign a posterior probability of expression to each probe. The prediction for each sequence of probe scores is then mapped onto the nucleotide sequence of the corresponding chromosome.
The tiling of the chromosomes only covers non-repeat sequence. Some predictions may span untiled regions. To avoid this, gaps in the tiling larger than ten bases are removed from the predictions.

Evaluation of performance
To establish the relative performance of ExpressHMM and Affymetrix's normalisation and transfrag method (ANTM) we test each method on the training set of probes covering RefSeq genes. To ensure that test and training material do not overlap, the performance of ExpressHMM is evaluated using six-fold cross-validation. The ANTM predictions are mapped onto the training set for comparison.
To acknowledge the most recent repeat predictions we remove, from all genomic predictions, all overlaps to the Human Repeat Masker annotation from the UCSC genome browser database (October 2005) [28]. Performance statistics are calculated using the Eval program [29]. Sensitivity is defined as true positives divided by the sum of true positives and false negatives but is calculated as the number of annotation objects (e.g. nucleotides or sequence blocks) overlapped by predictions, divided by the total number of annotation objects. Specificity is defined as true positives divided by the sum of true positives and false positives but is calculated as the number of prediction objects overlapping annotation objects divided by the total number of prediction objects. The statistics are calculated this way to ensure a sensible treatment of oneto-many overlaps. Since the motivation in tiling array analysis is to find un-annotated expression a high specificity based on known expression is not a goal in itself. For this reason we report 1 -specificity as "Novel Prediction".
To further evaluate the performance we compare the two sets of predictions to annotation of Human ESTs, mRNAs, and known genes (UCSC Oct. 2005). Only annotation of tiled genomic sequence is included. Before evaluation the expression annotation is pooled and all overlapping annotation is collapsed into a maximal set where each nucleotide is only represented once.
Within both sets of annotation only a subset are actually expressed in the SK-N-AS cell line, and a substantial fraction of expression remains un-annotated. Collapsing all overlapping annotation presents a similar problem especially affecting splice site statistics. Relying on this annotation will make estimates of sensitivity artificially low and the absolute values of this statistic should not be given too much weight. The relative values for ExpressHMM and ANTM, however, should accurately reflect relative performances.

Results
The performance statistics obtained from the cross-validation of our model is shown in Table 1. Statistics of the ExpressHMM genomic predictions on the ten chromosomes are shown in Table 2. The overlap of ExpressHMM predictions to ANTM predictions is shown in Figure 6. The distribution of predictions and annotation is exemplified by chromosome 6 shown in Figure 7.
In addition to the most likely categorisation of probes ExpressHMM also calculates a posterior probability of expression for each probe. When mapped to the nucleotide sequence this results in an expression probability curve parallel to the predictions allowing the user to evaluate the significance of predictions and to directly view the evidence of expression.
Two example predictions, together with the ANTM predictions, are shown in Figure 8 and 9. Note the correspondence between ExpressHMM predictions and conservation in human, chimp, mouse, rat, dog, chicken, fugu, and zebrafish [30].
Overlap between ExpressHMM and ANTM transfrags pre-dictions Figure 6 Overlap between ExpressHMM and ANTM transfrags predictions. The diagram shows the percentage of ExpressHMM predictions that are overlapped by one or more ANTM predictions and vice versa. Due to the one-to-many relationships of overlaps overlap percentages corresponding to each set of predictions are given.

ExpressHMM ANTM
All predictions can be downloaded and viewed from our web site [2].

Evaluation of normalisation strategies
In the process of developing our method we compared the performance of the LL and normalisation strategies on the training set. Six-fold cross-validation was used to ensure no over-training. The nucleotide-level sensitivities of LL and was approximately 38.8% and 44.6% respectively with corresponding specificities of 51.1% and 20.7%. Subsequent optimisation of the model resulted in further improvement of these values (see Table 2). The LL based normalisation is significantly more specific than the approach but at a significant cost to sensitivity. We have used the approach in this study because we view sensitivity as a more desirable trait than specificity for our predictor.
A tentative explanation for these rather surprising differences in performance can be found by considering Figure  10. Observe the heavy right tail of the LL distribution for probes tiling coding regions. For the method both the left and right tails of the distribution are heavy for probes tiling coding regions. The HMM only requires significantly different distributions between expressed and nonexpressed regions in order to perform well. The two heavy tails generated by the method seems to supply more discriminative information to the HMM than the single heavy tail generated by the LL method. The effect of quantile normalisation was tested. This did not improve performance, and is not used for the results presented. Relative performances on the training set As shown in Table 1, our probe sensitivity is higher than ANTM. Of the total number of probes annotated as expressed we predict 21% more than ANTM. Note that the upper limit to sensitivity is not 100% as only a fraction of the annotated sequence is actually expressed. For identification of exons ExpressHMM is equal to ANTM but better at predicting the borders of expressed regions correctly. The percentage of predicted probes that do not overlap RefSeq exons is 16% larger than ANTM. In contrast, the number of predicted transfrags not overlapping a RefSeq exon comprise only 27% of the ExpressHMM predictions whereas this number is 58% for ANTM. This means that even with the higher probe sensitivity ExpressHMM is more conservative than ANTM in predicting novel transfrags.
Predictions of the known ZSWIM4 gene and downstream expression including a miRNA cluster on the same strand Figure 9 Predictions of the known ZSWIM4 gene and downstream expression including a miRNA cluster on the same strand. Note the correspondence between the conservation track and the ExpressHMM curve. The miRNA track is extracted from the UCSC genome browser (October 2005). It is generated by aligning miRNAs from the miRNA Registry [33] to the genome. The remaining tracks are explained in Figure 8. Predictions of the known Dicer gene involved in the processing of pre-miRNA to miRNA Figure 8 Predictions of the known Dicer gene involved in the processing of pre-miRNA to miRNA. The conservation curve shows the conservation in human, chimp, mouse, rat, dog, chicken, fugu, and zebrafish, based on a phylogenetic hidden Markov model [30]. Note how peaks in the ExpressHMM curve identify expression that is not sufficiently characteristic to result in a prediction. All information including the Repeat Masker track is extracted from the UCSC genome browser (October 2005). The HMM in ExpressHMM models the signal variability within expressed and non-expressed regions. As a result the algorithm does not predict very short expressed or non-expressed regions unless the relevant probe scores are highly characteristic. ANTM's threshold approach allows for little variability in score along a region. For these reasons ANTM tends to split exons into smaller predictions whilst ExpressHMM tends to join proximal exons into a single prediction. These differences are reflected in the length statistics given in Table 2. In effect, the two methods complement each other.

Relative performances on the chromosome level
The comparison of genomic predictions to all Human EST, mRNA and known gene annotation within the tiled regions also reflects the good performance of ExpressHMM. The nucleotide sensitivity is 75% larger than ANTM's and our transfrag sensitivity is 23% larger. ExpressHMM identifies an 19% larger fraction of novel nucleotides. In contrast, however, the fraction of novel transfrags among ExpressHMM predictions is 24% smaller than among ANTM predictions.
ExpressHMM predicts roughly twice as many expressed nucleotides as ANTM. Still the total ExpressHMM predictions only amount to 1.5% of the genome. For comparison, recent cDNA based findings in Mouse indicates that 62% percent of the genome is transcribed and that the number of transcripts is at least an order of magnitude larger than the estimated number of genes [31].
As indicated by the nucleotide plot in Figure 7 the lower nucleotide specificity of ExpressHMM is not due to false positives uniformly distributed across the chromosomes. Rather, the distribution closely follows the density of both ANTM predictions and RefSeq exons. Note the excess of prediction in the sub-telomeric regions. This is most likely due to frequent duplication events in these regions as observed in the recent Chimpanzee-Human genome comparison [32]. These duplications are expected to cause extensive cross-hybridisation. The two large peaks around 4e+07 correspond to similar-sized peaks in the density of EST and mRNA annotation (data not shown).
The transfrag plot in Figure 7 shows the close correspondence between the distributions of ExpressHMM predictions and occurrence of RefSeq exons. shows the analogous distribution for . Note the heavy tails for both these statistics in the exonic distribution. Before settling on the HMM described above, a variety of models were tested. These included different order Markov models with and without modelling of transfrag borders. None of these offered apparent advantages over the one presented here. The order of the HMM has a pronounced effect on the results. For zeroth and first order emissions the nucleotide statistics are very similar. On the transfrag level, however, the zeroth order model is more sensitive (72% identified RefSeq exons) but less specific (55% novel transfrags predicted). We also experimented with excluding score sequences from the training set that did not seem unambiguously expressed. This improved the performance of some models but not that presented here.

Hidden Markov model
The training of the model parameters does not include any nucleotide sequence information. In addition, to avoid any discriminative length modelling characteristic of coding exons the trained intron and exon length distributions are geometric and estimated as one to make them identical. As a result, the training data does not bias ExpressHMM predictions towards the length distribution of coding exons.
The expression probability curve that accompanies the ExpressHMM prediction has a direct interpretation as the evidence of expression given the model. It is a valuable tool to visually assess the significance of each transfrag as well as to identify regions with evidence not sufficient to result in a transfrag. With slight modifications ExpressHMM can equally well be used for the analysis of tiling micro arrays with non-overlapping probes of various length.

Conclusion
In addition to performing better than the approach presented by Cheng et al. [1], the adaptive approach used by ExpressHMM avoids ad hoc thresholds for the analysis of signal data. Decisions regarding prediction are learned from the data in an automated fashion. In addition to predicting the most likely categorisation into expression and non-expression each genomic position queried by a probe is assigned a probability of expression. The resulting graph has a clearer interpretation than a smoothed fluorescence score. We expect that tiling arrays will play an increasing role in the investigation of genomic output.