Labtutorials.org

Archive for the ‘bioinformatika’ Category

A szekvencia motívumok szerepe a transzkripció szabályozásában

In bioinfo, bioinformatics, bioinformatika, DNA, molecular biology on June 8, 2019 at 8:18 pm

A magreceptorok motívumai

Szerző: Dr. Nagy Gergely

A magreceptorok csoportosítása

A magreceptorok beszédes névvel rendelkeznek, amely kifejezi, hogy tipikusan olyan fehérjékről van szó, amelyek jelmolekulákat ismernek fel és a sejtmagban fejtik ki hatásaikat. E fehérjék többsége sok más receptorhoz hasonlóan dimer formájában működik, viszont más receptorokkal ellentétben a magreceptorok nem membránkötöttek, hanem oldott formában a citoszólban, illetve a sejtmagban találhatóak. Ligandjaik a membránokon áthatolni képes lipid molekulák, beleértve a zsíroldékony hormonokat, vitamin-, szteroid- és zsírsavszármazékokat. A magreceptorok neve arra is utal, hogy nemcsak fehérje-fehérje kölcsönhatásokon keresztül képesek jelet továbbítani, hanem a sejtmagban, közvetlenül a DNS-hez kapcsolódva, mint transzkripciós faktorok szabályozzák a géneket. Azáltal, hogy ilyen rövid úton eljut a jel a célgénekhez, lényegesen lecsökken a sejtek adott körülményre adott válaszideje, nem úgy, mint a membránreceptoroktól induló, soklépéses jelátviteli útvonalak esetében.

A magreceptor szupercsalád emlősökben előforduló 19 családját 4 osztályba sorolják: a szteroid hormon receptorokra (I. osztály, 2 család), a retinoid X receptorral (RXR-rel) heterodimert alkotó ligandkötő receptorokra (II. osztály, 5 család), a dimerizáló árva receptorokra (III. osztály, 6 család) és a monomer árva receptorokra (IV. osztály, 6 család) (Mangelsdorf et al., Cell, 1995; Nuclear Receptor Nomenclature Committee, Cell, 1999; Evans and Mangeldorf, Cell, 2014). Az I. osztály tagjai homodimert alkotnak, és kizárólag szteroid hormonokat ismernek fel. A II. osztály tagjai a ligandok széles spektrumát képesek felismerni, mint például a tiroid hormont, az A- és D-vitamin, a zsírsavak, valamint a koleszterol származékait (Dawson and Xia, Biochim Biophys Acta., 2012). Az árva receptorok onnan kapták a nevüket, hogy eleinte nem sikerült a ligandjaikat azonosítani, később mégis kiderült, hogy a III. osztály fele képes valamilyen lipid természetű molekulát kötni. A valódi árva receptorok nem rendelkeznek működőképes ligandkötő doménnel, hanem mint más transzkripciós faktorok, fehérje-fehérje kölcsönhatások által vagy például foszforilációval szabályozódnak.

Magreceptor motívumok

A magreceptorok általában az AGGTCA motívumokat ismerik fel. Dimerek esetében ez a szekvencia kétszer szerepel egymás mellett, ezért magreceptor félhelynek is nevezik. Helytálló ez az elnevezés azért is, mert egy hatbázisos motívum, főleg, ha beleszámoljuk a lehetséges szekvencia variációkat, túl gyakran található meg a genomban (<46 vagy <45 = ~1000 bázisonként) és túl könnyen alakulhat ki véletlenszerű mutációk során ahhoz, hogy rendelkezzen a szükséges szelekciós erővel a génkifejeződés megfelelő szabályozásához. Hogy a magreceptor dimerek megtalálhassák az adott körülmények között szükséges szabályozó elemeiket, elsősorban a félhelyek egymáshoz viszonyított iránya és távolsága a felelős. Az I. osztály receptorai esetében például a félhellyel a tükörképe (például TGACCT) áll szemben, három bázissal elválasztva. Ezt a palindrom szekvenciát úgynevezett fordított ismétlődésnek vagy inverted repeat (IR)-nek nevezik, amit, mivel három, nagyjából véletlenszerű bázis van a közepén, IR3-ként emlegetnek. Ebben az osztályban az ösztrogén receptorok kivételesek az AGGTCA félhelyükkel, mivel az összes többi szteroid hormon receptor az AGAACA (illetve TGTTCT) szekvenciát preferálja.

A II-III. osztály dimerei ezzel szemben kivétel nélkül két, egymást azonos irányban követő magreceptor félhelyet, úgynevezett direct repeat (DR) elemet ismernek fel, ahol az elválasztó bázisok száma a leginkább meghatározó; és DR0-tól DR5-ig minden lehetőségre találunk specifikus dimereket (Umesono et al., Cell, 1991; Evans and Mangeldorf, Cell, 2014); de írtak már le működőképes DR8-at is. A DR0-t például GCNF homodimer, a DR1-et PPAR/RXR heterodimerek, valamint TR2/4 és HNF4 (homo)dimerek, a DR2-t RAR/RXR heterodimerek és REV-ERB (homo)dimerek, a DR3-at VDR/RXR heterodimerek, a DR4-et THR/RXR és LXR/RXR heterodimerek, a DR5-öt pedig RAR/RXR heterodimerek ismerik fel. Ezekben az osztályokban is vannak IR felismerő magreceptorok, illetve léteznek olyan dimerek is, amelyek, például a ligand minőségétől függően, különböző távolságra lévő félhelyeket kötnek. Az RAR/RXR heterodimerek az előbb említett DR5 és DR2 kötés mellett a DR1 elemeket is használhatják, a PXR/RXR heterodimerek esetében pedig leírták, hogy a pregnánszármazékok és másodlagos epesavak rugalmas kötése a konformációváltozás hatására különösen rugalmassá teszi a DR elemek felismerését is (Wu et al., Drug Discov Today, 2013; Frank et al., J Mol Biol., 2005).

Mivel a magreceptorok félhelye túlságosan gyakran fordul elő a genomban ahhoz, hogy specifikusan működhessen, a IV. osztály magreceptorai esetében a hat bázison felül általában további bázisok is hozzájárulnak az erős DNS-fehérje kölcsönhatáshoz. Ezek a bázisok minden érintett család esetén a félhelyek 5’ kiterjesztését jelentik. Az NR0B család kivételt képez ez alól, mert nem rendelkezik DNS-kötő doménnel (Ensembl). Az NR4A (NUR/NOR) fehérjék az AA-AGGTCA (Wilson, Milbrandt, Science, 1992), az NR3B (ESRR) és NR5A (SF-1, LRH1) családok tagjai a (T)CA-AGGTCA (Johnston, Mertz, Mol. Endocrinol., 1997; Lala, Parker, Mol. Endocrin., 1992; Laudet, Curr. Biol., 1995), az NR1F (ROR) fehérjék pedig az (A/T)AA(C/G)T-AGGTCA szekvenciákat ismerik fel (Giguere, Otulakowski, Genes Dev., 1994; IJpenberg, JBC, 1997). Ez utóbbi kiterjesztett félhely, az úgynevezett ROR válaszadó elem (RORE) azonban részét képezheti DR elemeknek is. Mind a PPAR/RXR, mind pedig a REV-ERB dimerek nagy affinitással kötik a kiterjesztett DR – DR1, illetve DR2 – elemeket, és ezeknek az elemeknek – a magreceptorok expressziós szintjének és az adott motívumokhoz való affinitásának függvényében – fontos szerepe van a sejtek napi ciklusának szabályozásában (Harding, Lazar, MCB, 1995; Duez, Stael, FEBS Letters, 2008; Zhang, Lazar, Science, 2015). Ez a kiterjesztés teheti specifikussá a DR1 elemek PPAR/RXR általi kötését a TR és HNF4 (homo)dimerekkel szemben, valamint a DR2 elemek REV-ERB általi kötését az RAR/RXR heterodimerekkel szemben.

Kiterjesztett magreceptor motívumok keresése

Az elmúlt három évtizedben lényegében négy olyan tényezőt azonosítottak, amely meghatározza a magreceptorok specifikus DNS kötését: a félhelyek szekvenciáját (AGGTCA vagy AGAACA), egymáshoz viszonyított irányát (IR vagy DR), egymástól való távolságát és 5’ kiterjesztését. Az alapszabályokkal ugyan tisztában vagyunk, de nem ismerjük minden magreceptor pontos szekvenciaigényeit. Ehhez az NGS módszerek, például a ChIP-seq vagy akár ATAC-seq és ezek elemző módszerei nagy segítséget nyújtanak (Heinz, Mol. Cell, 2010), mégsem mindig szembetűnő a különbség a különböző magreceptorok motívumai között. Egyszerre többféle DR vagy IR elem kiterjedt használata esetén, például az RXR cisztróm vizsgálatakor, megtörténhet a különböző motívumok teljes összekeveredése, összeolvadása („kiátlagolódása”), tehát akár egyetlen félhelyre redukálódása is (Dániel and Nagy, Genes. Dev, 2014). Mivel kisebb a kiterjesztett motívumok száma, mint azoké, amelyek nem rendelkeznek valamilyen 5’ kiterjesztéssel, a de novo motívumkeresések eredményeiben ezek általában nem hangsúlyosak vagy teljesen hiányoznak. Léteznek „trükkök” a motívumok szétválasztására a de novo motívumkeresés eszköztárában, ám ezek is szenvednek a módszernek attól az általános korlátjától, hogy csupán a bázisok gyakoriságát veszik figyelembe, ezekhez nem rendelik hozzá a fehérjekötés erősségét.

A de novo motívumokat kiegészítendő, kifejlesztettem egy motívum optimalizáló módszert, amely a motívumok bázisainak a fehérjekötéshez való hozzájárulását méri. Ennek segítségével lényegében egyetlen ChIP-seq minta alapján nagyon pontosan meghatározható volt a PPARg félhelyének a kiterjesztése. Ez a motívumkeresésen és -térképezésen alapuló módszer valójában bármely transzkripciós faktorra specifikus ChIP-seq adaton jól működhet, feltárva e fehérjéknek a gyakori motívumokon felüli szekvenciaigényeit. Bázisok kettőseit felhasználva több dimenzióban is tesztelhető a kettősök fehérjekötéshez való hozzájárulása, ezáltal akár különböző és átfedő motívumkiterjesztések, illetve távolabbi, úgynevezett szatellit elemek is azonosíthatóak. A PPARg mellett nagyszámú magreceptorra specifikus ChIP-seq adat érhető el nyilvánosan, például az NCBI SRA adatbázisában. Mivel elképzelhető, hogy a TR2/4 és HNF4 (homo)dimerek, valamint a THR/RXR és LXR/RXR heterodimerek DR1, illetve DR4 motívumaiban is található valamilyen eltérés, amely a specificitásukat adja, érdemes lehet e magreceptorok esetében is elvégezni a motívumok fehérjekötéssel kapcsolt optimalizálását; valamint feltételezhető, hogy az RAR, PNR és COUP-TF magreceptorok rugalmasabb DNS kötése mögött is van egy általános szabályszerűség. Kérdéses továbbá az is, hogy vajon minden, DNS-kötő doménnel rendelkező monomer árva receptor, beleértve az NR2E családot (TLX, PNR) is, vagy akár további dimerizáló receptorok is rendelkeznek-e kiterjesztett motívummal, illetve, hogy ezek a kiterjesztések mutatnak-e további specificitást.

Ha választ kapunk ezekre a kérdésekre, az közelebb visz a magreceptorok és motívumaik koevolúciós történéseinek a megismeréséhez is, amely egy sokkal teljesebb képet adhat a transzkripciós faktorok általi génszabályozásról és annak evolúciójáról.

Az Emberi Erőforrások Minisztériuma ÚNKP-18-4-DE-318 kódszámú
Új Nemzeti Kiválóság Programjának támogatásával készült.

Advertisements

GRO-seq

In bioinfo, bioinformatics, bioinformatika, molecular biology, RNA on May 24, 2019 at 7:57 pm

Szerző: Dr. Nagy Gergely

A módszer, amely elől nem bújhat el egyetlen RNS molekula sem

A teljes genom szintű (Global) Run-On (GRO) szekvenálás a naszcens transzkriptóm meghatározására alkalmas újgeneráriós szekvenálási (NGS) módszer. A transzkriptóm általánosságban a sejtek teljes RNS állományát jelenti, a GRO-seq lényege azonban éppen az, hogy csak az egy adott pillanatban átíródó RNS molekulákat, sőt azoknak is csak az éppen átíródó részét, tehát gyakorlatilag az átírást végző RNS polimeráz komplexek helyét mutatja meg a genomban. Ez úgy érhető el, hogy egy szarkozil nevű detergenssel (tisztítószerrel) meggátolják, hogy szabad polimerázok csatlakozzanak a DNS-hez, ellenben a már elkötelezett komplexek tovább tudnak működni. A run-on gyakorlatilag a polimerázok korlátozott „újraindítását” jelenti izolált sejtmagokban, jelölt nukleotid-trifoszfát szubsztrátok felhasználásával. Néhány tíz nukleotid felépítése elegendő ahhoz, hogy az RNS molekulák darabolása után a jelölés segítségével kifogják az új szakaszokat, és meghatározzák a bázissorrendjüket.

A GRO-seq eljárás során arra is ügyelnek, hogy az RNS molekulák bázissorrendjének az iránya is megismerhető legyen. Ehhez előbb a molekulák 5’, majd 3’ végéhez kapcsolnak végspecifikus adaptort. Az RNS molekulák töredékeinek az 5’ végén azonban nincs feltétlenül szabad foszfát csoport. A későbbi lépésekhez az mRNS-ek 5’ „sapkáját” el kell távolítani (TAP), valamint end-repair-rel mind az 5’, mind pedig a 3’ vég javítható (például foszforilálható, illetve defoszforilálható; PNK). A különböző adaptorokkal közrefogott RNS molekulákból reverz transzkripcióval DNS-t hoznak létre, majd ezt sokszorozzák (PCR) a szekvenáláshoz.

A rövid szekvencia-leolvasások tehát megmutatják a polimerázok általi RNS szintézis helyét és irányát, amely kiválóan használható bizonyos nyomon követéses kísérletek esetében. Ha kíváncsiak vagyunk, milyen gének kapcsolnak be vagy ki egy stimulus hatására, érdemes néhány vagy néhány tíz perces felbontásban mintát venni. Ilyen módon láthatóvá válik, hogy a polimeráz percenként 2,5-3 kilobázis távolságot halad a szabályozott géneken. Azonban vannak olyan hosszú gének is, amelyek átírásához órák kellenek, és ez idő alatt az mRNS teljes érése és fehérjére „fordítása” sem történhet meg. A rövidebb gének viszont hamar nagy mennyiségű fehérjeterméket eredményezhetnek, és amennyiben ezek képesek a transzkripciót szabályozni, például mint transzkripciós faktorok, megfigyelhetjük az általuk be-, illetve kikapcsolt gének egy újabb hullámát, amely szó szerint a gének lefedettségén is látható. Ha egy hosszú gén előbb indukálódik, majd nem sokkal később gátlódik, egy „csúcs” jelenik meg rajta, amely idővel (későbbi időpontokban) a gén vége felé „vándorol”. Kellően nagyszámú időpont vagy jól időzített időpontok használatával teljes transzkripciós kaszkádok térképezhetőek fel a módszer segítségével.

A GRO-seq-kel nyert génexpressziós adatokat azonban más okokból kifolyólag sem könnyű értelmezni. Például sokszor nincs egyszerű összefüggés a különböző RNS molekulák szintézisének gyakorisága és az érett RNS szintje között. Az érés sem feltétlenül egyszerű folyamat, de összességében talán az érett RNS molekulák stabilitása (féléletideje) a leginkább meghatározó tényező a génexpressziót tekintve. Csupán GRO-seq adatokból tehát nem sokat tudhatunk meg a génexpressziós szintekről, annál többet a génexpresszió kezdeti szabályozásáról. A polimerázok ugyanis nemcsak a géneken találhatóak meg, hanem transzkripciót mutatnak minden aktív szabályozóhelyen is, még ha nem is következik utána lánchosszabbítás (elongáció).

Polimerázok mindenütt

Bőven a GRO-seq előtt ismert volt, hogy nagyszámú polimeráz gyülekezik a promótereken, de ezeknek tipikusan csak töredéke tudja megkezdeni a génen való továbbhaladást, a többi csak vesztegel (pausing). Ez a GRO-seq adatok alapján úgy néz ki, hogy a gén kezdeti szakaszán van egy csúcs – rövid, úgynevezett abortált átiratokból –, ami többnyire jelentősen magasabb, mint a gén további szakaszának a – transzkriptumok elongációjából fakadó – lefedettsége. Előfordul az is, hogy a promóter jelentős aktivitást mutat, a génen pedig alig vagy egyáltalán nem detektálható transzkripció, valószínűleg valamilyen további aktiváló jel hiánya miatt. Magasan kifejeződő gének esetében viszont nem feltétlenül látható pausing, mivel közel minden megkezdett RNS molekula meghosszabbításra kerül. Ebben az esetben időegység alatt tovább is jutnak a polimerázok, mert gyorsabban tudnak haladni a tartósabban szétválasztott DNS-en.

Az átíródó szabályozó régiók alatt nemcsak a promótereket értjük, hanem az aktív enhanszereket (silencer-eket) is, melyek átírását ugyanúgy érintik a pozitív/negatív stimulusok, mint a fehérjekódoló génekét. Ezt kihasználva a promóterektől akár többszáz kilobázis távolságra elhelyezkedő, az adott stimulus hatására azonos expressziós mintázatot mutató szabályozó régiókat is a génekhez rendelhetjük, amely segíthet azt is megmondani, mely transzkripciós faktorok vesznek részt a szabályozásban. A promóterektől távol eső szabályozó helyeken általában nagy a pausing mértéke – tehát az abortált transzkriptumok aránya –, de ezeken a helyeken is történhet elongáció, melynek hosszú nem-kódoló RNS-ek lesznek a termékei. Elongáció hiányában egyszerűen enhanszer transzkripcióról beszélünk, amely tipikusan mindkét irányban megtörténik (divergens) a szabályozó régióhoz képest – valószínűleg azért, mert itt nincsenek olyan, a polimeráz aktivitás irányát meghatározó szabályozó, úgynevezett válaszadó elemek, mint a promóterek klasszikus elemei, például a TATA-box. Jóllehet, a legtöbb promóteren is jellemző divergens transzkripció, akár elongáció mindkét irányba; sőt többezer olyan fehérjekódoló génpár létezik, amely látszólag egyetlen promóteren osztozik.

Nem-kódoló RNS-ek

Ellentétben a génekkel, a hosszú nem-kódoló RNS termékek hossza a GRO-seq adatok alapján vélhetően nem azonos – minél távolabb jut a polimeráz, annál valószínűbb, hogy nem folytatja tovább az átírást. De mindig vannak kivételek: bizonyos hosszú nem-kódoló RNS-ek egy bizonyos pontig azonos lefedettséget, sőt akár a génekhez hasonló intronkivágódást is mutatnak. Például a „csak” mikroRNS-t kódoló „gének” is így viselkednek. Az a bizonyos pont, ameddig a gének és gén jelleget mutató hosszú nem-kódoló transzkriptumok nagyjából azonos polimeráz sűrűséggel bírnak, a transzkripció terminációs helye. (Csak első ránézésre) érdekes módon a terminációs helyet követően felerősödik a polimerázok jelenléte, majd a hosszú nem-kódoló termékekhez hasonlóan egyre kevesebb tovább hosszabbított terméket látunk. Ez a jelenség valószínűleg az RNS polimerázok lelassulásának tudható be, nem újabb komplexek csatlakozásának. A terminációs helyet követően a polimerázok nem válnak le rögtön a DNS-ről, de a sebességük lecsökken, így gyakrabban lehet detektálni a termékeiket; ez magyarázhatja a – magas expresszió esetén akár többtíz kilobázisos – továbbírást.

A polimerázok lassulása és gyorsulása valamennyire a géneken is érvényesül, attól függően például, hogy milyen a G/C bázisok aránya, milyen a kromatin szerkezete, vagy például van-e aktív szabályozó hely a génen. Főleg a promóterek közelében, de valójában bármelyik intronban lehet enhanszer transzkripciót látni, de az intronokon belül akár más gének promóterei is lehetnek aktívak, és bármelyik irányban keletkezhet róluk, akár elongált RNS termék. Nem könnyítik meg a transzkriptumok azonosítását az alternatív promóterekkel rendelkező gének sem. Referencia annotáció nélkül – illetve hiányos referencia annotáció esetén –, csak a lefedettség adatok alapján, sokszor nem lehet megállapítani, hogy egy hosszabb, alacsonyabb expressziójú transzkriptum variánst látunk-e, vagy egy eddig ismeretlen gént, amely ugyanazon a szálon található, és az ismert gén promótere előtt végződik. Az is előfordulhat, hogy az ismert promótertől downstream helyezkedik el egy eddig ismeretlen, intronikusnak látszó alternatív promóter, amelyet, ha alacsony expressziót mutat, könnyen enhanszernek nézhetünk. Az alternatív terminációs helyek nem gyakoriak, de még nehezebb kezelni őket.

A GRO-seq adatokban nemcsak ismert gének ismeretlen variánsait és sohasem látott enhanszer transzkriptumokat, hanem eddig teljesen ismeretlen, gén jelleget mutató transzkripciós eseményeket is lehet találni. Ilyen esetben meg lehet próbálni a nyitott olvasási keretek és exon-intron határok keresését, ami akár új gének felfedezését is eredményezheti. Az eddig említett transzkriptumok mellett természetesen megfigyelhető a kis sejtmagi és „magvacskai”, valamint a transzfer és riboszómális RNS molekulák expressziója is, bár ezek általában nem mutatnak jelentős időbeli változásokat.


GRO-seq adatok elemzése

A korábban bemutatott ChIP-seq elemző pipeline alkalmas GRO-seq adatok alapelemzésére is, a további elemzésekhez viszont az adatok összetettsége miatt komoly fantáziára is szükség lehet. 🙂

Az Emberi Erőforrások Minisztériuma ÚNKP-18-4-DE-318 kódszámú
Új Nemzeti Kiválóság Programjának támogatásával készült.

DNS szekvencia motívumok azonosítása II.

In bioinfo, bioinformatics, bioinformatika, DNA, molecular biology on January 14, 2019 at 1:15 pm

Szerző: Bojcsuk Dóra

Ahogyan az előző bejegyzést is indítottam, a két leggyakrabban feltett kérdés, amikor motívumanalízisről beszélünk, a következő:

  1. Milyen szekvencia motívumok „dúsulnak fel” a transzkripciós faktor kötőhelyek egy előre definiált csoportjában?
  2. Jelen van-e egy adott transzkripciós faktor motívuma a transzkripciós faktor kötőhelyek egy előre definiált csoportjában?

Arról, hogy mik is azok a motívumok, mit jelent maga a motívumfeldúsulás és milyen program segítségével lehet ezeket a feldúsulásokat azonosítani, a DNS szekvencia motívumok azonosítása I. bejegyzésben olvashattok, a következő néhány bekezdésben pedig arról lesz szó, hogy egy vizsgálni kívánt motívumról hogyan tudjuk eldönteni, hogy jelen van-e az általunk vizsgált régiókon belül.

Feltételezzük, hogy van egy 2000 transzkripciós faktor kötőhely pozícióit tartalmazó fájlunk (bed/txt kiterjesztésű) és szeretnénk csak azokat a kötőhelyeket, illetve a kötőhelyeken belül is csak azt a néhány bázispárnyi régiót visszanyerni, ahol megtalálható például az AP-1 fehérje motívuma. Ez a következő parancs begépelésével lehetséges:

annotatePeaks.pl peaks.bed hg19 -mbed output.bed -m AP1.motif -noann -nogene

Mi mit jelöl?

Az annotatePeaks.pl maga a program, amely a HOMER egyik nagyon hasznos eszköze. A peaks.bed a vizsgálni kívánt genomi régiókat tartalmazó bed fájl. A hg19 fogja a genomot megadni – amit a HOMER automatikusan felismer –, de ez a paraméter lehet mm10 (vagy mm9) is, amennyiben nem humán, hanem egér mintákkal dolgozunk. Minden egyéb modell organizmus használatakor meg kell adni a teljes elérési utat a vizsgálni kívánt genom FASTA fájljához. Eddig szinte minden ugyanúgy történik, ahogy a de novo motívumfeldúsulások keresésénél. Az -mbed paraméter után nevesítenünk kell egy bed fájlt, amely tartalmazni fogja a motívumtalálatok pozícióit (a fenti példában ennek az output.bed felel meg), az -m kapcsoló után pedig meg kell, hogy adjuk annak a motívumnak a mátrixát, amelynek jelenlétét vizsgálni szeretnénk az általunk megadott genomi pozíciókban. Végezetül egy kicsit gyorsíthatunk a motívumok keresésén a -noann és -nogene paraméterek megadásával; ezek használatával a vizsgált genomi pozíciók génekhez, ill. azok TSS-eihez történő annotálását a parancs nem fogja elvégezni.

Ezen felül az annotatePeaks-nél is működik és hasznos lehet a -size paraméter, mellyel a vizsgálni kívánt genomi régiók középpontjához viszonyítva megadhatjuk, milyen széles régión történjen a motívumkeresés.

Honnan szedjünk *.motif fájlt és mit érdemes a mátrixban változtatni?

Az előző bejegyzésben bemutattam, hogyan néz ki egy motívum mátrix és hogyan kell a benne található információkat értelmezni. Ezeket a mátrixokat a HOMER könyvtárunk homerResults vagy knownResults mappáiban találhatjuk, de akár készíthetünk újat, vagy paraméterezhetünk egy már meglévő mátrixot mi magunk is. Ezen felül létezik a HOMER-nek egy több, mint 400 ismert motívumot tartalmazó adatbázisa, melyet ide kattintva érhettek el: HomerMotifDB. A mátrixban a motívum score az, amit módosítani érdemes, annak függvényében, hogy mennyire szeretnénk szigorítani vagy lazítani a keresésen – bővebben erről is az előző bejegyzésben olvashattok.

No, de hogyan értelmezzük az eredményt?

Az output.bed kimeneti fájlunk 6 oszlopot fog tartalmazni. Az 1-3. oszlopok már nem az eredeti genomi pozíciókat fogják megadni, hanem pontosan azt a néhány bp-nyi régiót, ahol a keresett motívum megtalálható volt. A 4. oszlop a használt mátrix azonosítóját tartalmazza, amely a további munkálatok során nem releváns, az 5. oszlopban található score viszont annál inkább. A legalacsonyabb score legalább akkora lesz, mint a visszatérképezett mátrixban szereplő score-nál; annak lazításával a találatok száma növelhető.

A 6. oszlopban „+” vagy „–” jelöli, hogy a DNS-en pozitív vagy negatív irányban sikerült a motívumot azonosítani. Olyan fehérjék esetében, mint az AP-1, amely a TGAnTCA szekvenciához képes kötni, vagy a magreceptor szupercsalád bizonyos tagjai (például az ösztrogén recepor dimerek), melyek az AGGTCAnnnTGACCT szekvenciát preferálják, ha a reverz komplementerét vesszük a konszenzus szekvenciáiknak, mind a pozitív (+), mind a negatív (–) szálon olvasva ugyanazt a bázissorrendet láthatjuk. Ennek eredményeként előfordulhat, hogy a kimeneti fájl a 6. oszlopban eltérő irányultságot mutatva, de lényegében kétszer is tartalmazza ugyanazt a motívumot. A duplikátumok kiküszöbölése érdekében a kimeneti fájlunkat érdemes parancssorban merge-elni:

cat output.bed | sortBed | mergebed > output_v2.bed

 

Motívum score minden vizsgálandó genomi régióra? Lehetséges!

Az annotatePeaks-nek van még egy nagyon hasznos paramétere, mégpedig az -mscore. Ennek használatával a HOMER megkeresi a megadott mátrix által definiált motívumhoz legjobban hasonító szekvenciát minden egyes régióban, és kalkulál rájuk egy-egy motívum score-t. Ez az információ további szűréseket követően nagyon hasznos lehet abban az esetben, ha azt szeretnénk megvizsgálni, hogy egy adott motívum „erőssége” eltér-e különböző genomi régiók csoportjai között. Példaként, az alábbi ábra a TEAD, TCF, SIX, ERE, Fox és AP2 fehérjék motívumainak erősségét demonstrálják a „piros”, „lila” és a „kék” csoportok kötőhelyei alatt (Bojcsuk et al. bioRxiv, 2018):

boxes

Mivel az „erősebb”, tökéletesebb, vagy mondhatni kanonikus motívumok fehérje iránti affinitása sokkal nagyobb, az eltérések egyúttal utalhatnak a fehérjék kötésének meglétére vagy hiányára is.

 

A következő parancsot szükséges begépelnünk, hogy megkapjuk a motívum score-okat minden egyes kötőhelyre:

annotatePeaks.pl input.bed hg19 -m AP1.motif -mscore -noann -nogene -size 100 > output.txt

A kimeneti fájlban a következő oszlopok fognak szerepelni: PeakID, Chr, Start, End, Strand, Peak Score, Focus Ratio/Region Size, CpG%, GC%, Best Motif log-odds Score, melyből az utolsó oszlop lesz a meghatározott motívum score érték. Fontos, hogy az -m kapcsoló után nem csupán egy, hanem számos motívum mátrixát feltüntethetjük; például: -m AP1.motif AP2.motif ERE.motif TEAD.motif; ebben az esetben a kimeneti fájl utolsó oszlopai a mátrixok megadásának sorrendjében fogják a score-okat tartalmazni.

Az utóbbi parancs esetében mondhatni minden régióra „ráerőszakolunk” egy score-t, ezért előfordul az is, hogy negatív előjelű score-ral tér vissza az eredmény – emiatt is szükséges a további szűrés. Ha érdekel, én milyen feltételek alapján szűrtem a fenti ábra elkészítéséhez, olvassátok el a kéziratot! 🙂

EEM (1)

Az Emberi Erőforrások Minisztériuma ÚNKP-18-3-III-DE-253 kódszámú Új Nemzeti Kiválóság Programjának támogatásával készült.

DNS szekvencia motívumok azonosítása I.

In bioinfo, bioinformatics, bioinformatika, molecular biology on December 14, 2018 at 1:01 pm

Szerző: Bojcsuk Dóra

A két leggyakrabban feltett kérdés, amikor motívumanalízisről beszélünk:

  1. Milyen szekvencia motívumok „dúsulnak fel” a transzkripciós faktor kötőhelyek egy előre definiált csoportjában?
  2. Jelen van-e egy adott transzkripciós faktor motívuma a transzkripciós faktor kötőhelyek egy előre definiált csoportjában?

Ebben a bejegyzésben az első kérdés megválaszolására alkalmas bioinformatikai módszerről fogok írni – de előtte néhány gondolat arról, mik is azok a motívumok.

A motívumokról

A szekvencia motívumok viszonylag rövid, általában 6-20 bázispár hosszú, visszatérő mintázatok a DNS szekvenciájában, melyeket bizonyos DNS-kötő fehérjék, transzkripciós faktorok képesek felismerni. Ezekhez a motívumokhoz a transzkripciós faktorok sokkal nagyobb affinitással képesek kötődni, mint egy nem specifikus szakaszához a DNS-nek; ebből adódóan azonosításuk fontos a génszabályozás pontos megismerésének szempontjából.

A motívumok azonosítása régen és most

Elsőként 1975-ben David Pribnow azonosította azt a 6 nukleotid hosszúságú TATAAT motívumot (TATA-box), melyről kiderült, hogy mind eukariótákban, mind prokariótákban az egyik alapvető transzkripciós iniciációs helyet jelöli a gének promóter régiójában (10 bázispárra a kezdőponttól). Baktériumokban a TATA-boxon túl (35 bázispárra) megtalálható TTGACA motívum is fontos szereppel bír az RNS polimeráz enzim kiindulópontjainak kijelölésében.

A szabályozó szekvenciák azonosítása korábban az ún. footprint analízis (DNase footprinting) segítségével történt. Ennek során a tesztelni kívánt szekvenciát hordozó DNS darabokat radioaktív végjelöléssel látták el. Kontrollként egy olyan oldatot használtak, amely nem tartalmazta a vizsgálni kívánt fehérjét, csak a DNS-t, így minden, egy kiválasztott DNS-hasító enzim által létrehozott, jelölt fragmentum mérete láthatóvá vált gélen való futtatás és a gélkép előhívása után. A vizsgálni kívánt fehérjét tartalmazó oldatban viszont ott, ahol a DNS-fehérje interakció létrejött, az enzim nem volt képes hasítani, emiatt az érintett fragmentmérethez tartozó sáv nem volt látható a gélképen. A fehérje által kötött/megvédett szekvenciát nevezték footprint-nek (lábnyomnak). Végül az enzim által el nem hasított régiót a kontroll sávból visszanyerve DNS szekvenálással a fehérje által elfoglalt DNS szekvencia azonosíthatóvá vált.

Ma újgenerációs szekvenálási adatokból kiindulva és számítógépes programok segítségével sokkal egyszerűbb módon azonosíthatunk ismétlődő mintázatokat a DNS szekvenciájában – akár olyanokat is, amelyekről jelenleg nem is tudjuk, mely transzkripciós faktor(ok) kötheti(k).

A motívumok meghatározásához legtöbbször ChIP-seq, ATAC-seq vagy DNáz-seq adatokból indulunk ki, mert az ezekkel a kísérletekkel kapott csúcsok középső (~100-200 bázispárnyi) régiójáról feltételezhető, hogy a fehérjekötés középpontját, egyben a válaszadó elem helyét jelzik.

Feldúsult motívumok azonosítása

Az alábbiakban bemutatom, mely HOMER parancsok lehetnek segítségünkre a motívumok azonosításában és milyen paraméterekre érdemes figyelni.

Ahhoz, hogy megválaszoljuk, milyen szekvencia motívumok „dúsulnak fel” a transzkripciós faktor kötőhelyek egy előre definiált csoportjában, a HOMER findMotifsGenome.pl parancsát kell, hogy segítségül hívjuk. Ez alkalmas mind ismert, mind újonnan feldúsult, ún. de novo motívumok azonosítására is.

A használatához szükséges minimum paraméterek a következőek:

findMotifsGenome.pl peaks.bed hg19 output_dir -size 200 -len 8

Az általam használt egyéb paraméterekkel kibővítve:

findMotifsGenome.pl peaks.bed hg19 output_dir -len 8,10,12,14 -size 200 -dumpFasta -bits -preparse -homer2

Mi mit jelöl?

A findMotifsGenome.pl egy program a HOMER csomagból, amit meghívunk. A peaks.bed a vizsgálni kívánt genomi régiókat tartalmazó bed fájl. A hg19 fogja a genomot megadni – amit a HOMER automatikusan felismer –, de ez a paraméter lehet mm10 is, amennyiben nem humán, hanem egér mintákkal dolgozunk. Minden egyéb modell organizmus használatakor a vizsgálni kívánt genom FASTA fájljához meg kell adni a teljes elérési utat.

Az output_dir az eredmények helyét definiálja – ezt a könyvtárat előre létre kell, hogy hozzuk. A kimeneti mappában a homerResults könyvtár fogja tartalmazni a de novo találatokat, a knownResults könyvtár pedig azokat a feldúsult motívumokat, amelyek a HOMER adatbázisában (is) megtalálhatóak voltak.

A -len 8,10,12,14 kapcsoló segítségével mondhatjuk meg, milyen hosszúságú motívumokat keresünk. Bár a minimum motívum hossz a fenti példában 8 bázispár hosszúságú volt, ezzel a beállítással például a 6+1 bázispár hosszúságú AP-1 motívumot (TGAnTCA) is visszakaphatjuk (amennyiben fel van dúsulva).

A -size 200 paraméter segítségével definiálhatjuk, hogy a peaks.bed fájlban található genomi régiók középpontjához viszonyítva milyen széles régión szeretnénk a motívumokat azonosítani. A -size 200 paraméter esetében (amely egyébként az alapértelmezett beállításnak felel meg) a középponttól -100/+100 bázispáron belül eső régiót vesszük csak figyelembe, de a HOMER lehetőséget ad arra is, hogy a –/+ irányban eltérő hosszúságú régión keressünk. Például a -size 100,50 a középponttól -100/+50 bázispáron belül eső régiót veszi figyelembe, de kereshetünk akár a bed fájlban lévő teljes régiókon is a -size given paraméter megadásával.

A -dumpFasta kapcsoló használatával a HOMER kigyűjti két külön fájlba azon régiók szekvenciáit, melyeket a peaks.bed-ben megadtunk (target.fa), illetve az általa háttérszekvenciaként meghatározott régiókét is (background.fa). Ez abban az esetben hasznos, ha szeretnénk más programokat is bevonni a motívumanalízisbe; így ugyanazokhoz a háttérszekvenciákhoz hasonlíthatjuk a vizsgálni kívánt régióinkat.

Míg alapesetben méretarányos, ún. proporcionális motívum logókat rajzoltathatunk, ahol az A, T, C és G nukleotidok mérete annak megfelelően fog kirajzolódni minden pozícióban, hogy a motívumon belül mennyire volt gyakori az egyes nukleotidok előfordulása, a -bits kapcsoló használatával a nukleotidok ún. információtartalmával arányosan lesz súlyozott a karakterek magassága.

Proporcionális:

222

„Bits-es”:

          bits

 

Visszautalva a -size paraméterre, jól látható, hogy a fenti motívum 10 bázispár hosszúságú, de az AP-1 fehérje motívuma (TGAnTCA) ezzel a hosszal is szépen kirajzolódott, a többi pozícióban lévő nukleotidok pedig csak kisebb, kiegyenlítettebb valószínűséggel fordultak elő.

A -preparse használatával a HOMER minden keresés során új random szekvenciákat generál, a -homer2 pedig egyszerűen az új HOMER programot hívja meg a régi verzió helyett.

Hogyan értelmezzük az eredményt, mit jelent az, hogy motívumfeldúsulás?

A motívumkeresés eredményeit a kimeneti könyvtárban található homerResults.html fájl segítségével webböngészőben könnyen megjeleníthetjük, a talált motívumok logóit pedig a homerResults mappában kell keresni. Azt követően, hogy a homerResults mappát a html fájllal együtt letöltöttük a számítógépünkre, az eredményt ilyen formában láthatjuk:

2222

A táblázatból kiderül, hogy összesen 16188 genomi régióban kerestük a feldúsult motívumokat (Total target sequences = 16188), a HOMER pedig 33193 random genomi régió szekvenciáját használta kontrollként (Total background sequences = 33193). Ez azt jelenti, hogy a vizsgálni kívánt 16188 régióban feldúsult motívumokat a háttérként használt 33193 régióban is megkereste, mi pedig ebből már csak egy-egy százalékértéket látunk (% of Targets és % of Background), melyek azt jelölik, hogy a vizsgálni kívánt és a háttérként használt régiók hány százalékában fordult elő az adott motívum. Azonban könnyen félrevezethetjük magunkat, ha csupán azt vesszük figyelembe, hogy a vizsgálni kívánt régiók jelentős százalékában (pl. 32,45%) dúsult fel egy motívum, de a kontroll régiókra kapott %-ot (amely mutathat szintén jelentős, pl. 29,6%-os feldúsulást) figyelmen kívül hagyjuk. A HOMER minden találatra generál egy P-értéket (P-value) is, és ennek megfelelően rangsorolja a találatokat.

Szintén könnyen félrevezethetjük magunkat, ha elhisszük a HOMER-nek, hogy egy motívum valóban az, aminek ő nevezi. A fenti ábrán az 1. találat a BORIS transzkripciós faktor motívuma, azonban ha a More information hivatkozásra kattintunk, további lehetséges találatokról is tájékozódhatunk. Esetünkben az 1. motívumot valószínűleg csak az általánosan kifejeződő CTCF transzkripciós faktor tudja kötni, nem a paralógja – ennek megítélésére azonban valamilyen szinten ismernünk kell a modellrendszerünket; tudnunk kell, hogy mely fehérjék játszanak szerepet a vizsgált sejtben, ill. hogy adott fehérjecsaládból melyik fehérjék fejeződnek ki egyáltalán.

A további motívumtalálatokat is hasonló fenntartással kell, hogy kezeljük. A 2. találatot elegendő, ha C/EBP-nek nevezzük; az, hogy a fehérjecsalád mely tagja van jelen, szintén a vizsgált sejttípustól függ. A 3. találatot nevezhetjük AP-1-nek, mert az AP-1 fehérjecsoport tagjai képesek kötni, de nevezhetjük akár TRE-nek (TPA Reponse element) is, a TPA ligand válaszkészsége alapján. A 4. motívumot a promóter régiókra jellemző Sp1 fehérje képes kötni, de a motívum neve GC-box, amely pedig a szekvencia alapján kapta ezt a nevet, akárcsak a TATA-box, amit a bejegyzés elején említettem.

Bár az ábrán nem látszik, a piros csillaggal jelölt találatokat (*-possible false positive) az alacsony P-érték ellenére a HOMER lehetséges fals pozitívnak tekinti. Ha nagyon kevés régiót adunk meg, a HOMER szintén nem képes releváns feldúsulásokat eredményezni. Ebben az esetben azért a knownResults mappában található html fájlra is érdemes egy pillantást vetni, mert a HOMER visszatérképezi az adatbázisában található több mint 400 motívumot és a legnagyobb számban előforduló motívumokat kigyűjti (még darabszámot is megad), még ha a kevés kiindulási régió miatt azok nem is mutattak szignifikáns feldúsulást.

A fentiek helyes megítélése néha igényel egy kis kutakodást az interneten, de idővel könnyen rá lehet érezni. Fontos megjegyezni azt is, hogy az összes transzkripciós faktor felismerésére szolgáló motívum szekvenciája eltérhet egy-egy nukleotidban; ez a fehérje iránti affinitást nem feltétlenül, vagy csak kis mértékben befolyásolja.

Motívum mátrixok értelmezése

Minden motívum logóhoz tartozik egy motívum mátrix (ún. position weigth matrix) is, amiből készült, ill. amelyet a motif file (matrix) hivatkozásra kattintva nyithatunk meg, és az alábbi információkat tartalmazza:

>DRTTGCGHAA      3-DRTTGCGHAA,BestGuess:CEBPE/MA0837.1/Jaspar(0.925)  6.901336         -829.796349   0               T:1901.0(11.75%),B:1348.2(4.06%),P:1e-360

0.246      0.147      0.284      0.323
0.445      0.158      0.342      0.055
0.001      0.001      0.001      0.997
0.001      0.001      0.014      0.984
0.309      0.001      0.476      0.215
0.112      0.505      0.136      0.247
0.261      0.146      0.496      0.097
0.309      0.393      0.002      0.295
0.994      0.004      0.001      0.001
0.997      0.001      0.001      0.001

Rögtön az első sor a motívum nukleotidjait adja meg, melynek értelmezéséhez egy kis segítséget itt találhattok: https://www.bioinformatics.org/sms/iupac.html.

Szintén az első sor tartalmazza ugyanazokat az információkat, amelyek a táblázatos formában is szerepeltek. A motívum neve mellett sok esetben szerepel az is, hogy milyen korábbi ChIP-seq kísérletből származó motívumhoz hasonlít a legjobban. Szerepel még egy motívum score (érték) is (6.901336), amelyre még visszatérek, illetve itt is megjelennek a Target- és Background %-ok (T:1901.0(11.75%),B:1348.2(4.06%), illetve a P-érték (P:1e-360).

A bemutatott mátrix a fenti táblázatban szereplő C/EBP motívumot reprezentálja:

22

 

A mátrix minden sora a motívum egy bázisát írja le, balról jobbra haladva, az oszlopok pedig egy-egy nukleotidot jelölnek. Megfigyelhető, hogy az egy sorban szereplő számok összege 1-et tesz ki és azt írja le, hogy az adott pozícióban melyik nukleotid milyen valószínűséggel fordult elő. Ha megfigyeljük a vastagon kiemelt számokat a mátrixban, jól látható, hogy a 3. és a 4. pozícióban szereplő timinnek (T) a 4. oszlop felel meg, az utolsó két pozícióban szereplő adenint (A) pedig az 1. oszlop írja le. Végül pedig a 3. oszlop a guaninnak (G), a 2. oszlop pedig a citozinnak (C) felel meg:

   A             C             G            T
0.246      0.147      0.284      0.323
0.445      0.158      0.342      0.055
0.001      0.001      0.001      0.997
0.001      0.001      0.014      0.984
0.309      0.001      0.476      0.215
0.112      0.505      0.136      0.247
0.261      0.146      0.496      0.097
0.309      0.393      0.002      0.295
0.994      0.004      0.001      0.001
0.997      0.001      0.001      0.001

A mátrixban szereplő motívum score egyedi; függ a motívum hosszától és a mátrixban nagy gyakorisággal megjelenő bázisok számától is. Ezt az értéket akkor szoktuk módosítani, amikor adott genomi régiókon szeretnénk egy-egy motívum meglétét vizsgálni, de a keresés túl sok és nem specifikus találatot eredményezett. Ebben az esetben a score értéket megemeljük, így a mátrixban szereplő értékek bár ugyanazok maradnak, a keresés sokkal szigorúbb módon zajlik. Annak vizsgálatáról viszont, hogy jelen van-e egy adott transzkripciós faktor motívuma a transzkripciós faktor kötőhelyek egy előre definiált csoportja alatt, a következő bejegyzésben lesz szó bővebben.

EEM (1)

Az Emberi Erőforrások Minisztériuma ÚNKP-18-3-III-DE-253 kódszámú Új Nemzeti Kiválóság Programjának támogatásával készült.

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

Szerző: Bojcsuk Dóra

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.

 

 

 

 

 

 

 

 

 

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

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

Szerző: Bojcsuk Dóra

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

Szerző: Bojcsuk Dóra

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.

 

 

 

 

ChIP-seq_anal.sh

In bioinfo, bioinformatics, bioinformatika on June 9, 2018 at 4:08 pm

Szerző: Bojcsuk Dóra

Az újgenerációs szekvenálás (NGS) lényege, hogy több tízmillió (vagy akár százmillió), viszonylag rövid (50-200 nukleotid hosszú) DNS szekvenciát határozunk meg véletlenszerűen, és utólag rakjuk össze, mi rajzolódik ki belőlük. Ennek megvalósítására a mai napig fejlesztenek módszereket, így nem meglepő, hogy nincs egységes „pipeline” (csővezeték) az elemzésekre.

A funkcionális genomikai módszereket három nagyobb csoportra oszthatjuk az elemzés szempontjából: mRNS szekvenálás esetén az intronok kivágódásának köszönhetően a leolvasott szekvenciák többsége az egymást követő exonokra – akár három-négy exonra is – esik, így csak az exon-intron határok felismerésével történhet meg a szekvenciák referencia genomra illesztése. A kromatin 3D szerkezetének megismerését célzó módszerek esetében a meghatározott DNS fragmentumok – ideális esetben – többnyire két távoli genomi régióra esnek, amely e régiók fizikai közelségét jelenti, és amely egy teljesen másfajta technikai megközelítést igényel. A legegyszerűbb eset, amely a funkcionális genomikai módszerek többségére érvényes, hogy a meghatározott DNS szakaszok egy az egyben illeszkednek a referencia genom bizonyos szakaszaira. Ilyen módszerek például a kromatin nyitottságát vizsgáló DNase-, FAIRE-, Sono- és ATAC-seq, a különböző fehérjék ellenanyagát felhasználó ChIP-seq és az újonnan képződött RNS molekulák feltérképezésére szolgáló GRO-seq is. Az e módszerek alkalmazása során kapott szekvenciák feldolgozása tehát hasonlóképpen hajtható végre, és a szekvencia-feldúsulások jellegének megfelelően csak a végső lépésekben szükséges különbséget tenni az eredmények között.

A ChIP-seq_anal pipeline-unk (Barta, EMBnet.journal, 2011) az előbbieknek megfelelően nemcsak a kromatin csúcsszerű és szélesebb fehérjefeldúsulásait képes feltérképezni, hanem bármely olyan módszer szekvenciáinak a feldolgozására is alkalmas, amely a genomhoz viszonyítva folytonos leolvasásokat eredményez. Egy további nagy előnye, hogy ha megadjuk az ftp elérést, az NCBI újgenerációs szekvencia adatbázisából (SRA, Sequence Read Archive) közvetlenül letölti helyettünk a szekvenciafájlokat.

A legnagyobb előnye talán a teljes automatizáltság: ha létrehozunk egy-egy könyvtárat az ideiglenes (pl. FQ/fastq) és a feldolgozott fájloknak (pl. analysis), és megadunk egy listát arról, hogy mit és milyen néven szeretnénk elemezni, már indíthatjuk is a programot, minták számának függvényében pedig néhány nap alatt meg is kapjuk az eredményeket.

Hogyan kell kinéznie tehát egy ilyen listafájlnak?

Capture

A vi szövegszerkesztőt megnyitva, majd szerkesztési módba lépve, a fenti séma szerint kell elkészítenünk a listafájlunkat (samples.lst). A szövegfájl első oszlopa minden esetben tartalmazza a minta nevét, ahol fontos, hogy a genomok megkülönböztetése érdekében a mintanevek elején ’mm_’ ill. ’hs_’ karakterekkel kell a fajt megadni. Abban az esetben, ha a minta SRA adatbázisból történő letöltésére az elemzés során kerül sor, a beírt név a leendő mintanevet jelöli, ha viszont a FASTQ fájl már a szerverünkön, a fájljaink között található, az elemzésre váró minta pontos nevét kell feltüntetnünk. A második oszlop hivatkozik a minták SRA adatbázisban lévő elérhetőségére, az utolsó oszlop pedig a minta jellegétől függően felveheti az ‘input‘, ‘factor‘, ‘histone‘ vagy ‘groseq‘ címkéket. Ennek jelentőségéről részletesebben néhány bekezdéssel alább olvashattok. Fontos, hogy az oszlopok tartalmát egy-egy tabulátor választja el.

Ha párhuzamosan több listafájlt használunk, a fent említett futási idő a töredékére csökkenhet. Saját, vagy már letöltött adatok esetén nyilván nincs szükség az ftp cím megadására, csak a mintanevekre és a típusokra (ez esetben a listafájlban az ftp elérést egy plusz tabulátor helyettesíti). Ehhez azonban a szekvenciafájlt (FASTQ) vagy szimbolikus linkjét létre kell hozni a megfelelő könyvtárban (FQ/fastq).

 

A FASTQ szekvenciafájlok szekvenciánként 4 – tehát összesen n*4 – sorból állnak (ahol n a szekvenciák száma): egy ’@’ karakterrel kezdődő azonosító sorból, magából a szekvenciából, egy mindenképp ’+’ karakterrel kezdődő, opcionális „leírás” sorból és a szekvencia bázisonkénti minőségi mutatóinak sorából állnak. Az SRA-ban e hatalmas szövegfájlok tömörített verzióját, az ún. sra-lite fájlokat tárolják, amely a fastq-dump nevű programmal alakítható vissza szekvenciafájllá – természetesen ez is be van építve a pipeline-ba. UNIX parancssorban az ilyen mértetű fájlokat általában gzip-pel tömörítjük, amely zcat paranccsal gyorsan kibontható (a memóriába). A szekvenciák mm10 vagy hg19 referencia genomra illesztését a BWA programmal végezzük, amely végső soron egy olyan szövegfájlt hoz létre, amely mind a szekvencia adatokat, mind pedig a genomi koordinátákat tartalmazza (sequence/alignment map azaz SAM formátum), tárolásra ennek a binárisan tömörített (BAM) változatát használjuk, amely a SAMtools csomaggal hozható létre és indexelhető az adatok könnyebb elérése (pl. megjelenítése) érdekében.

Eddig a pontig nagyjából bármilyen szekvenciafájllal ugyanaz történik (kivéve persze az mRNS vagy 3D kromatin adatokat, amelyeknél már az illesztés is eltér), ezután lesz jelentősége a minta jellegének, amit a listafájlban kell definiálni. ChIP-seq esetén pl. a transzkripciós faktorokra jellemző a csúcsszerű feldúsulás az általuk kötött szabályozó hely körül, míg a különbözőképpen módosított hiszton fehérjék zöme szélesebb régiót fed le, és a szabályozó helyek pont kimaradnak, mert ott a hisztonokat más fehérjék váltják fel.

A polimeráz fehérje döntően géneket fed le, bár ez esetben is elkülöníthetőek csúcsok és szélesebb feldúsulások – a csúcsok a promóterek és enhanszerek környékén, az egyenletes lefedettség a gének testén jellemző. A kromatin nyitottságára specifikus módszerek is csúcsokat adnak, tehát hasonlóképpen kezelhetőek a faktor ChIP-seq adatokhoz. Az előbbiek szerint a ’factor’ és ’histone’ címkékkel lehet ellátni a mintákat, és az elemzés ennek megfelelően fog tovább folytatódni.

Ezek mellett kontroll mintákat is meg lehet adni – a listafájl(ok) elején ’control’ címkével, pl. aspecifikus IgG antitest esetében, vagy input minta esetében, amelyeknél nem várunk semmilyen feldúsulást, így alkalmasak a különböző minták hátterének és feldúsulásainak elkülönítésére. A GRO-seq, jellegénél fogva külön ’groseq’ címkét kapott, mivel szálspecifikusan az átíródó gének és enhanszerek határozhatóak meg vele, melyek pontos körülhatárolása és elkülönítése egy teljesen más kihívást jelent.

 

A szekvencia feldúsulások meghatározására a HOMER és MACS2 programcsomagokat használjuk. Ezek alkalmasak azoknak a régióknak a megtalálására, ahol a háttér lefedettségénél jelentősen több leolvasás található; ezek a régiók BED fájl formájában tárolódnak és tartalmazzák a csúcsok vagy szélesebb régiók genomi koordinátáit, de a teljes genom „domborzatát” is végigböngészhetjük ugyanakkor a HOMER és a MACS2 által létrehozott lefedettség (BEDGRAPH) fájlok segítségével. A HOMER alkalmas továbbá kiválasztott genomi régiók lefedettségének a kigyűjtésére – hisztogram vagy hőtérkép formájában –, valamint a csúcsok alatti szekvenciák feldúsulásainak a meghatározására és az így talált szekvenciamotívumok visszatérképezésére, nem mellesleg egyéb annotálási feladatok elvégzésére is…

…ezekről részleteiben viszont a további bejegyzésekben lesz szó.

 

Végezetül: milyen könyvtárakat kell definiálnunk a ChIP-seq pipeline indításakor és milyen sorrendben?

nohup sh <ChIP-seq_anal-v1_9mm10_fast.sh> <ChIP-seq/lists/samples.lst> <ChIP-seq/FQ/fastq> > <ChIP-seq/logs/samples.log> 2> <ChIP-seq/logs/samples.err> &

nohup: azért szükséges, hogy az elemzés a terminál bezárását követően is folytatódjon;

sh: parancs, amely jelzi, hogy egy shell szkriptet szeretnénk futtatni;

ChIP-seq_anal-v1_9mm10_fast.sh: ez maga a pipeline;

ChIP-seq/lists/samples.lst: a fent részletesen taglalt listafájl és elérési útja;

ChIP-seq/FQ/fastq: azt a mappát jelöli, amelyben a FASTQ fájlok megtalálhatóak, vagy ahová az SRA adatbázisból letöltést követően kerülni fognak;

a > jelöli a kimeneti fájlt, ebben az esetben egy log fájlt, amelynek kiterjesztése .log;

a 2> szintén egy kimeneti fájlt jelöl, viszont minden esetben hibafájlt, melynek kiterjesztése .err;

az & jellel a háttérbe küldjük a feladatot, így az éppen futó programok nem írják ki a képernyőnkre az aktuálisan futó vagy épp befejezett lépések automatikus üzeneteit, hanem azokat a *.log és/vagy *.err fájlokba gyűjti össze.

A fenti parancssort érdemes a ChIP-seq/analysis mappában indítani, így a pipeline a kimenetet az aktuális, analysis mappában fogja létrehozni.

 

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.

 

 

Mire jó a „cső” karakter?

In bioinfo, bioinformatics, bioinformatika on June 6, 2018 at 1:45 pm

Szerző: Bojcsuk Dóra

A „W” billentyű nem csupán egy betűt rejt, hanem egy másik, elég fontos karaktert is – ha az AltGr-rel együtt használjuk; a „|” („pipe” vagy „cső”) karakter ugyanis egy olyan, a UNIX-alapú operációs rendszerek egyik hasznos funkciójával rendelkezik, amely egyben rámutat e rendszerek kizárólagos előnyeire. A cső karakter egy parancs kimenetét egy újabb parancs bemenetévé alakítja, függetlenül attól, mekkora szövegtömbről is van szó (a UNIX rendszerekben nincs memóriakorlát, annyi memóriát használhatunk, amennyi csak van a számítógépünkben). Ezáltal egészen hosszú láncok, ún. „pipeline”-ok vagy „csővezetékek” építhetőek fel parancssorban, az adott feladatsornak megfelelően, és elkerülve a nagyszámú ideiglenes fájl létrehozását.

Ez jelenthet egészen egyszerű lépéseket is, pl. ha egy tömörített vagy átdolgozott szövegfájl első vagy utolsó sorait szeretnénk ellenőrizni, pl. zcat fajl.gz | head vagy sed ’s/ /\t/g’ fajl.txt | tail. Az előbbi példa egy kicsomagolt szövegfájl első 10 sorának parancssorban való megjelenítésére alkalmas, az utóbbi példa pedig egy szövegfájl „space” karaktereinek tabulátorokra való cseréje után az utolsó 10 sor megjelenítésére. Ezen a módon meglehetősen összetett dolgokat is ki lehet hozni az adatfájlokból: nemcsak karakterek cseréjére vagy bizonyos sorok kigyűjtésére alkalmas a parancsok egymás után fűzése, hanem gyakorlatilag bármilyen adattranszformáció vagy számítás megoldható ezen a módon. A csővezetékek (ill. a UNIX) további előnye, hogy az extrém nagy fájlokat is ugyanolyan hatékonyan kezelik, mint a kisebb fájlokat, és ciklus(ok)ba foglalva akár extrém sok azonos típusú fájlon is végrehajtható ugyanaz a feladat.

Emellett persze még ott van a párhuzamosítás lehetősége is, ha a fájlok csoportjait definiálva indítunk el párhuzamosan ciklusokat, amellyel megsokszorozhatjuk a feladat megoldásának sebességét, és tovább csökkenthetjük a várakozási és gépelési időt. Ha nemcsak bizonyos feltételeknek megfelelően strukturált fájlokat, hanem egységes fájlneveket és fájlrendszereket is használunk, tovább nő az automatizálhatóság lehetősége, amellyel nemcsak a gépelés mennyiségét, így pl. a gépelési hibákat küszöbölhetjük ki, hanem a munkaidőt tekintve is sokkal hatékonyabbak lehetünk. A parancsvégi „&” karakterrel háttérbe küldve a feladatokat új paramétereket tesztelhetünk, vagy új feladatokba is belekezdhetünk – nem lesz szükség a folyamatok felügyeletére, csak a végeredményt, a végső fájlokat kell ellenőriznünk.

Létezik ezen felül egy „nohup” (no hangup) parancs is, amelyet a csővezetékünk elé írva az internetkapcsolat bontásakor sem áll majd le a feladat; bár ez ciklusok esetében csak akkor alkalmazható, ha azokat Shell szkriptekbe építjük.

Így vagy úgy, de egy nagyon hasznos parancsról van szó, amely megengedi, hogy a szerver tovább dolgozzon akkor is, ha már nem vagyunk jelen. Ez különösen akkor fontos, ha napokig-hetekig fut valami, és nem rendelkezünk feladatütemezővel. Ami már egy másik történet…

 

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.

UNIX alapok

In bioinfo, bioinformatics, bioinformatika on June 2, 2018 at 1:36 pm

Szerző: Bojcsuk Dóra

A blognak bár nem célja programozásra tanítani az olvasót, úgy érzem, mégis ki kell térjek legalább egy alkalommal az általunk használt UNIX/LINUX (Shell) héjprogramozás alapszintű bemutatására és néhány alapvető parancs ismertetésére.

A UNIX (ejtsd: juniksz) egy közel 50 éve, 1969-ben létrehozott fejlesztői környezet, melyet Ken Thompson és Dennis Ritchei a Bell Laboratories-zal együttműködve, saját célra dolgozott ki. Bár sokak szerint furcsa dolog egy operációs rendszer nyelvét munkára bírni, az általa nyújtott lehetőségek pedig sok tekintetben elmaradhatnak az újabb, specializáltabb és magasabbszintű programnyelvek által nyújtott lehetőségektől – mint például a Python (Biopython), Perl vagy R programnyelvekétől, melyek 99%-ban lefedik a bioinformatikai programokat –, mégis nagyon hasznos eszközként kell tekintsünk rá. Nemcsak azért, mert könnyedén kiválthatjuk vele az Excel által kínált függvényeket, hanem azért is, mert ez egy szkriptnyelv, amely úgy működik, mintha csak szövegmanipulációra lett volna kitalálva. A bioinformatika pedig gyakorlatilag a nagy szövegfájlok (például szekvenciafájlok, táblázatok) feldolgozásának a tudománya. A UNIX által akár több millió sorból és/vagy oszlopból álló táblázatokat is hatékonyan kezelhetünk: alkalmazásával a feladatok végrehajtási ideje nagyságrendekkel lerövidülhet, és eltekinthetünk például attól a mindenki számára ismert problémától is, amikor a táblázatba foglalt értékeink dátummá alakulnak.

A UNIX fontos erénye, hogy kombinálható egyéb nyelvekkel. A SED (streamline editor) szövegszerkesztő és AWK „táblázatkezelő” például hamar a UNIX alapeszközeivé váltak, de egy R programnyelvben, kifejezetten bioinformatikai célra megírt elemző programot mint a DiffBind, is közvetlenül tudunk használni UNIX parancssorban (az Rscript parancs használatával). Egy UNIX platform egyszerre több felhasználónak képes munkafelületet biztosítani, és nagyszámú munkafolyamat egyidejű indítására is alkalmas. Az előző bejegyzésben bemutatott ngsdeb szerverünkön működő 236 (40+88+24+7×12) processzor lehetővé teszi azt is, hogy ezek a folyamatok párhuzamosan fussanak, más-más processzor igénybevételével, tovább gyorsítva így a felhasználók munkáját.

 A ngsdeb szerver több számítógép összekapcsolásából épül fel, melyeket node-oknak nevezünk. A bejelentkezést követően az ún. head (login) node-ra érkezünk, ahonnan munkánk megkezdése előtt – a házirend szabályai szerint – át kell lépjük valamelyik compute node-ra. Ehhez a következő karaktersort kell begépelni: ssh n001.

A korábbi, [username@admin ~]$ karaktersor helyett egy új sorban a [username@n001 ~]$ jelenik meg, amely továbbra is jelzi a felhasználónevünket, ill. hogy épp melyik node-ra jelentkeztünk át. A „~” karakter arra utal, hogy épp a /home könyvtárban tartózkodunk, a „$” (sorvég) karakter pedig arról árulkodik, hogy az előzőleg begépelt parancsunk befejeződött és újabb utasításokat adhatunk. Mivel a szervert használó munkacsoportok adatai más-más meghajtókon találhatóak, el kell jussunk a megfelelő könyvárig. No, de hogyan?

Azok számára, akik a terminál megnyitását követően már elveszettnek érzik magukat, nem tudják, hogyan tudnak lépegetni a mappák között, ill. hogyan tudnak fájlokat létrehozni, törölni vagy átmozgatni, összegyűjtöttem a legfontosabb UNIX parancsokat és azok magyarázatait:

pwd

pwd, azaz print working directory; kiírja a teljes elérési útját annak a könyvtárnak, ahol épp tartózkodunk; közvetlenül a bejelentkezést követően ez a /home/username könyvtárat fogja jelenteni

ls

kilistázza az adott mappa tartalmát (fájlokat és (akár al)könyvtárakat is egyaránt)

cd

cd, azaz change directory; begépelése után az adott könyvtárból elérhető alkönyvtár nevét megadva a kívánt könyvtárba léphetünk át (pl.: cd samples vagy cd samples/sample1 begépelésével egyből a samples könyvtárba, vagy a benne található sample1 könyvtárba juthatunk)

cd ..

visszatérés egy mappával „feljebb”; cd ../.. begépelésével pedig két mappával „feljebb”, stb.

mkdir

mkdir, azaz make directory;az mkdir all_samples begépelésével az aktuális könyvtárban létrehozhatunk egy új mappát all_samples néven

cp

cp, azaz copy paste; a cp sample1.bam ../. begépelésével a sample1.bam fájlt az eggyel „feljebb” található könyvtárba másolhatjuk át. A „../” után gépelt „.” karakter azt jelöli, hogy a megadott mappába szeretnénk másolni a fájlt. Amennyiben úgy szeretnénk áthelyezni, hogy új nevet is kapjon, a következő szerint kell eljárnunk: cp sample1.bam ../copy_of_sample1.bam

mv

mv, azaz move; a mv sample1.bam ../. begépelésével a sample1.bam fájlt egy könyvtárral „feljebb” mozgathatjuk, ugyanakkor a mv sample1.bam first_sample.bam paranccsal a sample1.bam fájlt az adott könyvtárban first_sample.bam-ra nevezhetjük át

rm

rm, azaz remove; rm sample1.bam begépelésével a sample1.bam törlésre kerül; egy könyvtár és annak teljes tartalmának törléséhez a -r kapcsolót is alkalmaznunk kell, pl.: rm -r samples

head

a head sample1.bed parancs kiírja a sample1.bed fájl első 10 sorát; head -25 sample1.bed pedig kiírja az első 25 sorát

tail

a tail sample1.bed parancs kiírja a sample1.bed fájl utolsó 10 sorát; tail -25 sample1.bed pedig kiírja az utolsó 25 sorát

cat

a cat sample1.bed kiírja a képernyőre a sample1.bed teljes tartalmát; ezt ebben a formában ritkán használjuk, inkább fájlok összefűzésére: cat sample1.bed sample2.bed sample3.bed > sample123.bed, amely a sample1-2-3.bed fájlok tartalmát egymás alá fűzi az általunk definiált sorrendben. A „>” karaktert követő fájlnevet tekinti a parancssor az újonnan létrehozott (vagy felülírt!) kimeneti fájl nevének

wc

wc, azaz word count; wc -l sample1.bed begépelésével a sample1.bed sorainak számát kapjuk meg

man

a man, azaz manual parancs használatával információt nyerhetünk arról, hogy egy adott parancs mire és milyen kapcsolókkal használható, pl. a man wc begépelésével megtudhatjuk, hogy a wc parancs a -w kapcsolóval együtt használva (man wc -w sample1.bed) kiírja, hogy a sample1.bed fájl hány szóból áll.

top

a top parancs kilistázza az aktuális node-on a felhasználók által elindított és épp futó job-okat

ps

a ps parancs begépelése információt nyújt csak az általunk elindított és futó job-okról az aktuális node-on

kill

a ps parancs által pl. „7698 pts/0”-ként kilistázott folyamatot, amely mondjuk a homer2 program futását jelöli, a kill 7698 pts/0 begépelésével a homer2 futását leállíthatjuk

exit

kilép az aktuális terminálból; ha valamelyik node-on tartózkodunk, akkor visszalép pl. a head node-ra

A fentiek alapján, ha szeretnénk megszámolni, hogy egy fájl első öt sora hány karakterből áll, megtehetjük azt, hogy a head -5 fájl.txt > fájl_v2.txt paranccsal létrehozunk egy új fájlt, amely csak a számunkra fontos 5 sort tartalmazza, majd a wc -m fájl_v2.txt paranccsal a képernyőre kiíratjuk az eredményt. Az átmeneti fájl (fájl_v2.txt) létrehozását azonban elkerülhetjük, munkánkat pedig meggyorsíthatjuk, ha a head és a wc parancsokat a megfelelő kapcsolókkal egy sorban alkalmazzuk: head -5 fájl_v2.txt | wc –m

Mit is csináltunk? Alkalmaztuk a „|”, azaz a „pipe” vagy „cső” karaktert, amely segítségével az első lépés kimenete lett a második bemenete.

Erről fog szólni a következő bejegyzés.

 

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.

Belépés az ngsdeb szerverre

In bioinfo, bioinformatics, bioinformatika on May 30, 2018 at 7:07 pm

Szerző: Bojcsuk Dóra

Munkánkhoz a Debreceni Egyetem Genomi Medicina és Bioinformatikai Szolgáltató Laboratóriumának ngsdeb szerverét használjuk, amely a közelmúltban jelentős bővítésen ment keresztül, és mint szuperszámítógép, 40+88+24+7×12 processzort, 3×128+256+6x48Gb memóriát és 40+6x11Tb merevlemezt tud magáénak.

Az új felhasználók számára a szerverre való belépéshez a rendszergazda biztosít egy bejelentkezési nevet és egy ahhoz tartozó jelszót, mely utóbbi az első belépést követően könnyen megváltoztatható. Bejelentkezés után egy UNIX (RedHot Linux) platformra érkezünk, amely UNIX operációs rendszerekről egy terminál megnyitásával közvetlenül elérhető, míg a Windows-t használóknak ugyanezt a UNIX felületet a PuTTY nevű, mindenki számára ingyenesen elérhető FTP kliens biztosítja.

A PuTTY letöltését és telepítését követően a szerverre való belépéshez a HostName mezőben a felhasználónevet és a szerver nevét kukaccal elválasztva kell megadni (username@ngsdeb.med.unideb.hu). A szerver titkosított hálózati kapcsolattal történő, távoli eléréséhez az SSH (Secure Shell) hálózati protokollt kell bejelölnünk, továbbá a Port beállítása is szükséges, melyet szintén a rendszergazda biztosít az új felhasználók számára. Amennyiben mindent begépeltünk, kattintsunk a megnyitás (Open) gombra.

1

A jelszavunk „vakon”, már a terminálablakba történő begépelése után egy Enter-rel tudjuk véglegesíteni a bejelentkezésünket, és máris a /home könyvtárban találjuk magunkat. Abban az esetben, ha a szervert nem az egyetem területéről próbáljuk elérni, az IP-címünket a rendszergazdának kell megadjuk ahhoz, hogy kivételként kezelhesse, és – akár otthonról is – bármikor hozzáférhessünk a szerverhez.

A Windows-t használóknak az asztali gépük és a szerver közötti fel- és letöltéshez a szintén ingyenes WinSCP FTP klienst is fel kell telepíteniük, amely szintén a kiszolgálót (ngsdeb.med.unideb.hu), a Port-ot, a felhasználónevünket és a jelszavunkat kéri a bejelentkezéshez.

2

Bejelentkezést követően a WinSCP két ablakra osztva jelenik meg; a bal oldali ablakban az asztali gépünk, jobb oldalon pedig a szerver mappái között lépegethetünk. Egyszerű áthúzással egy fájl/mappa, vagy többszörös kijelöléssel akár fájlok/mappák egyidejűleg is áthelyezhetőek.

3

Ha Windows operációs rendszert használunk, milyen programokat kell tehát telepítenünk ahhoz, hogy elkezdhessük a munkát?

PuTTY: https://www.chiark.greenend.org.uk/~sgtatham/putty/latest.html

WinSCP: https://winscp.net/eng/download.php

Mivel több felhasználó is dolgozhat egy időben a szerveren, használata mindenki számára azonos feltételekhez, szabályokhoz kötött.

A „házirend” a http://genbioinfo.med.unideb.hu/ honlapon érhető el.

 

 

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.

Bioinformatika – kell ez nekem?

In bioinfo, bioinformatics, bioinformatika on May 24, 2018 at 9:36 am

Szerző: Bojcsuk Dóra

Lassan húsz éve, hogy megjelentek az első közlemények az emberi genom teljes szekvenciájáról és a gének elhelyezkedéséről. E „térkép” megismeréséhez bár közel 15 évre volt szükség, nagyon fontos, mára már alapvető információt köszönhetünk a Humán Genom Projekt eredményeinek. Mindez magával hozta a DNS szekvenálás fejlődésének igényét, például az újgenerációs szekvenálási (NGS, Next-Generation Sequencing) módszerek megalkotását, finomítását, ezáltal a szekvenálás költséghatékonyabbá válását. A néhány száz genomi régióra irányuló chip-alapú technikákat felváltották a genomszintű vizsgálatok és specifikus antitestek használatával a fehérje-DNS kölcsönhatások genomi pozíciói is feltérképezhetővé váltak. Mindez hatalmas szekvenálási adathalmazt teremtett, tárolásuk pedig – jelenleg is – különféle adatbázisokban történik.

Manapság ha genomikáról, transzkriptomikáról vagy csak általánosságban molekuláris biológiáról beszélünk, már nem csupán a laboratóriumban végzett bonyolult kísérletekre gondolunk. A szekvenálási adatok önmagukban még nehezen értelmezhetetőek; az információ kinyeréséhez sok-sok lépcsőt kell megmásznunk, mely lépések bizonyos esetekben komolyabb informatikai és/vagy programozói ismereteket igényelnek.

A bioinformatikát mint kifejezést az 1970-es években Paulien Hogeweg és Ben Hesper holland biológusok alkalmazták először, ma pedig főként olyan folyamatok megnevezésére használjuk, mint a fehérjék 3D szerkezetének előrejelzése, biológiai hálózatok kutatása, vagy az említett szekvencia adatoknak az elemzése. Jelen van ugyanakkor az egészségügyben is, hiszen a képalkotó diagnosztikai módszerek, mint a röntgen, CT (komputertomográfia), MR (mágneses rezonancia) vagy a funkcionális MR képalkotás (fMRI) – amely egyenesen az agy működésének vizsgálatát teszi lehetővé – a fizika mellett túlnyomórészt az informatikára támaszkodnak. A bioinformatika interdiszciplináris tudományként tehát ötvözi a biológiát, informatikát, matematikát, fizikát, továbbá a statisztikát, így szakszerű alkalmazásához vagy a kísérletet végző biológusnak kell idejéből az informatika megismerésére áldozni, vagy az informatikusnak kell elmélyülnie a biológiában.

A genomi bioinformatikai módszereik alkalmazásában talán a legnagyobb nehézséget mégis az jelenti, hogy a szekvenálási adatok feldolgozásához nem létezik egységes protokoll. Csupán az alapelemzés az, ami automatizálható – minden további, specifikus kérdésre csak jól megtervezett és kiválasztott elemzési folyamat adhat pontos választ, melynek elengedhetetlen része a biológus és az informatikus közötti rendszeres és hatékony kommunikáció.

A kutatók körében kezd egyre inkább elfogadottá válni, hogy az adatbázisokban fellelhető nyers adatokat új információk reményében bárki felhasználhassa. A rohamosan növekvő NGS adatok – és sok esetben a redundáns kísérletek – lehetővé teszik, hogy a megfelelő mintagyűjtést követően akár saját laboratóriumi hozzájárulás nélkül – pusztán a bioinformatika alkalmazásával – is értékes funkcionális genomikai kutatásokat végezzünk. A bioinformatikai adatbázisokból történő információgyűjtés a biológusok munkáit is gyakran megalapozza; sőt, idejük legjava kezd átterelődni a számítógép elé. Nem kérdés tehát, hogy egy minimális bioinformatikai ismeret mindenképp a biológus javát szolgálja.

Magyarországon jelenleg leginkább a molekuláris biológusok és biotechnológusok nyerhetnek betekintést a bioinformatikába. Szakirányként választhatjuk a Szegedi Tudományegyetemen, de tantárgyként szerepel többek közt a Debreceni Egyetem molekuláris biológus, a Budapesti Műszaki és Gazdaságtudományi Egyetem villamosmérnök, ill. műszaki informatika, vagy az Eötvös Lóránd Tudományegyetem Természettudományi Karának képzésein is. Ugyanakkor, mint önálló szak, a felsőoktatásban még nem kapott helyet és az említett képzések tematikája is nagy varianciát mutat.

2017. január 1. óta Magyarország is tagja az ELIXIR (European Life-sciences Infrastructure for biological Information) programnak, melynek elsődleges célja a számítógépes biológiával foglalkozó európai kutatói közösség együttműködésének elősegítése. Támogatja régi és új adatbázisok, elemző módszerek fejlesztését, létrehozását, standardizálását, informatikai kurzusok szervezése révén pedig kezdők számára is lehetőséget biztosít egy-egy terület megismerésére.

Jelen sorozatban szeretném leginkább a biológiában jártas érdeklődők számára bemutatni, milyen NGS módszereket ismerünk, hogyan és milyen publikus adatbázisokból lehet szekvenálási adatokat letölteni, mik azok az alapvető lépések, amelyeket minden esetben alkalmaznunk kell a nyers szekvencia adatokon. Bemutatom az általánosan – és általunk is – használt programcsomagokat, majd specifikus elemző- és vizualizáló lépésekre is kitérek.

Remélem, bejegyzéseim sokak kedvét meghozza a bioinformatika tanulásához, de legrosszabb esetben is egy olyan alaptudást biztosít, melynek segítségével a biológus magabiztosan tud majd az informatikussal kommunikálni.

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: