Labtutorials.org

Genomi régiók lefedettségének ábrázolása hőtérképen

In bioinfo, bioinformatics, bioinformatika on June 22, 2018 at 3:02 pm

Hőtérképpel (heatmap) bizonyára már mindenki találkozott az internetes böngészése során, miközben egy időjárás-előrejelző oldalon rátekintett a Magyarország térképre és azt látta, ragyogó napsütésnek csupán a Balaton környékén van nyoma. Legtöbbször valóban a hőmérsékletet, csapadékmennyiséget jelző térképek esetében találkozhatunk ezzel a vizualizációs módszerrel, mely eltérő színek vagy színintenzitások használatával mindenki számára könnyen értelmezhető. Jócskán akad szerepe azonban a biológiában is, például génexpressziós vizsgálatok során.

Leggyakoribb a hideg-meleg érzést keltő zöld-piros, kék-piros vagy kék-sárga-piros, illetve ezen színek közötti átmenetek alkalmazása. Így például a nagyobb mértékben expresszálódó gének piros, az alacsony kifejeződést mutató gének pedig zöld négyzetként jelennek meg egy hőtérképeken. Az átmenetek értelmezéséhez azonban szükségünk van egy alsó-felső(-középső) értéket jelölő skálára is. Ebben a bejegyzésben azt fogom elmagyarázni, milyen módszerekkel tudjuk genomi régiók lefedettségét hőtérképen ábrázolni.

A genomi régiók lefedettségének megjelenítése kissé szofisztikáltabb a microarray vagy RNS szekvenálási módszerek által eredményezett génexpressziós értékek ábrázolásához képest. Utóbbi két módszer „jól meghatározott” értékeket szolgáltat, melyeket táblázatként kezelve akár a ClustVis weboldalon is pillanatok alatt ábrázolhatunk. A ClustVis bár számos R programcsomagot integrál (ggplot, pheatmap, RColorBrewer, stb.), egyetlen hátránya, hogy a megjeleníthető sorok száma maximalizálva van. Ezeket a programokat az ngsdeb szerverünkön keresztül az R programot meghívva is elérhetjük, aki pedig az R-hez készült fejlesztői környezetet, az RStudio-t használja, annak előbb telepítenie kell azokat. Mindhárom megközelítés a következő stílusú hőtérkép(ek)et eredményezné:

heatmap

Mit jelentenek a sorok és az oszlopok?

Minden sor egy gént reprezentál, az oszlopok pedig az adott sorban lévő gén 4 különböző expressziós értékének megfelelő színintenzitást veszik fel. Ezen az ábrán az A1 és A2, illetve a B1 és B2  oszlopok együtt értelmezendők, mivel azonos kondíciót (A1 és A2: kezeletlen; B1 és B2: kezelt) ábrázolnak. Ahhoz, hogy ne egy pepita képet lássunk, lehetőségünk van hasonlóság alapján klaszterezni az értékeket. Többféle klaszterezési mód létezik (hierarchikus, korreláció-alapú, Euklideszi távolság, Pearson, Spearman, Kendall-féle tau, stb.); a jobb oldali hőtérkép egy hierarchikus klaszteranalízis eredményét ábrázolja.

 

No, de hogyan lehet genomi régiók lefedettségét ábrázolni?

Egy korábbi bejegyzésben már volt szó a HOMER programcsomagról és a HOMER tag directory-jairól. Ezeket a tag könyvtárakat, amelyek a genomra térképezhető leolvasásokat tartalmazzák, a ChIP-seq_anal.sh pipeline-unk az alapanalízis során minden mintára egységesen létrehozza.

A leolvasások általi lefedettség értékek arányosak az adott pozícióban megvalósuló fehérjekötés gyakoriságával és erősségével, így tehát alkalmasak arra, hogy két, vagy akár több tíz mintában összehasonlíthassuk adott régiók egy bizonyos, vagy több fehérje általi lefedettségét.

Ehhez a következő parancssort kell begépelnünk:

annotatePeaks.pl peaks.bed hg19 -size 2000 -hist 50 -ghist -fragLength 150 -d ../tag_directory_of_the_sample/ > output_table.txt

Mi mit jelent?

Az annotatePeaks.pl maga a parancs;

peaks.bed fájl kell, hogy tartalmazza a régiókat, ahol a lefedettséget szeretnénk “megszámolni”;

a hg19 megadja a genomot (hg19: humán genom 19-es összállítása) és annak verzióját, de természetesen ez a fajtól és verziótól függően opcionális (hg18, mm9, mm10, stb.). Azt kell tudnunk, hogy a program ezen régiók határainak számtani közepét fogja középpontnak tekinteni, amely nem minden esetben jelenti a kötőhely valós csúcspontját. A nagymértékű eltolódások elkerülése végett legtöbbször a peak-ek summit pozícióját tartalmazó BED fájlt szoktuk megadni;

a -size 2000 paraméterrel pedig – jobb esetben – a summit pozíciókat egységesen kiterjesztjük -1000/+1000 bázispárnyi (bp-nyi) régióra. Így tehát összesen 2000 bp-nyi régión számolja a lefedettséget.

A -hist 50 megadja, hogy a 2000 bp-on belül hány bp lefedettségértékét átlagolja, azaz ebben az esetben 50 bp jelent egy ún. bin-t. Mivel a teljes régió 2000 bp lesz, mi pedig 50 bp-onként szeretnénk visszakapni egy értéket, azt jelenti, hogy a kimeneti fájlunkban 20+1+20 oszlop lesz (-1000 -950 -900 … 0 … +900 +950 +1000).

A -ghist kapcsoló régiónként kiszámolja a bin-ek lefedettségét;

a -fragLenght 150 paraméterrel az automatikus fragmenthossz számítása helyett azt tetszőlegesre állítom (opcionális). Megj.: A HOMER alapvetően a pozitív és negatív szálra térképeződő leolvasások egymáshoz viszonyított eloszlása alapján számolja ki az átlagos, elméleti fragmenthosszt.

Végezetül pedig a -d kapcsoló után kell megadnunk egy, vagy akár több minta HOMER tag könyvtárát.

A > karakter után a kimeneti fájl neve következik, melynek kiterjesztése .txt kell legyen, hogy a későbbiekben könnyen ábrázolni tudjuk..

 

Ábrázolás

A sok-sok sorból és oszlopból álló táblázatunkat ezek után már csak meg kell jelenítenünk. Ehhez mi a Java TreeView programot használjuk, ami bárki számára ingyen elérhető, letölthető a http://jtreeview.sourceforge.net/ weboldalról.

A program megnyitását követően tallózzuk be a táblázatunkat (File -> Open)! Automatikusan felajánlja, hogy az adott mappában csak a CDT vagy PCL kiterjesztésű fájlokat mutassa; itt ki kell válasszuk az All Files opciót. Innentől a táblázat megnyitása már csak annak méretén és a számítógépünkben található memória mennyiségétől függ. 🙂 Ha a betöltés leáll, próbáljuk meg egy jobb teljesítményű számítógépen ábrázolni, vagy osszuk két, vagy akár több részre a táblázatunkat!

Meg kell jegyezzem, hogy a TreeView-ban való megjelenítéshez a táblázat első két oszlopának egy-egy (megegyező) azonosítót kell tartalmaznia. Ez lehet akár a gének, a pozíciók neve, vagy akár egy sorszám. A táblázat első sorában is egy, az adott oszlopra vonatkozó azonosítónak kell szerepelnie. Ezek hiányában a táblázatban szereplő lefedettségértékek első két oszlopa és az első sora fogja ezt kiváltani, melyek emiatt nem lesznek ábrázolva.

Sikeres betöltést követően a következő kép fogad minket:

treeview1

 

Ezt követően a Settings -> Pixel Settings-re kattintva egy új ablak jelenik meg, ahol mind a négy skála esetében a Fixed Scale-ről át kell klikkelnünk a Fill-re:

 

treeview2

 

Szintén a Pixel Settings panelben tudjuk megváltoztatni az pozitív/negatív, a hiányzó és a nulla értéket felvevő bin-ek színét. Esetemben a “Positive”-at feketére, a “Zero”-t pedig fehérre változtatva, továbbá a skálát (Value) 3-ról 5-re átállítva a következő eredményt kaptam: 

 

treeview3.PNG

 

Az 5-ös érték azt jelenti, hogy egy bin-ben ha a lefedettségérték eléri az 5-öt, fekete színként jelenik meg. A skála csökkentésével kontrasztosabb képet kapunk, mivel egyre több olyan bin lesz, amely eléri a beállított értéket.

Az ábrán megjelenített kötőhelyek előzetesen egy másik minta középső 50 bp-nyi régiójának lefedettsége alapján sorba lettek rendezve, ennek köszönhető ez a mintázat.

Az elkészült ábrát az Export -> Export to Image menün keresztül tudjuk letölteni. Azonban itt is figyelni kell néhány dologra:

treeview4

Be tudjuk állítani az X- és az Y-tengely méretét, tehát az exportálni kívánt kép méretarányát. Amennyiben idő közben belekattintottunk a kirajzolódó képbe, a Selection Only-nál egy pipa jelenik meg – ezt ignorálnunk kell, ha az összes sort szeretnénk képként kinyerni. Végül pedig a bal oldalon kékkel kiemelt -1000 és Gene sorokra (amelyek esetemben az első sor és oszlop tartalmát jelzik) a Ctrl billentyű lenyomása mellett rá kell kattintani, hogy a kék kiemelés eltűnjön. Amennyiben ezt nem tesszük meg, a több tízezer soros táblázatunk minden sorának neve szerepelni fog az ábránkon. Természetesen néhány sor/oszlop megjelenítése esetén ez még akár jól is nézhet ki.

No, és mit rontottam el? Mivel -1000 van a név helyett, azt jelenti, hogy az első bin-emet elvesztettem az ábráról, mert nem duplikáltam meg az első oszlopot. 🙂

Utolsó lépésként mentsük el az ábránkat.

 

 EEM (1)
Az Emberi Erőforrások Minisztériuma ÚNKP-17-3-IV-DE-140 kódszámú Új Nemzeti Kiválóság Programjának támogatásával készült.

 

 

 

 

 

 

 

 

 

Advertisements
%d bloggers like this: