2 These last years the cost of sequencing genomes has been greatly
3 reduced, and thus more and more genomes are sequenced. Therefore
4 automatic annotation tools are required to deal with this continuously
5 increasing amount of genomical data. Moreover, a reliable and accurate
6 genome annotation process is needed in order to provide strong
7 indicators for the study of life\cite{Eisen2007}.
9 Various annotation tools (\emph{i.e.}, cost-effective sequencing
10 methods\cite{Bakke2009}) producing genomic annotations at many levels
11 of detail have been designed by different annotation centers. Among
12 the major annotation centers we can notice NCBI\cite{Sayers01012011},
13 Dogma \cite{RDogma}, cpBase \cite{de2002comparative},
14 CpGAVAS \cite{liu2012cpgavas}, and
15 CEGMA\cite{parra2007cegma}. Usually, previous studies used one out of
16 three methods for finding genes in annoted genomes using data from
17 these centers: \textit{alignment-based}, \textit{composition based},
18 or a combination of both~\cite{parra2007cegma}. The alignment-based
19 method is used when trying to predict a coding gene (\emph{i.e.}.
20 genes that produce proteins) by aligning a genomic DNA sequence with a
21 cDNA sequence coding an homologous protein \cite{parra2007cegma}.
22 This approach is also used in GeneWise\cite{birney2004genewise}. The
23 alternative method, the composition-based one (also known
24 as \textit{ab initio}) is based on a probabilistic model of gene
25 structure to find genes according to the gene value probability
26 (GeneID \cite{parra2000geneid}). Such annotated genomic data will be
27 used to overcome the limitation of the first method described in the
28 previous section. In fact, the second method we propose finds core
29 genes from large amount of chloroplast genomes through genomic
32 Figure~\ref{Fig1} presents an overview of the entire method pipeline.
33 More precisely, the second method consists of three
34 stages: \textit{Genome annotation}, \textit{Core extraction},
35 and \textit{Features Visualization} which highlights the
36 relationships. To understand the whole core extraction process, we
37 describe briefly each stage below. More details will be given in the
38 coming subsections. The method uses as starting point some sequence
39 database chosen among the many international databases storing
40 nucleotide sequences, like the GenBank at NBCI \cite{Sayers01012011},
41 the \textit{EMBL-Bank} \cite{apweiler1985swiss} in Europe
42 or \textit{DDBJ} \cite{sugawara2008ddbj} in Japan. Different
43 biological tools can analyze and annotate genomes by interacting with
44 these databases to align and extract sequences to predict genes. The
45 database in our method must be taken from any confident data source
46 that stores annotated and/or unannotated chloroplast genomes. We have
47 considered the GenBank-NCBI \cite{Sayers01012011} database as sequence
48 database: 99~genomes of chloroplasts were retrieved. These genomes
49 lie in the eleven type of chloroplast families and Table \ref{Tab2}
50 summarizes their distribution in our dataset.
54 \includegraphics[width=0.75\textwidth]{generalView}
55 \caption{A general overview of the annotation-based approach}\label{Fig1}
58 Annotation, which is the first stage, is an important task for
59 extracting gene features. Indeed, to extract good gene feature, a good
60 annotation tool is obviously required. To obtain relevant annotated
61 genomes, two annotation techniques from NCBI and Dogma are used. The
62 extraction of gene feature, the next stage, can be anything like gene
63 names, gene sequences, protein sequences, and so on. Our method
64 considers gene names, gene counts, and gene sequence for extracting
65 core genes and producing chloroplast evolutionary tree. The final
66 stage allows to visualize genomes and/or gene evolution in
67 chloroplast. Therefore we use representations like tables,
68 phylogenetic trees, graphs, etc. to organize and show genomes
69 relationships, and thus achieve the goal of representing gene
70 evolution. In addition, comparing these representations with ones
71 issued from another annotation tool dedicated to large population of
72 chloroplast genomes give us biological perspectives to the nature of
73 chloroplasts evolution. Notice that a local database linked with each
74 pipe stage is used to store all the informations produced during the
77 \input{population_Table}
79 \subsection{Genome annotation techniques}
81 For the first stage, genome annotation, many techniques have been
82 developed to annotate chloroplast genomes. These techniques differ
83 from each others in the number and type of predicted genes (for
84 example: \textit{Transfer RNA (tRNA)} and \textit{Ribosomal RNA
85 (rRNA)} genes). Two annotation techniques from NCBI and Dogma are
86 considered to analyze chloroplast genomes.
88 \subsubsection{Genome annotation from NCBI}
90 The objective is to generate sets of genes from each genome so that
91 genes are organized without any duplication. The input is a list of
92 chloroplast genomes annotated from NCBI. More precisely, all genomes
93 are stored as \textit{.fasta} files which consists in a collection of
94 protein coding genes\cite{parra2007cegma,RDogma} (gene that produce
95 proteins) organized in coding sequences. To be able build the set of
96 core genes, we need to preprocess these genomes
97 using \textit{BioPython} package \cite{chapman2000biopython}. This
98 step starts by converting each genome from FASTA file format to
99 GenVision \cite{geneVision} format from DNASTAR. Each genome is thus
100 converted in a list of genes, with gene names and gene counts. Gene
101 name duplications can be accumulated during the treatment of a genome.
102 These duplications come from gene fragments (\emph{e.g.} gene
103 fragments treated with NCBI) and from chloroplast DNA sequences. To
104 ensure that all the duplications are removed, each list of gene is
105 translated into a set of genes. Note that NCBI genome annotation
106 produces genes except \textit{Ribosomal (rRNA)} genes.
108 \subsubsection{Genome annotation from Dogma}
110 Dogma stands for \textit{Dual Organellar GenoMe Annotator}. It is an
111 annotation tool developed at University of Texas in 2004 for plant
112 chloroplast and animal mitochondrial genomes. This tool has its own
113 database for translating a genome in all six reading frames and
114 queries the amino acid sequence database using
115 BLAST \cite{altschul1990basic} (\emph{i.e.} Blastx) with various
116 parameters. Protein coding genes are identified in an input genome
117 using sequence similarity of genes in Dogma database. In addition in
118 comparison with NCBI annotation tool, Dogma can produce
119 both \textit{Transfer RNAs (tRNA)} and \textit{Ribosomal RNAs (rRNA)},
120 verify their start and end positions. Another difference is also that
121 there is no gene duplication with Dogma after solving gene
122 fragmentation. In fact, genome annotation with Dogma can be the key
123 difference when extracting core genes.
125 The Dogma annotation process is divided into two tasks. First, we
126 manually annotate chloroplast genomes using Dogma web tool. The output
127 of this step is supposed to be a collection of coding genes files for
128 each genome, organized in GeneVision file. The second task is to solve
129 the gene duplication problem and therefore we have use two
130 methods. The first method, based on gene name, translates each genome
131 into a set of genes without duplicates. The second method avoid gene
132 duplication through a defragment process. In each iteration, this
133 process starts by taking a gene from gene list, searches for gene
134 duplication, if a duplication is found, it looks on the orientation of
135 the fragment sequence. If it is positive it appends directly the
136 sequence to gene files. Otherwise reverse complement operations are
137 applied on the sequence, which is then also append to gene files.
138 Finally, a check for missing start and stop codons is performed. At
139 the end of the annotation process, all the genomes are fully
140 annotated, their genes are defragmented, and gene counts are
143 \subsection{Core genes extraction}
145 The goal of this stage is to extract maximum core genes from sets of
146 genes. To find core genes, the following methodology is applied.
148 \subsubsection{Preprocessing}
150 In order to extract core genomes in a suitable manner, the genomic
151 data are preprocessed with two methods: on the one hand a method based
152 on gene name and count, and on the other hand a method based on a
153 sequence quality control test.
155 In the first method, we extract a list of genes from each chloroplast
156 genome. Then we store this list of genes in the database under genome
157 nam and genes counts can be extracted by a specific length command.
158 The \textit{Intersection Core Matrix}, described in next subsection,
159 is then computed to extract the core genes. The problem with this
160 method can be stated as follows: how can we ensure that the gene which
161 is predicted in core genes is the same gene in leaf genomes? The
162 answer to this problem is that if the sequences of any gene in a
163 genome annotated from Dogma and NCBI are similar with respect to a
164 given threshold, then we do not have any problem with this
165 method. When the sequences are not similar we have a problem, because
166 we cannot decide which sequence belongs to a gene in core genes.
168 The second method is based on the underlying idea: we can predict the
169 the best annotated genome by merging the annotated genomes from NCBI
170 and Dogma according to a quality test on genes names and sequences. To
171 obtain all quality genes of each genome, we consider the following
172 hypothesis: any gene will appear in the predicted genome if and only
173 if the annotated genes in NCBI and Dogma pass a specific threshold
174 of \textit{quality control test}. In fact, the Needle-man Wunch
175 algorithm is applied to compare both sequences with respect to a
176 threshold. If the alignment score is above the threshold, then the
177 gene will be retained in the predicted genome, otherwise the gene is
178 ignored. Once the prediction of all genomes is done,
179 the \textit{Intersection Core Matrix} is computed on these new genomes
180 to extract core genes, as explained in Algorithm \ref{Alg3:thirdM}.
183 \caption{Extract new genome based on gene quality test}
186 \REQUIRE $Gname \leftarrow \text{Genome Name}, Threshold \leftarrow 65$
187 \ENSURE $geneList \leftarrow \text{Quality genes}$
188 \STATE $dir(NCBI\_Genes) \leftarrow \text{NCBI genes of Gname}$
189 \STATE $dir(Dogma\_Genes) \leftarrow \text{Dogma genes of Gname}$
190 \STATE $geneList=\text{empty list}$
191 \STATE $common=set(dir(NCBI\_Genes)) \cap set(dir(Dogma\_Genes))$
192 \FOR{$\text{gene in common}$}
193 \STATE $g1 \leftarrow open(NCBI\_Genes(gene)).read()$
194 \STATE $g2 \leftarrow open(Dogma\_Genes(gene)).read()$
195 \STATE $score \leftarrow geneChk(g1,g2)$
196 \IF {$score > Threshold$}
197 \STATE $geneList \leftarrow gene$
204 \textbf{geneChk} is a subroutine used to find the best similarity score between
205 two gene sequences after applying operations like \textit{reverse}, {\it complement},
206 and {\it reverse complement}. Algorithm~\ref{Alg3:genechk} gives the outline of
210 \caption{Find the Maximum Similarity Score between two sequences}
213 \REQUIRE $g1,g2 \leftarrow \text{NCBI gene sequence, Dogma gene sequence}$
214 \ENSURE $\text{Maximum similarity score}$
215 \STATE $score1 \leftarrow needle(g1,g2)$
216 \STATE $score2 \leftarrow needle(g1,Reverse(g2))$
217 \STATE $score3 \leftarrow needle(g1,Complement(g2))$
218 \STATE $score4 \leftarrow needle(g1,Reverse(Complement(g2)))$
219 \RETURN $max(score1,score2,score3,score4)$
223 \subsubsection{Intersection Core Matrix (\textit{ICM})}
225 To extract core genes, we iteratively collect the maximum number of
226 common genes between genomes and therefore during this stage
227 an \textit{Intersection Core Matrix} (ICM) is built. ICM is a two
228 dimensional symmetric matrix where each row and each column correspond
229 to one genome. Hence, an element of the matrix stores
230 the \textit{Intersection Score} (IS): the cardinality of the core
231 genes set obtained by intersecting one genome with another
232 one. Maximum cardinality results in selecting the two genomes having
233 the maximum score. Mathematically speaking, if we have $n$ genomes in
234 local database, the ICM is an $n \times n$ matrix whose elements
237 score_{ij}=\vert g_i \cap g_j\vert
240 \noindent where $1 \leq i \leq n$, $1 \leq j \leq n$, and $g_i, g_j$ are
241 genomes. The generation of a new core gene depends obviously on the
242 value of the intersection scores $score_{ij}$. More precisely, the
243 idea is to consider a pair of genomes such that their score is the
244 largest element in ICM. These two genomes are then removed from matrix
245 and the resulting new core genome is added for the next iteration.
246 The ICM is then updated to take into account the new core gene: new IS
247 values are computed for it. This process is repeated until no new core
248 gene can be obtained.
250 We can observe that the ICM is very large due to the amount of
251 data. As a consequence, the computation of the intersection scores is
252 both time and memory consuming. However, since ICM is a symetric
253 matrix we can reduce the computation overhead by considering only its
254 triangular upper part. The time complexity for this process after
255 enhancement is thus $O(\frac{n.(n-1)}{2})$. Algorithm ~\ref{Alg1:ICM}
256 illustrates the construction of the ICM matrix and the extraction of
257 the core genes, where \textit{GenomeList} represents the database
258 storing all genomes data. At each iteration, it computes the maximum
259 core genes with its two genomes parents.
261 % ALGORITHM HAS BEEN REWRITTEN
264 \caption{Extract Maximum Intersection Score}
267 \REQUIRE $L \leftarrow \text{genomes sets}$
268 \ENSURE $B1 \leftarrow \text{Max Core set}$
269 \FOR{$i \leftarrow 0:len(L)-1$}
270 \STATE $score \leftarrow 0$
271 \STATE $core1 \leftarrow set(GenomeList[L[i]])$
272 \STATE $g1 \leftarrow L[i]$
273 \FOR{$j \leftarrow i+1:len(L)$}
274 \STATE $core2 \leftarrow set(GenomeList[L[j]])$
275 \STATE $Core \leftarrow core1 \cap core2$
276 \IF{$len(Core) > score$}
277 \STATE $score \leftarrow len(Core)$
278 \STATE $g2 \leftarrow L[j]$
281 \STATE $B1[score] \leftarrow (g1,g2)$
287 \subsection{Features visualization}
289 The goal is to visualize results by building a tree of evolution. All
290 core genes generated represent an important information in the tree,
291 because they provide ancestor information of two or more
292 genomes. Each node in the tree represents one chloroplast genome or
293 one predicted core and labelled as \textit{(Genes count:Family name\_Scientific
294 names\_Accession number)}. While an edge is labelled with the number of
295 lost genes from a leaf genome or an intermediate core gene. Such
296 numbers are very interesting because they give an information about
297 the evolution: how many genes were lost between two species whether
298 they belong to the same family or not. By the principle of
299 classification, a small number of genes lost among species indicates
300 that those species are close to each other and belong to same family,
301 while a large lost means that we have an evolutionary relationship
302 between species from different families. To depict the links between
303 species clearly, we built a phylogenetic tree showing the
304 relationships based on the distances among genes sequences. Many tools
305 are available to obtain a such tree, for example:
306 PHYML\cite{guindon2005phyml},
307 RAxML{\cite{stamatakis2008raxml,stamatakis2005raxml}, BioNJ, and
308 TNT\cite{goloboff2008tnt}}. In this work, we chose to use
309 RAxML\cite{stamatakis2008raxml,stamatakis2005raxml} because it is
310 fast, accurate, and can build large trees when dealing with a large
311 number of genomic sequences.
313 The procedure used to built a phylogenetic tree is as follows:
315 \item For each gene in a core gene, extract its sequence and store it in the database.
316 \item Use multiple alignment tools such as (****to be write after see christophe****)
317 to align these sequences with each others.
318 \item we use an outer-group genome from cyanobacteria to calculate distances.
319 \item Submit the resulting aligned sequences to RAxML program to compute the distances and finally draw the phylogenetic tree.
323 \centering \includegraphics[width=0.75\textwidth]{Whole_system}
324 \caption{Overview of the pipeline}\label{wholesystem}
327 \section{Implementation}
328 We implemented the three algorithms to extract maximum core genes from large amount of chloroplast genomes. Table \ref{Etime}, show the annotation, execution time, and the number of core genes for each method:
333 \caption{Annotation, Execution Time, and core genes for each methodology. Annotation means the type of annotation algorithm used to annotate chloroplast genome, Features means the gene features which it is in two types: either gene name, gene sequence, or using the both. The execution time is represented in minute. The number of core genes in the super core is represented with NCBI and DOGMA. Bad genomes: are the number of genomes that can destroy core genes by low number of gene intersection}\label{Etime}
334 \begin{tabular}{ccccccccccc}
336 & \multicolumn{2}{c}{Annotation} & \multicolumn{2}{c}{Features} & \multicolumn{2}{c}{Exec Time} & \multicolumn{2}{c}{Core genes} & \multicolumn{2}{c}{Bad genomes} \\
337 ~ & NCBI & DOGMA & Name & Seq & NCBI & DOGMA & NCBI & DOGMA & NCBI & \multicolumn{1}{c}{DOGMA} \\
339 Gene prediction & $\surd$ & - & - & $\surd$ & ? & - & ? & - & 0 & -\\
340 Gene Features & $\surd$ & $\surd$ & $\surd$ & - & 4.98 & 1.52 & 28 & 10 & 1 & 0\\
341 Gene Quality & $\surd$ & $\surd$ & $\surd$ & $\surd$ & \multicolumn{2}{c}{1.29} & \multicolumn{2}{c}{4} & \multicolumn{2}{c}{1}
349 The second important factor is the amount of memory usage in each methodology. Table \ref{mem} show the amounts of memory consumption by each method.
354 \caption{Annotation, Execution Time, and core genes for each methodology. Annotation means the type of annotation algorithm used to annotate chloroplast genome, Features means the gene features which it is in two types: either gene name, gene sequence, or using the both. The execution time is represented in minute. The number of core genes in the super core is represented with NCBI and DOGMA. Bad genomes: are the number of genomes that can destroy core genes by low number of gene intersection}\label{Etime}
355 \begin{tabular}{cccccccc}
357 & & Load Genomes & T. genevision & Read genevision & ICM & Draw tree & Core Seq. \\
359 Gene prediction & ~ & ~ & ~ & ~ & ~ & ~ & ~\\
360 \multirow{2}{*}{Gene Features} & NCBI & 15.4 & 18.9 & 17.5 & 18 & 18 & 28.1\\
361 & DOGMA& 15.3 & 15.3 & 16.8 & 17.8 & 17.9 & 31.2\\
362 Gene Quality & ~ & 15.3 & >134 & 16.1 & 17 & 17.1 & 24.4
368 one algorithm used to extract core genes based on NCBI annotation, and the others based on NCBI and DOGMA annotation tool. Evolutionary tree generated as a result from each method implementation.
370 \subsection{Extract Core Genes based on Gene Contents}
372 \subsubsection{Core Genes based on NCBI Annotation}
373 The first idea to construct the core genome is based on the extraction of Genes names (as gene presence or absence). For instant, in this stage neither sequence comparison nor new annotation were made, we just want to extract all genes with counts stored in each chloroplast genome, then find the intersection core genes based on gene names. \\
374 The pipeline of extracting core genes can summarize in the following steps according to pre-processing method used:\\
377 \item We downloads already annotated chloroplast genomes in the form of fasta coding genes (\emph{i.e.} \textit{exons}).
378 \item Extract genes names and apply to solve gene duplication using first method.
379 \item Convert fasta file format to geneVision file format to generate ICM.
380 \item Calculate ICM matrix to find maximum core \textit{Score}. New core genes for two genomes will generate and a specific \textit{CoreId} will assign to it. This process continue until no elements remain in the matrix.
381 \item Evolutionary tree will take place by using all data generated from step 1 and 4. The tree will also display the amount of genes lost from each intersection iteration. A specific excel file will be generated that store all the data in local database.
384 There main drawback with this method is genes orthography (e.g two different genes sequences with same gene name). In this case, Gene lost is considered by solving gene duplication based on first method to solve gene duplication.
386 \subsubsection{Core Genes based on Dogma Annotation}
387 The main goal is to get as much as possible the core genes of maximum coding genes names. According to NCBI annotation problem based on \cite{Bakke2009}, annotation method like dogma can give us more reliable coding genes than NCBI. This is because NCBI annotation can carry some annotation and gene identification errors. The general overview of whole process of extraction illustrated in figure \ref{wholesystem}.
389 extracting core genes based on genes names and counts summarized in the following steps:\\
391 \item We apply the genome annotation manually using Dogma annotation tool.
392 \item Analysing genomes to store lists of code genes names (\textit{i.e. exons}). solve gene fragments is done by using first method in solve gene fragments. The output from annotation process with dogma is genomes files in GenVision file format. Sets of genes were stored in the database.
393 \item Generate ICM matrix to calculate maximum core genes.
394 \item Draw the evolutionary tree by extracted all genes sequences from each core. Then applying multiple alignment process on the sequences to calculate the distance among cores to draw a phylogenetic tree.
397 The main drawback from the method of extracting core genes based on gene names and counts is that we can not depending only on genes names because of three causes: first, the genome may have not totally named (This can be found in early versions of NCBI genomes), so we will have some lost sequences. Second, we may have two genes sharing the same name, while their sequences are different. Third, we need to annotate all the genomes.
399 \subsection{Extract Core Genes based on Gene Quality Control}
400 The main idea from this method is to focus on genes quality to predict maximum core genes. By comparing only genes names from one annotation tool is not enough. The question here, does the predicted gene from NCBI is the same gene predicted by Dogma based on gene name and gene sequence?. If yes, then we can predict new quality genomes based on quality control test with a specific threshold. Predicted Genomes comes from merging two annotation techniques. While if no, we can not depending neither on NCBI nor Dogma because of annotation error. Core genes can by predicted by using one of the
402 \subsubsection{Core genes based on NCBI and Dogma Annotation}
403 This method summarized in the following steps:\\
406 \item Retrieve the annotation of all genomes from NCBI and Dogma: in this step, we apply the annotation of all chloroplast genomes in the database using NCBI annotation and Dogma annotation tool.
407 \item Convert NCBI genomes to GeneVision file format, then apply the second method of gene defragmentation methods for NCBI and dogma genomes.
408 \item Predict quality genomes: the process is to pick a genome annotation from two sources, extracting all common genes based on genes names, then applying Needle-man wunch algorithm to align the two sequences based on a threshold equal to 65\%. If the alignment score pass the threshold, then this gene will removed from the competition and store it in quality genome by saving its name with the largest gene sequence with respect to start and end codons. All quality genomes will store in the form of GenVision file format.
409 \item Extract Core genes: from the above two steps, we will have new genomes with quality genes, ofcourse, we have some genes lost here, because dogma produced tRNA and rRNA genes while NCBI did not generate rRNA genes and vise-versa. Build ICM to extract core genes will be sufficient because we already check genes sequences.
410 \item Display tree: An evolution tree then will be display based on the intersections of quality genomes.
