Genome Annotation

---

Genome Annotation

David S. Wishart

Introduction

Thanks to rapid advances in DNA sequencing technology and DNA analysis software, genome

projects that used to take years and cost millions of dollars to finish can now be completed in

just weeks at a cost of a few thousand dollars. The typical workflow for a modern genome

sequencing project involves performing whole genome DNA sequencing of selected organ-

isms using a next generation DNA sequencer, running a variety of programs to assemble a

reference genome, and using software to locate and identify all of the protein-coding ribo-

somal RNA (rRNA) and transfer RNA (tRNA) genes within the genomic sequence. This last

process is called genome annotation and it is the primary subject of this chapter. Strictly speak-

ing, genome annotation is not genome prediction. Gene or genome prediction is a subfield of

genome annotation. In particular, gene prediction uses mathematical or probabilistic models

to analyze DNA sequences and to identify gene boundaries and gene structures. On the other

hand, genome annotation uses gene (and genome) prediction results along with other lines of

evidence such as gene expression data, protein expression data, sequence homology to other

annotated genomes, and even literature assessments to generate a set of genome annotations.

These annotations include not only the location of the genes on each chromosome but also

their names (based on homology), calculated properties (such as sequence length, amino acid

composition, and molecular weight), expression levels (if available), and probable functions.

Depending on the type of organism that has been sequenced, the task of genome annotation

can be either quite easy or quite difficult. Prokaryotes (including bacteria and archaea) have

relatively small genomes, typically no more than 5 million base pairs, consisting of one or

two circular chromosomes and perhaps one or two small plasmids. The gene structure for

prokaryotes is very simple, with each gene being a contiguous open reading frame (ORF).

Furthermore, the coding density for prokaryotes is very high, with at least 85–90% of their

DNA coding for proteins, tRNAs, and rRNAs (Hou and Lin 2009). This makes the identification

of genes in prokaryotes relatively simple. On the other hand, eukaryotic gene identification

is often quite difficult. This is because eukaryotes have very large genomes (often billions of

base pairs), with very low (often <2%) coding densities (Hou and Lin 2009). Eukaryotic gene

structure is also much more complex than prokaryotic gene structure. In particular, eukaryotic

genes are split into exons and introns, and most eukaryotic genes are separated by very large

stretches of non-coding DNA (called intergenic regions).

While the cellular machinery in eukaryotic cells is able to recognize and process gene sig-

nals with remarkable accuracy and precision, our understanding of the molecular mechanisms

by which eukaryotic sequence signals are recognized and processed remains incomplete. As a

result, currently available eukaryotic gene prediction methods are not very accurate. Therefore,

in the absence of additional experimental or extrinsic information (e.g. gene expression data),

one should assume that eukaryotic gene predictions are only approximate. Even with consid-

erable experimental data at hand, it is still quite difficult to fully annotate the best-studied

Bioinformatics, Fourth Edition. Edited by Andreas D. Baxevanis, Gary D. Bader, and David S. Wishart.

© 2020 John Wiley & Sons, Inc. Published 2020 by John Wiley & Sons, Inc.

Companion Website: www.wiley.com/go/baxevanis/Bioinformatics_4e

Genome Annotation

eukaryotic genomes. For instance, the DNA sequence for the human genome has been known

since 2001, but the actual number of genes encoded within our own genome has still not been

fully determined (Pennisi 2003; Ezkurdia et al. 2014).

This chapter briefly reviews some of the computational methods and algorithms underly-

ing computational gene prediction for both prokaryotes and eukaryotes. It also describes how

experimental evidence and database comparisons can be integrated into these gene prediction

tools to improve gene prediction performance and to ensure more complete genome annota-

tion. Methods for assessing the performance of computational gene finders are also described.

Finally, a number of genome annotation pipelines are highlighted, along with several tools for

visualizing the resulting annotations.

Gene Prediction Methods

Different methods for gene prediction have been developed separately for prokaryotes and

eukaryotes because of important differences in their overall gene organization. Gene-finding

programs, whether for prokaryotes or eukaryotes, fall into two general categories: intrinsic

(or ab initio) gene predictors and extrinsic (or evidence based) gene finders (Borodovsky et al.

1994).

Ab initio gene prediction approaches attempt to predict and annotate genes solely using

DNA sequence data as input and without direct comparison with other sequences or sequence

databases. Ab initio approaches involve searching for sequence signals that are potentially

involved in gene specification and/or looking for regions that show compositional bias that

has been correlated with coding regions. This combined approach to gene finding is called

searching by signal and searching by content. GeneMark (Borodovsky and McIninch 1993),

GLIMMER (Delcher et al. 1999, 2007), EasyGene (Larsen and Krogh 2003), and GENSCAN

(Burge and Karlin 1997) are well-known examples of intrinsic or ab initio gene-finding pro-

grams. In contrast, extrinsic gene-finding methods involve both homology-based and compar-

ative approaches, in which the gene structure is determined through comparison with other

sequences whose characteristics are already known. BLASTX is an example of an extrinsic

gene-finding program that has been frequently applied for gene identification in prokaryotic

genomes (Borodovsky et al. 1994). Extrinsic gene prediction methods depend on having exper-

imental evidence (such as messenger RNA (mRNA) or RNA-seq data) and/or a large body of

pre-existing experimental sequencing data to perform sequence comparisons and gene identi-

fications. We will discuss these extrinsic methods and their role in genome annotation a little

later in this chapter. To begin with, we will focus on the intrinsic or ab initio gene prediction

methods.

<!-- END LEGACY FRAGMENT: 01_Introduction_and_Gene_Prediction_Methods -->

---

Ab Initio Gene Prediction in Prokaryotic Genomes

A prokaryotic gene typically begins with a start codon (e.g. ATG), ends with one of three stop

codons (e.g. TAG, TAA, or TGA), and is usually at least 100 bases long (Figure 5.1). These

protein-coding genes are called ORFs. Most of the genes in prokaryotic genomes are organized

into operons, which are gene clusters consisting of more than one ORF that are under the con-

trol of a shared set of regulatory sequences. These regulatory sequences can include enhancers,

silencers, terminators, operators, or promoters. Regulatory sequences typically constitute the

10–15% of the prokaryotic genome that is not coding for protein sequences. A prokaryotic

gene promoter is a small segment of DNA that initiates transcription of a particular gene. Pro-

moters are located near the transcription start sites (TSSs) of genes, on the same strand and

upstream of the gene or ORF. In prokaryotes, the promoter contains two short sequence ele-

ments approximately 10 bases and 35 nucleotides upstream from the TSS. The element located

10 bases upstream is called the TATA box in archaea or the Pribnow (TATAAT) box in bacteria

(Pribnow 1975). These abbreviations or letters actually indicate the consensus DNA sequences

Ab Initio Gene Prediction in Prokaryotic Genomes

ATGACAGATTACAGA......TGCAGTTACAGGATAG

TATA box

Start codon

Stop codon

ORF

Figure 5.1 A simplified depiction of a prokaryotic gene or open reading frame (ORF) including the start

codon (or translation initiation site), the stop codon (TAG), and the TATA or Pribnow box.

seen for these regions. In addition to the TSS, almost all prokaryotic genes have a ribosome

binding site (RBS) that is 8–10 bases upstream of the start (ATG) codon. The start codon is

also called the translation initiation site (TIS). The RBS exhibits a specific nucleotide pattern

(AGGAGG) called a Shine–Dalgarno (SD) consensus sequence (Shine and Dalgarno 1975).

The SD sequence enables interactions between mRNA and the cell’s translational machin-

ery. In bacteria and archaea, translation initiation is generally thought to occur through the

base-pairing interaction between the 3′ tail of the 16S rRNA of the 30S ribosomal subunit and

the site in the 5′ untranslated region (UTR) of an mRNA that carries the SD consensus.

Consensus sequences, while providing a useful reminder or mnemonic, are never really

used in modern gene signal or gene site (i.e. TIS, RBS, TSS, and terminator) identification.

Instead, most gene signals can be identified by using positional weight matrices (PWMs) or

position-specific scoring matrices (PSSMs; see also Chapter 3). These scoring matrices are cal-

culated from carefully aligning a set of known functional signals and determining the adjusted

frequency with which specific bases may appear in certain positions. An example of how to

calculate a PSSM is given in Box 5.1. Once calculated for a given signal, signal-specific PSSMs

can be used to rapidly compute, along the length of a sequence of interest, the position and

likelihood of the selected gene signals. A simplified gene prediction protocol for prokaryotes

involves the following steps.

• Start at the beginning of the genome sequence at the 5′ end of one DNA strand and find an

ATG start codon that makes the longest ORF (minimum 150 bases), then move to the next

ATG downstream of the previously identified ORF and repeat the process for the rest of the

genome sequence.

• Repeat the above process for the opposite DNA strand.

• For all identified ORFs, score the quality of the TSS and RBS signals using site-specific

PSSMs to refine the ORF predictions and produce a final list of genes.

Box 5.1 Position-Specific Scoring Matrices

Position-specific scoring matrices (PSSMs), which are also called positional weight matri-

ces (PWMs) or positional specific weight matrices (PSWMs), are usually derived from a set

of aligned sequences that are believed to be functionally related. In this example, five dif-

ferent DNA sequences consisting of 10 bases each, which are believed to be functionally

related (as promoter regions), are aligned.

A T T T A G T A T C

G T T C T G T A A C

A T T T T G T A G C

A A G C T G T A A C

C A T T T G T A C A

From this alignment, a simple positional frequency matrix (PFM) can be generated.

In this matrix the frequency of the As, Cs, Gs, and Ts is tabulated (based on the above

alignment) for each of the 10 base positions. So in the first position there are three As,

one C, one G, and no Ts (see column 1). The PFM for the above alignment is:

(Continued)

Genome Annotation

Box 5.1 (Continued)

A 3 2 0 0 1 0 0 5 2 1

C 1 0 0 2 0 0 0 0 1 4

G 1 0 1 0 0 5 0 0 1 0

T 0 3 4 3 4 0 5 0 1 0

The PFM can now be converted to a positional probability matrix (PPM). A PPM is a

matrix consisting of a set of decimal values based on the percentage or frequency of

occurrences of each base in each position in the sequence alignment. In other words, we

must normalize the frequencies by dividing the nucleotide count at each position by the

number of sequences in the alignment. So if there are five sequences in the alignment

and three As in the first position, then the positional probability for A in the first position

is 3/5 = 0.6. Likewise, if there is one C in the first position, its positional probability is

1/5 = 0.2. One G corresponds to a positional probability of 0.2, and no Ts corresponds to a

positional probability of 0 (see column 1). Performing this same calculation across all 10

positions of the alignment, the full PPM would appear as follows:

A .6 .4

0 .2 0 0 1 .4 .2

C .2

0 .4

0 0 0 0 .2 .8

G .2

0 .2

0 1 0 0 .2

T

0 .6 .8 .6 .8 0 1 0 .2

The probabilities in the above PPM can be multiplied together to calculate the probabil-

ity that a given DNA sequence is closely related to the original five sequences. For instance,

if we wanted to know if the new sequence ATTTTGTATA is closely related, we could mul-

tiply the values for each sequence position to calculate that sequence’s probability:

p = 0.6 × 0.6 × 0.8 × 0.6 × 0.8 × 1 × 1 × 1 × 0.2 × 0.2 = 0.0055

Note that if we had performed this same calculation on an almost identical sequence

such as ACTTTGTATA (which differs by only one base) we would get p = 0. We get a 0 prob-

ability because C was not observed in the second position of our training set. Building a

PPM with only five sequences means you are very likely to underestimate (or overesti-

mate) the true fractional frequencies of each base, leading to problems in calculating

probabilities similar to what we just saw. To account for the small size of our multiple

sequence alignment (MSA) we should introduce pseudocounts. Pseudocounts are used to

avoid issues that result from matrix entries having a value of 0. Pseudocounting is equiv-

alent to multiplying each column of the PPM by a Dirichlet distribution, thereby allowing

the probability to be calculated for the “unseen” or unused sequences. A simple way of

doing this is to normalize the data to match the overall base composition of the genome(s)

being considered and to add a correction factor that scales as the square root of the num-

ber of sequences in the MSA. Hence, the following formula can be used to rescore each

base position in the PPM:

score (Xi) = (Qx + Px)∕(N + B)

where Qx is the number of counts of base type X at position i, Px is the number of

pseudocounts of base type X, which is equal to B × the frequency of base type X, N

is the total number of sequences in the MSA, and B is the number of pseudocounts

(assumed to be √N). For the genome or genomes of interest the frequency of As is 0.32,

Ts is 0.32, Cs is 0.18, and Gs is 0.18. Using this information the value for A in the first

position is (3 + (√5 × 0.32))/(5 + √5) = 0.51. The value for C in the second position is

(1 + (√5 × 0.18))/(5 + √5) = 0.19, and so on. The pseudocount corrected PPM is now:

A .51 .38 .09 .09 .24 .09 .09 .79 .38 .24

C .19 .06 .06 .33 .06 .06 .06 .06 .19 .61

G .19 .06 .19 .06 .06 .75 .06 .06 .19 .06

T .09 .51 .65 .51 .65 .09 .79 .09 .24 .09

Ab Initio Gene Prediction in Prokaryotic Genomes

Ideally each of the columns should sum to 1 but, because of rounding, the sums in

this example are sometimes slightly above or below 1. With this rescored matrix, you will

notice that there are now no zero entries. However, the calculation of probabilities through

multiplication is tedious (given the number of significant digits) and difficult. A simpler

way is to convert the PPM to a different type of matrix by taking the negative log10 of each

number in the PPM. This converts two-digit decimals to single-digit decimals and it also

allows one to add rather than multiply to calculate probabilities. If we take the −log10 of

the above PPM, we get:

A 0.3 0.4 1.0 1.0 0.6 1.0 1.0 0.1 0.4 0.6

C 0.7 1.2 1.2 0.5 1.2 1.2 1.2 1.2 0.7 0.2

G 0.7 1.2 0.7 1.2 1.2 0.1 1.2 1.2 0.7 1.2

T 1.1 0.3 0.2 0.3 0.2 1.0 0.1 1.0 0.6 1.0

This modified matrix is called a log likelihood scoring matrix or a PSSM. Using the above

PSSM, we can now calculate the score (or the log likelihood) for the query sequence

ATTTTGTATA: 0.3 + 0.3 + 0.2 + 0.3 + 0.2 + 0.1 + 0.1 + 0.1 + 0.6 + 0.6 = 2.8. The sequence

score gives an indication of how different the sequence is from a random sequence. The

higher the score, the more likely the sequence is a promoter/functional site and not a

random sequence. A score of 2.8 is very high. The sequence score can also be interpreted

in terms of the binding energy for that sequence.

However, such a simplified algorithm would only likely be 75–80% correct (Besemer et al.

2001). This is because prokaryotic genes are not always so simple to identify. For instance, the

ATG start codon is not always used for all bacterial genes. Among the 4284 genes identified

in Escherichia coli, 83% use ATG, 14% use GTG and 3% use TTG start codons (Blattner et al.

1997). Likewise, using a simple rule to identify only long ORFs may miss many short ORFs

or misidentify ORFs that have an unusual codon bias (indicating they are unlikely to code for

a gene). Indeed, the length distributions of ORFs known to code for proteins compared with

ORFs that occur by chance differ quite significantly. More specifically, coding ORFs have a

length distribution that resembles the gamma distribution (see Glossary), while non-coding

ORFs have a length distribution that resembles a simple exponential function (Lukashin and

Borodovsky 1998). In addition to these complications, it has recently been found that certain

prokaryotic genes have very unusual gene start signals because of a phenomenon called lead-

erless transcription (Slupska et al. 2001). In leaderless transcription, RNA transcripts have very

short 5′ UTRs, with a length < 6 bases. These regions are so short that they are unable to host

the RBS. This places the TSS at or very near to the TIS. In these cases, the promoter signal has

to be used for more accurate TIS identification.

Given the variations in the length and character of many prokaryotic gene signals, PSSMs are

not the most effective signal recognition tools available. More advanced methods of gene signal

recognition exist, such as Markov models (Box 5.2), hidden Markov models or HMMs (Box 5.3),

artificial neural networks, and support vector machines. These machine learning methods do a

far better job of handling variable lengths and conditional sequence dependencies that cannot

be captured with simple PSSMs.

Box 5.2 Markov Models

A Markov chain, model, or process refers to a series of observations in which the prob-

ability of an observation depends on a number of previous observations. The number of

observations defines the “order” of the chain. For example, in a first-order Markov model,

the probability of an observation depends only on the previous observation. In a Markov

chain of order 5, the probability of an observation depends on the five preceding obser-

vations. A DNA sequence can be considered to be an example of a Markov model because

the likelihood of observing a particular base at a given position may depend on the bases

(Continued)

Genome Annotation

Box 5.2 (Continued)

preceding it. In particular, in coding regions, it is well known that the probability of a

given base depends on the five preceding bases, reflecting observed codon biases and

dependencies between adjacent codons. In non-coding regions, such dependence is not

observed. When scanning an anonymous genomic region, one can compute how well the

local nucleotide sequence conforms to the fifth-order dependencies observed in coding

regions and assign appropriate coding likelihood scores.

Box 5.3 Hidden Markov Models in Gene Prediction

Hidden Markov models (HMMs) are used to provide a statistical representation of real

biological processes. They have found widespread use in many areas of bioinformatics,

including multiple sequence alignment, the characterization and classification of protein

families, the comparison of protein structures, and the prediction of gene structure.

In this chapter, all of the gene-finding methods that are described have two things in

common: they use a raw nucleotide sequence as their input and, for each position in the

sequence, they attempt to predict whether a given base is most likely found in an intron,

an exon, or within an intergenic region. In making these predictions, the algorithm applied

(HMM or otherwise) must take into account what is known about the structure of a gene,

showed in a simplified fashion in Figure 5.2.

Working from the 5′ to 3′ end of the gene, the method must take into account the

unique characteristics of promoter regions, transcription start sites, 5′ UTRs, start codons,

exons, splice donors, introns, splice acceptors, stop codons, 3′ UTRs, and polyA tails. In

addition to any conserved sequences or compositional bias that may characterize each of

these regions (Box 5.1), the method also needs to take into account that each of these

elements appears with a controlled syntax; for example, the promoter (and its TATA box)

must appear before the start codon, an initial exon must follow the start codon, introns

must follow exons, introns can only be followed by internal or terminal exons, stop codons

cannot interrupt the coding region, and polyA signals must appear after the stop codon.

Finally, an ORF must be maintained throughout to produce a protein once all is said and

done.

Each of the elements – exons, introns, and so forth – are referred to as states. The

sequence characteristics and syntactical constraints described above allow a transition

probability to be assigned, indicating how likely a change of state is as one moves through

Transcribed region

Exon 1

Exon 2

Exon 3

Intron 1

Intron 2

Start codon

5′ UTR

3′ UTR

Stop codon

Downstream

intergenic

region

Upstream

intergenic

region

Figure 5.2 A simplified depiction of a eukaryotic gene illustrating the multi-intron/exon structure,

the location of the start and stop codons, the untranslated regions (UTRs), and the intergenic

regions that surround the transcribed gene.

<!-- END LEGACY FRAGMENT: 02_Ab_Initio_Gene_Prediction_in_Prokaryotic_Genomes -->

<!-- BEGIN LEGACY FRAGMENT: 03_Ab_Initio_Gene_Prediction_in_Eukaryotic_Genomes -->

Source: Ch5 Genome Annotation / Ab Initio Gene Prediction in Eukaryotic Genomes

PDF Pages: 143-147 | Print Pages: 123-127

Boundary: starts at second true section title on PDF page 143; includes Figure 5.3-Figure 5.5 and stops before How Well Do Gene Predictors Work?

Note: PDF page 143 top contains the tail of Box 5.3 from previous section; excluded here.

---

Ab Initio Gene Prediction in Eukaryotic Genomes

A diagram of how eukaryotic genes are organized is shown in Figure 5.2. As can be seen from

this figure, eukaryotic genes are somewhat more complex than prokaryotic genes. In partic-

ular, the density of protein-coding regions for eukaryotic genomes (and especially vertebrate

genomes) is 90–100 times lower than it is for prokaryotic genomes. These sparse protein-coding

regions are separated by long stretches of intergenic DNA while their coding sequences (the

exons) are interrupted by large, non-coding introns. Genes are recognized and transcribed by

eukaryotic RNA polymerases, and the resulting long RNA transcripts are then cut by various

small ribonuclear proteins (snRNPs) to remove the introns (Will and Lührmann 2011). The

remaining exons are then spliced together to form the much smaller protein-coding transcript.

The snRNPs recognize specific cut sites at the exon/intron junctions to ensure that the splicing

is always performed precisely.

In the human genome, just 1.1% of the genome is composed of exons, 24% is composed of

introns, and 75% of the genome constitutes intergenic DNA. On average there are 5.48 exons

per gene with each exon encoding a peptide fragment of 30–36 amino acids (Sakharkar et al.

2002). The longest exon in the human genome is 11 555 bases, while the shortest is just two

bases long (Sakharkar et al. 2002). Not only are exons “rare,” they vary tremendously in length.

What is more, they can be alternately spliced to produce very different combinations of final

gene (transcript) products. This makes gene prediction significantly more difficult for eukary-

otes than for prokaryotes.

Computational gene prediction for eukaryotes essentially involves mimicking the biological

transcriptional and splicing process. In the biological process, various proteins and protein

complexes within the cell scan through the DNA sequence, recognize and bind to specific DNA

Genome Annotation

sites, transcribe the gene, and then cut and splice the transcript to form a final gene product.

In the computational process, the proteins are replaced with various algorithms that:

• identify and score suitable splice sites and start and stop signals along the query sequence

• determine the location of the candidate exons, as deduced through the detection of these

signals

• score and identify the best exons as a function of both the signals used to detect the exons

as well as the coding statistics computed on the putative exon sequence itself

• assemble (or “splice”) a subset of these exon candidates into a predicted gene structure. The

assembly is produced in a way that maximizes a particular scoring function that is dependent

on the score of each of the individual exon candidates.

The way in which each of these tasks is actually implemented varies from program to pro-

gram. Rather than discuss each program in detail, we will describe the three major processes

common to almost all ab initio eukaryotic gene prediction programs: predicting exon-defining

signals, predicting and scoring exons, and, finally, exon assembly.

Predicting Exon-Defining Signals

Just as prokaryotic genes have DNA signals, eukaryotic genes have distinct DNA signals

as well. Some of these elements are similar to prokaryotes, while others are quite different

(Figure 5.3). For instance, many eukaryotic genes have promoter elements that also share

some sequence similarity to prokaryotic genes. The most extensively studied core promoter

element in eukaryotes is known as the TATA box or the Goldberg–Hogness box (Lifton

et al. 1978), found 25–30 base pairs upstream from the TSS. The TATA box is also found

in archaea and bacteria and it appears to be a very ancient DNA signal. The TATA box in

eukaryotes has the consensus sequence TATA(A/T)A(A/T) and is often coupled to another

regulatory sequence called the CCAAT box (consensus: GGCCAATCT), located ∼150 base

pairs upstream of the TATA box. Only about 25–35% of mammalian genes contain TATA

boxes, while the rest contain other kinds of core promoter elements. Eukaryotic genes also

contain regulatory sequences beyond the core promoter, including enhancers, silencers,

and insulators. These regulatory sequences can be spread over a large genomic distance,

often hundreds of kilobases from the core promoters. In addition to having a wide variety of

promoter or enhancer signals, eukaryotic genes also have very specific DNA signals to define

the location of exons and introns.

More specifically, there are four basic DNA signals involved in defining exons: the TIS, the

5′ (or donor) splice site, the 3′ (or acceptor) splice site, and the translational stop codon. In

eukaryotes, the TIS is defined by the Kozak consensus sequence, often given as ACCATGG

(Kozak 1987), where the central ATG is the start codon. The 5′ donor splice site is typically

defined by a consensus sequence given as GG/GT, while the 3′ acceptor splice site has a con-

sensus sequence of CAG/G, where the slash indicates the cut sites for splicing (Figure 5.4). The

translational stop codons include the usual TAG, TAA, or TGA.

GC box

~200 bp

CCAAT box

~100 bp

TATA box

~30 bp

Gene

Transcription

start site

Exon

Exon

Intron

Figure 5.3 A schematic illustration of the upstream regions of a eukaryotic gene with the GC box located

∼200 bp upstream, the CCAAT box located ∼100 bp upstream, and the TATA box located ∼30 bp upstream

of the transcription start site.

Ab Initio Gene Prediction in Eukaryotic Genomes

Exon 1

Exon 2

Intron 1

Intron 2

Branchpoint site

5′ site

3′ site

AG/GT

CAG/NT

Figure 5.4 A schematic illustration of the splice site regions around exons and introns including the 5′

and 3′ splice sites and their consensus sequences.

The first methods used to identify exon-defining signals were simple PWMs or PSSMs.

These proved to be rather poor at identifying short DNA signals, such as splice sites. As a

result, these simple models have since given way to much more advanced pattern recognition

techniques such as HMMs (Box 5.3). These powerful pattern recognition approaches allow

very complex sequence patterns to be “learned” from large datasets consisting of well-known

or well-annotated exon-defining signals. An HMM is a statistical Markov model in which

the system being modeled is assumed to be a Markov process with unobserved (i.e. hidden)

states. HMMs are commonly used in many real-life applications such as speech, handwriting,

and gesture recognition. The application of HMMs in bioinformatics began in the early 1990s

(Krogh et al. 1994) and led to a significant advance in gene prediction accuracy. HMMs make it

possible to define highly complex patterns of variable lengths, including many exon-defining

signals such as protein-coding regions (discussed below), donor, acceptor, and lariat sites, as

well as translational start and end sites.

Predicting and Scoring Exons

In addition to the identification of exon-defining signals, the accurate prediction of exons also

depends on content-based features. Exons can be divided into three basic types:

• initial exons: ORFs delimited by a start site and a donor site

• internal exons: ORFs delimited by a 5′ (donor) site and a 3′ (acceptor) site

• terminal exons: ORFs delimited by a 3′ (acceptor) site and a stop codon.

Most transcribed genes are composed of one initial exon, multiple internal exons, and a

single terminal exon. Zhang (2002) provides a more comprehensive discussion of these types

of eukaryotic exons.

Exons, by definition, are protein-coding regions. Protein-coding regions are known to

exhibit characteristic compositional bias when compared with non-coding regions. These

include somewhat richer GC content and a distinctly non-random codon (triplet) frequency

preference. The observed codon bias results from the uneven distribution of amino acids

in proteins, the uneven use of synonymous codons, and natural selection for translational

optimization in coding regions. To discriminate protein-coding regions from non-coding

regions, a number of DNA content-based measures were developed in the 1990s (Fickett and

Tung 1992; Gelfand 1995; Guigó 1999). These content measures, which are also referred to

as coding statistics, reflect the likelihood that a given DNA sequence codes for a protein or

protein fragment. Many methods for the computation of content-based measures have been

published over the years. Some of the first methods measured patterns seen in codon triplet

frequencies. However, more information was found in the frequencies of pairs of triplets (i.e.

hexamers). As a result, hexamer frequencies, usually in the form of codon position-dependent

fifth-order Markov models (Box 5.2; Borodovsky and McIninch 1993), seem to offer the best

Genome Annotation

discriminative power to identify protein-coding regions in exons. Currently, these hexamer

frequencies lie at the core of all modern eukaryotic gene predictors.

Exon Assembly

Once the exons are predicted (using a combination of hexamer frequencies and HMMs to iden-

tify key gene signals and exon/intron boundaries), they need to be assembled into some sort of

multi-exon gene structure. The main difficulty in exon assembly lies in simple combinatorics:

the number of possible exon assemblies grows exponentially with the number of predicted

exons for any given gene. To address this problem, a number of dynamic programming tech-

niques have been developed. Dynamic programming is an optimization technique that allows

one to solve a complex problem by breaking it down into a collection of simpler subproblems.

Each of those subproblems is solved just once, and their solutions are stored. The next time

the same subproblem occurs, instead of recomputing its solution, one simply looks up the pre-

viously computed solution (Bellman 1957; see also Appendix 6.A for a detailed discussion).

For the optimal exon assembly problem, dynamic programming has been shown to find the

solution quite efficiently, without having to enumerate or consider each and every possible

combination of exons (Gelfand and Roytberg 1993). Nearly all modern eukaryotic gene predic-

tion tools now use some kind of dynamic programming method (called the Viterbi algorithm

by Markov modelers, but also known as the Needleman–Wunsch algorithm by most people

doing sequence alignment). By combining HMM-based exon signal identification with dif-

ferent HMM-derived scores for exons and then using dynamic programing to assemble the

exons, it is possible to generate robust eukaryotic gene predictions. Some early examples of

HMM-based gene prediction methods that use dynamic programming include GENIE (Kulp

et al. 1996) and HMMgene (Krogh 1997). Perhaps the most popular example of an HMM-based

eukaryotic gene predictor is GENSCAN (Burge and Karlin 1997), an ab initio gene predictor

that has been widely used to annotate hundreds of eukaryotic genomes.

Given the popularity of GENSCAN, it is perhaps worthwhile explaining how this program

works in a bit more detail and providing an example of how it can be used. For any given query

sequence, GENSCAN determines the most likely gene structure given an underlying HMM. To

model donor splice sites, GENSCAN introduced a method called maximal dependence decom-

position. In this method, a series of weight matrices (instead of just one) are used to capture

dependencies between positions in these splice sites. In addition, GENSCAN uses parameters

that account for many higher order properties of genomic sequences (e.g. typical gene den-

sity, typical number of exons per gene, and the distribution of exon sizes for different types of

exons). Separate sets of gene model parameters can be used to adjust for the differences in gene

density and G + C composition seen across genomes. Models have also been developed for use

with maize and Arabidopsis sequences. This leads to higher scores for exons exhibiting simi-

larity to known proteins, but decreased scores for predicted exons having little to no similarity

with known proteins.

A typical GENSCAN output is shown in Figure 5.5, using the human uroporphyrinogen

decarboxylase (URO-D) gene (U30787) as the query. Each exon in the prediction is shown in

a separate line. The columns, going from left to right, represent the gene and exon number

(Gn.Ex), the type of prediction (Type, either the exon type or an identified polyA signal),

the strand on which the prediction was made (+ or –), the beginning and endpoints for the

prediction, the length of the predicted exon, its reading frame, several scoring columns, and

a probability value (P). GENSCAN exons having a very high probability value (p > 0.99) are

97.7% accurate when the prediction matches a true, annotated exon. These high-probability

predictions can be used in the rational design of polymerase chain reaction primers for comple-

mentary DNA (cDNA) amplification, or for other purposes where extremely high confidence

is necessary. GENSCAN exons that have probabilities in the range of 0.50–0.99 are deemed to

be correct most of the time. The best-case accuracies for p values higher than 0.90 is on the

order of 88%. Any predictions having p < 0.50 should be deemed unreliable, and those data

How Well Do Gene Predictors Work?

Figure 5.5 Sample output from a GENSCAN analysis of the uroporphyrinogen decarboxylase gene. See the text for a more detailed descrip-

tion of the output.

are not given in the data table. The predicted amino acid sequence is given below the gene

predictions. In the example shown here, GENSCAN correctly predicted nine of the 10 exons

in URO-D. Only the initial exon was missed.

<!-- END LEGACY FRAGMENT: 03_Ab_Initio_Gene_Prediction_in_Eukaryotic_Genomes -->

---

How Well Do Gene Predictors Work?

The accuracy of gene prediction programs is usually determined using controlled, well-defined

datasets, where the actual gene structure has been determined experimentally. Accuracy can

be computed at either the nucleotide, exon, or gene level, and each provides different insights

into the accuracy of a predictive method. In the field of prokaryotic gene prediction, the results

are almost always reported at the gene level and given in terms of a percentage – that is, the

Genome Annotation

TP

FP

FN

FN

TN

TN

TP

Actual

Predicted

Sensitivity

Specificity

Sn = TP/(TP + FN)

Sp = TN/(TN + FP)

Figure 5.6 Schematic representation of measures of gene prediction accuracy at the nucleotide level.

The actual gene structure is illustrated at the top with confirmed exons identified with light blue bars and

confirmed introns in black lines. The predicted gene structure is illustrated at the bottom with predicted

exons identified with red bars and predicted introns in black lines. The four possible outcomes of a

prediction are shown: true positives (TP), true negatives (TN), false positives (FP), and false negatives

(FN). The equations for sensitivity and specificity are also shown using appropriate combinations of TP,

TN, FP, and FN.

number of correct gene predictions divided by the total number of known or validated genes in

the test set. In some cases, the number or percentage of over-predicted genes (false positives)

is also reported. In the field of eukaryotic gene prediction, performance reporting tends to be

somewhat more convoluted. This is because the evaluation problem is more complex and the

overall performance is often much worse. As a general rule, two basic measures are used: sen-

sitivity (or Sn), defined as the proportion of coding nucleotides, exons, or genes that have been

predicted correctly; and specificity (or Sp), defined as the proportion of coding and non-coding

nucleotides, exons, or genes that have been predicted correctly (i.e. the overall fraction of the

prediction that is correct). A more detailed explanation of sensitivity, specificity, and a number

of other evaluation metrics used in gene (and protein structure) prediction is given in Box 5.4.

Also introduced in this box are the concepts of true positives (TPs), true negatives (TNs), false

positives (FPs), and false negatives (FNs).

An example of a eukaryotic gene prediction with the four possible outcomes is shown in

Figure 5.6. This figure schematically illustrates the differences between a gene prediction and

the known (or observed) gene structure. Neither sensitivity nor specificity alone provides a

perfect measure of global accuracy, as high sensitivity can be achieved with little specificity

and vice versa. An easier to understand measure that combines the sensitivity and specificity

values is called the Matthews correlation coefficient (MCC or just CC), which is described

more formally in Box 5.4. The MCC ranges from −1 to 1, where a value of 1 corresponds to

a perfect prediction; a value of −1 indicates that every coding region has been predicted as

non-coding, and vice versa. Other accuracy measures are sometimes used as well; however,

the above-mentioned ones have been most commonly employed in the large assessment

projects on eukaryotic genome prediction such as the human ENCODE Genome Annotation

Assessment Project (EGASP; Guigó and Reese 2005), the RNA-seq Genome Annotation

Assessment Project (RGASP; Steijger et al. 2013) and the Nematode Genome Annotation

Assessment Project (nGASP; Coghlan et al. 2008).

Box 5.4 Evaluating Binary Classifications or Predictions in Bioinformatics

Many predictions in bioinformatics involve essentially binary or binomial (i.e. true/false)

classification problems. For instance, prokaryotic gene prediction can be framed as a

binary classification problem where one tries to distinguish open reading frames (ORFs)

from non-ORFs. Similarly, eukaryotic gene prediction can be posited as a binary classifica-

tion problem of predicting exons and non-exons (introns) or genes and intergenic regions.

Protein membrane helix prediction (discussed in Chapter 7) can be put in a similar binary

classification frame, where one distinguishes between membrane helices and non-helices

(or non-membrane regions). Binary classification problems can also be found in medicine,

where one tries to predict or diagnose sick patients versus healthy patients, or quality

control tasks in high-throughput manufacturing (pass versus fail).

How Well Do Gene Predictors Work?

The evaluation of binary classifiers or predictors normally follows a very standard prac-

tice with a common set of metrics and definitions. Unfortunately, this practice is not always

followed when bioinformaticians evaluate their own predictors or predictions. This is why

we have included this very important information box, one that is referred to frequently

throughout this book.

As shown in the diagram below, a binary classifier or predictor can have four combina-

tions of outcomes: true positives (TP or correct positive assignments), true negatives (TN or

correct negative assignments), false positives (FP or incorrect positive assignments), and

false negatives (FN or incorrect negative assignments). In statistics, the false positives are

called type I errors and the false negatives are called type II errors (see Chapter 18).

Observed state

positive

Predicted state

positive

True positive

(TP)

False positive

(FP)

Predicted state

negative

False negative

(FN)

True negative

(TN)

Observed state

negative

PREDICTION

OBSERVATION

Once a binary classifier has been run on a set of data, it is possible to calculate specific

numbers for each of these four outcomes using the above 2 × 2 contingency table. So

a gene predictor that predicted 1000 genes in a genome that had only 900 genes may

have 850 TPs, 200 TNs, 60 FPs, and 40 FNs. From this set of 4 outcomes it is possible to

calculate 8 ratios. These ratios can be obtained by dividing each of the four numbers

(TP, TN, FP, FN) by the sum of its row or column in the 2 × 2 contingency table. The most

important ratios and their names or abbreviations are listed (along with their formulae)

below. Also included are several other binary classifier evaluation metrics that are used

by certain subdisciplines in bioinformatics or statistics.

Name

Formula

Sensitivity (Sn)

Recall

True positive rate (TPR)

TP

TP + FN

Specificity (Sp)

True negative rate (TNR)

TN

TN + FP

Precision

Positive predictive value (PPV)

TP

TP + FP

False positive rate (FPR)

FP

FP + TN

False discovery rate (FDR)

FP

FP + TP

Negative predictive value (NPV)

TN

TN + FN

Accuracy (ACC), Q2

TP + TN

TP + FP + TN + FN

F1 score

F score

F measure

2TP

(2TP + FP + FN)

Matthews correlation

coefficient (MCC)

TP × TN −FP × FN

√

(TP + FP)(TP + FN)(TN + FP)(TN + FN)

(Continued)

Genome Annotation

Box 5.4 (Continued)

Sensitivity (Sn, recall, or TPR) measures the proportion of actual positives that are cor-

rectly identified as such, while specificity (Sp or TNR) measures the proportion of actual

negatives that are correctly identified as such. Precision (PPV) is the proportion of positive

results that are true positive results, while NPV is the proportion of negative results that

are true negative results. FDR is the binary (not the multiple testing) measure of false

positives divided by all positive predictions. Accuracy or ACC (for binary classification) is

defined as the number of correct predictions made divided by the total number of pre-

dictions made. ACC is one of the best ways of assessing binary test or predictor accuracy.

The F1 score is another measure of test accuracy and is defined as the harmonic average

of precision (PPV) and recall (Sn). MCC is a popular measure of test or predictor accuracy.

It is essentially a chi-squared statistic for a standard 2 × 2 contingency table. In effect, MCC

is the correlation coefficient between the observed and predicted binary classifications.

Different fields of science have different preferences for different metrics owing to dif-

ferent traditions or different objectives. In medicine and most fields of biology (including

bioinformatics), sensitivity and specificity are most often used to assess a binary classi-

fier, while in machine learning and information retrieval, precision and recall are usually

preferred. Likewise, different prediction tasks within bioinformatics tend to report perfor-

mance with different measures. Gene predictors generally report Sn, Sp, and ACC, while

protein structure predictors generally report ACC and MCC. The accuracy (ACC) score in

protein secondary structure prediction has also been termed Qn, where n is the number

of secondary structure classes (usually n = 3). In gene prediction, the ACC score is given

as Q2 since only two classes are identified (either exon/intron or ORF/non-ORF). Each of

the above ratios (except for MCC) can take on values from 0 to 1. For a perfect prediction,

Sn = 1, Sp = 1, PPV = 1, NPV = 1, ACC = 1, F1 = 1, MCC = 1, FPR = 0, or FDR = 0, whereas

a completely incorrect prediction would yield Sn = 0, Sp = 0, PPV = 0, NPV = 0, ACC = 0,

F1 = 0, MCC = −1, FPR = 1, or FDR = 1.

The performance of any binary predictor has to be assessed based on the existing bias

in the numbers within different classes (i.e. uneven class distribution). For example, an

ACC of 0.95 may seem excellent, but if 95% of the dataset belongs to just one class, then

the same ACC score could be easily achieved by simply predicting everything to be in that

class. This is the situation with many mammalian genomes, which have large intergenic

regions. Therefore predicting every nucleotide as being “intergenic” would easily create

a nucleotide-based gene predictor that would be >95% accurate. Such a predictor would,

of course, be completely useless.

Assessing the accuracy of gene prediction methods requires sets of reliably annotated genes

verified by experimental or computational evidence derived from complementary sources of

information. Experimental evidence can be provided by mass spectrometry-based proteomics

or by structural biology methods such as nuclear magnetic resonance spectroscopy or X-ray

crystallography (Chapter 12) that provide direct, visual confirmation of the protein sequences.

Computational evidence can appear in the form of similarity of the derived protein sequence

to the primary structures of proteins whose functions were verified experimentally. Extensive

gene prediction assessments have been done for both prokaryotic and eukaryotic organisms.

Assessing Prokaryotic Gene Predictors

The evaluation of prokaryotic gene predictors has been ongoing for many years, with each

publication that describes a new program (or a new version of an existing program) providing

a detailed performance assessment (Larsen and Krogh 2003; Delcher et al. 2007; Hyatt et al.

Assessing Eukaryotic Gene Predictors

2010; Borodovsky and Lomsadze 2011). One of the most recent and comprehensive assess-

ments of prokaryotic gene predictors was conducted by Hyatt et al. in 2010. In this paper the

authors compared five different programs: Prodigal 1.20 (Hyatt et al. 2010), GeneMarkHMM

2.6 (Borodovsky and Lomsadze 2011), GLIMMER 3.02 (Delcher et al. 2007), EasyGene 1.2

(Larsen and Krogh 2003), and MED 2.0 (Zhu et al. 2007) on two different tasks. The first task

involved predicting experimentally verified genes with experimentally verified TISs from 10

different bacterial and archaeal genomes. In this case, only 2443 (out of a possible 35 000+)

genes were considered experimentally verified. Hyatt et al. found that all five of the programs

were able to achieve a 98–99.8% level of accuracy (at the gene level) for the 3′ end of these ver-

ified genes and an 87–96.7% level of accuracy (at the gene level) for the complete genes (both

the 5′ and 3′ ends being correctly predicted). The second task involved predicting GenBank

(mostly manually) annotated genes from seven different bacterial genomes. In this case, a total

of 23 648 genes were evaluated. All of the programs were able to achieve a 95–99% level of accu-

racy (at the gene level) for the 3′ end of these genes. However, their performance for the full

gene prediction task (both 5′ and 3′ being correctly predicted) was much more variable, with

accuracy values ranging from 69% to 91% correct. The overall prediction average for all five

programs over all genes in this second task was about 80%. It is also notable that all five pro-

grams generally over-predicted the number of genes annotated in GenBank by about 4–5%,

with some programs (MED 2.0) over-predicting by as much as 40%.

Based on the data provided by Hyatt et al. (2010), the two best-performing prokaryotic gene

prediction programs were Prodigal and GeneMark, while the other three programs were only

marginally worse. Their results also show that the task of predicting the 3′ ends of prokaryotic

genes is essentially solved, while the challenge of predicting the 5′ ends of prokaryotic genes

needs more work. It is also evident that some prokaryotic genomes are harder to predict than

others, with the full-gene prediction performance on the E. coli genome often hovering about

90% while the performance for less studied genomes (such as Halobacterium salinarum) often

is around 70%. These results reflect the fact that ab initio gene predictors (both prokaryotic and

eukaryotic) require very extensive training on a large number of high-quality gene models.

Once trained, these tools can perform very impressively, especially in well-studied genomes

for which ample training data are available. However, the level of training required to reach

very high accuracy is often hard to achieve for newly assembled bacterial genomes.

Assessing Eukaryotic Gene Predictors

The assessment of eukaryotic gene predictors has been going on for more than 20 years. In

the early days, most eukaryotic gene prediction evaluations were conducted on single genes

whose exon/intron structure had been well characterized. This reflected the fact that very few

(if any) eukaryotic genomes had been fully sequenced and only a small number of eukaryotic

genes had their exon/intron structure fully determined. It also made the gene prediction tasks

much simpler, as the coding (exon) density is much higher (25–50%) than what would be found

over an entire genome (which is often <2%). This also led to overly optimistic performance

ratings. More recently, the field has evolved to assessing gene prediction performance over

entire genomes.

Burset and Guigó (1996) published one of the first systematic evaluations of eukaryotic gene

predictors. Their study evaluated seven programs, using a set of 570 vertebrate single-gene

sequences. The average CC at the nucleotide level for these programs ranged from 0.65 to 0.80.

Later, Rogic et al. (2001) performed a similar analysis of seven gene prediction programs, using

a set of 195 single-gene sequences from human and rodent species. The programs tested in the

Rogic et al. study showed substantially higher accuracy than those reported on in the Burset

and Guigó study, with the average CC at the nucleotide level ranging from 0.66 to 0.91. This

increase in the upper part of the range illustrates the significant advances that occurred in the

development of gene prediction methods over a relatively short period of time.

Genome Annotation

The early evaluations put forth by Burset and Guigó (1996), Rogic et al. (2001), and others all

suffered from the same limitation: the gene finders were tested using controlled datasets com-

prising short genomic sequences encoding a single gene with simple gene structures. These

datasets are obviously not representative of genomic sequences as a whole. Complete genome

sequences contain long stretches with low coding density, stretches coding for multiple or

incomplete genes (or both), and stretches having very complex or alternative gene structures.

As a result, two large-scale studies were conducted to assess the performance of ab initio

eukaryotic gene predictors on real-world mammalian genomic data. The first was based on an

analysis of human chromosome 22 (Parra et al. 2003) and the second was based on an analysis

of the human ENCODE regions (Guigó et al. 2006), covering about 1% of the human genome.

When human chromosome 22 was sequenced, it was subjected to very extensive manual

analyses, experimental confirmation, and detailed annotation by many experts (Dunham et al.

1999). This was done to provide a useful gold standard (at the time) for assessing genome

prediction and genome annotation tools. As a result, Parra et al. used the manually annotated

data for chromosome 22 to assess the performance of GENSCAN (Burge and Karlin 1997),

GenomeScan (Yeh et al. 2001), TBLASTX (Gish and States 1993), GeneID (Blanco et al. 2002),

and SGP-2 (Parra et al. 2003) at the nucleotide, exon, and whole gene/transcript level. The

results were quite disappointing. At the nucleotide level the programs had an average sensitiv-

ity/specificity ([Sp + Sn]/2) value ranging from 0.62 to 0.75 and a CC value ranging from 0.54 to

0.73. At the exon level, the programs had an average sensitivity/specificity value ranging from

0.54 to 0.62 and at the gene/transcript level the average sensitivity/specificity value ranged

from 0.05 to 0.11. The latter values are the numbers of greatest interest as they reflect the true

level of gene prediction performance. Interestingly, GENSCAN and GenomeScan performed

somewhat worse than GeneID and SGP-2. Indeed, SGP-2 consistently performed better than

all of the “pure” ab initio predictors as it also made use of comparative genomic data from

mouse chromosome 22. The inclusion of experimental sequence data technically made SGP-2

an extrinsic gene finder rather than a pure ab initio gene predictor.

A similar level of very high-quality manual annotation was achieved in 2005–2006 during

the first phase of the Encyclopedia of DNA Elements (ENCODE) project. The ENCODE project

is a long-term, multi-phase project that started in 2003 with the goal of identifying all of the

functional elements within the human genome sequence. During its pilot phase, a number of

regions from the human genome (approximately 1%) were selected for detailed investigation.

The availability of this “gold standard” dataset led to a second, much larger evaluation that

looked at the predictive performance of pure ab initio predictors, as well as gene finders that

used additional extrinsic data such as sequence homology and experimental sequencing data

(Guigó et al. 2006). For the Guigó et al. study, four ab initio predictors were tested: AUGUSTUS

(Hoff and Stanke 2013), GeneMark-A (Besemer and Borodovsky 2005), GeneMark-B (Besemer

and Borodovsky 2005), and GeneZilla (Allen et al. 2006). Once again, the results were quite dis-

appointing. At the nucleotide level, the programs had a CC value ranging from 0.53 to 0.76. At

the exon level the programs had an average sensitivity/specificity value ranging from 0.40 to

0.57, and at the gene or transcript level the average sensitivity/specificity value ranged from

0.05 to 0.14. Overall, AUGUSTUS performed significantly better than the other ab initio pro-

grams but not at a level that would allow one to use it to automatically annotate a eukaryotic

genome. However, the most important findings from this study were that significant improve-

ments (up to two times better at the exon level and up to four times better at the gene level)

could be made to the quality of eukaryotic gene annotations if comparative genomic data or

other experimental/extrinsic evidence were employed in the prediction process.

It was because of these studies that a significant change in the gene prediction community

occurred. In particular, the developers of gene predictors moved from reluctantly using experi-

mental or extrinsic data to wholeheartedly embracing experimental data. In other words, gene

prediction began to change to gene finding and genome prediction began to evolve toward

genome annotation. In doing so, genome analysis became a more holistic, evidence-based

process that combined ab initio gene prediction with extrinsic gene-finding methods. These

extrinsic gene-finding methods combined many other computational tools and other lines of

evidence including gene expression data, proteomic data, sequence homology to other anno-

tated genomes, and even literature-derived data.

<!-- END LEGACY FRAGMENT: 04_How_Well_Do_Gene_Predictors_Work -->

---

Evidence Generation for Genome Annotation

Genomic evidence is any information that can be used to identify or inform the structure

of a gene in an organism – be it a prokaryotic or a eukaryotic organism. Some of the most

useful evidence comes from experimental work such as transcriptional data (mRNA or

DNA data derived from RNA-seq experiments) or protein sequence data gathered about the

organism of interest or a closely related organism. Other kinds of evidence can be collected

through running various bioinformatic programs that identify genomic features such as

sequence repeats, tRNA and rRNA genes, pseudogenes, transcription factor binding sites,

retroviruses, prophages, and so on. In the following sections, we will briefly review some

of the evidence-generating approaches used for both extrinsic gene finding and genome

annotation.

<!-- END LEGACY FRAGMENT: 05_Evidence_Generation_for_Genome_Annotation -->

<!-- BEGIN LEGACY FRAGMENT: 06_Gene_Annotation_and_Evidence_Generation_Using_RNA_seq_Data -->

Source: Ch5 Genome Annotation / Gene Annotation and Evidence Generation Using RNA-seq Data

PDF Pages: 153-154 | Print Pages: 133-134

Boundary: starts at RNA-seq subsection heading on PDF page 153; stops before the next subsection heading on PDF page 154.

---

Gene Annotation and Evidence Generation Using RNA-seq Data

RNA sequencing (RNA-seq) is a next-generation DNA sequencing (NGS) technique that

involves converting RNA (mRNA, tRNA, and rRNA) transcripts into double-stranded cDNA

fragments, then sequencing them using low-cost NGS sequencing methods (Wang et al. 2009).

Over the last decade, RNA-seq has helped to revolutionize genome annotation methods for

both eukaryotes and prokaryotes (Trapnell et al. 2009; Sallet et al. 2014). A typical RNA-seq

experiment generates thousands of short DNA sequence reads corresponding to gene-coding

regions (also known as coding sequence or CDS segments). These sequences can then be

aligned to the reference genome sequence using gapped, short-read aligners to determine

which genome regions were being transcribed. Some of the more popular gapped short-read

aligners include TopHat2 (Kim et al. 2013), Stampy (Lunter and Goodson 2011), and GSNAP

(Wu et al. 2016). These alignments can be further processed into putative transcripts using

tools such as Cufflinks (Trapnell et al. 2012), StringTie (Pertea et al. 2015), or Trinity (Grabherr

et al. 2011). In this way, RNA-seq provides experimental evidence (through DNA sequencing)

regarding the location of gene-coding regions.

The improvement in gene-finding performance and gene annotation quality when RNA-seq

data are used is quite substantial. In RGASP (Steijger et al. 2013), a diverse set of 14 genome

annotation approaches were compared (including intrinsic/ab initio methods, extrinsic meth-

ods, and hybrid extrinsic/intrinsic methods). The gold standard for comparison was the refer-

ence human genome annotations from the GENCODE project, consisting of computationally,

manually, and experimentally determined gene annotations (Harrow et al. 2012). It turned

out that the best-performing programs for the task of identifying protein-coding genes were

gene annotation tools that used RNA-seq data. Examples of gene annotation programs that

incorporate RNA-seq data into their gene finders include AUGUSTUS (Hoff and Stanke 2013),

mGENE (Schweikert et al. 2009), Trembly, and Transomics (Sperisen et al. 2004).

As noted earlier, when processing RNA-seq data for genome annotation, one can either

splice-align the raw reads against the genome or, alternatively, transcript fragments can first

be assembled de novo and then aligned to the genome via BLASTN. This “mapping-first”

approach was shown to lead to more accurate annotations in the RGASP assessment and so it

is highly recommended. Spliced alignment can be done with tools such as GSNAP (Wu et al.

2016), Stampy (Lunter and Goodson 2011), TopHat2 (Kim et al. 2013), or STAR (Dobin et al.

2013). Integrating coverage information from RNA-seq data into a gene annotation tool can

typically be done by increasing the score of candidate exons that are covered by RNA-seq by

a certain factor that depends on the local coverage of each covered exonic region. Rewarding

Genome Annotation

individual splice sites, supported by RNA-seq evidence, is relatively easy for HMMs. Some

gene annotation tools integrate evidence for complete introns (i.e. splice site pairs) from

RNA-seq data.

Newer technologies that produce longer RNA-seq reads (10 000+ bases) have greatly

improved the ability to predict alternatively spliced transcripts in comparison with short reads

(100–400 bases), which mostly help to find local alternative splice variants. Long reads are

often near-complete transcripts and each spliced alignment gives the structure of a transcript,

albeit only approximately owing to the relatively high sequencing error rate. Gene finders

such as AUGUSTUS can integrate evidence from long-read alignments to further improve

their performance.

While RNA-seq has greatly improved the performance of many eukaryotic gene find-

ers, there is still a long way to go. According to the RGASP assessment (Steijger et al.

2013), the best-performing methods identified ∼59% of protein-coding transcripts from

the Caenorhabditis elegans genome (AUGUSTUS, mGene, and Transomics), 43% from the

Drosophila melanogaster genome (AUGUSTUS), and just 21% from the Homo sapiens genome

(Trembly). So, RNA-seq data have not (yet) been the key to “solving” the problem of accurate

and automatic eukaryotic genome annotation. Important issues still remain, including the

fact that a significant fraction of genes or splice forms may not be expressed in any RNA-seq

sample, that transcribed sequences may not be protein coding and, if they are, the correct

protein-coding ORFs remain to be identified, and that transcript assemblies and mapping

of the transcripts to the genome are notoriously error prone. These errors typically are seen

around exon boundaries, with the assemblies often extending into introns and, at times,

missing whole exons. Several programs have been developed to help address these mapping

problems, including Exonerate (Slater and Birney 2005) and GeneWise (Birney et al. 2004).

Both programs are “splice-aware” tools that can be used to polish BLAST alignments. These

polished alignments can then be used to improve the annotation of the exons, introns, splice

sites, and 5′ and 3′ UTRs.

<!-- END LEGACY FRAGMENT: 06_Gene_Annotation_and_Evidence_Generation_Using_RNA_seq_Data -->

<!-- BEGIN LEGACY FRAGMENT: 07_Gene_Annotation_and_Evidence_Generation_Using_Protein_Sequence_Databases -->

Source: Ch5 Genome Annotation / Gene Annotation and Evidence Generation Using Protein Sequence Databases

PDF Pages: 154-155 | Print Pages: 134-135

Boundary: starts at protein sequence databases subsection heading on PDF page 154; includes the page 155 tail despite the running header; stops before the true comparative gene prediction heading.

---

Gene Annotation and Evidence Generation Using Protein Sequence Databases

Just as RNA-seq data can be used as evidence for the existence of genes, so too can sequence

homology be used to locate or identify new genes in newly sequenced organisms. In

homology-based gene finding, the DNA sequence of the newly sequenced organism is trans-

lated into putative protein sequences and these putative sequences are then compared against

databases of known proteins. Homologous matches at the protein level can then be used to

annotate, identify, and locate the genes at the DNA level. A key advantage of homology-based

gene finding over ab initio gene prediction is that homology-based methods provide not only

the identification and location (as ab initio approaches do) but also the probable gene name

and probable gene function as inferred by the sequence similarity of the newly identified gene

to previously annotated proteins in the protein sequence databases.

Translated nucleotide searches such as BLASTX searches (Gish and States 1993) constitute

one of the simplest homology-based gene prediction approaches. These searches are partic-

ularly useful when comparing ORFs in prokaryotic genomes. However, when dealing with

the split nature of eukaryotic genes, BLASTX-like searches do not resolve exon splice bound-

aries particularly well. One useful approach is to use both the results of translated nucleotide

searches along with those produced through the use of ab initio methods. Examples of this

hybrid approach include programs such as GenomeScan (Yeh et al. 2001), GeneID (Blanco

et al. 2002), and AUGUSTUS (Hoff and Stanke 2013). GenomeScan is an extension of GEN-

SCAN that incorporates sequence similarity to known proteins using BLASTX.

A more sophisticated approach to eukaryotic gene prediction via sequence homology

involves aligning the genomic query against a protein target that is presumed to be homol-

ogous to the protein encoded in the genomic sequence that is being annotated. In these

alignments, often referred to as spliced alignments, large gaps corresponding to introns in the

query sequence are only allowed at “legal” splice junctions. Examples of programs using this

approach include PROCRUSTES (Gelfand et al. 1996), GeneWise (Birney and Durbin 1997),

Exonerate (Slater and Birney 2005), BLAT (Kent 2002), and GenomeThreader (Gremme et al.

2005).

The spliced alignment approach does not exploit all the information typically available for

homology-based gene prediction. In fact, for any given protein, a whole family of related pro-

teins is often available. Such an ensemble of sequences carries more information than just a

single protein. For instance, a well-constructed MSA shows which regions are well conserved

and which ones are prone to insertions or deletions. Using an MSA, it is possible to calculate

the probability that a certain amino acid occurs at a certain site. Using these data, it is possible

to calculate PWMs or PSSMs and to create what is called an MSA profile. While the task for

creating MSAs for prokaryotic genomes is relatively easy, the task of creating MSAs for eukary-

otic genomes is particularly challenging owing to the presence of repeats, as well as large-scale

genome rearrangements, duplications, and deletions.

So, rather than trying to find genes or exons with a set of individual protein sequences,

one may use MSAs of aligned protein families to do the job. These MSAs can be found in

orthology databases such as OrthoDB (Waterhouse et al. 2013). Several excellent software tools

have been developed to search the gene structures of members of a protein family given an

MSA profile representation of that family. These include GeneWise (Birney and Durbin 1997)

and AUGUSTUS-PPX (Keller et al. 2011), where PPX stands for Protein Profile eXtension.

AUGUSTUS-PPX has been shown to improve gene prediction accuracy over spliced align-

ment methods, especially when dealing with genes having large numbers of exons. However,

the MSA approach is limited to the availability of homologous families and by the degree of

sequence similarity. Therefore, MSA gene finding is best used in situations involving medium

to high sequence similarity.

More recently, this MSA concept has been extended to cover situations with more remote

sequence similarity through the development of BUSCO (Simão et al. 2015). BUSCO stands

for Benchmarking Universal Single-Copy Orthologs. These single-copy orthologs correspond

to a relatively small set of proteins that are highly conserved and that are found as single-copy

genes across many different phyla in the tree of life. The BUSCO dataset currently includes

3023 genes for vertebrates, 2675 for arthropods, 843 for metazoans, 1438 for fungi, 429 for

eukaryotes, and 40 universal marker genes for prokaryotes. Using HMMER (Eddy 2009), the

BUSCO gene set can be rapidly searched against any given query genome. The presence or

absence of these BUSCO genes in a given organism provides a good measure of the complete-

ness of a genome assembly. It also provides a good measure of the completeness of a given

genome annotation or a given genome prediction.

<!-- END LEGACY FRAGMENT: 07_Gene_Annotation_and_Evidence_Generation_Using_Protein_Sequence_Databases -->

<!-- BEGIN LEGACY FRAGMENT: 08_Gene_Annotation_and_Evidence_Generation_using_Comparative_Gene_Prediction -->

Source: Ch5 Genome Annotation / Gene Annotation and Evidence Generation using Comparative Gene Prediction

PDF Pages: 155-156 | Print Pages: 135-136

Boundary: starts at true comparative gene prediction subsection heading on PDF page 155; stops before the next non-protein-coding/foreign genes section on PDF page 156.

---

Gene Annotation and Evidence Generation using Comparative

Gene Prediction

Another approach to homology-based gene prediction exploits the fact that there is a large and

growing number of completely sequenced and well-annotated genomes now available. This

has given rise to a technique called comparative gene prediction. The rationale behind com-

parative gene prediction is that functional regions (the protein-coding regions) tend to be more

conserved than non-protein-coding regions. This observation provides the basis for identifying

protein-coding regions in newly sequenced genomes. Comparative gene prediction methods

exploit sequence homology but at a far more global scale than the protein sequence similar-

ity methods described above. In comparative gene prediction, the “known” and “unknown”

genomes are from different species, but the species are assumed to be so closely related that

their entire genomes can be aligned. Because these genomes are so long (millions to billions

of bases), the pairwise alignment or MSA is typically broken down into many local alignments

of syntenic (homologous) regions.

Genome Annotation

Early methods for comparative gene finding typically used just two genomic sequences as

input, such as DOUBLESCAN (Meyer and Durbin 2002), TWINSCAN (Korf et al. 2001), SLAM

(Alexandersson et al. 2003), or SGP-2 (Parra et al. 2003). SLAM is an HMM-based method

in which gene predictions and sequence alignments are performed simultaneously. TWIN-

SCAN and DOUBLESCAN are extensions of GENSCAN, whereas SGP-2 is an extension of

GeneID. Later on, comparative gene-finding methods were developed that could use more

than two genomic sequences to predict genes in a new genome, but they only did so for a sin-

gle target genome. These methods include programs such as N-SCAN (Gross and Brent 2006),

CONTRAST (Gross et al. 2007), and Mugsy-Annotator (Angiuoli et al. 2011). More recently, a

method called clade annotation has been developed and implemented in a version of AUGUS-

TUS known as “comparative AUGUSTUS” (König et al. 2016). Clade annotation allows the

simultaneous alignment and annotation of multiple target genomes. For example, compara-

tive AUGUSTUS can be used to simultaneously annotate the genomes of multiple (up to 20)

different mouse strains.

<!-- END LEGACY FRAGMENT: 08_Gene_Annotation_and_Evidence_Generation_using_Comparative_Gene_Prediction -->

<!-- BEGIN LEGACY FRAGMENT: 09_Evidence_Generation_for_Non_Protein_Coding_Non_Coding_or_Foreign_Genes -->

Source: Ch5 Genome Annotation / Evidence Generation for Non-Protein-Coding, Non-Coding, or Foreign Genes

PDF Pages: 156 | Print Page: 136

Boundary: starts at non-protein-coding/non-coding/foreign genes section heading; stops before tRNA and rRNA Gene Finding subsection.

---

Evidence Generation for Non-Protein-Coding, Non-Coding, or Foreign Genes

One of the best ways of determining the location of a protein-coding gene is to determine where

it is not. In other words, knowing that a DNA segment cannot possibly code for a protein allows

one to exclude it from the gene/protein-finding process. Prokaryotic genomes contain many

genes that do not code for proteins. These include tRNA and rRNA genes, along with many

foreign prophage genes that may or may not code for real phage proteins. Likewise, eukaryotic

genomes are filled with repeat regions, pseudogenes, retrotransposons, and retroviral genes,

along with an assortment of tRNA and rRNA genes. These non-protein-coding or non-coding

elements can account for 20–30% of a given prokaryotic genome (Casjens 2003) and more than

90% of a eukaryotic genome (Li et al. 2004).

<!-- END LEGACY FRAGMENT: 09_Evidence_Generation_for_Non_Protein_Coding_Non_Coding_or_Foreign_Genes -->

<!-- BEGIN LEGACY FRAGMENT: 10_tRNA_and_rRNA_Gene_Finding -->

Ch5 Genome Annotation — tRNA and rRNA Gene Finding

PDF page 156 下部 – page 157；印刷页码 136-137

Both prokaryotes and eukaryotes have a significant portion of their genome occupied by tRNA and rRNA genes. Prokaryotes (including bacteria and archaea) typically contain 70–80 copies each of tRNA genes and between three and 45 copies of rRNA genes. tRNA molecules are L-shaped adaptor RNA molecules, typically 76–90 nucleotides in length, that are essential for the translational process (Figure 5.7). In principle, a total of 61 tRNA genes are needed to permit the translation of all 61 coding (sense) codons. However, because of a phenomenon known as base wobble, many organisms are able to have a single tRNA serving two or more codons. As a result, most prokaryotes have 35–40 unique tRNA genes, but one or two copies of each.

rRNA molecules are the principal constituents (>60% by mass) of the ribosome, the translational engine of all cells. In prokaryotes the ribosome consists of two subunits, the small and the large subunit, which pair up to form the ribosome. The rRNAs present in prokaryotes are the 5S and 23S rRNAs in the large subunit and the 16S rRNA in the small subunit. Genes encoding these rRNAs are typically arranged into an operon (the rrn operon), with an internally transcribed spacer between the 16S and 23S rRNA genes. The number of rrn operons in prokaryotes ranges from one to 15 per genome.

tRNA and rRNA genes in eukaryotes share many similarities (in both structure and size) with those in prokaryotes. There are, however, some minor differences. For instance, eukaryotes typically have many more copies of tRNA genes than prokaryotes. There are 275 tRNA genes in Saccharomyces cerevisiae, 620 copies of tRNA genes in C. elegans, and 497 copies of tRNA genes in humans. All eukaryotes have 22 mitochondrial tRNA genes. Like prokaryotes, eukaryotic rRNA genes are also divided according to their location in the large or small subunit in the ribosome. However, instead of just two large subunit rRNAs, there are three rRNAs in the eukaryotic large subunit: 5S, 5.8S, and 28S. Just as with prokaryotes, eukaryotes have one rRNA gene (18S rRNA) for their small ribosomal subunit, but they also encode rRNA genes for their mitochondrial ribosomes (12S and 16S rRNA genes). Unlike prokaryotes, eukaryotes generally have many copies of the rRNA genes organized into tandem repeats. In humans, approximately 300–400 rRNA repeats are present in five clusters, located on five separate chromosomes. Unlike prokaryotic tRNA genes, some tRNA genes in eukaryotes are interrupted with introns.

The structure of tRNAs is highly conserved across all major kingdoms of life and there are a large number of tRNA sequences that are known for both prokaryotes and eukaryotes. As a result, most methods for identifying tRNA genes take advantage of common sequence motifs (recognizable via HMMs) and employ some kind of sequence homology or database comparison to identify tRNA genes. The best performing and most popular methods are RNAmmer (Lagesen et al. 2007), tRNAfinder (Kinouchi and Kogawa 2006), and tRNAscan-SE (Lowe and Eddy 1997). These programs are able to identify tRNA genes in both prokaryotes and eukaryotes with very high accuracy (>95%). In addition to these programs, there are several dedicated databases of tRNA sequences to assist with comparative tRNA identification approaches; these include tRNAdb and tRNA-DB-CE (Jühling et al. 2009; Abe et al. 2014). Currently the identification of tRNA genes is considered to be a "solved" problem.

Just like tRNA genes, rRNA genes also exhibit a very high level of sequence conservation, and there are a number of rRNA motifs that can be described by HMMs. These HMMs have been integrated into the program (and web server) called RNAmmer (Lagesen et al. 2007). RNAmmer is able to identify all rRNAs from prokaryotes and eukaryotes, with the exception of 5.8S rRNA. Owing to the complexity, size, and relatively poor annotation of rRNA genes, the performance for rRNA prediction is not yet at the same level as tRNA prediction. In addition to the RNAmmer predictor, there is also an RNA database called Rfam (Kalvari et al. 2018) that contains >2600 RNA families (including rRNA and tRNA sequence families). Each sequence family in Rfam is represented by an MSA, a consensus secondary structure, and a covariance model. Rfam can be used for rRNA (and other RNA gene) identification via sequence comparisons or MSAs. Regardless of their current shortcomings, the use of tRNA and rRNA gene identification tools invariably improves the accuracy of any protein-coding gene-finding or gene-prediction effort. It also enhances the quality of the overall genome annotation.

<!-- END LEGACY FRAGMENT: 10_tRNA_and_rRNA_Gene_Finding -->

<!-- BEGIN LEGACY FRAGMENT: 11_Prophage_Finding_in_Prokaryotes -->

Ch5 Genome Annotation — Prophage Finding in Prokaryotes

PDF page 157 下部 – page 158 中部；印刷页码 137-138

Prokaryotes are subject to constant attack by bacterial viruses called bacteriophages, which kill or disable susceptible bacteria. Bacteriophages are the most abundant biological entities on the planet and they play a major role in the bacterial ecosystem and in driving microbial genetic variation or genetic diversity. This genetic diversity is brought on through a particularly unique part of the bacteriophage lifecycle called lysogeny. Lysogeny involves the integration of the phage genome (often consisting of 10–20 genes) into the host bacterial chromosome at well-defined insertion points. The genetically integrated phages are called prophages. In some cases, prophages can become permanently embedded into the bacterial genome, becoming cryptic prophages (Little 2005). These cryptic prophages often serve as genetic “fodder” for future evolutionary changes of the host microbe (Bobay et al. 2014). Furthermore, prophages and cryptic prophages tend to introduce pathogenic elements or pathogenic islands that exhibit a very different base composition than the host genome. Prophages and cryptic prophages can account for up to 20% of the genetic material in some bacterial genomes (Casjens 2003), with some prophage genes coding for expressed proteins and others not.

Given their high abundance, the identification of these phage-specific genetic elements can be quite important, especially when it comes to annotating bacterial genomes. Prophage and cryptic prophage sequences exhibit certain sequence features (such as the presence of integrases and transposases, attachment sites, and an altered base composition) that can be used to distinguish them from “normal” bacterial genes. When combined with HMMs to improve the sequence feature recognition tasks, it is possible to identify prophage and cryptic prophage sequences with relatively good accuracy. The accuracy can be improved further if comparative genome analyses to databases of known phage sequences are performed. Several bacterial prophage-finding programs have been developed and deployed over the past decade, including Phage_Finder (Fouts 2006) and ProphageFinder (Bose and Barber 2006). More recently, phage finding has moved from stand-alone programs to web servers. In particular, two new web servers have been released that provide somewhat greater speed and improved accuracy for prophage finding over existing tools. These are known as PHAST (Zhou et al. 2011) and PHASTER (Arndt et al. 2016). Both web servers are between 85% and 95% accurate (depending on the test being conducted) and both provide rich graphical output, as well as detailed annotations of the prophage sequences and the surrounding bacterial genomic sequences (Figure 5.8).

Regardless of the method chosen, the annotation of prophage and cryptic prophage genes certainly enhances the quality of a prokaryotic genome annotation and it usually improves the accuracy of any prokaryotic gene predictions.

<!-- END LEGACY FRAGMENT: 11_Prophage_Finding_in_Prokaryotes -->

<!-- BEGIN LEGACY FRAGMENT: 12_Repetitive_Sequence_Finding_Masking_in_Eukaryotes -->

Ch5 Genome Annotation — Repetitive Sequence Finding/Masking in Eukaryotes

PDF page 158 中部 – page 161 上部；印刷页码 138-141

Unlike prokaryotes, eukaryotic genomes contain considerable quantities of repetitive DNA. These repetitive sequences include retrotransposons and DNA transposons, both of which are referred to as dispersed repeats, as well as highly repetitive sequences, typically called tandem repeats. The most abundant repeat sequences in eukaryotes are retrotransposons. Retrotransposons are genetic elements that can amplify themselves using a “copy-and-paste” mechanism similar to that used by retroviruses. To replicate and amplify, they are transcribed into RNA, then converted back into identical DNA sequences using reverse transcription and then inserted into the genome at specific target sites. In contrast to retrotransposons, DNA transposons do their copying and pasting without an RNA intermediate, instead using the protein called transposase. Approximately 52% of the human genome is made up of retrotransposons, while DNA transposons account for another 3% (Lander et al. 2001; Wheeler et al. 2013). In plants, retrotransposons are much more abundant, accounting for between 60% and 90% of the DNA in any given plant genome (Li et al. 2004).

Within the retrotransposon family are two subfamilies: long terminal repeat retrotransposons (LTR retrotransposons) and non-LTR retrotransposons. LTR retrotransposons are retrovirus-like sequences that contain LTRs that range from ∼100 bp to over 5 kb in length. In fact, a retrovirus can be transformed into an LTR retrotransposon simply through inactivation or deletion of certain genes (such as the envelope protein) that enable cell-to-cell viral transmission. Most LTR retrotransposons are non-functional endogenous retroviruses that are also called proviruses. In this regard, eukaryotic LTR retrotransposons can be thought of as the equivalent of prokaryotic prophages or cryptic prophages. Human endogenous retroviral sequences, all of which appear to be defective or non-replicative, account for about 8% of the human genome (Taruscio and Mantovani 2004).

Non-LTR retrotransposons consist of two subtypes: long interspersed nuclear elements (LINEs) and short interspersed nuclear elements (SINEs). LINEs are typically 7000 bp long and encode several genes to cover all the functions needed for retrotransposition. These include reverse transcriptase and endonuclease genes, as well as several genes needed to form a ribonucleoprotein particle. More than 850000 copies of LINEs exist in the human genome, covering 21% of all human DNA (Cordaux and Batzer 2009). However, more than 99% of LINEs are genetically “dead,” having lost their retrotransposition functions. In contrast to LINEs, SINEs are much smaller, typically consisting of DNA stretches spanning just 80–500 bp. SINEs are very abundant (in the millions of copies), accounting for about 10% of the DNA in the human genome. The most common SINEs in humans are the Alu repeats (Häsler and Strub 2006). Alu repeats are about 300 bp long. They are highly conserved in primates and are subject to frequent DNA methylation events.

In addition to transposable elements (or dispersed repeats), eukaryotes also contain large numbers of tandem repeats, including minisatellite DNA, microsatellite DNA (also known as short tandem repeats [STRs] or simple sequence repeats [SSRs]), and telomere repeats. Minisatellite DNA consists of repeats of 10–60 bp in length that stretch for about 2 kb and are scattered throughout the genome. Microsatellite DNA consists of repeats of 1–6 bp, extending for hundreds of kilobases, particularly around the centromeres. Telomere repeats consist of a highly conserved 6 bp sequence (TTAGGG) that is repeated 250–1000 times and found exclusively at the ends of eukaryotic chromosomes. Mini- and microsatellite DNA account for about 5% of the DNA in the human genome (Subramanian et al. 2003).

The fact that eukaryotes have so many repeat sequences, combined with the fact that these repeats account for such a large portion of their genomes (often >50%), has led to concerted efforts by genome annotators to identify, remove, or mask these sequences. This is because repeat sequences can seriously hinder gene identification and genome annotation activities. For instance, retrotransposons and DNA transposons can beeasily be mistaken as exons by ab initio gene predictors. Likewise, STRs can lead to spurious alignments when comparative genomics approaches are used in gene finding. STRs (which are also called low-complexity regions) can often be dealt with through two techniques: soft or hard masking. Soft masking is done by changing the case of the letters in the sequence file from uppercase to lowercase, while hard masking changes the offending sequence to Ns, thereby removing them completely from consideration. Soft masking prevents the masked region from seeding alignments but preserves sequence identity so that off-target alignments are minimized. Soft masking is routinely done by the programs SEG and DUST (Wootton and Federhen 1993), which are found in most versions of the BLAST sequence alignment suite.

While tandem repeats are relatively easy to deal with, repeat transposable elements (such as retrotransposons) are much harder to handle. This is because these sequences are far larger and far more complex. The Repbase (Jurka et al. 2005) database contains a comprehensive collection of repeat and transposable elements from a wide range of species. This resource is frequently used to identify repeat elements via comparative sequence analysis. However, if the transposon sequences are highly divergent from those found in RepBase, then other methods or other databases may need to be used. Dfam (Wheeler et al. 2013) is an example of a more advanced repetitive element database. In Dfam, the original Repbase sequences have been converted to HMMs. The use of these HMMs has allowed many more transposable elements to be identified (up to 54.5% vs. 44% in humans) with much improved accuracy (Wheeler et al. 2013).

In addition to Dfam (which is available as both a server and a downloadable resource), there are several stand-alone programs and web servers that have been developed to specifically identify retrotransposons, including RECON (Bao and Eddy 2002), RepeatScout (Price et al. 2005), RetroPred (Naik et al. 2008), LTR_FINDER (Xu and Wang 2007), LTRharvest (Ellinghaus et al. 2008), and MITE-Hunter (Han and Wessler 2010). These programs identify and label either LTR or non-LTR retrotransposons. While this information may be useful to some, many genome annotators are simply interested in removing retrotransposons from consideration. In this regard, RepeatMasker (Tarailo-Graovac and Chen 2009) has become the tool of choice as it simply hard masks (i.e. removes) all detectable retrotransposon sequences from the genome of interest.

As a general rule, hard masking of transposable elements is often the first step performed in a eukaryotic genome annotation. Hard masking with tools such as Dfam or programs such as RepeatMasker not only removes “uninteresting” genetic data, it also accelerates the gene identification process and improves annotation accuracy. Because coding exons tend not to overlap or to contain repetitive elements, ab initio gene prediction programs tend to predict fewer false-positive exons when using hard-masked sequences. For instance, when chromosome 22 was analyzed using different ab initio gene predictors, it was found that a significant reduction in false-positive gene predictions occurred (Parra et al. 2003). In particular, GENSCAN initially predicted 1128 protein-coding genes without using sequence masking, but when sequence masking was used the number of predicted genes dropped to 789. When GeneID was used, the number fell from 1119 to 730. The actual number of protein-coding genes in chromosome 22, according to the latest GENCODE annotation, is 489.

<!-- END LEGACY FRAGMENT: 12_Repetitive_Sequence_Finding_Masking_in_Eukaryotes -->

<!-- BEGIN LEGACY FRAGMENT: 13_Finding_and_Removing_Pseudogenes_in_Eukaryotes -->

Ch5 Genome Annotation — Finding and Removing Pseudogenes in Eukaryotes

PDF page 161 上部；印刷页码 141

A particular challenge with eukaryotic genome annotation is differentiating between predictions identifying “real” genes from those that correspond to non-functional pseudogenes. Database searches may not help to provide any clearer picture, as many pseudogenes are similar to functional, paralogous genes. The absence of an RNA transcript from an RNA-seq experiment cannot be used as a criterion either, because RNA transcripts do not always exist for actual genes because of variations in tissue expression or developmental stages. In general, intronless gene predictions for which multi-exon paralogous genes exist in the same genome are suspicious, as they may indicate sequences that have arisen through retrotransposition.

Multi-exon predictions, however, can also correspond to pseudogenes arising through a recent gene duplication event. If homologs in another organism exist, one solution is to compute the synonymous versus non-synonymous substitution rate (Ka/Ks; Fay and Wu 2003). Ka/Ks values approaching 1 are indicative of neutral evolution, suggesting a pseudogene. Support for multi-exon gene predictions can come from assessing the conservation of the overall gene structure in close homologs. For instance, the prediction or identification of homologous genes in two modestly related organisms (e.g. mouse and human) most likely indicates that the gene is real and is not a pseudogene (Guigó et al. 2003).

<!-- END LEGACY FRAGMENT: 13_Finding_and_Removing_Pseudogenes_in_Eukaryotes -->

In the early days of genome annotation, when it would often take years just to sequence a single organism, teams of researchers and bioinformaticians would gather and work together for many months, or even years, to assemble the genome, perform the initial ab initio gene predictions, manually collate the experimental or literature-derived evidence, conduct comparative sequence analysis, and then synthesize the data into a consensus genome annotation. This was routinely done for both bacterial and eukaryotic genomes (Lander et al. 2001; Winsor et al. 2005; Riley et al. 2006). Indeed, it is still being done for the GENCODE project, which has been preparing and updating the reference human genome annotation since 2003 (Harrow et al. 2012). However, these efforts have required (and continue to require) enormous resources and time. With the appearance of very high-throughput NGSs and the ability to routinely sequence an entire genome in a few days, these manual approaches to genome annotation have become unsustainable. Now, most genome annotations are done through automated pipelines that help users to synthesize multiple pieces of evidence and data to generate a consensus genome annotation.

The choice of a pipeline tool depends on the type of organism (eukaryote vs. prokaryote), the computational resources available, the available evidence (RNA-seq or no RNA-seq data), and the similarity of the organism to previously annotated organisms. For instance, if one is annotating a genome with a closely related, previously annotated species, a simple comparative analysis or sequence projection should be sufficient. If the organism of interest has no closely related annotated species, a pipeline that uses RNA-seq or experimentally acquired protein sequence data will generate more accurate annotations. The most advanced genome annotation pipelines require many programs and perform complex analyses that need supercomputers such as large multi-core machines or massive computing clusters (maintained locally or available via the Cloud). For example, to annotate the loblolly pine genome (which contains 22 billion bases – seven times more than the human genome) required 8640 central processing units (CPUs) running for 14.6 hours (Wegrzyn et al. 2014). In the following sections we will briefly describe some commonly used annotation pipelines for prokaryotes and eukaryotes.

<!-- END LEGACY FRAGMENT: 14_Genome_Annotation_Pipelines -->

<!-- BEGIN LEGACY FRAGMENT: 15_Prokaryotic_Genome_Annotation_Pipelines -->

Ch5 Genome Annotation — Prokaryotic Genome Annotation Pipelines

PDF page 162 中部 – page 163 图注；印刷页码 142-143

Annotation pipelines for prokaryotes typically do not require the same computational resources as for eukaryotes. Indeed, most bacterial genomes can be annotated in less than 30 minutes, whether on a web server or on a desktop computer. However, the recent shift toward metagenomics or community bacterial genomics is beginning to lead to significantly greater computational demands that will be discussed in more detail in Chapter 16. Some of the more popular publicly available prokaryotic genome annotation pipelines include Prokka (Seemann 2014), Rapid Annotation using Subsystem Technology (RAST; Overbeek et al. 2014), and the Bacterial Annotation System (BASys; Van Domselaar et al. 2005). Prokka is an open-source Perl program that runs with a command-line interface (on UNIX). Prokka can be used to annotate pre-assembled bacterial, archaeal, and viral sequences. With Prokka, a typical 4 million base pair bacterial genome can be fully annotated in less than 10 minutes on a quad-core computer. Prokka is also capable of producing standards-compliant output files for further analysis or viewing. Prokka’s appeal lies in its speed and ability to perform “private” annotations on local computers.

In contrast to Prokka, RAST and BASys are genome annotation web servers. Web servers are generally easier to use but they do not offer the privacy of a locally installed program. RAST is a registration-based web server that accepts standard, pre-assembled DNA sequence files and then identifies protein-encoding, rRNA, and tRNA genes, assigns functions to the genes, and finally uses this information to reconstruct a metabolic network for the organism. In contrast to RAST, BASys is an open access web server. BASys accepts pre-assembled FASTA-formatted DNA or protein files from bacteria, archaea, and viruses and performs many of the same annotation functions as RAST. However, BASys provides a much greater depth of annotation (covering more than 50 calculable properties) and produces colorful, easily viewed genome maps (Figure 5.9) using a program called CGView (Stothard and Wishart 2005).

Figure 5.9 A screenshot of a BASys bacterial genome annotation output for the bacterium Salmonella enterica. The BASys image can be interactively zoomed-in to reveal rich annotations for all of the genes in the genome.

<!-- END LEGACY FRAGMENT: 15_Prokaryotic_Genome_Annotation_Pipelines -->

<!-- BEGIN LEGACY FRAGMENT: 16_Eukaryotic_Genome_Annotation_Pipelines -->

Ch5 Genome Annotation — Eukaryotic Genome Annotation Pipelines

PDF page 163 中部 – page 164；印刷页码 143-144

Given the complexity of eukaryotic genomes, their corresponding annotation pipelines must do somewhat more than those used for prokaryotic genomes. In particular, eukaryotic genome annotation pipelines must combine not only the ab initio gene predictions (or multiple gene predictions from multiple sources) but also many other pieces of evidence, including experimental data. As a result, almost all modern eukaryotic genome annotation pipelines use a technique called "evidence clustering" to identify gene regions and then use the aligned RNA (from RNA-seq) and protein evidence to improve the accuracy of the gene predictors. Some pipelines go even further and make use of a "combiner" algorithm to select the combination of exons that are best supported by the evidence. Two combiner programs in particular are very good at this: JIGSAW (Allen and Salzberg 2005) and EVidence Modeler, or EVM (Haas et al. 2008). These programs assess different types of evidence based on known error profiles and various kinds of user input and then choose the best combination of exons to minimize the error. In particular, EVM combines aligned protein and RNA transcript evidence with ab initio predictions into weighted consensus gene models, while JIGSAW uses non-linear models or weighted linear combiners to choose a single best consensus gene model.

Among the most widely used eukaryotic genome annotation pipelines, all of which use some kind of combiner algorithm, are MAKER2 (Holt and Yandell 2010), Ensembl (Fernández-Suárez and Schuster 2010), the National Center for Biotechnology Information (NCBI) Eukaryotic Annotation Pipeline (Thibaud-Nissen et al. 2016), PASA (Haas et al. 2008), and BRAKER1 (Hoff et al. 2016). The MAKER2 annotation pipeline is a highly parallelizable, stand-alone program that aligns and polishes protein sequence and transcriptome (RNA-seq) data with BLAST; it also provides evidence-based hints to various gene predictors and it creates an evidence trail with various quality metrics for each annotation. Some of MAKER2's quality metrics include the number of splice sites confirmed by RNA-seq evidence, the number of exons confirmed by RNA-seq data, and the lengths of 5' and 3' UTRs. MAKER2 also uses a quality metric called the Annotation Edit Distance, or AED (Eilbeck et al. 2009). The AED value ranges between 0 and 1, with higher quality annotations being associated with lower AEDs. MAKER2 uses these AED values to choose the best gene predictions from which to build its final annotation. Like the MAKER2 pipeline, the Ensembl genome annotation pipeline builds its gene models from aligned and polished protein sequence- and RNA-seq-derived transcriptome data. To complete the annotation process, Ensembl merges identical transcripts, and a non-redundant set of transcripts is reported for each gene. Both MAKER2 and Ensembl supply hints to indicate intron/exon boundaries using protein and RNA-seq alignments to their internal gene predictors. This helps to generate gene models that better represent the aligned evidence. This approach also helps improve gene prediction accuracy for poorly (or insufficiently) trained gene finders. Like the Ensembl and MAKER2 pipelines, the NCBI Annotation Pipeline aligns and polishes protein and transcriptome data. It also generates gene predictions using the Gnomon gene-finding program (Souvorov et al. 2010). The NCBI system typically assigns higher weights to manually curated evidence over computationally derived models or computationally generated evidence. The PASA genome annotation pipeline is one of the oldest annotation pipelines and was one of the first to use a combiner or evidence-clustering algorithm (EVM). PASA aligns RNA transcripts to the reference genome using BLAT (Kent 2002) or GMAP (Wu et al. 2016). PASA is capable of generating annotations based on RNA transcriptome data, on pre-existing gene models, or on ab initio gene predictions. PASA, along with the MAKER2 and Ensembl annotation pipelines, is able to add UTRs to their genome annotations via RNA-seq data to further increase their accuracy. One of the latest additions to publicly available eukaryotic genome annotation pipelines is the BRAKER suite of programs (Hoff et al. 2016). BRAKER1 (and most recently BRAKER2) combines the strengths of GeneMark-ET with AUGUSTUS – both of which use RNA-seq data to improve their gene annotation accuracy. In the BRAKER pipeline, GeneMark-ET is used first to train and generate initial gene structures, then AUGUSTUS makes use of the initially predicted genes for further training and integrates RNA-seq data into the final gene predictions. BRAKER1 has been shown to be 10–20% more accurate than MAKER2 in terms of gene and exon sensitivity/specificity.

Even with an exon accuracy of >90% (rarely achieved by even the best eukaryotic genome annotation pipelines), most genes in a genome will have at least one incorrectly annotated exon. Incorrectly identified genes or mistaken gene annotations can have very serious consequences for experimentalists who are designing experiments to study gene functions. Indeed, many failed molecular biology or gene-cloning experiments can be traced back to incorrect gene annotations. Furthermore, incorrect annotations can propagate, leading to a cascade of errors that affect many other scientists. This happens when an incorrect annotation is innocently passed on to another genome project and then used as evidence in still more genome annotation efforts which eventually end up in public databases. To help prevent these errors or reduce the magnitude of these mistakes, most annotation pipelines include some kind of quality metrics which are attached to each and every gene annotation. Most of these metrics are based on a score that measures the agreement of a given gene annotation to an aligned RNA/protein sequence or on the basis of the homology and synteny of the gene to closely related species. Some pipelines use a simple star rating (ranging from zero to five). Zero stars correspond to an annotation where none of the exons is supported by aligned evidence, while a five-star rating corresponds to a situation where every exon is supported and every splice site is confirmed by a single full-length cDNA. Other pipelines use more sophisticated metrics, such as the AED score (mentioned above). Protein family domains can also be good indicators of annotation quality and annotation completeness. Certainly, any annotation that contains an identifiable protein domain is more likely to encode a functional protein than one that does not. Domain matching has been used to rescue a number of gene annotations that would have otherwise received a "failing" quality score owing to a poor sequence alignment. Both Ensembl and MAKER2 report the fraction of annotations containing a protein family domain as a quality measure. Interestingly, this fraction (0.69) appears to be quite constant across genomes; the closer a given genome is to this fraction, the more confidence one has in its quality. In addition to the domain-matching fraction, the presence or absence of BUSCO genes can also be used to provide a measure of the completeness of a genome annotation (Simão et al. 2015). Another excellent route to ensure good quality annotations is through manual inspection with genome visualization and editing software. This is discussed in more detail below.

<!-- END LEGACY FRAGMENT: 16_Eukaryotic_Genome_Annotation_Pipelines -->

<!-- BEGIN LEGACY FRAGMENT: 17_Visualization_and_Quality_Control -->

Ch5 Genome Annotation — Visualization and Quality Control

PDF page 165 上部；印刷页码 145

While automated or semi-automated pipelines for genome annotation have become the norm, there is still a need for a human factor in annotating genomes and assessing their quality. Having a knowledgeable biologist or some kind of "domain expert" carefully look through a genome annotation is essential to ensure that the annotations makes sense. This manual review process also allows one to catch and correct suspicious annotations or fill in missing annotations. However, to perform these manual reviews or curatorial tasks, it is necessary to visualize and interactively edit the annotations. Certainly two of the best known genome browsers are the University of California Santa Cruz Genome Browser (Casper et al. 2018) and Ensembl's Genome Browser (Fernández-Suárez and Schuster 2010), both of which have been thoroughly reviewed in Chapter 4. While these tools are excellent for visualizing genome annotations, there are also a number of other tools that support both visualizing and editing genome annotations, including Web Apollo (Lee et al. 2013), GenomeView (Abeel et al. 2012), and Artemis (Carver et al. 2012).

Web Apollo is both a visualization tool and a genome editor. More specifically, it is a web-based plug-in for JBrowse (Westesson et al. 2013) that provides an editable, user-created annotation track. All edits in Web Apollo are visible in real time to all members of the annotation team. This feature is particularly helpful when undertaking a community annotation project or when many investigators are involved in a particular genome analysis. GenomeView is an open-source, stand-alone genome viewer and editor that allows users to dynamically browse large volumes of aligned short-read data. It supports dynamic navigation and semantic zooming, from the whole genome level to the single nucleotide level. GenomeView is particularly noted for its ability to visualize whole genome alignments of dozens of genomes relative to a reference sequence. It also supports the visualization of synteny and multi-alignment data. Artemis is a genome browser and annotation tool that allows users to easily visualize, browse, and interpret large NGS datasets. It supports multiple sequence read views and variant displays, along with a comprehensive set of read alignment views and read alignment filters. It also has the ability to simultaneously display multiple different views of the same dataset to its users. Artemis can read EMBL and GENBANK database entries, FASTA sequence formats (indexed or raw), and other features in EMBL and GENBANK formats.

When reviewing an annotated genome (regardless of whether it is from a prokaryote or a eukaryote), it is always useful to randomly select a specific region and to use the chosen visualization/editing tools to carefully analyze the annotations together with the evidence provided. This evidence may include the ab initio predicted genes, the spliced RNA-seq alignments, or any homologous protein alignments. While browsing through the selected region, one may notice certain genes or clusters of genes that seem to contradict the displayed evidence. For instance, the RNA-seq data may support additional or different splice forms. Alternatively, certain cross-species proteins may map to genomic regions where no gene has been previously predicted. Visual inspection can also reveal certain systematic problems with the annotation process, such as a tendency to miss genes with known database homologs or the appearance of repeats that overlap or mask many protein-coding genes. These problems may be addressed by changing parameter settings on the genome annotation pipeline, performing the necessary edits manually, or by choosing another tool. Multiple, iterative rounds of manual reviewing and manual editing followed by automated pipeline annotation are often necessary to complete a full and thorough genome annotation.

<!-- END LEGACY FRAGMENT: 17_Visualization_and_Quality_Control -->

Genome annotation has evolved considerably over the past two decades. These changes have been driven, in part, by significant improvements in computational techniques (for gene prediction) and in part by a significant expansion in the number of known and annotated genomes from an ever-growing number of diverse species. The availability of improved gene prediction tools, along with significantly expanded databases of well-annotated genes, proteins, and genomes, has moved genome annotation away from pure gene prediction to a more integrated, holistic approach that combines multiple lines of evidence to locate, identify, and functionally annotate genes. When combined with experimental data such as RNA-seq data or protein sequence data (from structural proteomics or expression-based proteomics), it is possible to obtain remarkably accurate and impressively complete annotations. This comprehensive blending of evidence is the basis for many newly developed, semi-automated or automated genome annotation pipelines and to many of the newer genome browsers and editors. However, not all genome annotation efforts can yield the same quantity or quality of information. Certainly prokaryotic genome annotation is faster, easier, and much more accurate than eukaryotic genome annotation. Indeed, the challenge of prokaryotic genome annotation is essentially a “solved problem,” while the challenge of eukaryotic genome annotation has to be considered as a “work in progress.”

<!-- END LEGACY FRAGMENT: 18_Summary -->

<!-- BEGIN SUPPLEMENTAL BACK MATTER: PDF pages 166-173 -->

Supplemental source: Acknowledgments / Internet Resources / Further Reading / References

PDF pages: 166-173 | Print pages: 146-153

Acknowledgments

TheauthorthanksAndyBaxevanisandRodericGuigófortheirhelpfulcommentsandtheuse

ofmaterialfromprioreditionsofthisbook.

Internet Resources

AbInitioProkaryoticGenePredictors

EasyGene(server) www.cbs.dtu.dk/services/EasyGene

GeneMark.hmm(server) opal.biology.gatech.edu/GeneMark/gmhmmp.cgi

GeneMarkS(server) opal.biology.gatech.edu/GeneMark/genemarks.cgi

GLIMMER(program) www.cs.jhu.edu/~genomics/Glimmer

Prodigal(program) github.com/hyattpd/Prodigal

AbInitioEukaryoticGenePredictors

GeneID%28server%29%20genome.crg.es/geneid.html

GeneMark-ES(program) opal.biology.gatech.edu/GeneMark

GeneZilla(program) www.genezilla.org

GenomeScan%28server%29%20hollywood.mit.edu/genomescan.html

GENSCAN%28server%29%20hollywood.mit.edu/GENSCAN.html

HMMgene(server) www.cbs.dtu.dk/services/HMMgene

SNAP%28program%29%20korflab.ucdavis.edu/software.html

Hybrid/ExtrinsicEukaryoticGenomeFinders

AUGUSTUS(server) bioinf.uni-greifswald.de/augustus

AUGUSTUS-PPX(program) bioinf.uni-greifswald.de/augustus

CONTRAST(program) contra.stanford.edu/contrast

GeneID(server) genome.crg.es/software/geneid

GeneWise(server) www.ebi.ac.uk/Tools/psa/genewise

GenomeThreader(program) genomethreader.org

GSNAP(program) research-pub.gene.com/gmap

mGENE(program) www.mgene.org

Internet Resources 147

Hybrid/ExtrinsicEukaryoticGenomeFinders

Mugsy-Annotator(program) mugsy.sourceforge.net

SGP-2(program) genome.crg.es/software/sgp2

STAR(program) code.google.com/archive/p/rna-star

Transomics(program) linux1.softberry.com/berry.phtml?topic=transomics

tRNAandrRNAFinders

Rfam(server) rfam.xfam.org

RNAmmer(server) www.cbs.dtu.dk/services/RNAmmer

RNAMotif%28program%29%20casegroup.rutgers.edu/casegr-sh-2.5.html

tRNAdb(server) trnadb.bioinf.uni-leipzig.de/DataOutput/Welcome

tRNADB-CE(server) trna.ie.niigata-u.ac.jp/cgi-bin/trnadb/index.cgi

tRNAfinder(server) ei4web.yz.yamagata-u.ac.jp/~kinouchi/tRNAfinder

tRNAscan-SE(server) lowelab.ucsc.edu/tRNAscan-SE

Phage-FindingTools

Phage_Finder(program) phage-finder.sourceforge.net

PHAST(server) phast.wishartlab.com

PHASTER(server) phaster.ca

RepeatFinding/MaskingTools

Dfam(server) www.dfam.org

LTR_FINDER(server) tlife.fudan.edu.cn/tlife/ltr_finder

LTRharvest%28program%29%20genometools.org/index.html

MITE-Hunter%28program%29%20target.iplantcollaborative.org/mite_hunter.html

Repbase(server) www.girinst.org/repbase

RepeatMasker(program) www.repeatmasker.org

RepeatScout(program) bix.ucsd.edu/repeatscout

RetroPred%28program%29%20www.juit.ac.in/attachments/RetroPred/home.html

ProkaryoticGenomeAnnotationPipelines

BASys(server) www.basys.ca

Prokka(program) www.vicbioinformatics.com/software.prokka.shtml

RAST(server/program) rast.nmpdr.org

EukaryoticGenomeAnnotationPipelines

BRAKER1(program) bioinf.uni-greifswald.de/bioinf/braker

EVM(program) evidencemodeler.github.io

JIGSAW(program) www.cbcb.umd.edu/software/jigsaw

MAKER2%28program%29%20www.yandell-lab.org/software/maker.html

PASA(program) github.com/PASApipeline/PASApipeline/wiki

GenomeBrowsersand/orEditors

Artemis(program) www.sanger.ac.uk/science/tools/artemis

Ensembl%28program%29%20uswest.ensembl.org/downloads.html

GenomeView(program) genomeview.org

JBrowse(program) jbrowse.org

UCSCGenomeBrowser%20hgdownload.cse.ucsc.edu/downloads.html

WebApollo(program) genomearchitect.github.io

148 Genome Annotation

Further Reading

Hoff,K.J.andStanke,M.(2015).Currentmethodsforautomatedannotationofprotein-coding

genes.Curr.Opin.InsectSci. 7,8–14.Awell-writtenandup-to-datesummaryofsomeofthe

latestdevelopmentsingenomeannotationwithsomeverypracticaladviceaboutwhich

annotationtoolsshouldbeused.

Nielsen,P.andKrogh,A.(2005).Large-scaleprokaryoticgenepredictionandcomparisonto

genomeannotation. Bioinformatics.21,4322–4329.Averyreadableassessmentofprokaryotic

genepredictionandgenomeannotation.

Yandell,M.andEnce,D.(2012).Abeginner’sguidetoeukaryoticgenomeannotation. Nat.Rev.

Genet.13,329–342.Anice,easy-to-readintroductiontotheprocessesinvolvedineukaryotic

genomeannotationalongwithusefuldescriptionsoftheavailablecomputationaltoolsandbest

practices.

Yoon,B.(2009).HiddenMarkovmodelsandtheirapplicationsinbiologicalsequenceanalysis.

Curr.Genomics 10,402–415.AcomprehensivetutorialonHMMsthatprovidesmanyuseful

examplesandexplanationsofhowdifferentHMMsareconstructedandusedingeneprediction

andgenesequenceanalysis.

References

Abe,T.,Inokuchi,H.,Yamada,Y.etal.(2014).tRNADB-CE:tRNAgenedatabasewell-timedin

theeraofbigsequencedata. Front.Genet. 5:114.

Abeel,T.,VanParys,T.,Saeys,Y.etal.(2012).GenomeView:anext-generationgenomebrowser.

NucleicAcidsRes. 40(2):e12.

Alexandersson,M.,Cawley,S.,andPatcher,L.(2003).SLAM:cross-speciesgenefindingand

alignmentwithageneralizedpairhiddenMarkovmodel. GenomeRes. 13:496–502.

Allen,J.E.andSalzberg,S.L.(2005).JIGSAW:integrationofmultiplesourcesofevidenceforgene

prediction.Bioinformatics21:3596–3603.

Allen,J.E.,Majoros,W.H.,Pertea,M.,andSalzberg,S.L.(2006).JIGSAW,GeneZilla,and

GlimmerHMM:puzzlingoutthefeaturesofhumangenesintheENCODEregions. Genome

Biol.7(Suppl1,S9):1–13.

Angiuoli,S.V.,DunningHotopp,J.C.,Salzberg,S.L.,andTettelin,H.(2011).Improving

pan-genomeannotationusingwholegenomemultiplealignment. BMCBioinf 12:272.

Arndt,D.,Grant,J.R.,Marcu,A.etal.(2016).PHASTER:abetter,fasterversionofthePHAST

phagesearchtool. NucleicAcidsRes. 44(W1):W16–W21.

Bao,Z.andEddy,S.R.(2002).Automateddenovoidentificationofrepeatsequencefamiliesin

sequencedgenomes. GenomeRes. 12:1269–1276.

Bellman,R.E.(1957). DynamicProgramming.Princeton:PrincetonUniversityPress.

Besemer,J.andBorodovsky,M.(2005).GeneMark:websoftwareforgenefindinginprokaryotes,

eukaryotesandviruses. NucleicAcidsRes. 33(WebServer):W451–W454.

Besemer,J.,Lomsadze,A.,andBorodovsky,M.(2001).GeneMarkS:aself-trainingmethodfor

predictionofgenestartsinmicrobialgenomes.Implicationsforfindingsequencemotifsin

regulatoryregions. NucleicAcidsRes. 29:2607–2618.

Birney,E.andDurbin,R.(1997).Dynamite:aflexiblecodegeneratinglanguagefordynamic

programmingmethodsusedinsequencecomparison.In: ProceedingsoftheFifthInternational

ConferenceonIntelligentSystemsforMolecularBiology,Halkidiki,Greece(21–26June1997) ,vol.

5,56–64.MenloPark,CA:AAAIPress.

Birney,E.,Clamp,M.,andDurbin,R.(2004).GeneWiseandGenomewise. GenomeRes. 14:

988–995.

Blanco,E.,Parra,G.,andGuigó,R.(2002).Usinggeneidtoidentifygenes.In: CurrentProtocolsin

Bioinformatics,vol.1,unit4.3.NewYork:Wiley.

References 149

Blattner,F.R.,Plunkett,G.3rd,,Bloch,C.A.etal.(1997).Thecompletegenomesequenceof

Escherichiacoli K-12.Science277:1453–1462.

Bobay,L.-M.,Touchon,M.,andRocha,E.P.C.(2014).Pervasivedomesticationofdefective

prophagesbybacteria. Proc.NatlAcad.Sci.USA. 111:12127–12132.

Borodovsky,M.andLomsadze,A.(2011).Geneidentificationinprokaryoticgenomes,phages,

metagenomes,andESTsequenceswithGeneMarkSsuite. Curr.Protoc.Bioinformatics .Chapter

4,Unit4.5.1–17.

Borodovsky,M.andMcIninch,J.(1993).GeneMark:parallelgenerecognitionforbothDNA

strands.Comput.Chem. 17:123–133.

Borodovsky,M.,Rudd,K.E.,andKoonin,E.V.(1994).Intrinsicandextrinsicapproachesfor

detectinggenesinabacterialgenome. NucleicAcidsRes. 22:4756–4767.

Bose,M.andBarber,R.D.(2006).ProphageFinder:aprophagelocipredictiontoolforprokaryotic

genomesequences. InSilicoBiol.(Gedrukt) 6:223–227.

Burge,C.andKarlin,S.(1997).PredictionofcompletegenestructuresinhumangenomicDNA. J.

Mol.Biol. 268:78–94.

Burset,M.andGuigó,R.(1996).Evaluationofgenestructurepredictionprograms. Genomics.34:

353–357.

Carver,T.,Harris,S.R.,Berriman,M.etal.(2012).Artemis:anintegratedplatformforvisualization

andanalysisofhigh-throughputsequence-basedexperimentaldata. Bioinformatics28:464–469.

Casjens,S.(2003).Prophagesandbacterialgenomics:whathavewelearnedsofar? Mol.Microbiol.

49:277–300.

Casper,J.,Zweig,A.S.,Villarreal,C.etal.(2018).TheUCSCGenomeBrowserdatabase:2018

update.NucleicAcidsRes. 46(D1):D762–D769.

Coghlan,A.,Fiedler,T.J.,McKay,S.J.etal.,andnGASPConsortium.(2008).nGASP–the

nematodegenomeannotationassessmentproject. BMCBioinf 9:549.

Cordaux,R.andBatzer,M.A.(2009).Theimpactofretrotransposonsonhumangenome

evolution.Nat.Rev.Genet. 10:691–703.

Delcher,A.L.,Harmon,D.,Kasif,S.etal.(1999).Improvedmicrobialgeneidentificationwith

GLIMMER.NucleicAcidsRes. 27:4636–4641.

Delcher,A.L.,Bratke,K.A.,Powers,E.C.,andSalzberg,S.L.(2007).Identifyingbacterialgenesand

endosymbiontDNAwithGlimmer. Bioinformatics23:673–679.

Dobin,A.,Davis,C.A.,Schlesinger,F.etal.(2013).STAR:ultrafastuniversalRNA-seqaligner.

Bioinformatics29:15–21.

Dunham,I.,Shimizu,N.,Roe,B.A.etal.(1999).TheDNAsequenceofhumanchromosome22.

Nature402:489–495.

Eddy,S.R.(2009).Anewgenerationofhomologysearchtoolsbasedonprobabilisticinference.

GenomeInform. 23:205–211.

Eilbeck,K.,Moore,B.,Holt,C.,andYandell,M.(2009).Quantitativemeasuresforthe

managementandcomparisonofannotatedgenomes. BMCBioinf 10:67.

Ellinghaus,D.,Kurtz,S.,andWillhoeft,U.(2008).LTRharvest,anefficientandflexiblesoftware

fordenovodetectionofLTRretrotransposons. BMCBioinf 9:18.

Ezkurdia,I.,Juan,D.,Rodriguez,J.M.etal.(2014).Multipleevidencestrandssuggestthatthere

maybeasfewas19,000humanprotein-codinggenes. Hum.Mol.Genet. 23:5866–5878.

Fay,J.C.andWu,C.(2003).Sequencedivergence,functionalconstraint,andselectioninprotein

evolution.Annu.Rev.GenomicsHum.Genet. 4:213–235.

Fernández-Suárez,X.M.andSchuster,M.K.(2010).Usingtheensemblgenomeservertobrowse

genomicsequencedata. Curr.Protoc.Bioinformatics .Chapter1,Unit1.15.

Fickett,J.W.andTung,C.S.(1992).Anassessmentofproteincodingmeasures. NucleicAcidsRes.

20:6441–6450.

Fouts,D.E.(2006).Phage_Finder:automatedidentificationandclassificationofprophageregions

incompletebacterialgenomesequences. NucleicAcidsRes. 34:5839–5851.

Gelfand,M.S.(1995).PredictionoffunctioninDNAsequenceanalysis. J.Comput.Biol. 2:87–117.

150 Genome Annotation

Gelfand,M.S.andRoytberg,M.A.(1993).Predictionoftheexon-intronstructurebyadynamic

programmingapproach. Biosystems.30:173–182.

Gelfand,M.S.,Mironov,A.A.,andPevner,P.A.(1996).Generecognitionviasplicedsequence

alignment.Proc.Natl.Acad.Sci.USA. 93:9061–9066.

Gish,W.andStates,D.(1993).Identificationofproteincodingregionsbydatabasesimilarity

search.Nat.Genet. 3:266–272.

Grabherr,M.G.,Haas,B.J.,Yassour,M.etal.(2011).Full-lengthtranscriptomeassemblyfrom

RNA-seqdatawithoutareferencegenome. Nat.Biotechnol. 29:644–652.

Gremme,G.,Brendel,V.,Sparks,M.E.,andKurtz,S.(2005).Engineeringasoftwaretoolforgene

structurepredictioninhigherorganisms. Inf.SoftwareTechnol. 47:965–978.

Gross,S.S.andBrent,M.R.(2006).Usingmultiplealignmentstoimprovegeneprediction. J.

Comput.Biol. 13:379–393.

Gross,S.S.,Do,C.B.,Sirota,M.,andBatzoglou,S.(2007).CONTRAST:adiscriminative,

phylogeny-freeapproachtomultipleinformantdenovogeneprediction. GenomeBiol. 8:R269.

Guigó,R.(1999).DNAcomposition,codonusageandexonprediction.In: GeneticDatabases (ed.

M.Bishop)),53–80.Cambridge,MA:AcademicPress.

Guigó,R.andReese,M.G.(2005).EGASP:collaborationthroughcompetitiontofindhuman

genes.Nat.Methods 2:575–577.

Guigó,R.,Dermitzakis,E.T.,Agarwal,P.etal.(2003).Comparisonofmouseandhumangenomes

followedbyexperimentalverificationyieldsanestimated1,019additionalgenes. Proc.Natl.

Acad.Sci.USA. 100:1140–1145.

Guigó,R.,Flicek,P.,Abril,J.F.etal.(2006).EGASP:thehumanENCODEgenomeannotation

assessmentproject. GenomeBiol. 7(Suppl1):S2.1–S2.31.

Haas,B.J.,Salzberg,S.L.,Zhu,W.etal.(2008).Automatedeukaryoticgenestructureannotation

usingEVidenceModelerandtheprogramtoassemblesplicedalignments. GenomeBiol. 9:R7.

Han,Y.andWessler,S.R.(2010).MITE-Hunter:aprogramfordiscoveringminiature

inverted-repeattransposableelementsfromgenomicsequences. NucleicAcidsRes. 38:e199.

Harrow,J.,Frankish,A.,Gonzalez,J.M.etal.(2012).GENCODE:thereferencehumangenome

annotationforTheENCODEProject. GenomeRes. 22:1760–1774.

Häsler,J.andStrub,K.(2006).Aluelementsasregulatorsofgeneexpression. NucleicAcidsRes.

34:5491–5497.

Hoff,K.J.andStanke,M.(2013).WebAUGUSTUS–awebservicefortrainingAUGUSTUSand

predictinggenesineukaryotes. NucleicAcidsRes. 41(WebServerissue):W123–W128.

Hoff,K.J.,Lange,S.,Lomsadze,A.etal.(2016).BRAKER1:unsupervisedRNA-seq-basedgenome

annotationwithGeneMark-ETandAUGUSTUS. Bioinformatics32:767–769.

Holt,C.andYandell,M.(2010).MAKER2:anannotationpipelineandgenome-database

managementtoolforsecond-generationgenomeprojects. BMCBioinf 12:491.

Hou,Y.andLin,S.(2009).Distinctgenenumber–genomesizerelationshipsforeukaryotesand

non-eukaryotes:genecontentestimationfordinoflagellategenomes. PLoSOne 4(9):e6978.

Hyatt,D.,Chen,G.L.,Locascio,P.F.etal.(2010).Prodigal:prokaryoticgenerecognitionand

translationinitiationsiteidentification. BMCBioinf 11:119.

Jühling,F.,Mörl,M.,Hartmann,R.K.etal.(2009).tRNAdb2009:compilationoftRNAsequences

andtRNAgenes. NucleicAcidsRes. 37(Databaseissue):D159–D162.

Jurka,J.,Kapitonov,V.V.,Pavlicek,A.etal.(2005).Repbaseupdate,adatabaseofeukaryotic

repetitiveelements. Cytogenet.GenomeRes. 110(1–4):462–467.

Kalvari,I.,Argasinska,J.,Quinones-Olvera,N.etal.(2018).Rfam13.0:shiftingtoa

genome-centricresourcefornon-codingRNAfamilies. NucleicAcidsRes. 46(D1):

D335–D342.

Keller,O.,Kollmar,M.,Stanke,M.,andWaack,S.(2011).Anovelhybridgenepredictionmethod

employingproteinmultiplesequencealignments. Bioinformatics27:757–763.

Kent,W.J.(2002).BLAT–theBLAST-likealignmenttool. GenomeRes. 12:656–664.

Kim,D.,Pertea,G.,Trapnell,C.etal.(2013).TopHat2:accuratealignmentoftranscriptomesinthe

presenceofinsertions,deletionsandgenefusions. GenomeBiol. 14:R36.

References 151

Kinouchi,M.andKuoakawa,K.(2006).tRNAfinder:asoftwaresystemtofindalltRNAgenesin

theDNAsequencebasedonthecloverleafsecondarystructure. J.Comput.AidedChem. 7:

116–126.

König,S.,Romoth,L.W.,Gerischer,L.,andStanke,M.(2016).Simultaneousgenefindingin

multiplegenomes. Bioinformatics32:3388–3395.

Korf,I.,Flicek,P.,Duan,D.,andBrent,M.R.(2001).Integratinggenomichomologyintogene

structureprediction. Bioinformatics.17:S140–S148.

Kozak,M.(1987).Ananalysisof5 ′-noncodingsequencesfrom699vertebratemessengerRNAs.

NucleicAcidsRes. 15:8125–8148.

Krogh,A.(1997).TwomethodsforimprovingperformanceofaHMMandtheirapplicationfor

genefinding.In: ProceedingsoftheFifthInternationalConferenceonIntelligentSystemsfor

MolecularBiology,Halkidiki,Greece(21–26June1997) ,vol.5,179–186.MenloPark,CA:AAAI

Press.

Krogh,A.,Mian,I.S.,andHaussler,D.(1994).AhiddenMarkovmodelthatfindsgenesin E.coli

DNA.NucleicAcidsRes. 22:4768–4678.

Kulp,D.,Haussler,D.,Reese,M.G.,andEeckman,F.H.(1996).AgeneralizedhiddenMarkov

modelfortherecognitionofhumangenesinDNA.In: ProceedingsoftheFourthInternational

ConferenceonIntelligentSystemsforMolecularBiology ,vol.4,134–142,June12-15,1996,St.

Louis,MO.USA,AAAIPress,MenloPark,California.

Lagesen,K.,Hallin,P.,Rødland,E.A.etal.(2007).RNAmmer:consistentandrapidannotationof

ribosomalRNAgenes. NucleicAcidsRes. 35:3100–3108.

Lander,E.S.,Linton,L.M.,Birren,B.etal.(2001).Initialsequencingandanalysisofthehuman

genome.Nature409:860–921.

Larsen,T.S.andKrogh,A.(2003).EasyGene–aprokaryoticgenefinderthatranksORFsby

statisticalsignificance. BMCBioinf 4:21.

Lee,E.,Helt,G.A.,Reese,J.T.etal.(2013).WebApollo:aweb-basedgenomicannotationediting

platform.GenomeBiol. 14:R93.

Li,W.,Zhang,P.,Fellers,J.P.etal.(2004).Sequencecomposition,organization,andevolutionof

thecoreTriticeaegenome. PlantJ. 40:500–511.

Lifton,R.P.,Goldberg,M.L.,Karp,R.W.,andHogness,D.S.(1978).Theorganizationofthehistone

genesin Drosophilamelanogaster:functionalandevolutionaryimplications. ColdSpring

HarborSymp.Quant.Biol. 42:1047–1051.

Little,J.W.(2005).Lysogeny,prophageinduction,andlysogenicconversion.In: Phages:TheirRole

inBacterialPathogenesisandBiotechnology (eds.M.K.Waldor,D.I.FriedmanandS.L.Adhya),

37–54.Washington,DC:ASMPress.

Lowe,T.M.andEddy,S.R.(1997).tRNAscan-SE:aprogramforimproveddetectionoftransfer

RNAgenesingenomicsequence. NucleicAcidsRes. 25:955–964.

Lukashin,A.V.andBorodovsky,M.(1998).GeneMark.hmm:newsolutionsforgenefinding.

NucleicAcidsRes. 26:1107–1115.

Lunter,G.andGoodson,M.(2011).Stampy:astatisticalalgorithmforsensitiveandfastmapping

ofIlluminasequencereads. GenomeRes. 21:936–939.

Macke,T.J.,Ecker,D.J.,Gutell,R.R.etal.(2001).RNAMotif,anRNAsecondarystructure

definitionandsearchalgorithm. NucleicAcidsRes. 29:4724–4735.

Meyer,I.M.andDurbin,R.(2002).Comparativeabinitiopredictionofgenestructuresusingpair

HMMs.Bioinformatics18:1309–1318.

Naik,P.K.,Mittal,V.K.,andGupta,S.(2008).RetroPred:atoolforprediction,classificationand

extractionofnon-LTRretrotransposons(LINEs&SINEs)fromthegenomebyintegrating

PALS,PILER,MEMEandANN. Bioinformation2:263–270.

Overbeek,R.,Olson,R.,Pusch,G.D.etal.(2014).TheSEEDandtherapidannotationofmicrobial

genomesusingsubsystemstechnology(RAST). NucleicAcidsRes. 42(Databaseissue):

D206–D214.

Parra,G.,Agarwal,P.,Abril,J.F.etal.(2003).Comparativegenepredictioninhumanandmouse.

GenomeRes. 13:108–117.

152 Genome Annotation

Pennisi,E.(2003).Bioinformatics.Genecountersstruggletogettherightanswer. Science.301:

1040–1041.

Pertea,M.,Pertea,G.M.,Antonescu,C.M.etal.(2015).StringTieenablesimprovedreconstruction

ofatranscriptomefromRNA-seqreads. Nat.Biotechnol. 33:290–295.

Pribnow,D.(1975).NucleotidesequenceofanRNApolymerasebindingsiteatanearlyT7

promoter.Proc.Natl.Acad.Sci.USA. 72:784–788.

Price,A.L.,Jones,N.C.,andPevzner,P.A.(2005).Denovoidentificationofrepeatfamiliesinlarge

genomes.Bioinformatics21(Suppl1):i351–i358.

Riley,M.,Abe,T.,Arnaud,M.B.etal.(2006). Escherichiacoli K-12:acooperativelydeveloped

annotationsnapshot–2005. NucleicAcidsRes. 34:1–9.

Rogic,S.,Mackworth,A.K.,andOuellette,F.B.F.(2001).Evaluationofgene-findingprogramson

mammaliansequences. GenomeRes. 11:817–832.

Sakharkar,M.,Passetti,F.,deSouza,J.E.etal.(2002).ExInt:anexonintrondatabase. Nucleic

AcidsRes. 30:191–194.

Sallet,E.,Gouzy,J.,andSchiex,T.(2014).EuGene-PP:anext-generationautomatedannotation

pipelineforprokaryoticgenomes. Bioinformatics30:2659–2661.

Schweikert,G.,Behr,J.,Zien,A.etal.(2009).mGene.web:awebserviceforaccurate

computationalgenefinding. NucleicAcidsRes. 37(WebServerissue):W312–W316.

Seemann,T.(2014).Prokka:rapidprokaryoticgenomeannotation. Bioinformatics30:2068–2069.

Shine,J.andDalgarno,L.(1975).Determinantofcistronspecificityinbacterialribosomes. Nature

254:34–38.

Simão,F.A.,Waterhouse,R.M.,Ioannidis,P.etal.(2015).BUSCO:assessinggenomeassemblyand

annotationcompletenesswithsingle-copyorthologs. Bioinformatics.31:3210–3212.

Slater,G.S.andBirney,E.(2005).Automatedgenerationofheuristicsforbiologicalsequence

comparison.BMCBioinf 6:31.

Slupska,M.M.,King,A.G.,Fitz-Gibbon,S.etal.(2001).Leaderlesstranscriptsofthecrenarchaeal

hyperthermophilePyrobaculumaerophilum. J.Mol.Biol. 309:347–360.

Souvorov,A.,Kapustin,Y.,Kiryutin,B.etal.(2010).Gnomon–NCBIeukaryoticgeneprediction

tool.NatlCent.Biotechnol.Inf. 2010:1–24.

Sperisen,P.,Iseli,C.,Pagni,M.etal.(2004).Trome,trESTandtrGEN:databasesofpredicted

proteinsequences. NucleicAcidsRes. 32(Databaseissue):D509–D511.

Steijger,T.,Abril,J.F.,Engström,P.G.etal.,andRGASPConsortium(2013).Assessmentof

transcriptreconstructionmethodsforRNA-seq. Nat.Methods 10:1177–1184.

Stothard,P.andWishart,D.S.(2005).Circulargenomevisualizationandexplorationusing

CGView.Bioinformatics21:537–539.

Subramanian,S.,Mishra,R.K.,andSingh,L.(2003).Genome-wideanalysisofmicrosatellite

repeatsinhumans:theirabundanceanddensityinspecificgenomicregions. GenomeBiol. 4:

R13.

Tarailo-Graovac,M.andChen,N.(2009).UsingRepeatMaskertoidentifyrepetitiveelementsin

genomicsequences. Curr.ProtocBioinformatics .Chapter4,Unit4.10.

Taruscio,D.andMantovani,A.(2004).Factorsregulatingendogenousretroviralsequencesin

humanandmouse. Cytogenet.GenomeRes. 105:351–362.

Thibaud-Nissen,F.,DiCuccio,M.,Hlavina,W.etal.(2016).TheNCBIeukaryoticgenome

annotationpipeline. J.Anim.Sci. 94(Suppl4)):184.

Trapnell,C.,Pachter,L.,andSalzberg,S.L.(2009).TopHat:discoveringsplicejunctionswith

RNA-seq.Bioinformatics25:1105–1111.

Trapnell,C.,Roberts,A.,Goff,L.etal.(2012).Differentialgeneandtranscriptexpressionanalysis

ofRNA-seqexperimentswithTopHatandcufflinks. Nat.Protoc. 7:562–578.

VanDomselaar,G.H.,Stothard,P.,Shrivastava,S.etal.(2005).BASys:awebserverforautomated

bacterialgenomeannotation. NucleicAcidsRes. 33(WebServerissue):W455–W459.

Wang,Z.,Gerstein,M.,andSnyder,M.(2009).RNA-seq:arevolutionarytoolfortranscriptomics.

Nat.Rev.Genet. 10:57–63.

References 153

Waterhouse,R.M.,Tegenfeldt,F.,Li,J.etal.(2013).OrthoDB:ahierarchicalcatalogofanimal,

fungalandbacterialorthologs. NucleicAcidsRes. 41(Databaseissue):D358–D365.

Wegrzyn,J.L.,Liechty,J.D.,Stevens,K.A.etal.(2014).Uniquefeaturesoftheloblollypine( Pinus

taedaL.)megagenomerevealedthroughsequenceannotation. Genetics196:891–909.

Westesson,O.,Skinner,M.,andHolmes,I.(2013).Visualizingnext-generationsequencingdata

withJBrowse. BriefingsBioinf. 14:172–177.

Wheeler,T.J.,Clements,J.,Eddy,S.R.etal.(2013).Dfam:adatabaseofrepetitiveDNAbasedon

profilehiddenMarkovmodels. NucleicAcidsRes. 41(Databaseissue):D70–D82.

Will,C.L.andLührmann,R.(2011).Spliceosomestructureandfunction. ColdSpringHarbor

Perspect.Biol. 3(7),pii:a003707.

Winsor,G.L.,Lo,R.,HoSui,S.J.etal.(2005).Pseudomonasaeruginosagenomedatabaseand

PseudoCAP:facilitatingcommunity-based,continuallyupdated,genomeannotation. Nucleic

AcidsRes. 33(Databaseissue):D338–D343.

Wootton,J.C.andFederhen,S.(1993).Statisticsoflocalcomplexityinaminoacidsequencesand

sequencedatabases. Comput.Chem. 17:149–163.

Wu,T.D.,Reeder,J.,Lawrence,M.etal.(2016).GMAPandGSNAPforgenomicsequence

alignment:enhancementstospeed,accuracy,andfunctionality. MethodsMol.Biol. 1418:

283–334.

Xu,Z.andWang,H.(2007).LTR_FINDER:anefficienttoolforthepredictionoffull-lengthLTR

retrotransposons.NucleicAcidsRes. 35(WebServerissue):W265–W268.

Yeh,R.,Lim,L.P.,andBurge,C.(2001).Computationalinferenceofthehomologousgene

structuresinthehumangenome. GenomeRes. 11:803–816.

Zhang,M.Q.(2002).Computationalpredictionofeukaryoticproteincodinggenes. Nat.Rev.Genet.

3:698–709.

Zhou,Y.,Liang,Y.,Lynch,K.H.etal.(2011).PHAST:afastphagesearchtool. NucleicAcidsRes. 39

(WebServerissue):W347–W352.

Zhu,H.,Hu,G.,Yang,Y.etal.(2007).MED:anewnon-supervisedgenepredictionalgorithmfor

bacterialandarchaealgenomes. BMCBioinf 8:97.

<!-- END SUPPLEMENTAL BACK MATTER: PDF pages 166-173 -->