annotated.tex

   1
   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}.
   8
   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
  30 features extraction.
  31
  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.
  51
  52 \begin{figure}[h]
  53   \centering
  54     \includegraphics[width=0.75\textwidth]{generalView}
  55 \caption{A general overview of the annotation-based approach}\label{Fig1}
  56 \end{figure}
  57
  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. The extraction of gene feature, the  next stage, can be anything like gene names,  gene  sequences, protein  sequences,  and  so  on. Our  method considers gene  names, gene counts,  and gene sequence  for extracting core  genes and  producing  chloroplast evolutionary  tree. The  final stage   allows  to   visualize  genomes   and/or  gene   evolution  in chloroplast.    Therefore   we   use  representations   like   tables, phylogenetic  trees,  graphs,  etc.   to  organize  and  show  genomes relationships,  and  thus  achieve   the  goal  of  representing  gene
  61 evolution.   In addition,  comparing these  representations  with ones issued from  another annotation tool dedicated to  large population of chloroplast genomes  give us biological perspectives to  the nature of chloroplasts evolution. Notice that  a local database linked with each pipe stage is  used to store all the  informations produced during the process.
  62
  63 \input{population_Table}
  64
  65 \subsection{Genome annotation techniques}
  66
  67 To obtain  relevant annotated genomes, two annotation  techniques from NCBI and Dogma  are used. For  the first  stage, genome  annotation, many  techniques  have been developed  to annotate chloroplast  genomes.  These  techniques differ
  68 from  each others  in  the number  and  type of  predicted genes  (for example:  \textit{Transfer  RNA   (tRNA)}  and  \textit{Ribosomal  RNA (rRNA)}  genes). Two  annotation techniques  from NCBI  and  Dogma are considered to analyze chloroplast genomes.
  69
  70 \subsubsection{Genome annotation from NCBI}
  71
  72 The objective  is to generate sets  of genes from each  genome so that
  73 genes are organized  without any duplication.  The input  is a list of
  74 chloroplast genomes  annotated from NCBI. More  precisely, all genomes
  75 are stored as \textit{.fasta} files  which consists in a collection of
  76 protein  coding genes\cite{parra2007cegma,RDogma}  (gene  that produce
  77 proteins) organized in coding sequences.   To be able build the set of
  78 core    genes,     we    need    to     preprocess    these    genomes
  79 using  \textit{BioPython}  package \cite{chapman2000biopython}.   This
  80 step  starts by  converting  each  genome from  FASTA  file format  to
  81 GenVision \cite{geneVision}  format from DNASTAR. Each  genome is thus
  82 converted in  a list of genes,  with gene names and  gene counts. Gene
  83 name duplications can be accumulated during the treatment of a genome.
  84 These  duplications  come   from  gene  fragments  (\emph{e.g.}   gene
  85 fragments treated  with NCBI) and  from chloroplast DNA  sequences. To
  86 ensure that  all the  duplications are removed,  each list of  gene is
  87 translated  into a  set of  genes.  Note that  NCBI genome  annotation
  88 produces genes except \textit{Ribosomal (rRNA)} genes.
  89
  90 \subsubsection{Genome annotation from Dogma}
  91
  92 Dogma stands for \textit{Dual  Organellar GenoMe Annotator}.  It is an
  93 annotation tool  developed at  University of Texas  in 2004  for plant
  94 chloroplast and  animal mitochondrial genomes.  This tool  has its own
  95 database  for translating  a  genome  in all  six  reading frames  and
  96 queries     the     amino     acid     sequence     database     using
  97 BLAST  \cite{altschul1990basic}  (\emph{i.e.}   Blastx)  with  various
  98 parameters.  Protein  coding genes are  identified in an  input genome
  99 using sequence similarity of genes  in Dogma database.  In addition in
 100 comparison   with   NCBI    annotation   tool,   Dogma   can   produce
 101 both \textit{Transfer RNAs (tRNA)} and \textit{Ribosomal RNAs (rRNA)},
 102 verify their start and end  positions. further more, there is no gene duplication with gene annotations from Dogma after applying gene de-fragmentation process. In  fact, genome annotation with Dogma can be the key difference when extracting core genes.
 103
 104 The Dogma  annotation process  is divided into  two tasks.   First, we
 105 manually annotate chloroplast genomes using Dogma web tool. The output
 106 of this step is supposed to  be a collection of coding genes files for
 107 each genome, organized in GeneVision file. The second task is to solve
 108 the  gene   duplication  problem  and   therefore  we  have   used  two
 109 methods. The first method, based  on gene name, translates each genome
 110 into a set  of genes without duplicates. The  second method avoid gene
 111 duplication  through a  defragment  process. In  each iteration,  this
 112 process  starts by taking  a gene  from gene  list, searches  for gene
 113 duplication, if a duplication is found, it looks on the orientation of
 114 the  fragment sequence.   If it  is positive  it appends  directly the
 115 sequence to  gene files.  Otherwise reverse  complement operations are
 116 applied  on the sequence,  which is  then also  append to  gene files.
 117 Finally, a  check for missing start  and stop codons  is performed. At
 118 the  end  of  the  annotation  process,  all  the  genomes  are  fully
 119 annotated,  their   genes  are  defragmented,  and   gene  counts  are
 120 available.
 121
 122 \subsection{Core genes extraction}
 123
 124 The goal of  this stage is to extract maximum core  genes from sets of
 125 genes.  To find core genes, the following methodology is applied.
 126
 127 \subsubsection{Preprocessing}
 128
 129 In order  to extract  core genomes in  a suitable manner,  the genomic
 130 data are preprocessed with two methods: on the one hand a method based
 131 on gene  name and count,  and on  the other hand  a method based  on a
 132 sequence quality control test.
 133
 134 In the first method, we extract  a list of genes from each chloroplast
 135 genome.  Then we store this list of genes in the database under genome
 136 nam and  genes counts can be  extracted by a  specific length command.
 137 The \textit{Intersection  Core Matrix}, described  in next subsection,
 138 is then  computed to  extract the core  genes.  The problem  with this
 139 method can be stated as follows: how can we ensure that the gene which
 140 is  predicted in  core genes  is the  same gene  in leaf  genomes? The
 141 answer  to this problem  is that  if the  sequences of  any gene  in a
 142 genome annotated  from Dogma  and NCBI are  similar with respect  to a
 143 given  threshold,  the method is operational when the sequences are not similar. The problem of attribution of a sequence to a gene in the core genome come to light.
 144
 145 The second method is based on  the underlying idea that it is possible to predict the the best annotated  genome by merging the annotated  genomes from NCBI
 146 and Dogma according to a quality test on genes names and sequences. To
 147 obtain all  quality genes  of each genome,  we consider  the following
 148 hypothesis: any gene  will appear in the predicted  genome if and only
 149 if the  annotated genes  in NCBI and  Dogma pass a  specific threshold
 150 of  \textit{quality  control test}.   In  fact,  the Needle-man  Wunch
 151 algorithm  is applied  to compare  both  sequences with  respect to  a
 152 threshold. If  the alignment  score is above  the threshold,  then the
 153 gene will be  retained in the predicted genome,  otherwise the gene is
 154 ignored.   Once    the   prediction   of   all    genomes   is   done,
 155 the \textit{Intersection Core Matrix} is computed on these new genomes
 156 to extract core genes, as explained in Algorithm \ref{Alg3:thirdM}.
 157
 158 \begin{algorithm}[H]
 159 \caption{Extract new genome based on gene quality test}
 160 \label{Alg3:thirdM}
 161 \begin{algorithmic}
 162 \REQUIRE $Gname \leftarrow \text{Genome Name}, Threshold \leftarrow 65$
 163 \ENSURE $geneList \leftarrow \text{Quality genes}$
 164 \STATE $dir(NCBI\_Genes) \leftarrow \text{NCBI genes of Gname}$
 165 \STATE $dir(Dogma\_Genes) \leftarrow \text{Dogma genes of Gname}$
 166 \STATE $geneList=\text{empty list}$
 167 \STATE $common=set(dir(NCBI\_Genes)) \cap set(dir(Dogma\_Genes))$
 168 \FOR{$\text{gene in common}$}
 169         \STATE $g1 \leftarrow open(NCBI\_Genes(gene)).read()$
 170         \STATE $g2 \leftarrow open(Dogma\_Genes(gene)).read()$
 171         \STATE $score \leftarrow geneChk(g1,g2)$
 172         \IF {$score > Threshold$}
 173                 \STATE $geneList \leftarrow gene$
 174         \ENDIF
 175 \ENDFOR
 176 \RETURN $geneList$
 177 \end{algorithmic}
 178 \end{algorithm}
 179
 180 \textbf{geneChk} is a subroutine used to find the best similarity score between
 181 two gene sequences after applying operations like \textit{reverse}, {\it complement},
 182 and {\it reverse complement}. Algorithm~\ref{Alg3:genechk} gives the outline of
 183 geneChk subroutine.
 184
 185 \begin{algorithm}[H]
 186 \caption{Find the Maximum Similarity Score between two sequences}
 187 \label{Alg3:genechk}
 188 \begin{algorithmic}
 189 \REQUIRE $g1,g2 \leftarrow \text{NCBI gene sequence, Dogma gene sequence}$
 190 \ENSURE $\text{Maximum similarity score}$
 191 \STATE $score1 \leftarrow needle(g1,g2)$
 192 \STATE $score2 \leftarrow needle(g1,Reverse(g2))$
 193 \STATE $score3 \leftarrow needle(g1,Complement(g2))$
 194 \STATE $score4 \leftarrow needle(g1,Reverse(Complement(g2)))$
 195 \RETURN $max(score1,score2,score3,score4)$
 196 \end{algorithmic}
 197 \end{algorithm}
 198
 199 \subsubsection{Intersection Core Matrix (\textit{ICM})}
 200
 201 To extract  core genes, we  iteratively collect the maximum  number of
 202 common  genes   between  genomes  and  therefore   during  this  stage
 203 an \textit{Intersection  Core Matrix}  (ICM) is built.   ICM is  a two
 204 dimensional symmetric matrix where each row and each column correspond
 205 to   one   genome.   Hence,   an   element   of   the  matrix   stores
 206 the  \textit{Intersection Score}  (IS):  the cardinality  of the  core
 207 genes   set  obtained   by  intersecting   one  genome   with  another
 208 one. Maximum  cardinality results in selecting the  two genomes having
 209 the maximum score. Mathematically speaking, if we have $n$ genomes in
 210 local database, the ICM is an $n \times n$ matrix whose elements
 211 satisfy:
 212 \begin{equation}
 213 score_{ij}=\vert g_i \cap g_j\vert
 214 \label{Eq1}
 215 \end{equation}
 216 \noindent where $1 \leq i \leq n$, $1 \leq j \leq n$, and $g_i, g_j$ are
 217 genomes. The  generation of a new  core gene depends  obviously on the
 218 value  of the  intersection scores  $score_{ij}$. More  precisely, the
 219 idea is  to consider a  pair of genomes  such that their score  is the
 220 largest element in ICM. These two genomes are then removed from matrix
 221 and the  resulting new  core genome is  added for the  next iteration.
 222 The ICM is then updated to take into account the new core gene: new IS
 223 values are computed for it. This process is repeated until no new core
 224 gene can be obtained.
 225
 226 We  can observe  that  the ICM  is very  large  due to  the amount  of
 227 data. As a consequence, the  computation of the intersection scores is
 228 both  time and  memory consuming.  However,  since ICM  is a  symetric
 229 matrix we can reduce the  computation overhead by considering only its
 230 triangular  upper part.  The  time complexity  for this  process after
 231 enhancement is thus $O(\frac{n.(n-1)}{2})$.  Algorithm ~\ref{Alg1:ICM}
 232 illustrates the construction  of the ICM matrix and  the extraction of
 233 the  core  genes, where  \textit{GenomeList}  represents the  database
 234 storing all genomes  data. At each iteration, it  computes the maximum
 235 core genes with its two genomes parents.
 236
 237 % ALGORITHM HAS BEEN REWRITTEN
 238
 239 \begin{algorithm}[H]
 240 \caption{Extract Maximum Intersection Score}
 241 \label{Alg1:ICM}
 242 \begin{algorithmic}
 243 \REQUIRE $L \leftarrow \text{genomes sets}$
 244 \ENSURE $B1 \leftarrow \text{Max Core set}$
 245 \FOR{$i \leftarrow 1:len(L)-1$}
 246         \STATE $score \leftarrow 0$
 247         \STATE $core1 \leftarrow set(GenomeList[L[i]])$
 248         \STATE $g1 \leftarrow L[i]$
 249         \FOR{$j \leftarrow i+1:len(L)$}
 250                 \STATE $core2 \leftarrow set(GenomeList[L[j]])$
 251                 \STATE $core \leftarrow core1 \cap core2$
 252                 \IF{$len(core) > score$}
 253                   \STATE $score \leftarrow len(core)$
 254                   \STATE $g2 \leftarrow L[j]$
 255                 \ENDIF
 256         \ENDFOR
 257         \STATE $B1[score] \leftarrow (g1,g2)$
 258 \ENDFOR
 259 \RETURN $max(B1)$
 260 \end{algorithmic}
 261 \end{algorithm}
 262
 263 \subsection{Features visualization}
 264
 265 The goal is to visualize results  by building an evolutionary tree. All
 266 core  genes generated  represent  an important information  in the  tree,
 267 because they  provide ancestor information of two  or more
 268 genomes. Each  node in the  tree represents one chloroplast  genome or
 269 one predicted core and labelled as \textit{(Genes count:Family name\_Scientific
 270 names\_Accession number)}. While an edge is labelled with the number of
 271 lost  genes from  a leaf  genome or  an intermediate  core  gene. Such
 272 numbers are  very interesting because  they give an  information about
 273 evolution:  how many genes  were lost between two  species whether
 274 they  belong  to  the  same  lineage  or not. To depict  the links between
 275 species   clearly,  we   built   a  phylogenetic   tree  showing   the
 276 relationships based on the distances among genes sequences. Many tools
 277 are    available   to    obtain    a   such    tree,   for    example:
 278 PHYML\cite{guindon2005phyml},
 279 RAxML{\cite{stamatakis2008raxml,stamatakis2005raxml},    BioNJ,    and
 280 TNT\cite{goloboff2008tnt}}.    In   this  work,   we   chose  to   use
 281 RAxML\cite{stamatakis2008raxml,stamatakis2005raxml}   because   it  is
 282 fast, accurate,  and can build large  trees when dealing  with a large
 283 number of genomic sequences.
 284
 285 The procedure used to built a phylogenetic tree is as follows:
 286 \begin{enumerate}
 287 \item For each gene in a core gene, extract its sequence and store it in the database.
 288 \item Use multiple alignment tools such as (****to be write after see christophe****)
 289 to align these sequences with each others.
 290 \item Use an outer-group genome from cyanobacteria to calculate distances.
 291 \item Submit the resulting aligned sequences to RAxML program to compute
 292 the distances and finally draw the phylogenetic tree.
 293 \end{enumerate}
 294
 295 \begin{figure}[H]
 296   \centering \includegraphics[width=0.75\textwidth]{Whole_system} \caption{Overview
 297   of the pipeline}\label{wholesystem}
 298 \end{figure}
 299
 300 \section{Implementation}
 301
 302 All the different  algorithms have  been implemented using  Python on a personal computer running Ubuntu~12.04 with 6~GiB memory and
 303 a  quad-core Intel  core~i5~processor with  an operating  frequency of
 304 2.5~GHz. All the programs can be downloaded at \url{http://......} .
 305 genes  from large  amount of  chloroplast  genomes.
 306
 307 \begin{center}
 308 \begin{table}[H]
 309 \caption{Type of annotation, execution time, and core genes.}\label{Etime}
 310 {\scriptsize
 311 \begin{tabular}{p{2cm}p{0.5cm}p{0.25cm}p{0.5cm}p{0.25cm}p{0.5cm}p{0.25cm}p{0.5cm}p{0.25cm}p{0.5cm}p{0.2cm}}
 312 \hline\hline
 313  Method & \multicolumn{2}{c}{Annotation} & \multicolumn{2}{c}{Features} & \multicolumn{2}{c}{Exec. time (min.)} & \multicolumn{2}{c}{Core genes} & \multicolumn{2}{c}{Bad genomes} \\
 314 ~ & N & D & Name & Seq & N & D & N & D & N & D \\
 315 \hline
 316 Gene prediction & $\surd$ & - & - & $\surd$ & 1.7 & - & ? & - & 0 & -\\[0.5ex]
 317 Gene Features & $\surd$ & $\surd$ & $\surd$ & - & 4.98 & 1.52 & 28 & 10 & 1 & 0\\[0.5ex]
 318 Gene Quality & $\surd$ & $\surd$ & $\surd$ & $\surd$ & \multicolumn{2}{c}{$\simeq$3 days + 1.29} & \multicolumn{2}{c}{4} & \multicolumn{2}{c}{1}\\[1ex]
 319 \hline
 320 \end{tabular}
 321 }
 322 \end{table}
 323 \end{center}
 324
 325 \vspace{-1cm}
 326
 327 Table~\ref{Etime}  presents  for  each  method  the  annotation  type,
 328 execution time,  and the  number of core  genes. We use  the following
 329 notations:  \textbf{N}  denotes NCBI,  while  \textbf{D} means  DOGMA,
 330 and \textbf{Seq}  is for sequence. The first two {\it Annotation} columns
 331 represent the algorithm used to annotate chloroplast genomes. The next two ones {\it
 332 Features} columns mean  the kind  of gene feature used to extract core
 333 genes: gene name, gene sequence, or  both of them. It can be seen that
 334 almost all methods need low {\it Execution time} expended in minutes to extract core genes
 335 from the large set of chloroplast genomes. Only the gene quality method requires
 336 several days of computation (about 3-4 days) for sequence comparisons. However,
 337 once the quality genomes are well constructed, it only takes 1.29~minutes to
 338 extract core gene. Thanks to this low execution times that gave us a privilege to use these
 339 methods to extract core genes  on a personal computer rather than main
 340 frames or parallel computers. The lowest execution time: 1.52~minutes,
 341 is obtained with the second method using Dogma annotations. The number
 342 of {\it  Core genes} represents the  amount of genes in  the last core
 343 genome. The main goal is to  find the maximum core genes that simulate
 344 biological background of chloroplasts. With  NCBI we have 28 genes for
 345 96   genomes,   instead   of    10   genes   for   97   genomes   with
 346 Dogma. Unfortunately, the biological distribution of genomes with NCBI
 347 in core tree do not  reflect good biological perspective, whereas with
 348 DOGMA the  distribution of genomes is biologically  relevant. Some a few genomes maybe destroying core genes due to
 349 low  number  of  gene  intersection. More precisely, \textit{NC\_012568.1  Micromonas pusilla} is the only genome who destroyes the core genome with NCBI
 350 annotations for both gene features and gene quality methods.
 351
 352 The second important factor is the amount of memory nessecary in each
 353 methodology.   Table   \ref{mem}  shows  the  memory   usage  of  each
 354 method. In this table, the values are  presented in megabyte
 355 unit and \textit{gV} means  genevision~file~format. We can notice that
 356 the level  of memory which is  used is relatively low  for all methods
 357 and is available  on any personal computer. The  different values also
 358 show that the gene features  method based on Dogma annotations has the
 359 more   reasonable   memory   usage,   except  when   extracting   core
 360 sequences. The third method gives the lowest values if we already have
 361 the   quality   genomes,   otherwise   it  will   consume   far   more
 362 memory. Moreover, the  amount of memory, which is used by the third method also
 363 depends on the size of each genome.
 364
 365
 366 \begin{table}[H]
 367 \centering
 368 \caption{Memory usages in (MB) for each methodology}\label{mem}
 369 \tabcolsep=0.11cm
 370 {\scriptsize
 371 \begin{tabular}{p{2.5cm}@{\hskip 0.1mm}p{1.5cm}@{\hskip 0.1mm}p{1cm}@{\hskip 0.1mm}p{1cm}@{\hskip 0.1mm}p{1cm}@{\hskip 0.1mm}p{1cm}@{\hskip 0.1mm}p{1cm}@{\hskip 0.1mm}p{1cm}}
 372 \hline\hline
 373 Method& & Load Gen. & Conv. gV & Read gV & ICM & Core tree & Core Seq. \\
 374 \hline
 375 Gene prediction & NCBI & 108 & - & - & - & - & -\\
 376 \multirow{2}{*}{Gene Features} & NCBI & 15.4 & 18.9 & 17.5 & 18 & 18 & 28.1\\
 377               & DOGMA& 15.3 & 15.3 & 16.8 & 17.8 & 17.9 & 31.2\\
 378 Gene Quality  & ~ & 15.3 & $\le$3G & 16.1 & 17 & 17.1 & 24.4\\
 379 \hline
 380 \end{tabular}
 381 }
 382 \end{table}
 383
 384
 385
 386
 387
 388