\chapter{Praktická část}

\section{Testování softwarů}
V této části je popsáno testování softwarů, jež jsou uvedeny v kapitole 1.5,  na prokaryotickém genomu \textit{Escherichia coli} O111:H- str. 11128 DNA. Tento genom byl vybrán, jelikož na rozdíl od mnoha ostatních genomů v NCBI databázi, nebyl popsán pomocí automatických anotačních programů. V následujících kapitolách je nejprve popsáno ovládání jednotlivých programů s náhledem výstupních dat a poté je již vyhodnocena jejich účinnost.

\subsection{GeneMark.hmm}

Program GeneMark.hmm je dostupný na webových stránkách \url{http://exon.gatech.edu/gmhmmp.cgi}. Zde můžeme sekvenci, kterou chceme analyzovat, nahrát přímo pomocí textového pole nebo nahrát ze souboru. V obou případech musí být ve~FASTA formátu. Dále volíme ze seznamu organismů ten, který chceme analyzovat, což považuji za nevýhodu programu, protože pokud se náš organismus v seznamu nenalézá, nelze ho analyzovat. V dalším nastavení lze zvolit výstupní formát, dále, co vše chceme, aby bylo součástí výstupu a zda výsledky chceme poslat na e-mail. 

Výstupem je webová stránka, která je zobrazena na obr. \ref{vystup_genemark}. Zde najdeme FASTA hlavičku a dále predikované geny, jejich začátek, konec, délku, na jakém vlákně se nacházejí a třídu, která udává, zda se jedná o typický nebo atypický gen.

\begin{figure}[!h]
  \begin{center}
    \includegraphics[scale=0.6]{obrazky/genemarkhmm_vystup}
  \end{center}
  \caption{Výstupní obrazovka programu GeneMark.hmm.}
\label{vystup_genemark}
\end{figure}

\subsection{Glimmer}
Glimmer ve verzi 3.02 je přístupný na webových stránkách \url{http://www.ncbi.nlm.nih.gov/genomes/MICROBES/glimmer_3.cgi}. Zde opět můžeme použít k nahrání sekvence textové pole nebo lze sekvenci nahrát ze souboru. Sekvence stejně jako u~předcházejícího softwaru musí být ve formátu FASTA. Po nahrání vybereme, o~jaký genetický kód se jedná (Bacteria, Archea nebo Mycoplasma/Spiroplasma) a nakonec zvolíme topologii genomu (lineární nebo kruhovou).

Výstup nalezneme na webové stránce, kde se nachází FASTA hlavička sekvence, dále ID ORF, pozice start a stop kodonu, rámec a skóre, jak vidíme na obr. \ref{vystup_glimmer}. U~indexů pro začátky a konce genů si musíme dát pozor, jelikož na rozdíl od ostatních programů, jsou u komplementárního vlákna prohozené.

\begin{figure}[!h]
  \begin{center}
    \includegraphics[scale=0.5]{obrazky/glimmer_vystup}
  \end{center}
  \caption{Výstupní obrazovka programu Glimmer.}
\label{vystup_glimmer}
\end{figure}

\subsection{Prodigal}
Program je dostupný na webových stránkách \url{http://compbio.ornl.gov/prodigal/server.html}. Po otevření stránky musíme buď načíst sekvenci ve FASTA formátu ze~souboru nebo zadat znaky ručně. Dále volíme translační tabulka a to buď pro~standardní bakterie/Archea nebo pro Mycoplasma/Spiroplasma. Z možností vybereme, co chceme, aby bylo ve výstupu programu: pouze pozice genu, pozice genu s translací proteinu, pozice genu s detailními informacemi o startu nebo kombinaci všech uvedených možností.

Výstupem programu je jako v předchozích programech webová stránka, kterou vidíme na obr. \ref{vystup_prodigal}. Nalezneme zde hlavičku sekvence ve FASTA formátu, pozice jednotlivých genů a zda se nacházejí na vedoucím či komplementárním vlákně.
\clearpage
\begin{figure}[ht!]
  \begin{center}
    \includegraphics[scale=0.7]{obrazky/prodigal_vystup}
  \end{center}
  \caption{Výstupní obrazovka programu Prodigal.}
\label{vystup_prodigal}
\end{figure}

\subsection{Shrnutí výsledků}
Uvedené programy byly testovány na kompletním genomu  \textit{Escherichia coli} O111:H- str. 11128 DNA, který byl získán z databáze NCBI. Organismus \textit{E.coli} obsahuje 5264 anotovaných genů, přičemž v databázi jsou pro každý gen uvedeny indexy začátku a konce. Získané indexy jsme následně vyhledávali ve výsledcích programů pro predikci. Pro hledání pozic byla nastavena tolerance \begin{math}\pm 3\end{math} báze, což odpovídá jednomu kodonu. Výsledky predikce genů pro jednotlivé softwary jsou uvedeny v~tab. \ref{vysledky_softwary} a v tab. \ref{vysledky_softwary2}.

\begin{table}[h!]
\centering
\caption{Výsledky predikce u analyzovaných softwarů.}
\label{vysledky_softwary}
\begin{tabular}{|c|ccc|}
 \hline
 analyzovaný software & GeneMark.hmm & Glimmer & Prodigal \\ \hline
celkový počet nalezených genů       & 5327         & 5504    & 5217    \\
počet správně identifikovaných genů & 4542         & 4386    & 4727    \\
počet správně identifikovaných začátků & 4841 &4727 & 4904 \\
počet správně identifikovaných konců & 4820 & 4696 & 4919 \\
\hline 
\end{tabular}
\end{table}

\begin{table}[h!]
\centering
\caption{Úspěšnost analyzovaných softwarů.}
\label{vysledky_softwary2}
\begin{tabular}{|c|ccc|}
 \hline
 analyzovaný software & GeneMark.hmm & Glimmer & Prodigal \\ \hline
senzitivita [\%] & 86,28         & 83,32    & 89,80    \\
specificita [\%] & 85,26         & 79,69    & 90,61  \\
přesnost [\%] & 85,77         & 81,51    & 90,21    \\
\hline 
\end{tabular}
\end{table}

Senzitivita programů byla spočítána jako počet správně identifikovaných genů dělený počtem genů uvedených v databázi NCBI a specificita jako počet správně identifikovaných genů dělený celkovým počtem nalezených genů pro daný genom. Nakonec byla určena přesnost jednotlivých programů jako součet senzitivity a specificity dělený dvěma. Senzitivita, specificita i přesnost byly nakonec vynásobeny stem, aby výsledné hodnoty vyšly v procentech. Vzorce pro výpočet byly převzaty z~\cite{15}.

Z celkového počtu 5264 genů, nalezly první dva zmíněné programy více genů, než ve skutečnosti genom obsahuje. Tato skutečnost může být dána tím, že ne každý start kodon musí automaticky znamenat start, ale může být i uprostřed genu a zde kódovat aminokyselinu.

Z hlediska přesnosti dosáhl nejvyšší úspěšnosti Prodigal, jehož přesnost vyšla 90,21 \% a poté GeneMark.hmm s 85,77 \%. Přesnost Glimmeru dosáhla nejmenší hodnoty a to 81,51 \%.

\section{Program Gene\_finder}
Program \textit{Gene\_finder} byl vytvořen v prostředí MATLAB R2015a. Při jeho tvorbě byl využit i Bioinformatics toolbox. Úkolem programu je načíst FASTA soubor, jež obsahuje sekvenci prokaryotického organismu a v dané sekvenci najít jednotlivé geny. Výsledkem je tabulka nalezených začátků, konců a délky predikovaných genů a dále, ve kterém čtecím rámci se nacházejí. Pokud geny leží na vedoucím vlákně, je hodnota čtecího rámce kladná, pokud leží na komplementárním vlákně, je hodnota záporná. Výstup programu je ukládán ve formátu 'results\_of\_gene\_finder.xls'. Náhled na výstupní tabulku je zobrazen na obr. \ref{vystup_gene_finder}.

Samotný program se skládá z několika funkcí, jejichž princip je teoreticky popsán v následujících sekcích. Vývojový diagram navrženého algoritmu je uveden v~příloze~\ref{vyvojak}.

\begin{figure}[!h]
  \begin{center}
    \includegraphics[scale=0.8]{obrazky/gene_finder_vystup}
  \end{center}
  \caption{Výstup programu v Excelu.}
\label{vystup_gene_finder}
\end{figure}

\subsection{Načtení sekvence}
V prvním kroku dochází k načtení sekvence pomocí funkce \textit{readfasta}, do které vstupuje název souboru s analyzovanou sekvencí ve formátu 'nazev.fasta'. Výstupem funkce je samotná sekvence a její hlavička.

Následně je pomocí funkce \textit{leading\_and\_complement} vytvořen reverzní komplement sekvence, který je uložen do jedné proměnné společně s hlavní sekvencí. Program \textit{Gene\_finder} probíhá nejprve pro vedoucí vlákno a následně pro reverzní komplement.

\subsection{Hledání stop kodonů}
Do funkce \textit{find\_ends\_in\_frame}, která vyhledává stop kodony, vstupuje sekvence a čtecí rámec, ve kterém chceme hledat. Nejprve provedeme kontrolu, zda je zadán správný čtecí rámec, tj. zda rámec nabývá hodnot 1, 2 nebo 3. Následně provedeme samotné hledání stop kodonů. Začínáme od začátku sekvence. Vždy vezmeme triplet bází a porovnáme jej se zadanými stop kodony (TAA, TAG, TGA). Pokud nalezneme shodu, dojde k uložení pozice první báze stop kodonu. Hledání kodonů končí, když projdeme celou sekvenci. Zároveň je zde ošetřen případ, kdy bychom brali poslední triplet, který by již nebyl úplný (chyběla by jedna nebo dvě báze). 

Výstupem funkce je seznam pozic, které náleží první bázi v nalezených stop kodonech.

\subsection{Hledání start kodonů}
V dalším kroku programu provádíme hledání start kodonů pomocí dvou funkcí a to \textit{find\_gene\_with\_start} a \textit{find\_gene\_in\_frame}. Vstupem první funkce je sekvence, start kodon a pozice dvou stop kodonů, které jsou vedle sebe. Nejprve se provede kontrola, zda zadané start a stop kodony jsou opravdu start a stop kodony. Následně se provádí již samotné vyhledávání start kodonu, které probíhá od konce sekvence. Vezmeme pozici posledního a jemu předcházejícího stop kodonu a mezi nimi provedeme vyhledávání start kodonu. Zde může nastat několik možností, které jsou rozebrány níže.

Pokud žádný start nenalezneme, přesouváme se po vlákně dopředu, tudíž z předcházejícího stop kodonu se stane koncový a poté hledáme mezi ním a jemu předcházejícím stopem.

Pokud nalezneme jeden start kodon, provedeme kontrolu, zda leží ve stejném rámci jako stop kodon. Pokud ne, přesouváme se mezi další dvojici stop kodonů, ale pokud ano, uložíme si pozici první báze start a poslední báze stop kodonu potencionálního genu.

Jestliže se mezi dvěma stop kodony nachází více možných startů, provádíme kontrolu rámce od nejvzdálenějšího a pokud ji nesplňuje, pokračujeme dále směrem blíže ke stop kodonu.

Tímto způsobem projdeme celé vlákno, přičemž když se dostaneme k pozici prvního stop kodonu, start hledáme mezi ním a první bází sekvence. Na výstupu funkce tedy získáme seznam pozic začátků a konců potencionálních genů.

Funkce \textit{find\_gene\_in\_frame} má za vstupy sekvenci, pozice stop kodonů a čtecí rámec. Tato funkce volá funkci předchozí a zároveň provádí hledání startů pro všechny možné start kodony. Nejprve se mezi dvěmi stopy hledá start kodon ATG. Pokud není nalezen, proběhne hledání GTG a pokud se ani ten v daném úseku nenachází, vyhledá se kodon TTG. Tato posloupnost není náhodná, je zvolená tak, aby odpovídala reálné situaci, kde nejčastěji vyskytující se start kodon je právě ATG, poté GTG a nakonec TTG.  Start kodonem mohou být i triplety ATT a CTG, které se ovšem vyskytují velmi vzácně, a proto je ve vyhledávání neuvažujeme.

Na konci funkce dochází k volání funkce \textit{gene\_length\_ok}, která odstraňuje ze seznam potencionálních genů ty, které jsou kratší 200 bází, jelikož existuje pravděpodobnost, že delší otevřené čtecí rámce kódují skutečné geny.

\subsection{Frekvence kodonů}
Následně provedeme výpočet frekvence kodonu skrz funkci \textit{codon\_frequency\_map}, kam vstupuje sekvence a čtecí rámec. V každém čtecím rámci sekvence spočítáme, kolikrát se vyskytuje který kodon. Získaný počet vydělíme počtem všech kodonů a vynásobíme tisícem, čím dosáhneme hodnoty vztažené na tisíc kodonů. Vypočtené frekvence pro jednotlivé kodony použijeme v následujícím kroku.

\subsection{Skórování potencionálních genů}
Abychom mohli z potencionálních genů vybrat finální geny, ohodnotíme je pomocí skórovacího systému, k čemuž využijeme funkci \textit{get\_scoring\_windows} a \textit{gene\_score}.  Potencionální geny převedeme do textových řetězců pomocí funkce \textit{genes\_to\_string}. Abychom zabránili tomu, že delší geny budou mít nižší skóre než kratší geny, provádíme hodnocení v okně. Okno počítáme pomocí první zmíněné funkce, kam vstupují jednotlivé textové řetězce. Délku okna jsme zvolili 48. Tato hodnota byla vybrána tak, aby byla dělitelná třema a dvěma. Také z experimentálních pokusů vycházeli hodnoty nejlépe pro výše zvolenou velikost. Okno se posouvá o půlku své délky. Pomocí funkce vypočítáme, kde okno začíná a jak bude dlouhé. Pokud délka poklesne pod 48, je jasné, že jsme narazili na konec genu a dál již okno posouvat nebudeme. Výstupem funkce je matice, kde v prvním sloupci jsou hodnoty, kde okno začíná a ve druhém jak bude dlouhé.
 
Poté již přistupujeme k samotnému skórování pomocí druhé zmíněné funkce. Sem nám vstupuje textový řetězec, matice s okny a frekvenční mapa. Postupně bereme jednotlivé kodony v rámci okna. Vezmeme první kodon, vyhledáme ho pomocí frekvenční mapy a přiřadíme k němu danou frekvenci. Posuneme se na další kodon, opět najdeme jeho frekvenci a vynásobíme ji s předcházející hodnotou. Takto pokračujeme do konce okna. Okno posuneme o polovinu jeho délky a znovu provádíme jednotlivé výpočty skór. Po projití celého genu určíme celkové skóre a to tak, že vezmeme hodnoty pro jednotlivá okna, která nejprve sečteme a následně je vydělíme celkovým počtem oken pro daný gen.

\subsection{Výběr finálních genů}
Pomocí funkce \textit{apply\_scoring} provedeme výběr genů, které označíme jako finální. Vstupují nám sem potencionální geny a jejich skóre. Pokud je skóre potencionálního genu vyšší, než je průměrná hodnota všech skór, označíme gen jako finální. Při tomto způsobu výběru genů se geny mohou překrývat v rámci čtecích rámců, což odpovídá jejich skutečnému umístění.

\subsection{Uložení finálních genů}
Poslední funkcí \textit{save\_as\_table} spojíme výsledky pro vedoucí a reverzní komplementární vlákno dohromady a výsledek převedeme do tabulky, kterou poté uložíme jako 'results\_of\_gene\_finder.xls'. Zároveň zde probíhá převedení pozic genů na reverzním komplementu, a to do takové podoby, aby měly výsledky stejnou strukturu, jakou mají zmiňované programy.

\section{Testování program Gene\_finder}
Program pro vyhledávání genů v prokaryotických organismech byl testován celkem na 5 genomech: \textit{Escherichia coli} str. K-12 substr. MC4100, \textit{Escherichia coli} strain JEONG-1266, \textit{Escherichia coli} strain MRE600, \textit{Finegoldia magna} ATCC 29328 DNA a \textit{Escherichia coli} O111H-str.11128 DNA, na kterém byly testovány i programy uvedené v předchozích kapitolách. První čtyři zmíněné genomy byly anotovány automaticky, tudíž popis jednotlivých genů nemusí být stoprocentně správný. Všechny uvedené genomy jsou volně dostupné v databázi GenBank pod následujícími přístupovými čísly: HG738867.1, CP014314.1, CP014197.1, NC\_010376.1, AP010960.1.

Jednotlivé genomy s počtem genů, které obsahují a počtem genů, které našel testovaný program, jsou uvedeny v tab. \ref{genomy}.

\begin{table}[h!]
\centering
\caption{Přehled testovaných genomů.}
\label{genomy}
\begin{tabular}{|c|cc|}
 \hline
analyzovaný genom & počet genů&počet nalezených genů\\ \hline
\textit{Escherichia coli} str. K-12 substr. MC4100 &4026& 4999  \\
\textit{Escherichia coli} strain JEONG-1266 &5481&5896 \\
\textit{Escherichia coli} strain MRE600 &4814  &5429\\
\textit{Escherichia coli} O111H-str.11128 DNA &5264&5886  \\
\textit{Finegoldia magna} ATCC 29328 DNA &1645&694  \\
\hline 
\end{tabular}
\end{table}

Vyhledávání bylo prováděno s tolerancí \begin{math}\pm 3\end{math}, \begin{math}\pm 30\end{math}, \begin{math}\pm 60\end{math}, \begin{math}\pm 90\end{math} bází. Pro každý genom byly určeny následující hodnoty: počet celkem nalezených genů, počet správně nalezených genů, počet správně nalezených začátků a konců genů, senzitivita, specificita a přesnost. Na obr. \ref{presnost3}, obr. \ref{presnost30}, obr. \ref{presnost60} a obr. \ref{presnost90} jsou uvedeny grafy pro specificitu, senzitivitu a přesnost predikce pro jednotlivých genomy a to pro všechny tolerance. Kompletní získané výsledky jsou uvedeny v tab. \ref{kompletvysledky1} a tab. \ref{kompletvysledky2}.

Nejlépe si program Gene\_finder vedl u genomu \textit{Escherichia coli} str. K-12 substr. MC410, kde dosáhl nejvyšší přesnosti. Z celkového počtu 4026 genů nalezl celkem 4999, přičemž správně identifikovaných jich bylo při nejvyšší toleranci 2438 a při nejnižší 2042. Přesnost predikce se tedy v tomto případě pohybuje od 46 \% do 55~\%.

U genomu \textit{Escherichia coli} O111H-str.11128 DNA, na kterém byly testovány předešlé zmíněné programy, bylo identifikováno 5886 genů. Ve skutečnosti genom obsahuje 5264 genů. Správně určeno bylo při toleranci \begin{math}\pm 3\end{math} báze 2468 genů a při \begin{math}\pm 90\end{math} bází 2938 genů. Výsledná přesnost se v tomto případě pohybuje o 2 \% níže než v předcházejícím případě.

Nejhůře dopadlo testování \textit{Escherichia coli} strain MRE600. Program nalezl celkem 5429 genů, z toho správně identifikovaných při nejnižší toleranci bylo pouhých 2117 genů. Přesnost zde dosáhla pouhých 41,5 \%.

\begin{figure}[!h]
 \centering
    \includegraphics[scale=0.53]{obrazky/presnost_3}
  \caption{Úspěšnost predikce pro $\pm 3$ báze.}
  \label{presnost3}

    \includegraphics[scale=0.53]{obrazky/presnost_30}
  \caption{Úspěšnost predikce pro $\pm 30$ bází.}
  \label{presnost30}
\end{figure}
\clearpage

\begin{figure}[!h]
 \centering
    \includegraphics[scale=0.53]{obrazky/presnost_60}
  \caption{Úspěšnost predikce pro $\pm 60$ bází.}
  \label{presnost60}

    \includegraphics[scale=0.53]{obrazky/presnost_90}
  \caption{Úspěšnost predikce pro $\pm 90$ bází.}
  \label{presnost90}
\end{figure}
\clearpage

Z uvedených výsledků je patrné, že program najde větší množství genů než se ve skutečnosti v genomu objevuje. To může být zapříčeněno tím, že pokud jsme nenašli start kodon ATG, automaticky jsme začali hledat alternativní start kodony a neuvažovali jsme, že by se mohlo jednat o nekódující oblast. Dále z výsledků vidíme, že nachází správně skoro stejný počet začátků a konců genů. Ovšem tyto hodnoty jsou nižší než je celkový počet genů v organismu, a proto by bylo zapotřebí dosáhnout jejich zvýšení.

\begin{table}[h!]
\centering
\caption{Srovnání analyzovaných a vytvořeného softwaru.}
\label{srovnani}
\begin{tabular}{|c|cccc|}
 \hline
 analyzovaný software & GeneMark.hmm & Glimmer & Prodigal &Gene\_finder \\ \hline
přesnost [\%] & 85,77         & 81,51    & 90,21   & 44,41\\
\hline 
\end{tabular}
\end{table}

Jak vidíme v tab. \ref{srovnani}, za pomoci základního námi vytvořeného programu lze dosáhnout 40 - 50 \% úspěšnosti, z čehož ve srovnání s ostatními plyne, že pro lepší výsledky je potřeba využít mnohem pokročilejších algoritmů a metod, které přesahují obsah této práce.













