1. Úvod Buňka je základní stavební jednotkou všech živých organismů (mimo virů) a bez ohledu na to, že buňky mohou nabývat nejrůznější tvary a velikosti, se při jejich členění omezujeme na dvě základní třídy. První třídou jsou prokaryotické buňky. Tyto malé buňky jsou základem jednobuněčných (prokaryotických) organismů a na rozdíl od eukaryotických buněk, jež tvoří třídu druhou, jsou relativně primitivní. Kromě mnohých rozdílů ve struktuře, které nejsou předmětem této práce, se eukaryotické buňky od prokaryotických odlišují na mnohem základnější, genetické bázi, přesněji v genové expresi. Geny eukaryotických organismů obsahují kromě kódujících sekvencí (exonů) i části nekódující (introny). Tyto introny jsou při transkripčních úpravách vystříhávány (děje se tak při operaci zvané „splicing“ - sestřih), čímž dojde ke spojení exonů a vzniku funkční mRNA, která je následně využita při translaci ke vzniku proteinů. Exony tedy obsahují části sekvence, které kódují nějakou specifickou část určitého proteinu. Během evoluce dochází k tzv. genetické rekombinaci, při které se pohybují bloky DNA v rámci genomu organismu. Spekuluje se, že při těchto procesech dochází ke kombinaci domén a jiných funkčních částí v rámci proteinů, což by mohl být následek pohybu exonů mezi jednotlivými geny. Tento mechanismus pro pohyb domén byl poprvé představen v roce 1978 Walterem Gilbertem v rámci jeho teorie o pohybech exonů [1]. Tato teorie naznačovala, že by každý exon mohl kódovat v proteinu jednu doménu, což by znamenalo, že kombinací exonů z různých genů by se mohlo dosáhnout vytvoření široké palety proteinů, jež by v sobě obsahovaly libovolnou kombinaci funkčních vlastností a domén z existujících proteinů. V roce 2005 byla publikována studie, která zkoumala možnou korelaci mezi doménami proteinů a exony. Jejím výsledkem bylo, že některé druhy domén s hranicemi exonů korelovaly velice silně, kdežto jiné typy jako např. immunoglobulin tuto vlastnost nevykazoval [2]. Cílem mé bakalářské práce bylo vytvořit automatizovaný nástroj, který by mapoval hranice exonů na sekvence proteinů a tyto hranice následně vyhodnotil oproti hranicím domén a sekundárních struktur v rámci jednotlivých proteinů a poskytl odpověď na otázku, nakolik exon-intron struktura genů koreluje s organizací domén a sekundárních struktur v proteinech. K dosažení tohoto cíle program využívá záznamy z veřejně dostupných databází genů, mRNA sekvencí a proteinů, jejichž popisem se zabývá první část této práce. Druhá část se zaobírá návrhem a implementací programu. Zde představím jednotlivé módy, ve kterých program pracuje, nastíním několik problémů a popíši mnou zvolená řešení. Ve třetí a závěrečné části píši o výsledcích, jež program dosahuje, a o jeho přínosu k řešení dané problematiky.
1
2. Veřejně dostupná biologická data Bioinformatika je poměrně mladá vědní disciplína; samotný její název vznikl na konci sedmdesátých let. I když její primární cíl zůstává stále stejný, tedy lepší pochopení biologických procesů, je nutno dodat, že za 30 let její existence se tvář informatiky značně změnila a nemalou zásluhu na tom má rapidní rozvoj a rozšíření Internetu a také velice rychlý růst (v některých případech až exponenciální) dat obsažených v biologických a bioinformatických databázích. Pro nalezení optimálního řešení mnoha problémů, se kterými se potýkají např. molekulární biologové, jsou tyto databáze doslova nepostradatelné. Jmenujme například návrh nových léčiv nebo predikci struktury proteinů, kde je téměř nutnost vycházet ze znalosti již známých struktur proteinů. Biologické databáze jsou ve své podstatě velké knihovny záznamů, které byly shromážděny z nejrůznějších zdrojů, jakými jsou vědecké pokusy, publikace a odborné články nebo také automatické predikční nástroje. Protože současných biologických a bioinformatických dat je nepřeberné množství, jsou tato data rozdělena do specializovaných veřejných databází, které jsou často vzájemně propojeny. Takové databáze poté můžeme rozdělit do několika logických skupin, např. databáze genomů (Ensembl, ERIC, CAMERA), databáze sekvencí proteinů (PDB, Swiss-Prot, UniProt), databáze struktury proteinů (CATH, SCOP), databáze metabolických cest (KEGG, MANET) a další. Existují také vysoce specializované databáze, jako je například databáze PhenCode, jež sbírá data pro modelování a lepší pochopení vztahů mezi genotypem a fenotypem u lidí [3]. Speciálním typem databází jsou tzv. metadatabáze, které v sobě ukrývají hned několik různých databází, a záznamy z těchto databází zpřístupňují v nějaké jednotné formě. Typickým zástupcem takových databází je Entrez.
2.1 Entrez Jak už bylo řečeno, Entrez je metadatabáze, tedy databáze, která zpřístupňuje vetší množství na sobě nezávislých databází. I když se o Entrezu zcela běžně vyjadřujeme jako o databázi, není tato definice zcela korektní. Mnohem přesnější by bylo mluvit o Entrezu jako o vyhledávači a softwarovém rozhraní pro získávání informací z biologických databází. Toto rozhraní bylo vytvořeno Národním centrem pro biotechnologické informace (dále jen NCBI), jež je součástí Národního zdravotního institutu (NIH). NCBI bylo založeno v roce 1988 a od svého vzniku se zaměřuje na shromažďování výsledků sekvenací DNA. Kromě toho NCBI vytváří veřejně dostupné databáze, provádí výzkum v oblasti výpočetní biologie a vytváří různé softwarové nástroje pro zkoumání biologických dat. Entrez je právě jedním z těchto druhů nástrojů. Poskytuje rozhraní k vyhledávání dat ve více než 25 databázích. Mezi tyto databáze patří například databáze Gene poskytující informace ohledně genů, databáze Protein a Nucleotide, které obsahují data týkající se sekvencí nebo CDD, což je databáze konzervovaných domén, a mnohé další. Myšlenka Entrezu se zakládá na existenci logických vztahů mezi jednotlivými záznamy z různých veřejných databází [4]. Například předpokládejme, že existuje odborný článek, jenž je evidován v databázi PubMed, který pojednává o určitém specifickém genu, jehož záznam se nachází v databázi Gene a jeho sekvence v databázi Nucleotide. Námi zvolený gen dále může kódovat nějaký protein a informace o sekvenci tohoto proteinu je přístupná z databáze Protein; tento protein může obsahovat konzervovaní domény a informace o těchto doménách můžeme vyhledat
2
za použití databáze CDD atd. K získání všech těchto informací nám stačí pouze znalost výše zmíněného článku. Jak je již z mého popisu patrné, Entrez je opravdu silným nástrojem, což je také důvod jeho značné oblíbenosti mezi bioinformatiky a velmi častého využívání v rámci nejrůznějších výzkumů a projektů. Kromě webového rozhraní k Entrezu, jež je dostupné přímo z domovské stránky NCBI na adrese http://www.ncbi.nlm.nih.gov/, NCBI též poskytuje tzv. Entrez Programming Utillities, což je nástroj, díky kterému se dá přistoupit k databázím pomocí zaslání speciálně naformátovaného URL dotazu NCBI serveru. V tomto URL dotazu uživatel specifikuje, k jaké databázi chce přistoupit (v současné době je jich přes tyto dotazy přístupných 14). Jako další údaj musí být uvedeno primární ID záznamu, který chceme, aby nám byl zaslán. URL dotaz obsahuje také několik volitelných položek jako je například email.adresa, přes kterou nás může NCBI kontaktovat, pokud by se vyskytl problém s našimi dotazy. Velice užitečná je možnost vybrat si v URL dotazu formát, v jakém nám bude námi vyžádaný záznam zaslán (např. Fasta, GenBank formát, INSDSeq strukturovaný soubor a další). Typický URL dotaz může pak vypadat následovně: http://www.ncbi.nlm.nih.gov/entrez/eutilsefetch.fcgi?db=protein&id=NP_001079&rettype=gb V tomto dotazu žádám o zaslání záznamu proteinu s identifikátorem NP_001079 z databáze Protein, navíc specifikuji, že by tento záznam měl být poslán ve formátu GenBank. U podobně řešených služeb vždy existuje riziko přetížení serveru, a to jak neúmyslné, tak úmyslně zapříčiněné. Z tohoto důvodu NCBI omezuje počet dotazů od jednoho uživatele na 3 za 1 sekundu [5].
2.2 CDD The Conserved Domain Database (dále jen „CDD“) je jednou z databází vyvinutou NCBI, jež je také zpřístupněna přes rozhraní Entrez. CDD obsahuje informace o konzervovaných doménách, což jsou samostatné jednotky identifikovatelné v terciální struktuře proteinu. Tyto jednotky jsou často spojeny s určitou specifickou funkcí v proteinu a jejich identifikace v sekvenci proteinu může být prvním krokem k pochopení funkce proteinu jako celku [6]. Protein, mimo jiné i v závislosti na své délce, může obsahovat až několik desítek takových konzervovaných domén, ale také nemusí obsahovat žádnou. Databáze CDD získává data o doménách z více různých zdrojů. Hlavním zdrojem je databáze Pfam, což je největší a stále se rozrůstající databáze domén, které rozděluje do rodin podle sekvenční podobnosti. Pfam se skládá ze dvou částí. Pfam-A je manuálně upravovaná část databáze obsahující více než 9000 záznamů. Pro každý záznam má Pfam-A uloženo zarovnání a Markovův model. Druhá část nese označení Pfam-B a obsahuje velké množství automaticky generovaných malých rodin domén. Podle zaznamenaných statistik pro 74 % sekvencí proteinů existuje alespoň 1 záznam v databázi Pfam [7]. Další zdroje, ze kterých CDD získává data jsou databáze a kolekce domén SMART, COGs, Protein Clusters a PRK [8]. Nevýhodou při využívání více databází je jejich značná redundance, proto CDD všechny získané záznamy shlukuje do velkých superrodin, čímž se redundanci v záznamech vyhýbá. Kromě samotné databáze CDD obsahuje interaktivní službu CD-Search (Conserved Domain Search). Tato služba je webový nástroj pro detekci konzervovaných domén v sekvencích proteinů. Uživatel musí zadat sekvenci ve formátu FASTA, případně identifikátor proteinu z databáze Protein
3
a vybere si, jakou databázi chce využít k hledání domén. Standardní volbou je databáze CDD, ale použita může být kterákoli databáze, z které CDD čerpá. Pro prohledání databáze a nalezení příslušných domén použije CD-Search program BLAST a výsledky zobrazí v grafické podobě, kde znázorní každou doménu a její pozici na proteinu. Pro každý nalezený výsledek navíc uvede E-value.
2.3 PDB Při studiu proteinů je znalost primární struktury, tedy sekvence aminokyselin, pouze prvním, a nejlehčím krokem. Proteiny se skládají do sekundárních struktur jako je alfa šroubovice nebo beta skládaný list. Několik takových sekundárních struktur dále tvoří tzv. motivy, které se dále skládají do terciální (dále jen „3D“) struktury. Jedna z největších výzev a úkolů pro současné bioinformatiky je nalezení spolehlivé metody pro odvození této 3D struktury proteinů ze znalosti jejich primární struktury. Znalost přesného tvaru proteinu by byla velice užitečná pro pochopení jeho funkce nebo např. pro návrh nových léků. Databáze PDB shromažďuje a je v současnosti hlavním zdrojem informací o 3D struktuře makromolekul (především proteinů a nukleových kyselin). Převážná část těchto dat (asi 86 %) byla získána za použití rentgenové krystalografie, ostatní metody zahrnují především NMR spektroskopii nebo elektronovou mikroskopii. Každá struktura, která je součástí databáze, je v PDB representována čtyřmístným identifikátorem. Pro vytváření těchto identifikátorů platí určitá pravidla. První znak musí být numerický, ostatní tři znaky alfanumerické (např. 3ULA nebo 1K5L). Databáze PDB vznikla v roce 1971 a na začátku obsahovala strukturní informace o sedmi proteinech. Na konci roku 1992 obsahovala již 887 struktur. Od roku 1993 začal počet záznamů v PDB velice rychle narůstat a k datu 18. 5. 2009 obsahuje již 57 558 záznamů. Na konci roku 2009 by počet záznamů měl přesáhnout 60 000 [9]. Všechny tyto záznamy jsou volně přístupné prostřednictvím webových stránek na adrese http://www.rcsb.org/pdb/home/home.do. Druhou možností je stáhnout příslušné záznamy v PDB formátu prostřednictvím FTP serveru, k tomuto úkonu stačí pouze znalost identifikátoru struktury. Kromě samotné primární struktury a informací o pozicích všech atomů příslušného proteinu obsahuje PDB databáze i data o sekundární struktuře, jež je odvozena pomocí algoritmu DSSP. Tento algoritmus je standardním způsobem zjišťování sekundární struktury u proteinů, u nichž známe souřadnice všech atomů a rozlišuje celkem osm typů sekundární struktury. PDB dále poskytuje informace o strukturálních doménách, které získává z databáze SCOP.
4
3. Návrh a implementace programu 3.1 Popis programu Jak už jsem shora uvedl, cílem mnou napsaného programu bylo zjistit hranice exonů příslušného genu nebo genů, tyto hranice převést do sekvence proteinu (případně proteinů), které tento gen kóduje, a následně je také vyhodnotit vůči hranicím konzervovaných domén a sekundární struktuře. Je tedy zřejmé, že pro splnění tohoto cíle budou zapotřebí informace z biologických databází jak o genu, jenž bude předmětem analýzy, tak o tímto genem kódovaném proteinu. Otázkou zůstává, jakým způsobem bude program záznamy obsahující příslušné informace získávat. Nabízejí se dvě základní varianty. První přístup by bylo využití záznamů o genech a proteinech, které se nacházejí na lokálním počítači, ze kterého je program spouštěn. Druhou možností by bylo využít biologických databází a možnosti přístupu k nim přes Internet. Obě řešení mají své výhody i nevýhody. Za zmínku jistě stojí největší výhoda prvního přístupu, a tou je především jeho rychlost. Data na lokálním počítači jsou velice rychle a snadno přístupná, kdykoli je program potřebuje. Tím odpadá závislost na třetí straně. Jak jsem se ale již zmínil v předchozí kapitole, bioinformatické a biologické databáze jsou značně rozsáhlé, obsahují desetitisíce záznamů, které jsou pravidelně rozšiřovány o záznamy nové. Záznamy staré jsou občas upravovány, rušeny nebo jiným způsobem měněny. Tedy vyvstává jednak problém s fyzickou pamětí, jež by byla nutná pro uložení všech záznamů, a další obtíží by se staly pravidelné aktualizace záznamů, které by byly nutné pro dosažení správných výsledků. Takový program, pomineme-li jeho možnou serverovou podobu, by byl v podstatě nepoužitelný pro někoho, kdo by nebyl ochotný zároveň uchovávat spolu s ním na disku svého počítače kompletní databázi, resp. databáze Druhý přístup spoléhá na využití databází volně přístupných v rámci Internetu. Uživatel pro správné užívání programu potřebuje pouze přístup k Internetu, což je v dnešní době považováno téměř za samozřejmost. V této souvislosti mluvíme o nesporné výhodě oproti předchozímu způsobu, ale ani v tomto případě se nejedná o dokonalé řešení. Především si musíme uvědomit, že se program stává pro své správné fungování na databázích závislý a v tomto kontextu existuje řada neovlivnitelných faktorů, jež mohou způsobit komplikace. V první řadě je to fakt, že databázové servery nemusí vždy fungovat, případně mohou být momentálně přetížené, což může v lepším případě způsobit, že si na odpověď od serveru počkáme o chvíli déle a v horším pak, že se jí nedočkáme vůbec. Ani připojení uživatele nemusí být vždy bezchybné a může se stát, že uprostřed analýzy dojde k výpadku. Na zmíněné shora musí program samozřejmě reagovat a vypořádat se s tím, případně alespoň informovat uživatele, kde se stala chyba. I kdyby jsme však předpokládali nulovou chybovost databází a zanedbali možné potíže s připojením k Internetu, stále zůstává největší nevýhoda spojená s využitím databází, a tou je čekání na data od serveru, které i v tom nejlepším případě bude nezanedbatelně dlouhé. Tato doba se především projeví u analýzy velkého množství genů. Zejména proto, že program byl od začátku určen k distribuci mezi uživatele, jsem zvolil řešení druhé, tedy to, v němž jsou využívány databáze, konkrétně se jedná o databáze Entrez a PDB. Uživatel po spuštění programu zadá identifikátor genu nebo případně genů, které chce analyzovat. Program za využití Entrez Programming Utillities přistoupí k databázi Gene a obdrží záznamy příslušného genu (příp. genů) a v těchto záznamech vyhledá identifikátory mRNA a proteinů. Tyto identifikátory obratem použije při přístupu k databázím Nucleotide a Protein. Z těchto záznamů se dají získat
5
informace o sekvencích, ale také hranice exonů, nebo v případě proteinů, konzervovaných domén. Pokud abstrahujeme od samotné analýzy, pak posledním krokem by byl přístup k databázi PDB, kde se program pokusí získat informace o sekundární struktuře proteinu. Tedy za předpokladu, že uživatel zadal smysluplný (jak v rámci dané databáze, tak v rámci programu) identifikátor genu, program musí vykonat nejméně čtyři přístupy k databázím. Později uvidíme, že tento počet je často ještě vyšší. Na přiloženém schématu (obrázek 3.1) je shrnuto, s kterými entitami program během své činnosti přichází do styku.
Obrázek 3.1
Program je schopen pracovat ve třech různých módech. V prvním (standardním) módu uživatel zadává pouze jeden identifikátor genu, na kterém chce provést analýzu. Výstupem jsou textové informace popisující sekvenci proteinu (nebo proteinů) kódovaného genem spolu s hranicemi konzervovaných domén a jeho sekundární strukturou, pokud je známá. Druhý mód je módem statistickým. Zde uživatel zadává cestu k textovému souboru obsahujícímu libovolné množství genových identifikátorů. Výstupem tohoto módu je jednak textový soubor obsahující zprávu o průběhu analýzy a statistické informace pro všechny zpracované geny, ale především také několik grafů, které tyto statistické informace znázorňují graficky. Poslední mód se nazývá módem uživatelským a podobá se módu prvnímu. Tento mód zbavuje uživatele omezení analyzovat hranice exonů pouze proti hranicím domén a sekundární struktuře a poskytuje možnost rozdělit sekvenci proteinu na části podle svého uvážení a hranice těchto částí poté analyzovat proti hranicím exonů. Na vstup je zadána cesta ke speciálně naformátovanému textovému souboru, jenž obsahuje identifikátor genu a hranice nových segmentů v genem kódovaném proteinu. Výstup tohoto módu je obdobný výstupu módu prvního. I když se módy navzájem odlišují, v programu samotném používají buď stejné nebo minimálně se lišící úseky kódu. Práce všech tří módů se pak dá logicky rozdělit na několik podčástí. V první části
6
se získávají záznamy o jednotlivých genech a z těchto záznamů se extrahují informace o mRNA a proteinech, jež jsou s tímto genem spojeny. Druhá část se skládá ze získání záznamů o příslušných mRNA sekvencích a proteinech, z extrakce pro program důležitých informací (jako např. hranice domén, hranice exonů), a z jejich zpracování, především z namapování hranic exonů na sekvenci proteinu. Ve třetí části program získá data o sekundární struktuře proteinu za využití databáze PDB. Tato část je společná pouze pro první dva módy, protože pro poslední mód jsou data o sekundární struktuře zcela irelevantní. Čtvrtá část je částí, kde se vyhodnotí do té doby získaná data, tedy např. dochází k analýze hranic domén proti hranicím exonů. V poslední části jsou vypsány výsledky, resp. vytvořeny grafy, pokud se nacházíme ve statistickém módu. V následujících podkapitolách se budu věnovat bližšímu popisu těchto jednotlivých částí s tím, že občas upozorním na drobné rozdíly mezi jednotlivými módy.
3.2 Práce s genovými záznamy Ihned po spuštění se program uživatele zeptá, ve kterém ze tří různých módů chce pracovat. V případě, že uživatel vybere standardní mód, bude vyzván, aby na standardní vstup vložil identifikátor genu, který chce analyzovat. V opačném případě bude požádán, aby zadal cestu k textovému souboru, ve kterém se nachází buď seznam identifikátorů nebo identifikátor s hranicemi segmentů proteinu, v závislosti na zvoleném módu. Pokud se uživatel rozhodne pro práci v módu statistickém, je požadováno, aby byl textový soubor obsahující genové identifikátory náležitě strukturován. Na jednom řádku se pro správné fungování programu musí vyskytovat právě jeden identifikátor. Tento identifikátor se může nacházet na řádku na libovolné pozici, soubor může dokonce obsahovat prázdné řádky, ale více než jeden identifikátor na jednom řádku pravděpodobně způsobí zkreslení statistických výsledků. Pokud uživatel zvolil práci v uživatelském módu, požaduje se rovněž, aby byl textový soubor speciálně naformátován. První řádek musí začínat znakem „>‟ bezprostředně následovaný identifikátorem genu. Další řádky obsahují jednotlivé segmenty a opět platí pravidlo ne více než jednoho segmentu na jeden řádek. Prázdné řádky jsou povoleny, ale ty neprázdné musí začínat dvěma čísly oddělenými pomlčkou, následovanými jednou mezerou a názvem segmentu. Název segmentu je povinný údaj, nesmí začínat mezerou a musí obsahovat alespoň jeden nebílý znak. V případě nedodržení těchto pravidel program skončí a sdělí uživateli, na kterém řádku souboru se vyskytla chyba. Následující text je příkladem správně naformátovaného souboru:. >651 1-120 Segment 1 124-530 Segment 2 454-651 Segment 3 Už několikrát jsem v předešlém textu použil výraz identifikátor genu, ale ještě jsem tento pojem nedefinoval. V kontextu mého programu se identifikátorem genu myslí identifikátor používaný databází Entrez, což je řetězec numerických znaků, jinými slovy nezáporné celé číslo. Program sám o sobě nijak nekontroluje, zda-li jsou zadané identifikátory v nějakém daném rozmezí, ani jestli obsahují nenumerické znaky. Tuto práci nechává na Etrezu, jenž v případě zaslání dotazu s chybným identifikátorem zašle upozornění.
7
Pro další běh programu předpokládejme, že uživatel úspěšně zadal identifikátor (ať již prostřednictvím standardního vstupu nebo souboru). Program nyní vyšle URL žádost o zaslání záznamu tohoto genu z databáze Gene. K tomu účelu je využita služba efetch, která je součástí Entrez Programming Utillities. Na tomto místě se mohou v principu stát pouze dvě věci. Buď program od databáze obdrží data nebo ne. Druhý případ může znamenat výpadek Internetu, přetížení serveru databáze nebo případně ztracení dotazu nebo odpovědi. Program za těchto okolností vypíše na standardní chybový výstup zprávu pro uživatele a ukončí svou činnost. Po obdržení dat z databáze Gene začne testovací fáze. V první řadě data, která přišla, nemusí být vůbec záznamem z databáze samotné, protože program před zasláním URL dotazu netestuje validitu identifikátoru. V případě neexistujícího identifikátoru genu zašle Entrez upozornění, že nastala chyba. V takovém případě program opět ukončí svou činnost. V opačném případě program obdrží požadovaný záznam, a to ve formátu ASN1. Na tomto místě by bylo dobré si uvědomit, že zdaleka ne všechny geny, které se v databázi Gene vyskytují, jsou pro účely programu použitelné. Především databáze obsahuje celkem deset různých druhů genů. Druh genu se v záznamu nachází pod položkou „gene type“ a možné hodnoty jsou: tRNA, rRNA, snRNA, scRNA, snoRNA, miscRNA, proteincoding (kódující protein), pseudo, other (jiný) a unknown (neznámý). Nebudu zde všechny druhy popisovat, nicméně pro náš záměr jsou vhodné pouze geny s typem protein-coding, které tvoří většinu genů v databázi. Zadruhé existuje možnost, že záznam genu byl z nějakého důvodu skončený („discontinued“) a tedy již není aktuální. Ani tento záznam není pro další analýzu vhodný. Zatřetí se může stát, že příslušný záznam byl nahrazen záznamem jiným. A konečně začtvrté je nutné, aby se v záznamu vyskytoval alespoň jeden pár identifikátorů pro mRNA a protein svázaný s genem (viz dále). Kdyby tyto údaje v záznamu chyběly, znamenalo by to, že buď nejsou známy, nebo že je gen prokaryotický. Všechny zde popsané případy opět způsobí konec programu (spolu s vypsáním zprávy na standardní chybový výstup). V posledních dvou odstavcích jsem popsal všechny možné okolnosti spojené s genovým identifikátorem, které se projeví ukončením programu. Toto je žádoucí chování u standardního a uživatelského módu, kde program analyzuje jediný gen, ale už je méně vhodné u módu statistického. Představme si situaci, kdy má uživatel v úmyslu zanalyzovat geny s identifikátorem 0-10 000, což zabere poměrně netriviální kvantum času. Po úspěšné analýze 9000 genů dojde ke chvilkovému přetížení databázového serveru, následkem čehož program neobdrží žádnou odpověď na svůj URL dotaz. Jistě by za těchto okolností nebylo žádoucí, aby program skončil s chybou bez toho, aniž by poskytl uživateli jakékoli výsledky. Z tohoto důvodu, pokud se program nachází ve statistickém módu a dojde k jednomu z výše zmíněných případů, si program poznamená identifikátor genu (pro pozdější zprávu uživateli), u kterého se vyskytla chyba, a pokračuje analýzou dalšího genu na seznamu. Na tomto místě již vycházíme z toho, že program obdržel záznam genu, který je vhodný pro další zpracování. Poslední činností, kterou program v této části vykonává, je extrakce identifikátorů pro mRNA a pro proteiny, jež gen kóduje. Tyto identifikátory se skládají ze dvou velkých písmen, následovaných podtržítkem a řetězcem číselných znaků. Podle počátečních dvou písmen se dá odvodit druh mRNA nebo proteinu, čehož program využívá. mRNA může mít celkem dva různé druhy identifikátorů. První začíná písmeny NM, druhý písmeny XM. mRNA, jejíž identifikátor začíná písmeny XM je pouze modelová, tedy teoretická mRNA nalezená při analýze sekvence příslušného chromosomu, jejíž existence nebyla potvrzena. Taková mRNA není pro program ničím zajímavá, protože v jejím záznamu nejsou uvedeny žádné informace o exonech. Program se v tuto chvíli také musí vypořádat s dalším fenoménem eukaryotických genů, alternativním splicingem, tedy s faktem, že jeden eukaryotický gen může kódovat více než jeden
8
protein. Děje se tak například z důvodu, kdy jeden (nebo více exonů) může být vynechán při transkripci, což se následně projeví i při translaci a na výsledném proteinu. Ve standardním a uživatelském módu program analyzuje a později podá zprávu o všech transkriptech a proteinech, jejichž identifikátory se nachází v záznamu genu. Zde se opět drobně odlišuje mód statistický. Statistický mód by mohl taktéž analyzovat všechny produkty alternativního splicingu, ale problémem je existence genů (např. gen s identifikátorem 51), jejichž alternativní traskripty se liší pouze v částech, které neovlivní výsledný protein. Tedy jeden gen může kódovat několik různých proteinů, které jsou sice v databázi uloženy pod jinými identifikátory, ale jinak jsou naprosto identické. Program by ve statistickém módu analyzoval jeden a ten samý protein několikrát, což by ve svém důsledku znamenalo nežádané ovlivnění konečných statistik. Program tedy musí ke každému genu vybrat pouze jeden transkript a k němu náležící protein, na kterém provede analýzu. Do úvahy by přicházelo vybrat protein, který je nejdelší nebo naopak nejkratší, ale tento přístup se ukazuje jako nepříliš vhodný. Důvodem je, že v záznamu genu se naneštěstí nenacházejí informace o délce proteinů, pouze jejich identifikátory. Pro zjištění délky proteinu je nutné vyžádat si jeho záznam od databáze, což představuje přílišné zdržení. Proto jsem se při psaní programu rozhodl, že statistický mód bude tento výběr provádět čistě náhodně. Samotná extrakce informací o mRNA a proteinech ze záznamu genu je založená na klíčové vlastnosti Entrezu, tedy na tom, že Entrez hledá a eviduje logické vztahy mezi záznamy z jednotlivých databází. Je to ovšem dvousečná zbraň. Díky této vlastnosti záznamy o genech obsahují identifikátory jejich transkriptů a proteinů, ale také mimo jiné obsahují odkazy na velké množství dalších záznamů, které mají s genem souvislost, ale pro program jsou bezvýznamné. Proto je nutné si dát při extrakci dat ze záznamů velký pozor a dbát na pečlivost při zadávání parametrů pro jejich vyhledávání.
3.3 Práce se záznamy mRNA a proteinů V předchozí podkapitole jsem se pokusil vysvětlit, jakým způsobem program získává záznamy o genech a jak je využívá k vyhledání identifikátorů mRNA a proteinů. V této podkapitole shrnu, jakým způsobem jsou tyto identifikátory využívány v dalších částech programu. Prvním krokem je získání databázových záznamů o mRNA a proteinech. Princip je stejný jako u genů. Opět je využita služba efetch, jenom v tomto případě se odkazujeme na databáze Nucleotide a Protein. Odlišuje se i formát, ve kterém si program nechává záznamy posílat. V tomto případě je jisté, že použité identifikátory jsou správné, takže programu odpadá povinnost testování, zda došlá data opravdu představují databázový záznam, nicméně program stále musí testovat, jestli se nevyskytne chyba zapříčiněná připojením nebo přetížením serveru. Reakce na takovou chybu je opět stejná jako v případě, kdy program získával genový záznam (včetně odlišné reakce pokud se program nachází ve statistickém módu). Když program posílá URL dotazy databázi Nucleotide nebo Protein, požaduje, aby mu byly poslány záznamy ve formátu GenBank, což je v podstatě pouze strukturovaný ASCII text, jenž je velice dobře čitelný pro lidi (narozdíl například od formátu ASN1). NCBI tento formát používá pro ukládání záznamů pro databázi nesoucí stejné jméno jako zmíněný formát. Záznam ve formátu GenBank je rozdělen na několik částí. První část se označuje slovem LOCUS a obsahuje informace, jako je jméno proteinu (příp. mRNA nebo dalších), jemuž je záznam věnován, délku sekvence nebo čas poslední modifikace záznamu. Další části jsou například SOURCE s informacemi o organismu, ve kterém byl protein lokalizován, REFERENCE, kde jsou odkazy na články a publikace od autorů záznamu vztahující se k danému proteinu, nebo KEYWORDS obsahující klíčová slova. Poslední dvě
9
části jsou pro program obzvláště důležité. Část, jež nese jméno FEATURES, obsahuje informace o biologicky důležitých regionech, jako jsou vazná místa (binding sites) nebo pro nás velmi důležité exony (v případě mRNA). Poslední část má označení ORIGIN a je v ní zaznamenána sekvence, tedy řetězec aminokyselin nebo nukleotidů [10]. Pro účely programu jsou podstatné především části LOCUS, FEATURES a ORIGIN. Z části LOCUS program vyextrahuje a uloží označení mRNA nebo proteinu pro pozdější použití při vypsání výsledků analýzy. V případě mRNA záznamu obsahuje část FEATURES informace o hranicích exonů a pole se jménem CDS (kódující sekvence), což je souvislý region nukleotidů v mRNA, který odpovídá sekvenci aminokyselin v proteinu. Je důležité, že CDS nemusí a taky většinou nekončí/nezačíná tam, kde začíná/končí mRNA. Z toho mimo jiné vyplývá, že ne všechny exony, ze kterých se mRNA skládá, se přímo účastní translace (viz dále). Pokud jde o záznamy proteinů, je velice příjemné, že se v části FEATURES dají nalézt hranice a jména konzervovaných domén, které se v proteinu nachází, což je další příklad logického propojení databází zajišťované Entrezem. Kdyby tyto informace v záznamech proteinů nebyly, musel by program využít službu CD-Search (viz kapitola 2.2). Je jenom škoda, že se v záznamech proteinů nenachází i informace o sekundární struktuře.1 Ty se musí získávat poněkud obtížnějším způsobem. Z části ORIGIN si program uloží sekvenci mRNA/proteinu. Hranice domén jsou důležité pouze u prvních dvou módů. Mód uživatelský se jimi nezabývá, protože jeho cílem je výhradně analýza exonů vůči uživatelem definovaným segmentům proteinu. Po získání všech důležitých dat ze záznamů čeká program v této části poslední netriviální problém. V současné chvíli má program uložené hranice exonů tak, jak se nacházejí v mRNA. Tyto hranice musí být převedeny na sekvenci proteinu, tedy aby začátek a konec exonu tvořila čísla reprezentující pozice aminokyselin, na kterých exon začíná a končí. K tomu bude program potřebovat kromě starých hranic také informace z pole CDS. První komplikací, na kterou během převádění hranic program narazí, je fakt, že ne všechny exony v mRNA musí být kódující. Ve skutečnosti lze exony rozdělit až do dvanácti různých kategorií v závislosti na jejich transkripčních a translačních hranicích [11]. Jádrem celého problému je fakt, že se mRNA neskládá pouze z CDS, ale obsahuje i další části (viz. obrázek 3.2).
Obrázek 3.2
5‟ Cap je pouze modifikovaný guanine, Poly-A tail (ocásek) je dlouhá sekvence nukleotidů adenin. Popis jejich funkcí přesahuje rámec této práce, ale důležité je, že spolu s UTR (nepřekládané regiony) nemají žádný přímý vliv na aminokyselinové složení proteinu, který po translaci vznikne. Důsledkem, se kterým se musíme vypořádat je, že některé exony mají levou hranici v 5‟ UTR regionu, 1
Pro upřesnění: Informace o sekundární struktuře se v záznamech proteinů nacházejí, ale jedná se pouze o proteiny převzaté z databází PDB nebo UniProt/Swiss-Prot. S těmito záznamy program během své činnosti do styku nepřijde.
10
kdežto pravá hranice je součástí CDS. Podobně tomu je s 3‟ UTR regionem. Některé exony nemusí do CDS vůbec zasahovat. Otázka tedy zní: Které exony budou zahrnuty do analýzy? Odpověď asi není moc překvapivá. Analyzovány budou všechny exony, které alespoň jednou aminokyselinou přispívají do sekvence proteinu, což ve svém důsledku znamená, že po převedení hranic exonů, bude mít levá hranice prvního exonu s velkou pravděpodobností zápornou hodnotu a naopak pravá hranice posledního exonu bude přesahovat za konec proteinu. Exony, které se na tvorbě proteinu nepodílejí, program zapomene a nebude je uvádět ani v závěrečné zprávě. Další, tou závažnější komplikací je, že začátky a konce exonů nekorespondují se začátky a konci kodónů, což následně způsobuje nejednoznačnosti při mapování hranic exonů na protein. Kodóny jsou trojice nukleotidů, které určují složení proteinů, i když přesněji řečeno to, co určuje složení proteinů je spíše genetický kód, který mapuje jednotlivé kodóny na aminokyseliny. V ideálním případě by exony začínaly na první pozici v kodónu a končily na pozici poslední. Tomu tak ale vždy není. Konec jednoho exonu a začátek následujícího exonu často sdílejí jednu aminokyselinu a program musí rozhodnout o tom, který z nich má na tuto aminokyselinu nárok. Jedno z řešení, které se nabízí, je povolit, aby v již převedených hranicích končil jeden exon na té samé aminokyselině, na které začíná exon další, což by ovšem mělo velmi špatný vliv při pozdějším vyhodnocení hranic, kde se samostatně analyzují začátky a konce exonů. Jiný přístup je preferování buď začátků nebo konců hranic, ale i toto řešení by způsobilo pozdější zkreslení statistik. Proto jsem se rozhodl pro řešení, které je vůči začátků, a koncům exonů spravedlivé a zároveň nejméně ovlivňující statistické vyhodnocení. Využil jsem toho, že kodón je tvořen lichým počtem nukleotidů, proto aminokyselina bude přiřčena tomu exonu, který v jejím kodónu dosáhne majority. Jinými slovy exon potřebuje v kodónu nejméně dva nukleotidy ze tří.
3.4 Získání sekundární struktury Zmínil jsem se již o tom, že v záznamech proteinů, které program stahuje přes rozhraní Entrez, scházejí informace o sekundární struktuře. Proto, pokud bychom chtěli dělat analýzu hranic exonů nejen vůči doménám, musíme nalézt nějaký alternativní zdroj informací, ze kterého by se dala data o sekundární struktuře proteinů získat. Mluvím samozřejmě o jiných databázích proteinů, které kromě sekvence obsahují i informace ohledně sekundární struktury. Mezi nejvýznamnější databáze tohoto druhu patří např. PDB nebo Swiss-Prot. Já jsem se rozhodl pro využití databáze PDB, a to ze dvou důvodů. Především databáze Swiss-Prot obsahuje kromě proteinů se známou sekundární strukturou i spoustu proteinů, u nichž sekundární struktura známá není, což by jistě ztížilo hledání.2 Druhý důvod se spíše týká jednoduchosti výsledného kódu a usnadnění si práce, protože pro práci s PDB databází existuje jedna velice šikovná knihovna, která většinu práce udělá za nás. Tato knihovna je volně k dispozici na adrese http://search.cpan.org/~miorel/WWW-PDB-0.00_03/lib/WWW/PDB.pm. Třetím důvodem by mohla být služba, kterou PDB poskytuje, která pro libovolný řetězec aminokyselin ve formátu FASTA prohledá a provede zarovnání se všemi záznamy proteinů z databáze programem BLAST a která nalezne všechny podobné sekvence. Tu stejnou službu ovšem nabízí i databáze SwissProt. Knihovna WWW::PDB disponuje velkým množstvím užitečných funkcí pro práci s databází PDB. Poskytuje přístup k informacím jako je sekvence, délka sekvence, status proteinu, sekundární 2
I databáze PDB obsahuje proteiny, u kterých chybí informace o sekundárních struktuře, ale ve srovnání s databází Swiss-Prot to je zanedbatelné množství.
11
struktura apod. Vše, co je potřeba, je znalost čtyřmístného identifikátoru proteinu, o který se zajímáme. Kromě již vyjmenovaných funkcí může uživatel knihovny také použít funkci se jménem blast, která zprostředkovává vzdálený přístup k výše zmíněné PDB službě využívající program BLAST. Program pro získání sekundární struktury postupuje následovně. Na začátku přidělí celému proteinu sekundární strukturu: neznámá. Existuje totiž nemalá šance, že se buď nepovede odvodit sekundární strukturu, nebo že sekundární strukturu odvodíme jen pro určité segmenty z celkové sekvence proteinu. Program se dále snaží nalézt PDB proteiny s podobnou sekvencí aminokyselin. Pokud budeme mít štěstí, bude programem analyzovaný protein zcela odpovídat nějakému proteinu z PDB. To se ovšem příliš často nestává. Mnohem pravděpodobnější je, že v PDB existuje několik proteinů, které mají se zkoumaným proteinem společný segment několika desítek nebo stovek aminokyselin. V takovém případě tomuto segmentu přiřadíme zjištěnou sekundární strukturu a u všech ostatních aminokyselin proteinu ponecháme označení: neznámá. K hledání PDB proteinů, které by mohly pomoci k odvození sekundární struktury, využívá program již zmíněnou funkci knihovny WWW::PDB blast. Tato funkce (alespoň ve formě, v jaké ji využívá program) přijímá 4 parametry. Prvním parametrem je samotná sekvence, pro kterou chceme nalézt podobné proteiny. Druhý a třetí parametr se týkají nastavení BLASTu, resp. filtrování potenciálních výsledků. Je to minimální E-hodnota, kterou musí výsledná zarovnání splňovat. E-hodnota je jednoduše řečeno číslo, charakterizující jedno výsledné zarovnání, určující pravděpodobnost, že existuje zarovnání lepší než to nalezené. Tedy čím menší hodnotu nastavíme, tím méně výsledků bychom měli dostat. Třetím parametrem je matice, která bude použitá programem BLAST při určování skóre při zarovnávání. Tyto matice jsou konstruovány tak, aby při zarovnávání braly do úvahy pravděpodobné mutace. Nejtypičtějšími zástupci jsou matice PAM a BLOSUM. Poslední parametr funkce blast je formát, ve kterém budou vráceny výsledky. Můj program používá při volání funkce následující parametry: WWW::PDB->blast($sequence, 1e-17, 'BLOSUM80', 'HTML') Funkce blast typicky vrátí několik málo výsledků seřazených podle jejich E-hodnoty. Každý výsledek representuje jedno zarovnání a obsahuje krom jiných následující informace: identifikátor PDB proteinu, se kterým se provádělo zarovnání, počet identických aminokyselin v zarovnaných segmentech, počet positivních shod a počet mezer obsažených v zarovnání. Pod těmito údaji se nacházejí oba zarovnané segmenty včetně číselných hranic. Program z dosažených zarovnání vybírá pouze ta, která mají počet identických aminokyselin větší nebo roven 75 % a která neobsahují žádné mezery. Bohužel ani v této části se neobejdeme bez problémů. Program sekvenčně prochází jednotlivá zarovnání podle jejich dosažené E-hodnoty. Pokud některé zarovnání splňuje kladené požadavky, program si poznamená hranice obou zarovnaných segmentů pro pozdější přiřazení sekundární struktury. Co ale dělat, pokud se vyskytnou dvě nebo více různých zarovnání, pro která platí, že zarovnané segmenty programem analyzovaného proteinu mají nenulový průnik? Mělo by se upřednostňovat delší zarovnání nebo to s lepším skóre? Měly by se vůbec brát do úvahy zarovnání, která mají nenulový průnik s již přijatými zarovnáními? To jsou otázky, kterými se při odvozování sekundární struktury musíme zabývat. Řešení nastíněného problému existuje několik. Nejjednodušší by asi bylo ignorovat jakékoli zarovnání, které má nenulový průnik s nějakým již akceptovaným zarovnáním. Toto řešení není ve své podstatě špatné, protože statisticky vzato se většina nejlepších zarovnání týká stejného
12
segmentu proteinu, pouze se liší v několika málo aminokyselinách. Může ovšem nastat i případ, kdy zarovnání bude mít s nějakým již přijatým průnik o délce jedné aminokyseliny. Ignorováním tohoto zarovnání bychom ztratili velké množství potenciálně jinak nenahraditelných informací. Druhým extrémem by bylo u každého zarovnání splňující požadavky pouze useknout průnik a zbytek přidat k již akceptovaným zarovnáním. Toto řešení trpí jiným nedostatkem. Tím je případ, kdy výsledkem useknutí průniku od zarovnání bude jedna nebo obecně několik málo aminokyselin. Toto ve skutečnosti nastává velice často. Já jsem zvolil střední cestu, kdy program postupně zkoumá jednotlivá zarovnání, a pokud nalezne PDB protein, jehož zarovnání neobsahuje žádnou část obsaženou v zarovnání jiném, pak je jisté, že zarovnaný segment v tu chvíli analyzovaného proteinu dostane přiřazenu sekundární strukturu právě tohoto PDB proteinu a nic to už nemůže změnit. Pokud se objeví zarovnání, které má s nějakým jiným již přijatým nenulový průnik, tak se spočítá velikost tohoto průniku (myšleno průnik se všemi již přijatými zarovnáními). Pokud je tento průnik menší než 50 % celkové délky zarovnání, je toto zarovnání přijato, ale části náležící do průniku budou useknuty. Příklad můžeme vidět na obrázku 3.3. Zarovnání 1 a 2 měla postupně nejlepší E-hodnoty a byla přijata. Zarovnání 3 má s předchozími dvěma nenulový průnik, ale ten není již od pouhého oka větší než 50 % jeho délky, takže je po nezbytných úpravách následně také přijato.
Obrázek 3.3
Dalším krokem programu je získání sekundární struktury pro každý segment, jenž byl součástí nějakého zarovnání. K tomuto poslouží funkce WWW::PDB->get_kabsch_sander. Této funkci stačí pouze předat jako parametr PDK identifikátor proteinu a ona vrátí řetězec s informacemi o sekundární struktuře. Z tohoto řetězce pak stačí vyjmout příslušnou část a přiřadit ji segmentu, kterému patří. Poslední akcí, kterou program provádí, je úprava sekundární struktury. Pro připomenutí PDB proteiny mají sekundární strukturu odvozenou algoritmem DSSP. Tento algoritmus přiřazuje jednotlivým aminokyselinám osm různých typů sekundární struktury (s přiřazenými jednopísmenovými zkratkami G, H, I, E, B, T, S, L). Program těchto osm typů převede na mnohem častější členění. Místo G, H a I je aminokyselinám přiřazena alfa šroubovice (alfa helix), E a B jsou nahrazeny beta skládaným listem (beta sheet) a ostatní jsou nahrazeny smyčkou (loop). Ještě by se asi hodilo připomenout, že ne každá aminokyselina musí mít přiřazenou sekundární strukturu. V každém
13
proteinu existují krátké segmenty, u kterých je sekundární struktura nepřiřazená (nezaměňovat s neznámou sekundární strukturou). Tato část programu je v uživatelském módu vynechána, a to ze stejných důvodů proč byla vynechána analýza domén.
3.5 Vyhodnocení dat Shrňme si, jaká data program do této doby získal (alespoň, pokud vše probíhalo bez komplikací). Z druhé části program má informace o hranicích exonů, které převedl, aby odpovídaly jejich umístění na proteinu. Ze stejné části pocházejí i hranice domén. Poslední část se zabývala odvozením sekundární struktury. Pokud se sekundární struktura aminokyseliny proteinu odvodit nedá, bude označena za neznámou. Program tedy má všechna potřebná data pro vyhodnocení hranic exonů. Celkově budou vykonány dvě různé analýzy, proti sekundární struktuře proteinu a proti hranicím domén. Obě tyto analýzy budou mít rozdílný průběh i rozdílný výsledek. Je tomu tak proto, že se dá předpokládat, že domén je v proteinu téměř vždy méně než exonů, což znamená, že zdaleka ne všechny exony (resp. jejich hranice) budou mít ve své blízkosti nějakou konzervovanou doménu. Jinak je tomu u sekundární struktury. Zde se dá naopak předpokládat, že mezi každým začátkem a koncem exonu se vyskytuje alespoň jedna aminokyselina s přiřazenou sekundární strukturou (do této množiny je započítána i neznámá sekundární struktura). Z předešlého textu je asi jasné, že v uživatelském módu se hranice exonů vůči sekundární struktuře analyzovat nebudou. To samé by se dalo říci o doménách, nicméně analýza segmentů proti hranicím exonů, což je jediná analýza, která se bude v uživatelském módu provádět, má naprosto stejný průběh jako analýza proti hranicím domén. Z tohoto důvodu se popisem analýzy segmentů v dalším textu nebudu zabývat. Nejprve popíši analýzu hranic exonů vůči sekundární struktuře. Důležité je, že program zachází jinak s začátky a jinak s konci exonů. Analyzuje je samostatně, což se projeví především na závěrečných grafech statistického módu. Analýza sekundární struktury je velmi jednoduchá. Pro každou hranici exonu program kontroluje sekvenci v obou směrech a hledá nejbližší sekundární strukturu. Pokud ji nalezne, zapamatuje si typ sekundární struktury a vzdálenost od hranice. Tato vzdálenost může být i záporné číslo, aby se v závěrečné zprávě nebo grafu odlišily jednotlivé směry na sekvenci. V případě, že je nalezena sekundární struktura v obou směrech od hranice a ve stejné vzdálenosti, program si zapamatuje a do pozdějšího vyhodnocení zahrne obě dvě. Aminokyseliny, které mají nepřiřazenou sekundární strukturu, jsou ignorovány. Jinak tomu ale je u aminokyselin s neznámou sekundární strukturou. Tyto přeskakovány nejsou, právě naopak, narazí-li program na neznámou sekundární strukturu, dál již nehledá a přechází k další hranici. Kdyby tomu tak nebylo, došlo by k výraznému zkreslení konečných statistik ve statistickém módu, které se týkají vzdálenosti sekundární struktury od hranic exonů. Zde popsaný proces je společný všem hranicím exonů kromě začátku prvního exonu a konce posledního exonu. U těchto se nejprve musí testovat, zda leží mimo hranice proteinu, a pokud ano (což je velice pravděpodobné), pak jsou z této analýzy vynechány. U analýzy hranic exonů vůči hranicím domén je situace o poznání složitější. Stejně jako u sekundárních struktur hledáme v obou směrech od hranice, ale v tomto případě, pokud je hranicí začátek exonu, hledáme pouze začátky domén, a stejně u konců exonů hledáme pouze konce domén. Principiálně bychom mohli postupovat obdobně jako u předešlé analýzy, musíme mít ale na paměti pár věcí. Především u sekundárních struktur jsme vycházeli z předpokladu, že předtím, než v našem pátrání po sekvenci dosáhneme další nebo předešlé hranice exonu, nalezneme
14
buď aminokyselinu s přiřazenou sekundární strukturou, nebo aminokyselinu, jejíž sekundární struktura bude neznámá. Podobný předpoklad si na tomto místě dovolit nemůžeme, a to kvůli důvodům již zmíněným na začátku podkapitoly. Je rozhodně nežádoucí, už jenom kvůli možným důsledkům ve statistickém módu, aby program pro každou hranici prohledával sekvenci až do té doby, dokud nenalezne začátek/konec domény. Kdyby totiž protein obsahoval jenom jednu doménu, byl by její začátek přiřazen všem začátkům exonů a stejně tak i její konec. Naším cílem by mělo být, aby každý začátek/konec domény byl přiřazen nejvýše jedné hranici exonu. Pro každou hranici exonu musíme tedy nějakým způsobem omezit prostor, ve kterém bude program pátrat po hranicích domén. Omezení tohoto prostoru pevně zvolenou hodnotou je zcela jistě nejsnadnějším řešením, ale opět zdaleka ne ideálním. Mnohem elegantnější je řešení, ve kterém při hledání hranic domén nesmí být překročena hranice (myšleno hranice exonu) stejného druhu, jinými slovy při hledání začátků domén pro začátek exonu nesmíme překročit začátek předcházejícího a začátek dalšího exonu. Musíme samozřejmě dát pozor na odlišné chování prvního a posledního exonu. To ale samo o sobě nestačí. Pořád existuje možnost, že jeden začátek domény bude přiřazen právě dvěma začátkům exonů, ale to se dá řešit. Stačí, když přidělíme začátek/konec domény té bližší ze dvou hranic. Ani teď ale nemůžeme být zcela spokojeni. Uvažme situaci na obrázku 3.4.
Obrázek 3.4
Protein je rozdělen šesti exony a obsahuje dvě domény. Hranice těchto domén jsou přiřazeny hranicím exonů, jak je znázorněno na obrázku. Začátek domény B je tak přiřazen začátku exonu č. 4, a to i přesto, že vzdálenost k začátku exonu č. 5 je několikanásobně kratší. Začátek exonu č. 5 už ale má svou doménu přiřazenou. To se může zdát jako zcela korektní, ale jenom do té doby, než začneme uvažovat nad důsledky pro statistický mód. Jedna ze statistik tohoto módu se zabývá vzdáleností začátků/konců domén od začátků/konců exonů a případy tohoto typu by se na ní mohly neblaze podepsat. Řešení se nachází v dalším omezení prostoru, kde program pátrá po doménách, a to přesně o polovinu. Na počátku si každý začátek/konec exonu zjistí vzdálenosti k začátku/konci předchozího a následujícího exonu a v následném hledání domén nepřekročí jejich poloviny v jednom nebo druhém směru. Takže například v předchozím příkladě začátek domény B nebude přiřazen žádnému začátku exonu (stejně jako její konec). Nakonec bych ještě rád upozornil, že na rozdíl od analýzy sekundární struktury se v analýze domén pracuje i s těmi hranicemi exonů, které zasahují mimo protein. Ve statistickém módu program vykonává ještě jednu jednoduchou analýzu a tou je zjištění, jaký druh aminokyseliny se nachází na pozici, kde má exon své hranice. V ostatních módech
15
by tyto informace neměly smysl z důvodu prezentace jejich výsledků, ze které se tyto údaje dají velice jednoduše zjistit.
3.6 Statistika a prezentace výsledků V tento moment má program všechna potřebná data, stačí je pouze vhodně prezentovat a v případě druhého módu statisticky vyhodnotit. Standardní mód a uživatelský mód prezentují výsledky své činnosti textovou zprávou, kterou pošlou na standardní výstup. Uživatel tuto zprávu může přesměrovat do svého souboru, pokud mu tento způsob nevyhovuje. Standardní mód vznese ještě před začátkem analýzy dotaz na složku, ve které by měly být výsledky uloženy. Pokud tato složka neexistuje, je vytvořena. Výsledné zprávy standardního a uživatelského módu jsou velice podobně formátovány. V hlavičce zprávy je napsáno jaký gen byl zpracováván a kolik různých transkriptů bylo z databáze pro daný gen získáno a následně zanalyzováno. Podle počtu transkriptů je zbytek zprávy rozdělen na části, z nichž každá obsahuje podrobnou analýzu daného transkriptu a k němu patřícímu proteinu. V záhlaví jsou vždy uvedeny základní informace jako název mRNA a proteinu a jejich ID, pod kterými je možné je vyhledat přes Entrez. Součástí názvu bývá jméno organismu, ze kterého mRNA a protein pocházejí. V následující části je vypsána sekvence proteinu včetně sekundární struktury ke každé aminokyselině. Tato část je speciálně strukturována pro dobrou orientaci v sekvenci. Na jednom řádku je vypsaných maximálně šedesát aminokyselin včetně pozic, kde se příslušná část v proteinu nachází. Tyto informace jsou následovány výčtem exonů a konzervovaných domén, včetně jejich hranic v proteinu. V případě uživatelského módu jsou namísto domén vypsány uživatelem definované segmenty. Na konci jsou vypsány výsledky analýz, a to samostatně pro začátky a konce exonů. Součástí příloh jsou textové dokumenty obsahující příklady zpráv obou módů. Výsledky práce statistického módu se skládají z jednoho textového dokumentu pojmenovaného „report”, jenž obsahuje sumární informace o analýze a osm sloupcových grafů. Samotné grafy jsou vytvořeny prostřednictvím funkcí v knihovně GD::Graph. V dokumentu report jsou postupně uvedeny informace jako počet analyzovaných genů, kolik tyto geny obsahovaly celkově exonů, kolik genů nemohlo být zanalyzováno včetně jejich identifikátorů, kolika exonům se povedlo přiřadit sekundární strukturu/doménu nebo z jakých aminokyselin se skládaly proteiny (vyjádřeno procenty). Tato zpráva obsahuje také procentuálně vyjádřené informace, které jsou součástí grafů, jako je např. typ sekundární struktury na začátku nebo konci exonu nebo druh aminokyseliny v místě, kde začínají/končí exony. Poslední významné údaje obsažené ve zprávě jsou výčty názvů domén a průměrné vzdálenosti jejich hranic od hranic exonů včetně směrodatné odchylky. Zmíněné grafy obsahují informace o čtyřech různých statistikách (vždy jeden graf pro začátky exonů a jeden pro konce exonů). První statistika se týká vzdálenosti hranic exonů a k nim nejbližších hranic domén. Tyto vzdálenosti jsou v grafu vyjádřeny pomocí jedenácti intervalů, každému intervalu odpovídá jeden sloupec v grafu. Druhou statistikou je vzdálenost nejbližší sekundární struktury k hranicím exonů. Zde mohou vzdálenosti spadnout do sedmi různých intervalů. Poslední dvě statistiky jsou druh aminokyselin na pozici hranic exonů a druh nejbližší sekundární struktury. Grafy ve skutečnosti obsahují sloupce dvou různých barev. Červené sloupce representují již zmíněné statistiky, které byly spočítány pro zadané geny a reálné exony v těchto genech zaznamenané. Zelené sloupce obsahují statistiky pro stejnou množinu genů, ale náhodné hranice exonů. Program ve statistickém módu tedy dělá všechny analýzy a statistiky dvakrát. Jednou pro reálné hranice exonů a jednou pro hranice, které byly náhodně zvoleny. Celkový počet exonů
16
v genech se ale nikdy nemění. Tyto dodatečné informace jsou uvedeny pro porovnání a mohou napomoci lepšímu vyhodnocení důležitosti výsledků. Pro představu zde prezentuji dva příklady grafů (viz obrázek 3.5 a 3.6), které byly získány po analýze genů s identifikátory 0-100. Příklady všech typů grafů jsou obsaženy v přílohách.
Obrázek 3.5
17
Obrázek 3.6
3.7 Implementace Program byl vyvíjen v jazyku Perl a pro svou práci využívá knihovny LWP, GD::Graph, WWW-PDB a POSIX. Testování probíhalo pod systémem Windows Vista s dual-procesorem Intel, 2.00 GHz. Všechny proměnné, komentáře, výsledné zprávy a statistiky programu jsou v anglickém jazyce z důvodu snadnější rozšiřitelnosti. Samotný program dostal název „exon examiner“. Myslím, že i bez komplexního testování by mělo být na první pohled patrné, kde se nachází největší problém celého programu. Je jím statistický mód, resp. jeho nepříliš uspokojivá rychlost. Pro získání opravdu směrodatných a vypovídajících statistik je nutné zpracovat velké množství genů, což znamená spoustu přístupů k internetovým databázím a velké množství spotřebovaného času. Jedním ze způsobů, jak by se tomu dalo předejít, je kód paralelizovat a zpracovávat tímto způsobem několik genů současně. Bohužel to není v praxi možné, protože NCBI omezuje počet dotazů na Entrez od jednoho uživatele na 3 za 1 sekundu.
18
4. Dosažené výsledky V této chvíli máme k dispozici plně funkční nástroj, který by nám při správném použití měl vyhledat informace o exonech, doménách a sekundární struktuře a tato data vyhodnotit. Je jasné, že pro získání dostatečně relevantních a vypovídajících statistik je nutné otestovat velké množství genů. Na tomto místě bych rád uvedl výsledky testu, kdy program dostal za úkol zanalyzovat všechny geny databáze Entrez s identifikátorem 1 – 10 000 (výsledná zpráva a grafy jsou součástí doplňkových materiálů). Z výsledné zprávy vyplývá, že z celkového objemu 10 000 genů, pouze 6597 bylo skutečně předmětem analýzy a přispěly tak k výsledným hodnotám. 3403 genů otestováno nebylo. Mohlo se tak stát z nejrůznějších důvodů, z nichž ty nejpodstatnější byly popsány ve třetí kapitole. Je jen velice těžké odhadnout, kolik genů nemohlo být zpracováno z důvodu chyb nebo přetížení databází, ale je velice pravděpodobné, že by pro stejná vstupní data program při opakovaném běhu vyprodukoval trochu odlišné výsledky. Tyto rozdíly by ale neměly být pro výsledné statistiky nijak podstatné. Program našel a zpracoval 71 670 exonů, z tohoto počtu pouze 18 % exonů dostalo přidělenou doménu (hodnoty se nepatrně liší pro začátky a konce exonů), což potvrzuje, že konzervovaných domén se v proteinu obyčejně nachází méně než exonů. Typ sekundární struktury v místech, kde má exon své hranice, byl zjištěn u přibližně 19,5 % exonů. Nyní se podrobněji zaměřme na jednotlivé grafy, začneme se sekundární strukturou. Pokud jde o vzdálenost nejbližší aminokyseliny se známým typem sekundární struktury od hranice exonu (ať už začátku nebo konce), dopadlo vše podle předpokladů. Ve většině případů je typ sekundární struktury známý již přímo pro aminokyselinu na hranici nebo pro aminokyselinu 1-2 pozice vzdálenou v obou směrech od hranice. Mnohem zajímavější data poskytl druhý typ grafu. Z jeho výsledků, které jsou opět téměř totožné pro začátky i konce exonů, vyplývá, že nejčastějším typem sekundární struktury nejblíže k místu hranice je smyčka. V těsném závěsu je alfa šroubovice a nejméně zastoupený typ aminokyseliny je beta složený list (viz tabulka 4.1).
Typ sekundární struktury
Typ sekundární struktury nejblíže k začátku exonu
Typ sekundární struktury nejblíže ke konci exonu
Celkové zastoupení v analyzovaných proteinech
Alfa šroubovice
35,01 %
36,08 %
45,38 %
Beta složený list
26,81 %
27,06 %
27,41 %
Smyčka
38,18 %
36,86 %
27,21 %
Tabulka 4.1
Údaje o nejbližším typu sekundární struktury se stanou zajímavé v okamžiku, kdy je porovnáme s procentuálním zastoupením jednotlivých typů v proteinech, které byly v průběhu programu podrobeny analýze. Čistě statisticky by tyto hodnoty měly být velmi podobné, ne-li téměř totožné s hodnotami získanými z grafů. Jak je tedy možné, že smyčka dosáhla tak vysokých hodnot? Jeden z možných důvodů, který se na první pohled nabízí, je, že typ sekundární struktury smyčka i
19
přes své nižší zastoupení v celkové sekundární struktuře proteinu prokazuje vysokou korelaci s hranicemi exonů. Existuje však i jiné vysvětlení. Pokud bychom zkoumali sekundární strukturu jednotlivých proteinů, brzy bychom si všimli, že typy alfa šroubovice a beta složený list většinou pokrývají delší, spojité segmenty proteinů. Naopak typ smyčka jen velmi vzácně pokrývá souvislé úseky delší než pět aminokyselin. Když se zamyslíme nad důsledky, vyplyne nám, že pokud začátek/konec exonu nespadne do úseku proteinu se sekundární strukturou alfa šroubovice nebo beta složený list a my musíme nejbližší typ sekundární struktury hledat dále od místa hranice, je větší šance, že jako první narazíme právě na smyčku. Je velice těžké určit, jakou měrou je výsledek právě tímto ovlivněn. Nápovědou by mohla být část grafu, která se vztahuje k náhodným hranicím exonů. I zde totiž typ smyčka dosáhl poměrně vysokých hodnot, což jenom podporuje správnost druhého vysvětlení. Nyní rozeberme grafy znázorňující nejčastěji zastoupené aminokyseliny v místě hranice. V místě, kde exon začíná dosahuje největších hodnot glycin (18.76 %) a valin (12.84 %). Pokud tyto hodnoty porovnáme s hodnotami pro zastoupení jednotlivých aminokyselin v celých sekvencích proteinů, nemůžeme si nevšimnou výrazných rozdílů (6.85 % pro glycin a 6.08 % pro valin). Stejně tomu je i na místě, kde exony končí. Zde mají největší zastoupení lysin (17.14 %), glutamin (13.49 %), arginin (11.51 %) a kyselina glutamová (10.26 %) a stejně jako v předešlém případě jsou tyto hodnoty dvakrát až třikrát vyšší než by statisticky měly být. I když se na první pohled dosažená čísla jeví zajímavě, jsou pouhým důsledkem nerovnoměrného rozložení nukleotidů v místech hranic mezi exony a introny (toto rozložení je dobře graficky znázorněné na stránce http://weblogo.berkeley.edu/examples.html). Nakonec se věnujme analýze hranic domén vůči hranicím exonů, která byla hlavní motivací pro vznik celého programu. Z grafů je patrné, že hranice domén se poměrně zřídka nachází dále než dvacet pět aminokyselin od místa hranice, což by naznačovalo jistou korelaci. Abychom však tato data mohli zasadit do správného kontextu a správně je vyhodnotit, potřebujeme ještě znát průměrnou délku exonů, což v našem případě je 91,92 aminokyselin. Kdyby všechny exony měly tuto délku, pak by hranice domény od hranice exonu nemohla být nikdy vzdálená více než čtyřicet pět aminokyselin, což uvádí námi získaná data na pravou míru. I přesto stojí za povšimnutí, že pro začátky exonů byly nejčastěji naměřeny hodnoty v rozmezí -1 až -10 a pro konce exonů v rozmezí 1 až 10, a také případů, kdy se hranice exonů a domén sešly na stejné aminokyselině, bylo více, než by se dalo statisticky předpokládat. Na tomto místě je už korelace jasně patrná a dává za pravdu teorii, že případný pohyb exonů v rámci genomu organismu by mohl způsobovat znovupoužití domén v jiných proteinech. Za zmínku také určitě stojí seznam domén, který se dá nalézt na konci výsledné zprávy. Pro každou doménu je uveden počet výskytů v rámci analýzy, průměrná vzdálenost od začátků/konců exonů a směrodatná odchylka. Můžeme teoretizovat, že domény s nejvyšším výskytem by měly být ty, které se účastnily genetické rekombinace a mělo by se to patřičně odrazit v jejich naměřených hodnotách. A opravdu jisté druhy domén vykazovaly vzhledem k počtu výskytů nízké hodnoty průměrné vzdálenosti, ale co je důležitější i nízké hodnoty u směrodatné odchylky. Naopak u jiných vysoce zastoupených domén naměřené hodnoty odpovídaly spíše nahodilému rozložení v rámci proteinů. Některé příklady obou extrémů jsou uvedeny v tabulce 4.2 (v závorkách jsou uvedeny směrodatné odchylky).
20
Název domény
Počet výskytů
Prům. vzdálenost od začátku exonu
Prům. vzdálenost od konce exonu
CCP
137
3,12 (6,49)
-0,37 (8.95)
EGF_CA
207
3,97 (10,05)
-2,40 (13,42)
SPEC
118
-0,99 (38,14)
-0,42 (34,77)
SH3
120
5,29 (38,86)
4,07 (31,23)
Tabulka 4.2
21
5.
Závěr
Program prezentovaný v této práci prošel během svého vývoje několika fázemi, dokud nevykrystaloval do současné podoby. Nicméně hlavní požadavky byly od počátku jasné. Program měl být automatickým, snadno použitelným, ale přesto užitečným nástrojem pro analýzu biologických dat, který měl poskytnout výsledky v dobře čitelné, snadno pochopitelné podobě, čehož bylo doufám dosaženo. Korektnost výstupu programu se pokládá za samozřejmost. Uživatel dostává na výběr z několika módů práce, včetně možnosti vlastní segmentace analyzovaného proteinu. Všechny módy mají své využití a mohou při správném využití poskytnout hodnotné výsledky. I když jsem s fungováním programu spokojen, uvědomuji si, že je vždy co vylepšovat a existují věci, na kterých je možno do budoucna zapracovat. Příkladem může být například nízká hodnota u počtu aminokyselin s přiřazenou sekundární strukturou. Ideální by bylo kdyby program mohl spolupracovat s účinným predikčním nástrojem, který by mohl doplnit zbývající části. Takových nástrojů je nepřeberné množství, bohužel ty opravdu kvalitní z nich pracují formou, jež je pro využití programem zcela nevhodná. Jedná se o proces, kdy jsou výsledky obdrženy e-mailem na dohodnuté adrese. Jinou možností je využití dalších databází nebo případné rozšíření hledání v PDB databázi i na proteiny s menší mírou podobnosti. Dalším vylepšením by se mohlo jevit spojení uživatelského a statistického módu, jinými slovy zbavení uživatele omezení analyzovat pouze jeden gen, jehož protein si sám rozdělí na segmenty a výsledky následně vyhodnotit statisticky. Během psaní programu se několikrát vyskytla nutnost provést určitá rozhodnutí, jako byl způsob mapování hranic exonů na protein nebo přiřazování jednotlivých hranic domén k hranicím exonů apod. Tato rozhodnutí měla vliv na celý program, především na výsledky statistického módu. Snažil jsem se k těmto problémům přistupovat zodpovědně, případně se o jejich řešení poradit se zkušenějšími lidmi než jsem já sám, ale hlavně jsem se snažil nalézt co nejlepší řešení a tato řešení v této práci představit, a to včetně důvodů pro jejich výběr.
22
6. Použitá literatura [1] Gilbert, W.: Why genes in pieces? Nature, 271 (5645), 1978, s. 501 [2] Mingyi Liu, Shaoping Wu, Heiko Walch and Andei Grigoriev: Exon-domain correlation and its corollaries. Bioinformatics, 21, 2005, s. 3213-3216 [3] Giardine, B. aj.: PhenCode: connecting ENCODE data with mutations and phenotype. Human Mutation, 28, 2007, s. 554-562 [4] Baxevanis, D. Andreas, Quellette Francis B. F.: A practical guide to the analysis of genes and proteins. John Willey & Sons, 2001 [5] Oficiální stránky NCBI – Entrez Utilities http://www.ncbi.nlm.nih.gov/entrez/query/static/eutils_help.html [6] Marchler-Bauer, A. aj.: CDD: a conserved domain database for interactive domain family analysis. Nucleic Acids Research, 35, 2007, s. 237-240 [7] Stránky Wikipedie, http://en.wikipedia.org/wiki/Pfam [8] Marchler-Bauer, A. aj.: CDD: specific functional annotation with the Conserved Domain Database. Nucleic Acids Research 37, 2009, s. 205-210 [9] Oficiální stránky RCSB (PDB) http://www.rcsb.org/pdb/statistics/contentGrowthChart.do?content=total&seqid=100 [10] Oficiální stránky NCBI – GenBank record http://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html [11] Zhang, M. Q.: Statistical features of human exons and their flanking regions. Human Molecular Genetics, 7, 1998, s. 919-932
Obsah přiloženého DVD: zdrojový kód programu, Readme soubor s instrukcemi pro uživatele, knihovny WWW::PDB a GD::Graph, textové dokumenty obsahující příklady zpráv pro standardní a uživatelský mód, kompletní výstupní data statistického módu pro příklad z kapitoly 4, příklady vstupních dat pro statistický a uživatelský mód.
23
