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

Adatok megjelenítése genom böngészőben

In bioinfo, bioinformatics, bioinformatika on June 16, 2018 at 5:10 pm

Az alapanalízis során “készülnek” olyan fájlok, amelyek alkalmasak genom böngészőkben történő megjelenítésre. Ilyen például a BAM fájl, amely pontosan megmutatja, hogy egy adott minta esetében a genom mely szakaszaira és mennyi read (leolvasás) térképeződött, illetve jelzi a leolvasásokban azokat a nukleotidokat is, amelyek a referencia genomhoz képest eltérést mutatnak.

De megjelenítésre alkalmasak a különféle BEDGRAPH-ok is, melyek a genom teljes “hosszában” Gauss-i eloszlást mutató csúcsokat rajzolnak ki; ebben az esetben a csúcsok magassága és szélessége attól függ, adott pozícióra hány read térképeződött. Megjeleníthetjük ugyanakkor a BED fájlokat is, önmagukban vagy BEDGRAPH-okkal együtt; utóbbi abban az esetben lehet hasznos, ha szeretnénk látni, egy kívánt lókuszon mely csúcsok érték el azt a küszöbértéket, amely alapján a MACS2 vagy a HOMER programok őket a kötőhelyek csoportjába sorolta, vagy sem.

Talán az IGV (Integrative Genomics Viewer) a légszélesebb körben használt genom böngésző. Számos hasznos tulajdonsággal rendelkezik; az egyszerű ide-oda húzgálás vagy a vizsgálni kívánt genomi régió célzott megjelenítése mellett egy listát létrehozva egyszerre több régió is megjeleníthető, amely funkció különböző ábrák készítésekor nagyon hasznos szereppel bír. A betöltött track-ek, ún. mintasávok sorrendje, elnevezése, a sáv mérete vagy épp színe egy-két klikkeléssel könnyen megváltoztatható.

Honnan tölthető le az IGV?

A https://software.broadinstitute.org/software/igv/ weboldalon a Downloads panel alatti nyílra kattintva megjelenik egy újabb oldal, ahol  kiválaszthatjuk az általunk használt operációs rendszerrel kompatibilis verziót, ill. ha nem rendelkezünk Java-val, közvetlenül innen azt is letölthetjük – a Java ugyanis, mint kiegészítő modul (plug-in), elengedhetetlen a IGV futtatásához.

java1

Letöltést és indítást követően egy “üres ablak” fogad bennünket; a fontosabb paneleket funkcióit az alábbi ábrán feliratozva találjátok:

igvvvvvvv_2

Végezetül pedig nézzük meg, hogyan néznek ki a fent említett fájltípusok az IGV-ben:

ddd

Első mező (coverage): BEDGRAPH fájlt, ami a transzkripciós faktor kötőhelyeket (ún. peak-eket) mutatja.

Második mező (peaks): a kék négyzetek a BED fájl tartalmát jelenítik meg, azaz jelölik a peak-ek pontos kezdő- és végpozícióit. Jól mutatja, mely feldúsulások valódi kötőhelyek és melyek nem.

Harmadik mező (BAM): a BAM fájl tartalmát jeleníti meg, amely az adott genomi pozícióra térképeződött leolvasásokat mutatja.

 

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.

 

A szuper-enhanszerekről és azonosításukról

In bioinfo, bioinformatics, bioinformatika on June 12, 2018 at 11:53 am

Az enhanszerek a génektől távol elhelyezkedő DNS szakaszok, mely szakaszokon belül enhanszer elemek találhatóak. Az enhanszer elemeket kötő fehérjék transzkripciós faktorokként működnek, melyek kötődésük révén – nevükből adódóan is – a transzkripciós szabályozásban vesznek részt. Egy tipikus-enhanszer régióhoz képest a szuper-enhanszer első sorban abban különbözik, hogy egy viszonylag nagy DNS szakaszon belül sok, akár több tíz enhanszer elemét is magába foglalhatja és  az egyedi, kifejezetten aktív enhanszer kötések közötti távolság nem haladja meg a 12.5 kb-nyi távolságot  (1. ábra)

ábra1

1. ábra.  Az Oct4 és Sox2 fehérjék egy (tipikus-)enhanszer és egy szuper-enhanszer régiója genom böngészőben megjelenítve. Megjelenített adat: ChIP-seq. (Whyte et al. Cell, 2013)

Jellemző a szuper-enhanszerekre a Mediátor komplex egy kiemelt tagjának, a MED1 fehérjének a jelenléte, a nyitott kromatin régióra jellemző H3K27ac és DNáz I hiperszenzitivitás, továbbá az aktív, korábban aktív, vagy aktiválható régiót jelölő H3K4me1 markáns jelintenzitása is (2. ábra).

 

ábra2

2. ábra. Az aktív kromatin- és enhanszer markerek (Med1, H3K27ac, H3K4me1 és DNaseI) kötéserősségben becsült különbsége a tipikus- és a szuper-enhanszereken. (Whyte et al. Cell, 2013)

Úgy tartják, hogy a szuper-enhanszerek a sejtek identitásáért felelős géneket szabályozzák, tehát egy tumoros sejtben többnyire onkogének közelében lelhetőek fel. Meg kell azért jegyeznem, bármilyen enhanszerről is beszélünk, csupán a pozíciójuk alapján nem tudhatjuk, pontosan mely gént vagy géneket szabályozzák, hiszen a genom térbeli szerkezete DNS-hurkok kialakulása révén lehetővé teszi két, vagy akár több távoli DNS szakasz térben egymás mellé kerülését (3. ábra). A fentebb említett Mediátor komplexnek is hasonló szerepe van: képes a promóter és az enhanszer régiók között távoli interakció révén kapcsolatot teremteni. De bármennyire is tűnik mindez bonyolultnak, az utóbbi években rengeteg, a kromatin interakció megismerésére alkalmas módszert dolgoztak ki (3C-seq, Hi-C, Capture-C, stb.). Egy tipikus-enhanszerhez képest a gének kifejeződésében sokkal markánsabb hatást képesek elérni (már csupán az additív hatásukat figyelembe véve is), tehát a szuper-enhanszerek feltérképezése és az általuk szabályozott gének azonosítása semmiképp sem hátrány (3. ábra)

 

ábra3

3. ábra. Egy tipikus-és egy szuper-enhanszer hatása az általuk szabályozott gének kifejeződésének tekintetében. (Evan et al., Clin Cancer Res, 2017)

 

Hogyan azonosíthatóak a szuper-enhanszerek?

Az előző bejegyzésben már volt szó a ChIP-seq módszerről, amely egy viszonylag könnyen kivitelezhető, mára már jól optimalizált módja az enhanszer régiók genom-szintű feltérképezésének. Általánosan elfogadott, hogy H3K27ac jel alapján határozzák meg a vizsgált sejt szuper-enhanszereit, de akadnak tanulmányok, amelyekben egy, a sejt számára meghatározó szereppel bíró fehérje kötése révén definiálták őket. Utóbbi szerint jártunk el mi is, amikor kifejezetten az ösztrogén receptor alfa (ERα)-vezérelt szuper-enhanszerekre voltunk kíváncsiak. Ahhoz tehát, hogy szuper-enhanszereket tudjunk azonosítani, vagy H3K27ac, vagy egy fontos fehérje ChIP-seq adataira van szükségünk. Szóba jöhet még az ATAC-seq (Assay for Transposase Accessible Chromatin with high-throughput sequencing) módszer által nyert adat is, mely egy sejt összes nyitott kromatin régióját képes a behasító transzpozáz enzim segítségével feltérképezni. Ebben az esetben annyi hátrányunk lehet, hogy ha csupán az ATAC-seq adatok állnak rendelkezésünkre, nem tudjuk elkülöníteni egy nyitottnak vélt régióról, hogy aktivátor vagy represszor fehérje köti-e. Mindezt a bizonytalanságot egy hiszton ChIP-seq, vagy egy, az aktív enhanszerekre jellemző bármilyen fehérje (pl. p300) ChIP-seq adataival összevetve elsimíthatjuk.

Volt szó az előző bejegyzésben a HOMER programcsomag pipeline-ba integrált eszközeiről is. A szuper-enhanszerek predikcióját a pipeline bár nem tartalmazza, de egy egyszerű parancssorral elvégezhető és ehhez csupán egy jól megválasztott ChIP-seq adat HOMER tag directory-jára van szükségünk. A tag directory-k (tag könyvtárak) a BAM fájlhoz hasonlóan sort-olva, azaz genomi pozíció szerint rendezve, lefedettségértékeket tartalmaznak, melyeket kromoszómánként egy-egy *.tags.tsv fájl tárol az alapanalízis során. A szuper-enhanszer régiók meghatározásához a következő parancsot kell begépelnünk:

findPeaks <tag directory> -i <input tag directory> -style super -o auto

A findPeaks maga a parancs, a tag directory-t, pontosabban a *.tags.tsv fájlokat tartalmazó mappa nevét a teljes elérési úttal a vizsgálni kívánt mintára vonatkozóan kell megadni. A -i kapcsolót követően egy inputként szolgáló minta tag könyvtárára hivatkozhatunk, de ettől, ha ilyen nem áll rendelkezésünkre, eltekinthetünk. A -style super kapcsolóval jelezzük a findPeaks parancsnak, hogy szuper-enhanszer régiókat keressek, a -o kapcsolót követően pedig csupán egy kimeneti könyvtárat kell megadnunk.

Még több infó (angol nyelven) a szuper-enhanszerek HOMER programcsomag segítségével történő meghatározásáról itt.

Egyéb, sokak által használt, a szuper-enhanszerek prediktálására szintén alkalmas program még a ROSE (Rank Ordering of Super-Enhancers), illetve létezik egy dbSUPER nevű adatbázis is, ahol megtalálhatjuk számos primer sejt vagy sejtvonal szuper-enhanszer régióinak gyűjteményét és azok egyéb jellemzőit.

 

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.

 

 

 

 

%d bloggers like this: