Feny? D. (Ed.) Computational Biology

Подождите немного. Документ загружается.

30 Miklós

The authors used numerical approaches to ﬁnd the tree topol-

ogy and edge lengths that maximize the probability of generating

the multiple alignment. Once the Maximum Likelihood tree has

been found, the most likely generation by the SCFG was calcu-

lated using the CYK algorithm. The CYK algorithm is also a

dynamic programming algorithm that ﬁnds for each substring

and nonterminal the most likely generation of the substring, start-

ing with the nonterminal. The running time of the CYK algo-

rithm is O(L

), where L is the length of the RNA string and

M is the number of nonterminal symbols. Since the number of

nonterminal symbols is ﬁxed, this algorithm – just like the

Nussinov algorithm and the Zuker–Sankoff algorithm – runs in

cubic time with the sequence length.

As also can be seen in Fig. 6, any generation of the Knudsen–

Hein grammar deﬁnes a secondary structure in an unequivocal

way, and this secondary structure is the prediction for the com-

mon secondary structure. It is also possible to calculate posterior

probabilities that a nucleic acid is single stranded or base-paired

with a particular partner nucleic acid. The posterior probability

for being single stranded means the conditional probability that

the character (or alignment column) was generated by the L → s

rewriting rule, with the condition that the SCFG generated the

sequence (or alignment). Similarly, the posterior probability that

the nucleic acids in positions i and j are base-paired is the condi-

tional probability that the nonterminal symbol F generated the

substring

+−

…

11ij

, with the condition that the SCFG generated

the sequence. Indeed, whenever an F appears on the right hand

side of the rewriting rules, it comes together with two ds, which

go to positions i and j if F generates the substring

+−

…

11ij

Knudsen and Hein showed that this posterior probability corre-

lates with the probability that the secondary structure prediction

is correct for that particular position (25).

The Knudsen–Hein grammar is very simple. Though it can

generate any pseudoknot-free secondary structure, the distribu-

tion of the structures it generates is far from the distribution that

we can ﬁnd in biological databases. Indeed, the grammar can

generate a run of single-stranded nucleic acids with a geometric

distribution, disregarding where these single-stranded nucleotides

are, in a hairpin, an internal loop, a bulge, or multiloop. A better

distribution can be achieved by increasing the number of non-

terminals. The different nonterminals are used for generating dif-

ferent structural elements such as hairpin loops, internal loops,

bulges, and multiloops. Nebel introduced such an extended

grammar and showed that such an extended grammar indeed

improves the goodness of predictions (28).

Covarion Models can be considered as the CFG equivalent of

proﬁle Hidden Markov Models (29–31). While proﬁle-HMMs

3.2. Covarion Models

RNA Structure Prediction

describe the proﬁle of a multiple sequence alignment of protein

sequences, Covarion Models describe the proﬁle of a multiple

sequence alignment of RNA sequences that belong to a fold fam-

ily. They cannot generate arbitrary secondary structures, and they

are speciﬁc to the common structure of the RNA sequences in the

multiple alignment. The most likely generations (also called most

likely parsings) align the RNA sequences together via their

Covarion Model: if two characters in two RNA sequences are

generated by the same nonterminal, then they are also aligned

together in the multiple alignment. Similarly, if two pairs of nucle-

otides are generated by the same nonterminal, the two pairs are

predicted to form base pairs, and both the left hand sides and the

right hand sides are aligned together.

Since Covarion Models are speciﬁc to an RNA family, they

are used to decide if a sequence belongs to a fold family. The dif-

ference between the energies of mfe secondary structures of a

randomly drawn RNA sequence and a functional RNA sequence

is not statistically sufﬁcient to ﬁnd novel RNA genes in genomic

sequences (32, 33). On the other hand, Covarion Models are

very successful in ﬁnding speciﬁc RNA genes. One of the most

successful applications is the tRNAScan-SE (34), which can detect

~99% of eukaryotic nuclear or prokaryotic tRNA genes, with a

false positive rate of less than one per 15 Gb, and with a search

speed of about 30 kb/s. This high-throughput method ﬁrst

quickly select a small part of the genome that contains basically all

tRNA genes using a simple Hidden Markov Model. Then, this

smaller part is analyzed further with a Covarion Model speciﬁc for

tRNA sequences.

The main drawback of the Knudsen–Hein method is that it needs

an accurate alignment to predict the common secondary struc-

ture of RNA sequences. Indeed, structure prediction for mis-

aligned parts is very hard, and it is impossible if base-paired

homologous nucleic acids are not aligned together. Covarion

Models give both an alignment and a predicted structure for each

sequence, however, they need a large set of homologous sequences

for an accurate prediction.

The Knudsen–Hein grammar and the Covarion Models can

predict only pseudoknot-free structures. Witwer introduced a

method that can predict planar pseudoknotted structures for

aligned RNA sequences (35). Her method is based on the

Maximum Weighted Matching (MWM) algorithm (36). The

MWM algorithm ﬁnds the set of base pairings that maximizes

the sum of the weights of base pairings, without any constraint.

Although it can predict arbitrary pseudoknots, the outcome

might also contain a subset of base pairings that clearly cannot be

formed due to steric constraints, see Fig. 7. Witwer suggested to

start with an MWM algorithm, and then remove base pairs to get

3.3. Joint Prediction

of Alignments

and Secondary

Structures

32 Miklós

a planar pseudoknot or pseudoknot-free structure; In a ﬁnal step,

the helices are extended, if possible. This algorithm also needs a

good initial alignment.

When the number of homologous sequences is small and the

sequences are hard to align without knowing their secondary

structure, no good initial alignment is available. An ideal approach

would coestimate the alignment and the secondary structure.

Sankoff suggested the problem to align and estimate the com-

mon structure of two RNA sequences, and he gave an algorithm

that runs in O(L

) running time, where L is the geometric mean

of the sequence lengths (37). Unfortunately, this algorithm is too

slow to be used in practice. Holmes and Rubin introduced the

SCFG approach for pairwise RNA structure comparison (38),

and Holmes introduced an evolutionary model for evolving RNA

structures (39). Although this latter method has been accelerated

(40), the acceleration is just a constant factor, and no theoretical

breakthrough has been achieved.

Meyer and Miklós introduced a Markov chain Monte Carlo

approach for the joint estimation of multiple alignments, evolu-

tionary trees, and RNA secondary structures including pseudo-

knots (41). The drawback of the approach is that the convergence

of the Markov chain might be quite slow, and it might take a lot

of computational time to draw a sufﬁcient number of samples

Fig. 7. A nonsense secondary structure that might be easily the outcome of a MWM

algorithm. The MWM algorithm maximizes the sum of weights of base pairs, disregard-

ing whether or not all base pairs can be formed. Clearly, the showed structure cannot

exist in real life due to steric constraints.

RNA Structure Prediction

from the Markov chain. The authors also provide software for the

analysis of the samples from the Markov chain. It is possible to

highlight the consensus structure of the sampled structures and

to give posterior probabilities for each base pair.

Finally, we mention the CARNAC method (42, 43) that tries

to ﬁnd a set of conserved helices in a family of homologous RNA

sequences. The set of conserved helices are not allowed to form a

pseudoknot. The approach is based on the k-clique problem,

which is known to be NP-hard. CARNAC has an implementation

that is reasonably fast on small datasets, and hence it provides a

practical approach.

There are both experimental (44–46) and statistical evidences

(47) that RNA sequences start folding during their translation.

Gultyaev published the ﬁrst work on simulating RNA folding

pathways (48, 49). Isambert and his colleagues implemented a

software package called KineFold that also simulates RNA folding

(50, 51). They also considered the possibility of cotranscriptional

folding, namely, the RNA sequence is being folded during its

transcription in the computer simulation. The authors also pro-

vide a graphical interface for visualizing the folding dynamics.

KineFold is also available on a webserver. Pseudoknotted struc-

tures are allowed in Kinefold.

So far, all implemented methods simulate the folding of a

single RNA sequence. Comparative approaches do not consider

directly the folding dynamics, but they try to ﬁnd common local

structures, for example, in mRNA sequences (52, 53).

In this paper, we have given an overview of approaches and

techniques for predicting RNA secondary structures. There are

three main approaches for predicting RNA secondary structures:

combinatorial optimization methods, comparative methods, and

folding simulations. The speed of algorithms based on combi-

natorial optimization depends on whether or not pseudoknots

are allowed, and if so, what kind of pseudoknots might be con-

sidered. Although pseudoknots are common in biological RNA

sequences (54], the general pseudoknot prediction problem is

hard. Planar pseudoknots can be predicted in polynomial time;

however, the O(L

) or O(L

) running time makes these algorithms

impractical.

4. Folding Kinetics

5. Discussion

34 Miklós

Comparative approaches are also common, and if there are

more than one sequence within the same secondary structure,

there is at least a hope that the common structure can be pre-

dicted with a better accuracy due to the increased amount of

information represented in the set of homologous sequences.

Kinetic approaches try to simulate the folding dynamics of

RNA sequences. Although these simulations are very impressive,

so far, no large-scale analysis has been published about the accu-

racy of such methods. Nevertheless, the kinetic approach has not

yet been combined with comparative approaches. We even do not

know how well the folding dynamics is conserved during

evolution.

Acknowledgements

The author thanks the following grants: Bolyai postdoctoral

fellowship, OTKA grant F61730. Caddy and Kanazawa are

thanked for their support.

References

1. Woese CR (1967) The Genetic Code. New

York, Evanston and London: Harper and

Row.

2. Crick FHC (1968) The origin of the genetic

code. J Mol Biol 38: 367.

3. Orgel LE (1968) Evolution of the genetic

apparatus. J Mol Biol 38: 381.

4. Kruger K, Grabowski PJ, Zaug AJ, Sands J,

Gottschling DE, Cech TR (1982) Self-splicing

RNA: Autoexcision and autocyclization of the

ribosomal RNA intervening sequence of

Tetrahymena. Cell 31: 147–157.

5. Nussinov R, Jacobson A (1980) Fast algo-

rithm for predicting the secondary structure

of single-stranded RNA. Proc Natl Acad Sci U

S A 77: 6309–6313.

6. Tinoco I Jr, Uhlenbeck OC, Levine MD

(1971) Estimation of secondary structure in

ribonucleic acids. Nature 230: 362–367.

7. Tinoco I Jr, Borer PN, Dengler B, Levine

MD, Uhlenbeck OC, et al (1973) Improved

estimation of secondary structure in ribo-

nucleic acids. Nature New Biol 246:

40–41.

8. Zuker M, Sankoff D (1984) RNA secondary

structures and their prediction. Bull Math

Biol 46: 591–621.

9. Wuchty S, Fontana W, Hofacker I, Schuster P

(1999) Complete suboptimal folding of RNA

and the stability of secondary structures.

Biopolymers 49: 145–165.

10. Zuker M (2003) Mfold web server for nucleic

acid folding and hybridization prediction.

Nucleic Acids Res 31: 3406–3415.

11. Miklós I, Meyer IM, Nagy B (2005) Moments

of the Boltzmann distribution for RNA sec-

ondary structures. Bull Math Biol 67(5):

1031–1047.

12. Mathews D, Sabina J, Zuker M, Turner D

(1999) Expanded sequence dependence of

thermodynamic parameters improves predic-

tion of RNA secondary structure. J Mol Biol

288: 911–940.

13. Hofacker IL (2003) Vienna RNA secondary

structure server. Nucleic Acids Res 31:

3429–3431.

14. Lyngsø R, Zuker M, Pedersen C (1999) Fast

evaluation of internal loops in RNA secondary

structure prediction. Bioinformatics 15:

440–445.

15. McCaskill JS (1990) The equilibrium parti-

tion function and base pair binding probabili-

ties for RNA secondary structure. Biopolymers

29: 1105–1119.

16. Rivas E, Eddy SR (1999) A dynamic program-

ming algorithm for RNA structure prediction

including pseudoknots. J Mol Biol 285:

2053–2068.

RNA Structure Prediction

17. Rivas E, Eddy SR (2000) The language of

RNA: A formal grammar that includes pseudo-

knots. Bioinformatics 16: 334–340.

18. Akutsu T (2000) Dynamic programming

algorithms for RNA secondary prediction

with pseudoknots. Discrete Appl Math 104:

45–62.

19. Lyngsø R, Pedersen C (2000) RNA pseudo-

knot prediction in energy based models.

J Comput Biol 7: 409–428.

20. Lyngsø R, Pedersen C (2000) Pseudoknots in

RNA secondary structures. In: Shamir R,

Miyano S, Istrail S, Pevzner P, Waterman M,

editors. Proceedings of the Fourth Annual

International Conference on Computational

Molecular Virology. New York: ACM Press.

pp. 201–209.

21. Dirks RM, Pierce NA (2003) A partition

function algorithm for nucleic acid secondary

structure including pseudoknots. J Comput

Chem 24: 1664–1677.

22. Reeder J, Giegerich R (2004) Design,

implementation and evaluation of a practi-

cal pseudoknot folding algorithm based on

thermodynamics. BMC Bioinformatics 5:

104.

23. Lyngsø R (2004) Complexity of pseudoknot

prediction in simple models. In: Diaz J,

Karhumäki J, Lepistö A, Sannella D, editors.

Proceedings of the 31st International

Colloquium on Automata, Languages, and

Programming (ICALP), 12–16 July 2004,

Turku, Finland. pp. 919–931.

24. Knudsen B, Hein J (1999) RNA secondary

structure prediction using stochastic context-

free grammars and evolutionary history.

Bioinformatics 15: 446–454.

25. Knudsen B, Hein J (2003) Pfold: RNA sec-

ondary structure prediction using stochastic

context-free grammars. Nucleic Acids Res 31:

3423–3428.

26. Durbin R, Eddy S, Krogh A, Mitchison G

(1998) Biological sequence analysis:

Probabilistic models of proteins and nucleic

acids. Cambridge: Cambridge University

Press. p. 356.

27. Felsenstein J (1981) Evolutionary trees from

DNA sequences: A maximum likelihood

approach. J Mol Evol 17(6): 368–376.

28. Nebel M (2004) Identifying good predictions

of RNA secondary structure. Proc Pac Symp

Biocomput 9: 423–434.

29. Eddy SR, Durbin R (1994) RNA sequence

analysis using covariance models. Nucleic

Acids Res 22: 2079–2088.

30. Hofacker IL, Fontana W, Stadler PF,

Bonhoeffer S, Tacker M, et al (1994) Fast

folding and comparison of RNA secondary

structures. Monatsh Chem 125: 167–188.

31. Sakakibara Y, Brown M, Underwood R, Mian

IS, Haussler D (1994) Stochastic context-free

grammars for modeling RNA. In: Proceedings

of the 27th Hawaii International Conference

on System Sciences. Honolulu: IEEE

Computer Society Press. pp. 283–284.

32. Rivas E, Eddy SR (2000) Secondary structure

alone is generally not statistically signiﬁcant

for the detection of noncoding RNAs.

Bioinformatics 16(7): 583–605.

33. Workman C, Krogh A (1999) No evidence

that mRNAs have lower folding free energies

than random sequences with the same dinu-

cleotide distribution. Nucleic Acids Res

27(24): 4816–4822.

34. Lowe T, Eddy S (1997) tRNAscan-SE: A pro-

gram for improved detection of transfer RNA

genes in genomic sequence. Nucleic Acids

Res 25: 955–964.

35. Witwer C (2003) Prediction of conserved and

consensus RNA structures [dissertation].

Vienna: Universität Wien. p. 187.

36. Tabaska J, Cary R, Gabow H, Stormo G

(1998) An RNA folding method capable of

identifying pseudoknots and base triples.

Bioinformatics 14: 691–699.

37. Sankoff D (1985) Simultaneous solution of

the RNA folding, alignment and protose-

quence problems. SIAM J Appl Math 45:

810–825.

38. Holmes I, Rubin G (2002) Pairwise RNA

structure comparison with stochastic context-

free grammars. Pac Symp Biocomput 2002:

163–174.

39. Holmes I (2004) A probabilistic model for

the evolution of RNA structure. BMC

Bioinformatics 5: 166.

40. Holmes I (2005) Accelerated probabilistic

inference of RNA structure evolution. BMC

Bioinformatics 6: 73.

41. Miklós I, Meyer IM (2007) SimulFold:

Simultaneously inferring RNA structures

including pseudoknots, alignments, and trees

using a Bayesian MCMC framework. PLoS

Comput Biol 3(8): e149.

42. Perriquet O, Touzet H, Dauchet M (2003)

Finding the common structure shared by two

homologous RNAs. Bioinformatics 19:

108–116.

43. Touzet H, Perriquet O (2004) CARNAC:

Folding families of related RNAs. Nucleic

Acids Res 32: W142–W145.

44. Boyle J, Robillard G, Kim S (1980) Sequential

folding of transfer RNA. A nuclear magnetic

resonance study of successively longer tRNA

36 Miklós

fragments with a common 59 end. J Mol Biol

139: 601–625.

45. Morgan SR, Higgs PG (1996) Evidence for

kinetic effects in the folding of large RNA

molecules. J Chem Phys 105: 7152–7157.

46. Heilmann-Miller SL, Woodson SA (2003)

Effect of transcription on folding of the

Tetrahymena ribozyme. RNA 9: 722–733.

47. Meyer IM, Miklós I (2004) Co-transcriptional

folding is encoded within RNA genes. BMC

Mol Biol 5: 10.

48. Gultyaev A (1991) The computer-simulation

of RNA folding involving pseudoknot forma-

tion. Nucleic Acids Res 19: 2489–2493.

49. Gultyaev A, von Batenburg F, Pleij C (1995)

The computer-simulation of RNA folding

pathways using a genetic algorithm. J Mol

Biol 250: 37–51.

50. Isambert H, Siggia E (2000) Modeling RNA

folding paths with pseudoknots: Application

to hepatitis delta virus ribozyme. Proc Natl

Acad Sci U S A 97: 6515–6520.

51. Xayaphoummine A, Bucher T, Thalmann F,

Isambert H (2003) Prediction and statistics of

pseudoknots in RNA structures using exactly

clustered stochastic simulations. Proc Natl

Acad Sci U S A 100: 15310–15315.

52. Pedersen JS, Forsberg R, Meyer IM, Hein J

(2004) An evolutionary model for protein-

coding regions with conserved RNA struc-

ture. Mol Biol Evol 21: 1913–1922.

53. Pedersen JS, Meyer IM, Forsberg R,

Simmonds P, Hein J (2004) A comparative

method for ﬁnding and folding RNA second-

ary structures within protein-coding regions.

Nucleic Acids Res 32: 4925–4936.

54. Staple DW, Butcher SE (2005) Pseudoknots:

RNA structures with diverse functions. PLoS

Biol 3: e213. doi:10.1371/journal.

pbio.0030213.

David Fenyö (ed.), Computational Biology, Methods in Molecular Biology, vol. 673,

Chapter 3

Normalization of Gene-Expression Microarray Data

Stefano Calza and Yudi Pawitan

Abstract

Expression microarrays are designed to quantify the amount of mRNA in a speciﬁc sample. However, this

can only be done indirectly through quantifying the color intensities returned by labeled mRNA mole-

cules bound to the array surface. Translating pixel intensities into transcript expression requires a series

of computations, generically known as preprocessing and normalization steps. In this chapter, we intro-

duce the basic concepts and methods, and illustrate them using data from three commonly used com-

mercial platforms.

Key words: Affymetrix, Agilent, Illumina, mRNA expression, Oligonucleotide arrays

Along with the intrinsic biological variability, microarray data will

show nonbiological or technical variability due to many potential

sources: slide fabrication, biological material extraction, quanti-

ties of mRNA hybridized, sample labeling, dye afﬁnity, scanner

settings, image acquisition conditions, spatial anomalies, etc. The

normalization step is the process that attempts to remove or

reduce these systematic technical biases among the samples. This

is a crucial step that can substantially affect the results of down-

stream statistical analyses.

In this chapter, we introduce some basic concepts, followed

by several normalization methods. The list is not exhaustive, but

rather represents a description of the most commonly-used

procedures. For practical illustrations, we use real data from a

comparison study (1) that employed the three currently most

commonly-used commercial platforms, namely, Affymetrix (2),

Agilent (3), and Illumina (4). Multiple use of technology plat-

forms on the same samples is ideal for illustration, since each has

1. Introduction

38 Calza and Pawitan

its own unique probe design and issues in the normalization step.

The usage of software packages available under the R statistical

environment (5) and the Bioconductor framework (6) will be

demonstrated.

The single slide or chip is ﬁrst read by a scanner, typically followed

by some image-analysis steps: gridding, segmentation, and inten-

sity extraction. Gridding is the process of assigning coordinates to

each of the probes/spots; segmentation classiﬁes the pixels as

foreground (belonging to a speciﬁc probe) or as background

(outside the probe); during extraction, foreground and back-

ground intensities are summarized for each probe. There could

be other image-analysis steps speciﬁc to each technology; we refer

to the manufacturer’s manual for details.

Background noise may derive from many sources such as target

binding to slide substrate, processing effects, and cross-hybridiza-

tion or nonspeciﬁc binding. Automated hybridization reduces

background noise, thanks to tighter control of temperature and

chemical conditions, as well as robot intervention. Nevertheless,

a certain amount of baseline noise is inevitable, and it must be

accounted for in the preprocessing.

When both background and foreground values are available

(e.g., in Agilent arrays), background correction can be performed

via subtraction of the background from the foreground values,

based on a simplistic assumption of an additive background effect.

As an undesired side effect, this method might generate a lot of

negative values that are problematic since most downstream anal-

yses are done in log scale. According to a recent comparison paper

(7), the best method is the so-called normexp, which is a modiﬁ-

cation of the Affymetrix background correction implemented in

the RMA method below.

Because of the tight probe arrangement on the chip, the

Affymetrix oligonucleotide array cannot use the pixels surround-

ing the probe/spot to estimate background levels. Therefore,

probe intensities have to be used to estimate both foreground

and background signals. Several procedures are available, though

there are two that are most commonly used. The ﬁrst one is

implemented in the so-called MAS5 algorithm (8). It works at

single-array basis and uses both perfect-match (PM) and mis-

match (MM) probe intensities. Two corrections are applied: (1) a

location-speciﬁc background adjustment, aimed at removing

overall background noise, (2) PM correction based on the MM

2. Preprocessing

2.1. Image Processing

2.2. Background

Correction

Normalization of Gene-Expression Microarray Data

value. For (1), each array is divided in 16 equally sized regions

arranged in 4-by-4 grids. Within region k, a background value b

is estimated using the mean of the lowest 2% of probe values, and

a noise level n

based on their standard deviation. Each probe

intensity P(x, y) is then adjusted using a weighted average for the

background values B(x, y), as well as for the noise values N(x, y),

as follows

(1)

where a is a tuning parameter (defaults to 0.5), and

(

)

(

)

∑

Bxy w xyb

, and the weight w

is a function of the

Euclidean distance of the probe location (x, y) and the centroid of

grid k.

Perfect-match correction based on MM value was originally

This turned out to be problematic as around 30% of MM values

are bigger than the corresponding PM values (9), thus causing

negative values. The remedy proposed by Affymetrix is based on

the computation of the ideal mismatch, a quantity that is smaller

than PM (10). It should be noted that the estimate of both the

foreground intensity and its local background will inevitably be

affected by measurement errors, so that subtracting one noisy

estimate from another is going to increase the overall variability

of the signal.

The second procedure is that used by the Robust Multichip

Analysis (RMA) algorithm (11, 12). It is based on the assumption

that the observed PM intensity X is the sum of a Gaussian com-

ponent for the real probe signal S and an exponential component

for the background B. In contrast to the MAS5 algorithm, the

MM intensities are not used. First, the mode of the distribution

of PM values is estimated. Then, points above the mode are used

to estimate the parameters of the exponential, while a half-Gaussian

is ﬁtted to those below the mode to estimate the Gaussian param-

eters. The ﬁnal background-corrected intensities are computed as

the conditional mean E(S|X ), which is the best predicted value of

the signal S given the observed data X.

Most normalization procedures rely on rather strong assumptions:

(1) the great majority of genes are unaltered between arrays,

and (2) genes are expected to be roughly symmetrically distrib-

uted among the upregulated and the downregulated. These

assumptions seem reasonable in many lab studies, where an

intervention is not expected to have extensive changes in global

=−( , ) max{ ( , ) ( , ); ( , )}.Pxy Pxy Bxy Nxya

3. Different

Methods of

Normalization