Az élet kódjai

Nyuszikát és egyéb családi ünnepet hátrahagyva bepattantam Archee kocsijába és elhúztam Németországba az idei Revisionre. Az esemény a tervezési fázisában még jó ötletnek tűnt (kihagyni az értelmetlen zabálásokat, hallgatni, ahogy megvető megjegyzéseket tesznek az ajándékokra, vagy bármi másra, amit az ember tesz). De ahogy itt ülök és nem találom a helyemet, már nem olyan vicces.

Az utazás elég hosszúra sikerült. Mivel Archee-ék péntek reggel korán akartak indulni, én már csütörtökön Pestre utaztam, hogy le ne késsem az indulást. Kislányomnak elmagyaráztam, mi fog történni, amit ő így összegzett:

- El akar menni az apukád. - el az a hülye barom.

Utoljára még ágyban aludtam, mert az elkövetkező három napban csak a föld lesz a társam. Reggel Archee Kolumbiai barátnője finom rántottát készített, majd bepakoltunk és indultunk.

Mivel a lány viszonylag keveset tudott angolul, főleg spanyolul ment a társalgás, amibe én nem tudtam bekapcsolódni. A tíz órás út, sok bóbiskolás mellett hamar elröppent. Rendkívüli esemény nem történt.

Ez az első party, hogy nem hozok semmi release-t. A mini csapat egyik tagja sem készült semmivel.

A magyar kontingens egy helyre települt, de mikor megérkeztünk, már csak három szabad hely volt azon a részen. Egy ideig néztem, ahogy kétségbeesetten igyekeznek a három székből négyet varázsolni, majd úgy döntöttem, maga oldom meg a problémát. Ez nem is volt olyan egyszerű, mert minden csapat foglalta a haveroknak a helyet. Végül a Mercury mellett volt egy magányos hely, oda telepedtem.

A Revisionnek mindig van egy témája, ami köré csoportosul a party. Idén ez egy utóposztikus technokrata világ, ahol a demoscene az egyetlen lázadó erő a mindent és mindenkit uralni akaró gonosz multinacionális cégekkel szemben. Az ötlet egész jó, nekem már az invitációnál bejött.

Megtartották a Meteoriks díj-átadót, ahol az elmúlt év legjobb alkotásait díjazták. A Conspiracy hideglelős 64K-ja három díjat is nyert (közönségdíj, legjobb high-end alkotás, legjobb történet). Ez úton is gratulálunk nekik.

A lineáris regresszió igen hasznos és sokoldalú, de nem minden esetben van normális eloszlású mintánk, sőt nem minden esetben dolgozunk folytonos adatokkal.

Diszkrét adatok esetén a negatív binomiális regressziót használhatjuk. A név kicsit félrevezető, mert a függő változó igazából pozitív értékeket vesz fel.

Használata R-be nagyon hasonló a lineáris regresszióhoz. Nézzünk meg egy példát:

library(MASS)
model <- glm.nb(y ~ x, data = d)

 Vicces módon a lineáris regressziónál bemutatott példát itt is használhatjuk, mert a génekben található mutációk száma csak egész lehet. Most nincs szükség az adatok transzformálására sem.

genes <- read.table("genes.tsv")
genes.log <- log(genes+1)
nb.model <- glm.nb(snps ~ length, data = genes)
lin.model <- lm(snps ~ length, data = genes.log)

A glm.nb() fügvény futtatásánál kapunk néhány figyelmeztető üzenetet, de ezzel most nem törődünk. Valószínű az adat nem teljesen tiszta, vagy túl sok olyan gén van, amiben a mutációk száma nulla.

Tehát van egy lineáris és egy negatív binomiális modellünk. Melyik a jobb? Amelyiknél a reziduálisok négyzetösszege kisebb:

sum(lin.model$residuals**2)
[1] 34722.41
sum(nb.model$residuals**2)
[1] 28060.83

Mint látjuk, a negatív binomiális modell jobban illeszkedik az adatokra.

Mire használhatjuk még a negatív binomiális regressziót? Genomi töréspontok meghatározására. Ebben az esetben a töréspontok várható számát becsülték meg egy regressziós modellel, majd vizsgálták, mikor találnak az elvártnál magasabb értékeket. Hasonló logika húzódik meg a rákot okozó gének keresésében is. A különbség, hogy a várható értéket a szinoním mutációk alapján határozták meg, és nézték, hogy az aminosav cserével járó mutációk száma szignifikánsan eltér-e vagy sem.

Használhatjuk még RNA-seq-ben a változó expressziót mutató gének felderítésére, és így vissza is kanyarodunk az RNA-seq sorozatunkhoz. A kvantifikálás során ugyanis egy egész értékekkel feltöltött mátrixot kapunk. Az egy génre (vagy transzkriptre) jutó szekvenciák száma az összes mintára nézve, negatív binomiális eloszlást fog követni. A regresszió segítségével meghatározhatjuk a szórást, a szórás segítségével pedig eldönthetjük, hogy a változás, amit látunk csupán zaj, vagy valódi expressziós különbség. A helyzet természetesen bonyolultabb, mert a különböző minták szekvenálása során eltérő mennyiségű readet kapunk, de a mögötte álló logikán ez mit sem változtat.

A negatív binomiális regressziónak van egy gyenge pontja, amiről még nem beszéltünk. Ezek pedig a "túl sok nulla". Nem tudom egészen pontosan megfogalmazni, mit jelent a túl sok, de azt tudom, hogy emiatt ez a fajta modell nem használható egysejtes RNS szekvenálások feldolgozásához. Ilyen esetben a "zero-inflated negative binomial regression"-t használják, aminek nem tudom, mi a magyar megfelelője.

Akit részletesebben érdekel a téma, annak ezt a könyvet ajánlom. Nem egyszerű olvasmány, de gyakorlatilag a teljes koncepciót levezeti.

Néha egyszerűen semmi nem megy. Vannak olyan helyzetek, amikor a legjobb kikapcsolni a számítógépet és semmihez nem érni, ami kicsit is bonyolultabb egy arithmometernel. Ha a vállunkra ül a Bitmumus, semmit nem tehetünk, mással kell foglalkozni, hogy elunja magát és odébb álljon. A legrosszabb, ha még nagyobb energia befektetéssel dolgozni kezdünk.

Egy publikáció ábrájához kellett színskálát készítenem. A clusterProfiler remek csomag, sok szép képet lehet vele készíteni. Csupán egy baj van vele: semmilyen jelölést nem ad az ábrához. Nem tudni, hogy az enrichMapen a piros árnyalatok milyen p-értékekhez tartoznak, a körök mérete hány gént jelöl, stb.

A forráskód böngészése után arra gondoltam, elkészítem magam. Bitmumus ekkor mászott a vállamra és csendben figyelni kezdte, mit csinálok.

R-ben a színátmenetekért a colorRampPalette felelős. Segítségével létrehozhatunk tetszőleges színek között átmeneteket. Mivel gyorsan akartam eredményt elérni, ezért szokás szerint egy karakteres változó nevet használtam. Bitmumus a c-t javasolta és alig hallhatóan kuncogni kezdett.

Az R ezután furcsán kezdett viselkedni. Látszólag értelmes kódok hibaüzenetet kezdtek dobálni. Azt még érdemes hozzátenni, hogy Jupyter alatt tevékenykedtem. Miután láttam, hogy nem vagyok képes a színátmenettel düllőre jutni, elkezdtem mással foglalkozni, természetesen ugyan abban a munkafolyamatban. Ott is egyre-másra kaptam a furcsa hibaüzeneteket.

Kb. fél óra elteltével rájöttem, hogy semmi nem működik, amiben a c() függvényt használom. Bitmumus ekkor már a földön fetrengett a röhögéstől. Nekem még mindig nem esett le. Miért nem megy az egyik legalapvetőbb R függvény? Hozzá sem értem...Kivéve, mikor létrehoztam egy c változót. Biiip. Rossz válasz! A colorRampPalette függvényt ad vissza, nem változót. Sikerült felül definiálnom egy beépített függvényt!

Oké, legalább a hiba megvan. Már csak ki kell javítani. Kernel újra indítás? Bitmumus rekeszizmai begörcsöltek, csak kezével jelezte: csináld, te hülye. Az újra indítást követően még mindig nem kaptam vissza az eredeti függvényt.

Végül bezártam az egész rendszert, kinyírtam minden futó folyamatot. Még a böngészőt sem kíméltem. Visszakaptam az eredeti, értintetlen rendszert. Az órára néztem. Egy óra telt el, de még mindig nem volt színátmenetem. Szépen felálltam és kimentem a szabad levegőre. Bitmumus légszomjal küzdve, négykézláb mászva új áldozat után nézett.

Jessica Chang: Transforming gene discovery by radically open data sharing

A humán gének 51%-áról nem tudni, milyen hatással van a fenotípusra. Az újgenerációs szekvenálásoknak és a bioinformatika fejlődésének hála a felfedezések száma ugrasszerűen megnőtt. A következő 15 évben valószínűleg megugrik a feldolgozásra váró adatok mennyisége, de az elemzési, validálási kapacitás nem fog emelkedni, mert a kísérletek áteresztő képessége még mindig alacsony. Az adatok korai megosztása nem jellemző a kutatókra, mert mindenki magas impaktú lapba akar publikálni, ráadásul a páciens adatok nyilvánossá tétele etikai aggályokat is felvet. A megoldás a MyGene2 lenne, ahol a páciens oszthatja meg a saját vizsgálati adatait. Ez a Web2 alapokon nyugvó ötlet azért is előnyös, mert a kutatók a ritka betegségek kapcsán, ha megírták a publikációt, nem térnek vissza a betegekhez. Itt viszont a hasonló betegségben szenvedő páciensek kapcsolatban lehetnének egymással és információt kaphatnának az újabb terápiás lehetőségekről. Kórkép, fenotípus és a VCF feltölthető lesz, amit bárki nézegethet.

Konrad Karczenski: Prediction of splice variants and transcript-level effects improves the identification of gene-disrupting variants

A funkcionális régiókban a variációk száma lecsökken, de a szekvenálási hibák aránya ugyan az marad. Ha a variáció az intron kivágódás helyén keletkezik, az nem jelent feltétlenül a transzkript megszűnését, ha az intron tartalmaz GT vagy AG szakaszt az eredeti hely közelében, akkor a transzkript átíródhat. Ezen hibák felderítésére a Loftee VEP plug-in alkalmas. A Broad Institute Exac adatbázisa pedig összegyűjti ezeket. A potenciális kivágódási helyek megtalálására a MaxEstScan a legjobb program. A programból lesz új verzió is, amibe már integrált gépi tanulásos algoritmusok segítik a felhasználót.

Martin Kircher: Advancing massively parallel reporter assays for interpreting regulatory variants

A ritka variációk alacsony mintaszáma miatt nehéz szignifikáns eredményeket kihozni, ezért a megtalált variációk klinikai hatását a rcare nevű programmal lehet megbecsülni. A nem kódoló szakaszban előforduló variációk is fontosak, ezek felderítésére egy lenti vírus MPRA riporter assey-t kell használni. Ez beépül a genomba, de a beépülés helye nem lehet tetszőleges. A módszerük elég jól reprodukálható, de nem eléggé prediktív az eredmény, ezért a módszer további finomítására szorul.

Jennifer Harrow: Improving the speed of genome interpretation via network collaboration and knowledge sharing

Az Illumina felhő alapú genom elemző megoldása, a BaseSpace, több pipeline-t is tartalmaz. A szolgáltatáshoz tartozó AppStore-ban a felhasználó kiválasztja, milyen eszközökkel akar elemezni, majd fizet és megkapja az eredményeket. Nukleotid variációk keresésére még mindig a BWA+GATK a legjobb, struktúrális variációk felderítésére viszont a Delly-nél jobb a Manta. A megtalált variációk biológiai hatásának értelmezésére is van eszközük, ami állítólag 15 perc alatt lefut egy humán VCF-re. Ugyancsak van páciens fenotípus predikció is. A végső riport készítő alkalmazásuk pedig gyönyörű, publikáció kész összesítéseket bocsát a felhasználó rendelkezésére.

Jared Simpson: Analysis methods for nanopore sequencing data

Nanopore szekvenálás egy membrán két oldalán feszültség különbséget mér. A membránon keresztül vezet egy protein csatorna, amin ha átmegy a DNS, megváltoztatja a feszültséget. 1D-nek nevezik a kapott szekvenciát, ha egy szál halad csak át a csatornán, 2D, ha a reverz komplementer is átmegy. A szignál korrelál az átvándorló DNS mintázatával. Egy eseményből (frekvenia változás) 5 vagy 6 méretű k-mer szekvenciájára lehet következtetni (bár a homopolimer hiba magasabb 6-mernél). A hivatalos base caller a Metrichor, de a közösség létrehozta saját base callerjeit. A Nanocall egy rejtett Markov lánc alapú megoldás, csak úgy, mint a Metrichor, de a DeepNano már recurrent neurális hálózatot használ. A szekvenálás hatékonysága is sokat javult az új flowcelleknek hála. Az R7-nél még 81%, de az R9-nél már 92%, ami az új pólus fehérjének köszönhető. De novo összeszerelésnél hasznos eszköz a Nanopolish, ami egy rejtett Markov lánc alapú konszenzus készítő. Segítségével tovább csökkenthető a hibás nukleotidok száma. A szerkezet képes minden féle bonyolult könyvtárkészítési procedura nélkül metilált nukleotidok azonosítására. A base caller ebben az esetben kiterjesztett nukleotid készletet kap, hatékonysága 87% és a módszer hibája úgy tűnik kontextus függő.

 Kerstin Howe: Curation of multiple genome assemblies of the same species and utilisation for reference improvement

Az előadó a GRC referencia genomok felügyeletétől jött. A referencia ellenőrzésére a gEVAL, Gbench programokat használják, de vannak saját curátor eszközök is. Az információt több helyről kapják. Vannak kapilláris szekvenálási eredmények, alternatív módon összeillesztett referenciák (például CHM1, HuRef embernél, MGSCv3 egérnél). Kérdés, mennyire jók a szekvenciák, mennyire lehet javítani a genomokat? Embernél az elsődleges cél a chr11 javítása, amihez optikai illesztést (optical mapping) és  BAC klónokat, PacBio szekvenálásokat használnak. A kapott konszenzus szekvenciákat összehasonlítják. Meglepő módon még mindig vannak hiányzó és hiányos gének, mert a repetitív részeknél fragmentálódik a szekvencia. A kromoszómára nem illeszkedő részeket igyekeznek bevonni az összeszerelésbe, átrendezik a scaffoldokat, eltávolítják a lyukakat és hibás duplikátumokat. Munkájuknak hála 1300 problémát oldottak meg és 3%-al csökkentették az N-ek számát. Aki részletesebben is érdeklődik munkájuk iránt, nézze meg az oldalukat.

Frank Nothaft: Processing genomics data at scale with ADAM and Toil

A szekvenálás igen olcsóvá vált, de a számítógépek még nem tartanak ott, hogy lépést tartsanak az adat áradattal. A processzor mag szám növekszik, a I/O keresztmetszetét horizontálisan skálázódó eszközökkel igyekeznek növelni. A felhő sem jelent megoldást mindenre, mert a performancia sok esetben rosszabb, mint a HPC-k esetében. Egy lehetséges megoldás az Apache Spark. Ez a memória alapú adatbázis iteratív feladatokra a legjobb, erős Python/R/SQL támogatással. Hogyan lehet ezt genomi adat feldolgozásra használni? Először is más logika alapján kell szervezni a munkát, hogy jobban kihasználhassuk a párhuzamos feldolgozás adta előnyöket. A régi fájlformátumok (BAM/VCF) nem támogatják ezt a fajta feldolgozást, mert szekvenciális módon tárolják az információt, valamint az egymásra épülő feldolgozási lépések bizonyos feltételekkel hajtódnak csak végre (legyen a BAM sorba rendezve). Erre a Toil lehet a megoldás, ami egy munkafolyamat építő eszköz, Python alapokon. ADAM-al együtt használva 30-szor gyorsabb sebességet érhetünk el, miközben az eredmények nem változnak. Természetesen nem minden munkafolyamat ültethető át erre a logikára. Haplotípus azonosítás már nem skálázódik hasonló hatékonysággal.

Andrew Farrel: RUFUS: Reference free variant detection improves accuracy

A Rufus egy de-novo mutáció kereső program, ami illesztés nélkül találja meg az eltéréseket. Bemeneti adatai trió adatok. K-merekre bontja a readeket, hashbe gyűjti őket, majd de-novo assemblyt készít az egyedi szekvenciákból. Mivel csak a szekvenciákat nézi, minden élőlényen működik, nem szükséges paramétereket állítgatni. GATK-hoz nagyon hasonló eredményt ad ember esetén, de például a férfi X kromoszómán pontosabbak a találatok, mert a GATK feltételezi, hogy diploid szekvenciákkal dolgozik. https://github.com/jandrewrfarrell/RUFUS

Alicia Oshlack: SuperTranscript: a linear sequence assembled for the visualization and analysis of transcriptome data

Először a readeket összeszerelik a transzkriptekké, majd génekké klaszterezik őket a Corset nevű programmal. A kérdés, hogyan érdemes vizualizálni az egyes izoformákat? Ők úgynevezett szuper transzkripteket készítenek, ami minden lehetséges szekvenciát tartalmaz. Az így kapott mesterséges szekvencia akár differenciál expressziós vizsgálatokra is használható. Előnye, hogy könnyen áttekinthető az exon lefedettség, nem kell aggódni a több helyre illeszkedő readek miatt. A Lace nevű program képes a beadott klaszterekből szuper-transzkripteket építeni. Nem model organizmusokhoz is használható és akár variációk detektálására is használható. Az előadáson bemutatott egyik ROC görbe, ami az egyes transzkriptek expresszióját hasonlította össze, nem volt teljesen meggyőző, de egy másik példában olyan új csirke szekvenciát találtak, ami nincs benne a genomban. A hagyományos differenciál expresszióval való összehasonlítás eltérést mutat a két módszer között, de az előadó szerint a szuper-transzkript ennek ellenére is jobb módszer.

 Jouni Siren: GCSA2: A scalable approarch to indexing population variation graphs

Hogyan indexeljünk útvonalakat gráfokban? Ez egy kompromisszum az index mérete, részletessége és a keresés hatékonysága között. A gráfok számítógépes reprezentálásá folyamatosan fejlődött. Először deBruijn gráfokat használtak, később ezt felváltotta az FM index. A legújabb adatstruktúra a GCSA2. A gráfoknak lehetnek redundáns algráfjai, amiket ha összevonunk, kisebb memória méretet érhetünk el. Egy LF térképezésnek nevezett eljárás segítségével hosszú, egzakt egyezéseket derítenek fel a gráfokban. Az implementáció az 1000 genom variációit 7 óra alatt 387 GB-ra tömöríti és 4 mikroszekundum alatt megtalálja a keresett szekvenciát. Sebességben a BWA szintjén jár.

Birte Kehr: Unknown and non-repetitive sequence in the Icelanders genomes

Az izlandi emberek szekvenálása során olyan szekvenciákra bukkantak, amelyek nincsenek benne a referencia genomban, nem repetitívek és nem transzpozonok. Tizenötezer ember genomjából a nem illeszkedő readeket kiszedték, kiszűrték a mikrobiális genomokat, a maradékot Velvet segítségével összeszerelték. Csoportokba rendezték a kontigokat és variációkat kerestek bennük. A fenotípus még nem teljesen beazonosított, de struktúrális variációkkal találtak asszociációkat.

Katie Pollard: Decoding enhancer function with machine learning

A távoli enhanszerek beazonosítása, a szabályozott gének megtalálása és az enhanszerek mutációinak hatásvizsgálatára gépi tanulási módszereket vetettek be. Sok esetben a Chip-seq nem az enhanszereket azonosítja, a peakhez legközelebbi gén nem minden esetben a szabályozott gén. A szabályozás ennél komplexebb. Az EnhancerFinder képes elkülöníteni az aktív enhanszereket, egy másik eszköz, a MotifDiverge segítségével pedig megállapítható, hogy a transzkripciós faktorban azonosított mutáció a funkció elvesztésével jár-e. Chip-seq adatok mellett publikus Hi-C adatokat is feldolgoztak és azt vették észre, hogy a peakhez közeli gének sok esetben nem szabályozottak. A TargetFinder segítségével a komplex interakciók is azonosíthatóak. Azt a hipotézist vették alapul, hogy ha mutáció keletkezik a szabályozó régióban, akkor az a DNS alakjának megváltozásával jár, mert a transzkripciós faktorok térbeli alakokat ismernek fel. DNAshape program 5mereket rendel struktúrálist formákhoz. Egy hipergeometrikus teszttel megnézték, hogy ezekhez a strukturális formákhoz lehet-e szekvenciát rendelni. A peakek 25%-hoz sikeresen rendeltek összetartozó térbeli formát és szekvenciát. Az eredményeket CRISPR/Cas rendszeren ellenőrzik.

Chris Probert: DeepNuc: A deep learning model that accurately predicts genome-wide nucleosome positioning from ATAC-seq

ATAC-seq-el a nyitott kromatin régiókat lehet meghatározni. V-plotot használnak az eredmények megtekintéséhez, de a nukleoszóma kötőhelyek megtalálása így is nagy kihívás. Nagy mélységű neurális hálót tanítanak be, ahol az első réteg egy konvolúciós mátrix, majd 6 rejtett réteget használnak és különböző szűrőket használnak a különböző jellemzők azonosítására. A kimenet a MNáz szignál. Az eredmények összehasonlításából az derül ki, hogy bár a random forest módszerek is egészen jó eredményeket adnak, az előadó módszere még annál is jobb, split-ATAC-seq-el egyenesen tökéletes.

Michael Hoffman: Transcription factor expression and its effects on binding site occupancy and motif preference

A transzkripciós faktor kötése sejtvonal függő, ezért RNA-seq-et és Chip-seq-et csináltak ugyan abból a szövetből, hogy meghatározzák a kötés preferenciáját. Azt találták, hogy a peakek erőssége korrelál a gén expresszióval. Azt is észrevették, hogy különböző motivumok más-más szövetekben aktívak. A motivumokat klaszterezése után azt találták, hogy a motivumok sejttípus függőek, de biológiai validálást még nem végeztek.

Anshul Kundaje: Deep learning transcription factor binding sites and regulatory sequence grammars in diverse cell types and lineages

A gépi tanuláson alapuló módszerükhöz a tréning adata az 1kb hosszú genomi régiók, amelyek átfednek a transzkripciós faktor peakekkel. A negatív adatszett azokat a régiókat tartalmazta, ahol nincs átfedés. Ha elég sok konvolúciós réteget halmoznak egymásra, motivum kombinációkat is képes megtalálni a rendszer. Habár különböző bemeneti adatokkal tanították a rendszert (ATAC-seq, MNase-seq), a rendszer nem tűnik elég kiforrottnak. Például a Homerrel történő összehasonlítás során az eredmények inkább az utóbbi felé hajlottak. Valószínűleg a negatív adatszettel lehetett valami probléma. Link, link.

Vera Kaiser: Mutational biases drive elevated rates of substitution at regulatory sites across cancer types

A nem kódoló variációkat nem kutatják olyan intenzíven. A transzkripciós adatbázisok átnézése közben azt vették észre, hogy kicsi az átfedés közöttük. Khi-négyzet probával összehasonlították a rákhoz köthető gének szabályozó régióit egy kontrol régióval, nézték a mutációs spektrumot és vizsgálták a variációkat. Azt vették észre, hogy a CTCF kötőhelyeknél a mutációk aránya eltolódott.  Regressziós analízissel igyekeztek meghatározni, mely hatások a legjelentősebbek (pl. replikációs idő), de csak azt találták, hogy mindegyk fontos.

Kim Pruitt: NCBI graphical display tools – Hidden in plain sight

Az NCBI-nál egy új genom nézegetőt készítettek. Érdekessége a többi, hasonló kezdeményezéssel szemben, hogy ez tetszőleges weboldalba integrálható, majd saját annotációt lehet hozzárendelni. A terv az, hogy a GEO-ból is elérhető legyen. Bemutatták a Genome Workbench legújabb verzióját is, a tervek szerint ez fogja leváltani a SeqIn-t.

David Powell: RNA-seq visualization using Degust

A Degust, ahogy az előadás címe is mondja, egy kollaboráció központú RNA-seq vizualizációs program. Mivel az RNA-seq analízsre számtalan megoldás létezik, ezért ezeket építették be a Degustba és inkább a kutatók együttműködésére helyezték a hangsúlyt. Számtalan interaktív eszköz áll a kutatók rendelkezésére, hogy a differenciál expressziós eredményeket nézegessék. A bemenet egy nyers mátrix, amit webszerver alapú alkalmazás ezután feldolgoz. Magas szinten működik, de a hozzáértőknek lehetőségük van, hogy akár az ábrákat generáló R kódot is megtekintsék. A rendszert lehet az eredmények tárolására is használni. Érdekessége, hogy eredetileg Haskelben írták, de végül NodeJS-ben implementálták újra: https://degust.org

Jonathan Manning: ShinyNGS: Interactive data mining with R and Shiny

A Shiny egy webes keretrendszer R-hez, interaktív weboldalak készítéséhez. Ezen alapul a ShinyNGS, ami az analízis végén helyezkedik el és a már analizált eredmények megjelenítésével foglalkozik. Az előadásban példákat láttunk a használatára. A szerzők szeretnék kiterjeszteni egy sejtes RNA-seq-re is. https://github.com/pinin4fjords/shinyngs

Davis McCarthy: Using single-cell RNA-seq data for inference on gene regulation

Az egy sejtes RNA-seq népszerű, mert egyszerű dimenzió redukciós módszerekkel a sejt differenciálódás nyomon követhető. Az expresszió a differenciált sejteknél alacsonyabb a pszeudó idő skálán. A jó eredményekhez természetesen jó bemeneti adatok is kellenek. A scater csomag segítségével a QC könnyen megállapítható. A normalizálásra még nincs általánosan elfogadott módszer, de az scran elfogadható eredményt ad. A különböző batch effektek eltávolítására még mindig az ismétlésszám növelése a legjobb módszer. A hagyományos RNA-seq módszerek nem alkalmazhatóak.

Hagen Tilgner: Comprehensive transcriptome analysis using synthetic long read sequencing reveals molecular co-association and conservation of distant splicing events

A transzkriptek száma magas, hosszúságuk eltérő, de legtöbbször 400-700 bázispár hosszúak. lncRNS-eknek több izoformájuk van protein kódoló társaikénál, míg utóbbiak hosszabbak. Hosszú PacBio readekkel akarták azonosítani az egyes izoformákat. Sikerült új izoformákat azonosítaniuk, allélspecifikus variációkat is találtak, sőt egyes esetekben még új exononokat is felfedeztek. A kivágódások száma a polyA véghez közelebb megnő. Azt is észrevették, hogy összefüggés van a promóterek és az alternatív splicing között.

Simon Hardwick: Spliced synthetic genes as internal controls in RNA sequencing experiments

Az RNA-seq szekvenálások minőségéről az úgynevezett spike-in kontrol használata ad felvilágosítást. Ezeket a mintáktól jól el kell tudni különíteni. A szintetikus szekvenciák 78 gént tartalmaznak, 164 izoformát, egytől 36 exonig. A méret eloszlásuk 250 nukleotidtól 7kb-ig terjed. A hozzáadott kontrolok koncentrációja fontos, mert ha túl alacsony, akkor az FPKM és a tényleges koncentráció között nagyon nagy lesz a szórás. (Az előadó szerint 5.17 FPKM alatt a model nem használható.) www.sequin.xyz

Pali Meisted: Fast fusion detection, assembly, and quantification using kallisto

A fúziós géne fontos szerepet töltenek be a rák diagnosztikában. Az alap Kallisto k-merek segítségével azonosítja a transzkripteket és ha inkonzisztenciát talál a referenciával, eldobja azt. Ezért az alap program nem jó fúziós gének keresésére. A módosított program ellenben visszaadja azokat az eredményeket, amelyek nem fednek át ismert transzkriptekkel. A fals pozitívok száma magas, de a futási idő csupán 5 perccel növekedett meg. A fals pozitívok egy része a rossz annotációból ered.

A Sangerben általában hatékonyan mennek a dolgok. Természetesen itt sokan panaszkodnak, hogy a rendelési rendszer bonyolult, az állatházban fintorognak, ha egyedi kérés van, stb, de összességében szerintem tervszerűen haladnak.

De nem hiba nélkül. Ha pedig hibáznak, akkor az nagyot szól. Összességében kétféle hibát lehet megkülönböztetni. Az első kategória a "hagyományos" hiba, amikor valamit elterveznek, de nem úgy sikerül végrehajtani. A második kategória, amikor eltervezik, végrehajtják, de onnan nincs továbblépés, mert azt már nem tervezték el.

Az első bemutatandó példa Katerina esete. Közel 760 PDX mintát kellett volna feldolgoznia. A PDX (patient-derived xenograft) alapja, hogy a humán rákos sejteket egérbe oltják és ott növesztik. A szekvenálás során tehát a minta meghatározatlan mennyiségben egér DNS-t is tartalmazott. Ezt így tervezték, csak épp a bioinformatikusoknak nem szóltak előtte.

Mivel a pipeline-t nem hibridek feldolgozására tervezték, az IT nem nyújtott segítséget. Katerina talált egy R szkriptet, ami megoldja a problémát, de mivel csak neki volt rá szüksége, az IT nem telepítette azt a központi gépekre. Pont helpdeskes voltam, amikor felkeresett, hogy segítsek neki telepíteni a modult.

Biztos fantasztikus ötlet húzódik a program hátterében, de a megvalósítás hagy némi kívánni valót maga után. Először is, a program saját maga akarta elkészíteni a BAM indexeket és ezért a szimbolikus linkeket is visszakövette eredeti helyükre. Ami az iRODS szerveren van, halandók számára csak olvasási joggal.

Mikor erre rájöttem, módosítottam a modul kódját és onnantól már használható is volt.

A másik eset az én csoportomban történt. Közel 1400 mintát szekvenáltak meg és a DNS izolálás több tálcán történt. Kutatónk, hogy gyorsítsa a munkát, nem feliratozta a csöveket, hanem a tálca-pozíció alapján azonosította a mintákat. Az egyik lépésnél, amikor pipettázni kellett egyik csőből a másikba, a monoton munka hatására kimaradt egy cső, amit csak az egész folyamat végén vett észre.

Megvolt a három technikai ismétlés, de csak a Nedves Labor Nagy Szelleme tudta, hogy melyek azok. Az egyetlen dolog, amit tehettem, hogy mind az 1400 mintát hierarchikusan klasztereztem, hogy megtudjam, mik tartoznak össze. Ez nem volt olyan egyszerű, mert igencsak az elérhető memória (1,5 Tb) felső határán járt a szkript és ott futott több, mint négy napig. Kétszer is le kellett futtatnom, mert a fát tartalmazó PDF nem volt elég széles, hogy olvashatóak legyenek a feliratok.

Egyébként a legszaftosabb szekvenálási eredmények a Kísérletes Szekvenálási csoport produkálja. Ami tőlük jön, az olyan, mintha az állatorvosi lovat megfőznék a boszorkány konyhában. Csakhogy néhány példát említsek: 1 nukleotid hosszú adaptor, több száz sejtmagot tartalmazó sejtek (egysejtes) szekvenálása. Ezek egyike sem futhat a pipeline-ban. Gyakran a jól bevált programokat sem tudjuk használni. BWA például úgy nyelte el azokat a rövid adaptorokat, hogy öröm volt nézni. Aztán csak néznek nagy boci szemekkel: nem lehetne valamit mégis csinálni? Természetesen lehet, csak had indítsak el egy vim-et.

Ha az ember önéletrajzot készít, hajlamos gyakran összeegyeztethetetlen dolgokat is beleírni:

- Tud Ön egyedül dolgozni?
- Igen, nincs szükségem senkire.
- Ön csapatjátékos?
- Természetesen, akkor érzem igazán elememben magam, ha sok ember vesz körül és közös erővel dolgozunk.

De legutóbb arra jöttem rá, hogy nem vagyok csapatjátékos.

Mit csinál egy csapatjátékos? Egy jó csapatjátékos. Természetesen elfogadja, hogy ő egy kis fogaskerék, megkapja a maga kis inputját, hozzáteszi, amit tud, majd továbbadja a következő fogaskeréknek. Mindenki egy valamivel foglalkozik, nem akar senki több lenni a másiknál, egyformán mosolygunk és kis hangya szívünket melegség tölti el, ha arra gondolunk, hogy a nagy gépezet olajozottan működik.

Ez mind rendben is van. Semmi bajom vele. De vannak pillanatok, amikor a kis fogaskerék arra lesz figyelmes, hogy nincs input. Nincs min csámcsogni és nem lehet mit tovább adni. Természetesen a jó csapatjátékos ilyenkor hátradől és várja, hogy megjavuljanak a dolgok, esetleg egy kicsit zsörtölődik félhangosan, csak hogy lássák, saját hibáján kívül kénytelen Facebookozni.

Ez az, ami nekem nem megy.

Itt a Sangerben a DNS hosszú utat tesz meg az élőlényből a képernyőig. A kutató izolálja a DNS-t, átadja a Központnak, ahol elkészítik a könyvtárakat, betöltik a szekvenáló masinákba. A FASTQ fájlokat illesztik, az illesztést feltöltik egy fájlszerverre. A Rák genom program (illetve idéntől Rák, öregedés és szomatikus mutációk program) ügyintézői levelet kapnak, hogy új szekvenciák érhetők el. Létrehozzák az adatbázis bejegyzéseket (projektet, mintákat) egy weboldalon keresztül, majd a szkriptek működésbe lépnek és megindul az előfeldolgozás. Az előfeldolgozás során megnézik, hogy a szekvenálás milyen minőségű és az eredménytől függően az ügyintézők elfogadják vagy eldobják azokat.

A kutató végül kap egy levelet, hogy kész a projektje, majd a weboldalon keresztül kijelöli, hogy milyen vizsgálatokat kíván végrehajtani és a farm processzorai felizzanak. Ha bármilyen probléma adódna, a program bioinformatikusai (helpdesk) megoldják. Aztán már nincs is más hátra, mint táblázatokat generálni.

Az ügyintézők hárman vannak. Véletlen szerűen rendelik őket egy-egy szekvenáláshoz, de ha hozzárendelték őket, akkor az adott folyamattal kapcsolatos összes levelezést ők kapják. Mivel én tagja vagyok a genom program bioinformatikai csapatának is, nyomon tudom követni a szekvenciák útját, de nincs jogosultságom minden tevékenységhez.

Mi történik, ha például az egyik ügyintéző (vagy kettő) nyáron elmegy két hét szabadságra? Az egész rendszer megreked. Hiába van ott a harmadik, nem kap semmilyen értesítést a már futó projektekről. Az új igényeket megpróbálja kiszolgálni és ezzel el is megy minden ideje.

Én láttam, hogy kész minden szekvenálás, de hiába vártam, hogy elkészítsék a projektet. Végül felmentem az ügyintézőkhöz, de meglepetésemre az egész iroda sötét volt. Visszamentem a gépemhez és létrehoztam a projekteket. Ez nem volt nehéz, még dokumentáció is volt hozzájuk. Eltelt még egy nap és meg kellett volna nézni a minőségi mutatókat. Ehhez már nem volt jogosultságom.

Gondoltam egyet, elmentem az egyik principális bioinformatikushoz és megkérdeztem, kaphatok-e jogosultságot. Azt mondta, ha csak a saját projektjeimet birizgálom, kaphatok. Becs'szóra megígértem. Mire az ügyintéző visszajött, már kész volt minden. (Aztán egy másik alkalommal leszúrtak, mert létrehoztam egy projektet, de a pipeline nem tudta használni azt. Ugyanis nem tudtam, hogy bizonyos referencia szekvencia, analízis kombinációkat nem állíthatok be. Szépen elmagyarázták nekem, hogy a projekt létrehozás nem az én dolgom.)

A másik példa deviáns viselkedésemre és a szabályok figyelmen kívül hagyására Katerina esete. Majd részletesebben is írok róla, most csak annyit elöljáróban, hogy volt kb. 700 szekvenciája, ami vegyesen tartalmazott egér és ember readeket. A pipeline-t úgy tervezték, hogy csak egy referenciával dolgozzon, tehát szerencsétlennek mindent kézzel kellett volna csinálni. Talált valami béna szkriptet, ami a leírás alapján képes megoldani a problémáját, de képtelen volt működésre bírni. A bioinformatikus csoport (a nevünkben benne van, hogy támogató, de ez alkalommal nem támogatták) kétszer hajtotta el, mindenféle segítségnyújtás nélkül. Először azért, mert a rendszer nem kezel hibrideket, másodszor, mert a csoport által nem támogatott szoftvert akart használni. Pont helpdeskes voltam, amikor sokadszorra keresett valakit, aki megoldja a problémáit.

Nekem is egyszerű lett volna a szabályokra hivatkozva lepattintani (a feladatkörben benne van, hogy ha egy munka várhatóan 2 óránál hosszabb időt vesz igénybe, nem kell megcsinálni, hanem akkor a főnökség kirendel valakit a megoldására, ha prioritást élvez), mégsem tettem. A kollégák a szabályokra hivatkozva később is mondták, hogy ne segítsek neki, de mivel nekem nem főnökeim, nem kell hallgatnom rájuk.

Ha ezt jelenti csapatjátékosnak lenni, akkor köszönöm, nem kérek belőle.

A lineáris regresszió az egyik legegyszerűbb gépi tanulás módszer. Tudat alatt mi is sokszor használjuk, ezért a mögötte álló logika könnyen elsajátítható. Az alapkoncepció ugyanis az, ha egy mennyiséget növelünk, egy másik érték is növekedni fog. Egy hétköznapi példával élve, ha a kutyánknak több kaját adunk vacsorára, a súlya is növekedni fog. Itt az étel mennyisége és a kutya súlya áll kapcsolatban.

Ha tovább gondoljuk a példát, akkor érezzük, azért csupán a vacsora mennyisége, habár befolyásolja a testtömeget, nem elég, hogy pontosan meghatározzuk kutyánk tömegét. Tegyük fel, hogy rossz gazdik vagyunk és megengedjük állatunknak, hogy a szemétben is kotorásszon a séták idején. Kedvencünk súlyát tehát a szemétből kihalászott finomságok és a vacsora együttesen jobban meghatározzák, mint a vacsora egymaga.

Modellünknek immár két függő változója van (szemét + vacsi) és egy független változója (testsúly). Természetesen összetettebb modellel is próbálkozhatunk (megivott víz mennyisége, szőr között megbúvó kullancsok száma, stb.), de intuitíve érezzük, a súly meghatározása nem lesz sokkal pontosabb.

Ha eldicsekszünk másoknak is modellünkkel és ők is kipróbálják saját kutyájukon, csalódni fognak. Hamar ki fog derülni, hogy bár nekünk jól működik, más kutyák testsúlyát valami miatt rosszul határozza meg. A modellünk ugyanis nem elég általános, a szomszéd kutyája jónevelt, csak suttyomban eszik a kukák mellől, egy másik ismerősünk rendszeresen fut a kutyával, aki emiatt nehezebben hízik, stb. Jó lenne, ha ezeket a nehezen mérhető különbségeket is megjósolná a modellünk. Ha ezeket az adatokat is bevonjuk a modell generálás folyamatába, mi fog ezután történni? Valamelyik kutya súlyát pontosabban megjósolja, másokét kevésbé pontosan.

Úgy is mondhatnánk, egy idealizált kutya paramétereit kapjuk meg. Egyes kutyák jobban hasonlítanak majd erre a nem létező kutyára, mások (köztük természetesen a sajátunk is), kevésbé. Azt, hogy mennyire hasonlítanak az egyes kutyák erre az elméleti állatra, a modell által megjósolt testsúly és a ténylegesen mért érték különbségéből lehet kiszámolni (egészen pontosan a különbségek négyzetösszegéből).

Lássuk mindezt a gyakorlatban, bioinformatikához is köthető adatokkal. Töltsük le a Drosophila variációkat és genom annotációt. Készítsünk egy táblázatot, ami a gének nevét, hosszát és a bennük található SNP-k számát tartalmazza. Erre a következő szkriptet használhatjuk:

#!/usr/bin/python

import sys
import gzip
import re

class Gene:
    def __init__(self, start, stop, name):
        self.start = start
        self.stop  = stop
        self.name  = name
        self.count = 0

def findGene(genes, pos, starti, stopi, vnum):
    if genes[starti].start < pos and genes[starti].stop > pos:
        genes[starti].count += vnum
    if genes[stopi].start < pos and genes[stopi].stop > pos:
        genes[stopi].count += vnum
    if stopi - starti < 2:
        return
    midi = starti + (stopi - starti) / 2
    if pos > genes[starti].start and pos < genes[midi].stop:
        return findGene(genes, pos, starti, midi, vnum)
    if pos > genes[midi].start and pos < genes[stopi].stop:
        return findGene(genes, pos, midi, stopi, vnum)
    return

geneidpatt = re.compile('Name=([^;]+);')

genes = dict()

gff = gzip.open(sys.argv[1])
for line in gff:
    if line.startswith('#'):
        continue
    fields = line.rstrip().split()
    chrx   = fields[0]
    ftype  = fields[2]
    start  = int(fields[3])
    stop   = int(fields[4])
    match  = geneidpatt.search(fields[8])
    if ftype != 'gene' or not match:
        continue
    genename = match.group(1)
    if chrx not in genes:
        genes[chrx] = list()
    actual = Gene(start, stop, genename)
    genes[chrx].append(actual)
gff.close()

for c in genes:
    genes[c] = sorted(genes[c], key = lambda x:x.start)

vcf = gzip.open(sys.argv[2])
for line in vcf:
    fields = line.rstrip().split()
    if line.startswith('#'):
        continue
    chrx = fields[0]
    pos  = int(fields[1])
    vnum = len(fields[4].split(','))
    findGene(genes[chrx], pos, 0, len(genes[chrx]) - 1, vnum)
vcf.close()

print "length\tsnps"
for chrx in genes:
    for entry in genes[chrx]:
        print "%s\t%d\t%d" % (entry.name, entry.stop - entry.start, entry.count)

Töltsük be R-be a táblázatot és nézzük meg a változóink eloszlását:

data <- read.table("genes.tsv")
hist(data$length)
hist(data$snps)

Ez bizony baj. Kicsit sem hasonlítanak a normál eloszlásra, márpedig ez feltétele a lineáris regresszió használatának. Szerencsére egy egyszerű transzformáció segítségével közel normál eloszlásúvá alakíthatjuk őket. Ezután viszont tartsuk észben, hogy az összefüggésünk nem a gén hossza és a variációk száma között fog fennállni, hanem a transzformált értékek között!

data$length <- log(data$length + 1)
data$snps <- log(data$snps + 1)

R-ben a modell alkotás egyszerű folyamat. Meg kell adni a független változót (snps), jön egy furcsa karakter '~', majd felsoroljuk a függő változókat (jelen esetben ez csak a length). Az összes többi feladatot az R megoldja helyettünk. Érdemes alaposan összeismerkedni a modellekkel (az R dokumentáció formula néven emlegeti), mert teljesen új dimenziókat nyit meg előttünk.

model <- lm(snps ~ length, data = data)

Mit tudunk kezdeni ezzel a modellel? Jósolhatunk. Mennyi SNP-t várunk egy 23926 hosszú génben? Csak jusson eszünkbe, hogy a modellünk logaritmikus skálán működik!

a <- log(23926)
b <- predict(model, data.frame(length = a))

A válasz 7.100061. Ami 1212 mutációt jelent. Ha most megnézzük a 23926 hosszú géneket, akkor azt látjuk, hogy az stnB gén 1407 SNP-t tartalmaz, míg az stnA egyetlen egyet sem. Itt el is kezdhetünk gondolkodni, hogy mindez mit is jelent, de az túlmutat jelen poszt témáján.

Aki további részleteket akar megtudni a lineáris modellekről, annak ajánlom ezt a könyvet. Annyira szájbarágós, hogy még én is megértettem matematikai előképzettség nélkül. A Courseran is több kurzus értinti a témát. (ez, ez és még ez is, csak azokat említve, amiket láttam is.)

Itt a Sangerben hihetetlen figurákkal lehet összefutni. Az egyikük Tony, aki az EBI-nál a szekvenciák tárolásával foglalkozott. Előtte 15 évig a CERN-ben volt, hogy a részecske ütköztető által generált töménytelen adat tárolását megoldja. Elég jó rálátása van mindkét tudományterületre, a jelen poszt a vele folytatott beszélgetés kivonata.

A két terület közel azonos mennyiségű adatot tárol. Korábban a részecske fizikusok vezettek, de már a szerves vonalon is kezdünk annyi információt felhalmozni, hogy nincs okunk szégyenkezni. De még mi mindent központosítunk (NCBI, EBI, DDJB), addig ők elosztott adatbázist használnak. Nem félnek attól, hogy adat-paraziták menő cikkeket írnak mások keservesen összekuporgatott eredményeiből.

Érdekes módon a detektor által gyűjtött anyag nagy részét eldobják. Ez egy biológus számára felfoghatatlan, mert míg mi nem tudjuk, mikor lesz később szükség valamire, addig a fizikusok pontosan tudják, hogy mi az, ami soha az életben nem fog kelleni. Nem is tárolják. Pontos arányokra nem emlékszem, de szinte csak töredék információt tárolnak. (Ami így is elég nagy, hogy álmatlan éjszakákat okozzon Tony-nak és csapatának.)

Ez a tervezett hozzáállás egyébként az adatfeldolgozás minden lépésére jellemző. Mint megtudtam, a LHC építése előtt már megvoltak a fájl formátumok, programok, minden, ami a kiértékeléshez kell. A részecskefizikusok 15 éves formátumokkal dolgoznak, változtatás nélkül. Nálunk ez egy kicsit kaotikusabb: Juj, kéne valami, amiben tároljuk az illesztést. Kész a SAM. Ja, ez helypazarló, csináljunk BAM-ot. Most jut eszembe, hát ott a referencia, minek tároljunk redundáns információt. Nesze, itt a CRAM. De még a jó öreg FASTA is mennyit változott az idők során. És akkor a különböző kutatók által kidolgozott dialektusokról ne is beszéljünk.

Ugyan ez áll a programozási nyelvekre is. A fizikusok sokáig Fortranban írták a kódot és Tony állítása szerint egy kollektív döntés volt, hogy áttérnek C++-ra. Én már annak is örülök, ha egy bioinformatikai program csak nyelven íródott (IRAP,  Trinity).

Az alkalmazott módszerek viszont gyorsabban változnak nálunk. Amit ma napi szinten űzünk, nem biztos, hogy öt év múlva is csinálni fogjuk.

Még tavaly nyáron elmentünk az egyik angliai élményparkba. Elég sok ilyen van, közös jellemzőjük, hogy kiakakítanak egy egyéni stílust, majd telezsúfolják ugráló várakkal és szurikátákkal. De komolyan, néha már az az érzésem, hogy ezek az állatok innen származnak.

Ez a látogatás nem sikerült zökkenőmentesen. Kicsit felületesen néztem meg az odavezető utat, aminek az lett a következménye, hogy a sokadik busz, amivel aznap mentünk, kitett minket egy város széli megállóban és én nem tudtam, merre menjünk tovább.

De nem voltunk egyedül. Egy másik, népesebb indiai család is oda igyekezett, amit a sofőrrel folytatott beszélgetésükből tudtam meg. Beszélgetni kezdtünk és együtt tettük tettük meg az a két kilómétert, ami még hátra volt. A végére már járda sem volt, csak az autóúton bandukoltunk, libasorban. Ők még babakocsit is vittek. Biztos remekül mutattunk, de az autósok remekül tolerálták kis karavánunkat.

Kiderült, hogy a családfő a Thompson Reutersnek dolgozik, valami rákkutatással kapcsolatos részlegben. Hiába, itt Cambridge környékén köpni sem lehet, hogy el ne találnánk egy molekuláris biológust. Minden áron beszélni akart velem a bioinformatikai elemzésekről.

Megbeszéltünk egy találkozót, amolyan munkaebédet, ahol megtudtam, ők készítik a Mindentcsinálót. A program nevére nem emlékszem, csak arra, hogy differenciál expressziót számol, anyagcsere utakat elemez, asszociációkat határoz meg és potenciális rákterapeutikai célpontokat is visszaad. Természetesen, ahogy frissen szerzett barátom elmondta, rajtam kívül mindenki ezt használja.

Az eszköz egy weboldalon keresztül üzemel, de van hozzá REST API, ami a hozzám hasonló vén csontoknak szánt csali, hogy ráharapjanak. A mögötte megbúvó adatbázist manuálisan kurálják és naponta jelenik meg hozzá frissítés. Tehát kétszer nem lehet úgy elvégezni egy elemzést, hogy ugyan azt az eredményt kapjunk.

A legjobb dolog persze hátravolt: mindezt a kísérletet végző kutató tudja futtatni, bioinformatikus segítsége nélkül.

- Mit gondolsz, mit tenne a főnököd, ha ilyen szoftvere lenne? - kérdezte tőlem a fickó.
- Kirúgna engem - mondtam rövid gondolkodás után. Erre a válaszra nem számított.

Azt tudni kell minden üzletkötőről, hogy csak pozitív üzenetet szeretnek sugározni, amit a hasonló ügyfelekkel történő beszélgetések során jól be is gyakorolnak. Ha ez a jól bevált séma léket kap, összezavarodnak. Most is ez történt. Elkezdte mondani, hogy akkor én mást csinálhatok. Megkérdeztem, ugyan mit? Ettől még jobban összezavarodott. Végül megsajnáltam és azt mondtam neki, csak vicceltem. Ekkor megnyugodott és visszatért a beszélgetés az általa irányított mederbe.

De a helyzet az, hogy tényleg nem tudom, mi lenne a feladatom egy Mindentcsináló mellett. Talán az adatok Mindentcsináló-kompatíbilis formára alakítása? Vagy csak olyan elemzéseket végzek, amire a program nincs felkészítve?

A fickó azzal próbált még győzködni, hogy a bioinformatikai elemzés szokott lenni a szűk keresztmetszet. Ezzel a kijelentéssel vitatkoztam. Én úgy tapasztaltam eddig, hogy nem az elemzés volt a publikációk gátja, hanem a szekvenálás sebessége. Itt a Sangerben kb. 4-6 hónapig terjed, mire a DNS mintákból szekvencia adatok lesznek. Nekem utána egy hét elég, hogy az elsődleges eredményeket tálaljam (ha szabvány elemzés kell. Egyéni kísérleti elrendezésnél persze minden sokkal tovább tart.).

Természetesen az ábrakészítés lassan megy, mert a minták sorrendje rossz, a színek nem illeszkednek a publikáció színvilágába, vagy egyeszerűen csak, miért használok logaritmikus skálát? Ezen egyetlen program sem fog segíteni.

Ami még lassan szokott menni, az az, amit én csak képernyő bámulásnak hívok. Tudjátok, amikor ott van előtted minden információ és fogalmad sincs, hogy ez mit is jelent. Nem nyomsz le egyetlen gombot sem. Az egér mozdulatlanul pihen a kezed alatt. Az agyad meg max fordulatszámon pörög (vagy rossz esetben sokkhatás alatt áll).

Ilyenkor megy a szöszmötölés: ennek a génnek furcsa neve van, mit csinál? Miért nem találom a kutató kedvenc génjét a listában? Az adatok két harmada miért viselkedik ellentétesen, mint ahogy várom? A képernyő bámulást még a Mindentcsináló sem tudja felgyorsítani.

Búcsúzóul megkért, hogy a főnökömnen feltétlenül meséljek a Mindentcsinálóról én meg azt hazudtam, hogy természetesen. Eljátszottam a gondolattal, hogy írok neki egy levelet: A főnök úgy döntött, megtart engem és nem használja az Önök programját. Indoklás: a fizetésem kevesebb, mint a szoftver havi díja.

A klaszter adminisztráció nem egyszerű feladat. Emlékszem, még az egyetem alatt az egyik tanár mesélte az új labor felállításáról, hogy a telepítést úgy oldották meg,  hogy egy gépre elvégezték, majd a dd paranccsal a többi gép merevlemezére másolták azt. Nyilván, ha egyforma gépeink vannak, azokat ugyan olyan feladatot látnak el, ez működik is.

A mai klaszterek nem homogének. Változatos konfigurációban állnak rendelkezésre és van közöttük webszerver, tárolásra kialakított gép és számításokat végző egységek. A helyzetet bonyolítja, hogy ezek a szerepek az idő folyamán változhatnak.

Miként biztosítsuk a flexibilitást és a könnyű kezelhetőséget? A Sangernél az OpenStack használatával. Minden egységen egy virtuális gép fut, és ezen gépekre távolról rakják fel a képfájlokat. Ha valamelyik gépnek új funkció kell, csak a megfelelő képfájlt kell ráküldeni és már kész is. Gyakorlatilag ezt a módszert alkalmazzák az Amazonnál is. Csakhogy nekik nem kell törődni a szoftver környezettel, mert az a felhasználó dolga.

Itt viszont ha frissíteni kell egy programot, akkor új képfájlt kell létrehozni, ami lehetőleg hibamentes (tehát tesztelni kell), lehetőleg automatikusan készüljün el (ne kelljen a nulláról installálni mindent, csak mert megjelent egy új verzió az awk-ból) és az sem árt, ha a képfájl elkészítésének folyamatát egy verzió követő rendszerben megőrizzük.

Természetesen mindegyikre van megoldás. Jöhet a szakzsargon zsonglőrködés.

A verziókövető rendszer a Git, gondolom nem kell bemutatni. A stabil verzió a fő ágon van, minden újítást egy másik ágon kezdenek. Ha elkészültek a módosítással, letesztelték azt, mennek vissza a fő ágra. Akár csak egy szoftver fejlesztési projektnél. De a Git alapvetően kódot (szöveges állományt) tárol. Milyen kód kerüljön bele?

A virtuális képfájlokat a Packerrel készítik. Ez egy konfigurációs állományt állítja elő kódból. Támogat számtalan virtuális környezetet, tehát ha valami miatt az OpenStack nem válik be és egy másik technológiát kell követni, a Packer kód változatlan maradhat.

Az elkészült képfájl fordítását ellenőrizni kell, mert több száz program telepítése során igen sok hiba fordulhat elő. Erre a GitLab CI-t (continuous integration)  használják. Ha a CI nem talál hibát, a képfájl elkészülhet.

Ezután a kész képfájlt kell ellenőrizni, hogy a számtalan szoftverkomponens nem ütközik-e. Itt bizony nincs mese, a képfájlt telepíteni kell egy tesztgépre és különböző tesztprogramokat kell futtatni. A Test Kitchen célja pontosan ez. Csupán egyetlen konfigurációs fájlt kell készíteni és a rendszer lefuttatja a benne definiált teszteket.

Magukat a teszteseteket nem a Kitchen definiálja, az csupán egy keretrendszer. A tényleges tesztek Bats-ban futnak. A szoftver komponenseken kívül a szerver állapotát is tesztelni kell (felcsatolta-e a meghajtókat, müködik-e a hálózati csatoló, stb.), amire ServerSpec-t használnak.

Tehát ha új szoftver kell a klaszterre, először frissítik a Packer képet, a régi, stabil verzió felépítése megtalálható a Git-en. Megírják a teszteket, lefuttatják azokat, majd küldik ki a képfájlokat a csomópontokra. Igen, kicsit összetettebb, mint a dd.

Kijött az idei első cikk! Az elemzés egyik lépéséről már készült egy korábbi poszt. A vizsgált régióra még a tesztelésre megkapott egyik MinION-t is ráengedtük, de nem volt elég nagy a lefedettség, hogy bármi érdemi eredményt kapjunk.

Amilyen gyors volt a bioinformatikai elemzés, a kísérletes validálás annál lassabban ment, de ez nem is meglepő, mert elég sok fajtán kellett letesztelni az eredményeket és nincs mindenhol pipettázó robot. A mélység tehát nem csak a szekvenálásban, de a kísérletes oldalon is megvolt.

Sok pletykálni való nincs, a vizsgálat szempontjából érdekes variációk nagyon hamar a feltűntek, de akkor még nem tudtuk, mit is jelent ez. Mivel a kódoló régióban vártuk az érdekes mutációkat, meglepett minket, hogy a promóterben volt. Minden esetre a későbbi vizsgálatok alátámasztották a mutáció jelentőségét.

Ezt a projektet is nagyon szerettem, mert előtte még soha nem dolgoztam poolozott mintákkal. A témát alaposan körbejártuk (vagy inkább járták, mert abban már nem vettem részt) és a cikk is szép kerek egészet alkot.

Hétfő óta nincs víz a lakásban. Ha megnyitom a csapot, csak valami furcsa hang jön ki rajta: MagyarországcsatlakozottazElixirhezMagyarországcsatlakozottazElixirhez...

Előző évben egy buszos beszélgetés alkalmával kiderült, hogy az Elixir egyik adminisztrátorával szoktam utazni. Mikor megtudta, hogy honnan származom, mindig csepegtetett egy kis információt az események aktuális állásáról. Mivel a csatlakozás folyamatában személyesen nem vett részt, túl sok érdemi dolgot nem tudott említeni, de mondta, hogy épp melyik stádiumban van a kérelem. De pont ebben a hónapban nem sikerült beszélnünk, ezért én is más csatornán keresztül értesültem a sikeres csatlakozásról.

Korábban egy félnapos Elixir bemutatót is szerveztek itt, a Campuson, amire nem mentem el, de akik ott voltak, csak annyit tudtak mondani róla, a kaja jó volt. A legtöbb leírás elég semmitmondó, tele olyan marketing szagú szövegekkel, mint biztonság, innováció, cloud, fenntartható fejlődés, stb.

Oké, mihez is csatlakoztunk? Az Elixir egy elosztott infrastruktúra lesz. Míg korábban mindent integrálni akartak egy kézbe, mint az EBI, most a feladatokat szétosztják az Elixir tagjai között. Magyarország szerepe a fehérje szekvenciával és felépítéssel kapcsolatos eszközök, adatbázisok fejlesztése lesz.

Ezen kívül az Elixir szárnyai alatt szerveznek képzéseket, biztosítanak számítógépes kapacitást, és segítik a kollaboráció kiépítését.

De azért nem véletlen, hogy kevés konkrétumot hallani. Ha megnézzük, a terveken kívül nem sok kézzel fogható dolgot látni. Most még hurrá fázisban vagyunk, mint a gyerekek a játszótéren, akik összegyűjtik a haverokat a nagy homokvár építéséhez. Remélem kitart a lelkesedés akkor is, ha piszkos lesz a kezünk.

Az illesztést követően meg kell határozni, hogy a vizsgált régiókra mennyi read esik. A vizsgált régió lehet gén, transzkript, de akár exon is. Annak ellenére, hogy ez első látásra egy egyszerű feladat, hamar érdekes kérdésekre kell választ adnunk. A régiók ugyanis átfedhetnek. Antiszensz transzkriptek esetén az átfedő régióból keletkező readek megfeleltetése nem egyszerű feladat. Az itt bemutatott programok végtelenül egyszerűen oldják meg ez a problémát: eldobják a readeket. Ezért is mondtam, hiába van szuper illesztőprogramunk, ha utána nem vesszük figyelembe a bonyolult esetekből származó eredményeket.

De ha egy génről több transzkript képződik, a közösen használt exonok miatt akkor is nehéz megmondani, melyik read honnan is származik. (És akkor még nem is beszéltünk alternatív UTR-ekről, intron megtartásról és az alternatív splicing többi fajtájáról.)

HTSeq

Ez a szkript egyszerű, mint a faék. Minden primitívsége ellenére sokáig egyeduralkodó volt. Nem csinál semmit, csak összeszámolja, hogy génünk, transzkriptünk mennyi readet tartalmaz és eldob minden olyan esetet, ahol a read származási helye nem egyértelmű. Annyira egyszerű, hogy sokáig még pozíció alapján sorba rendezett SAM fájlokkal sem boldogult (de ha tehetjük, ne használjuk ezt a funkciót most se). A dokumentáció szerint olvas BAM fájlokat, de a gyakorlatban én még ezzel nem találkoztam.

samtools view input.bam  | htseq-count -q -m intersection-nonempty --stranded=yes --idattr=gene_id -r pos  - reference.gtf >count_table.tsv

A fenti kód az esetek 80%-ban működik. Ha a "Maximum alignment buffer size exceeded while pairing SAM alignments" hibaüzenet jelenik meg, kénytelenek leszünk a BAM fájlokat név szerint rendezni. A másik rejtélyes esetem a programmal, amikor az egyik bioinformatikusunk optimalizáció végett kromoszómánként akarta elkészíteni a végső táblázatot. Az eredmény eltért attól, amit akkor kaptunk, ha egyben futtattuk a teljes genomon. Egy ideig nyomoztuk az esetet, de végül félbe maradt az ügy. Ha mindezek nem tántorítottak el senkit a használatától, megemlítem, hogy még lassú is.

FeatureCount

Ez a program méltó utódja az előbbinek. Gyors, hála az optimalizált kódnak, a változatos paramétereknek köszönhetően kellően flexibilis és képes az összes BAM fájlt beolvasva a teljes mátrix legyártására (bár ez sokkal több időbe kerül, mintha BAM-onként futtatnánk és utólag fűzzük egybe az eredményt.) Egyetlen dolog nem tetszik csak: rengeteg járulékos, felesleges információ van a kimenetben. A következő opciókkal a HTSeq-hez elég hasonló eredményt kapunk.

featureCounts -a ../annotation.gtf -t gene -f -s 1 -o count.tsv input.bam

A sorozatot kicsit felfüggesztem, és egy látszólag teljesen független témával fogom folytatni, de akár csak a sok szereplős, a szálak között ide-oda ugráló romantikus filmeknél, az összefüggésekre végül fény derül.

A cím kicsit félrevezető. A fájl nélküli (filess) nem azt jelenti, hogy nincsenek fájlok. A szerver nélküli (serverless) esetén se gondoljunk arra, hogy a szolgáltatások a semmin futnak. A két koncepció csupán annyit jelent, hogy a fájlok és szerverek, mint elsődleges célok megszűnnek és helyüket átadják egy magasabb absztakciós szintnek. Mik is ezek a magasabb szintek?

Fájltalan

Erről a koncepcióról először a Broad Institute egyik előadásán hallottam, ahol az intézet infrastruktúrális hátteréről beszéltek. Náluk nincs hatalmas számítógép park, amit lehet, kereskedelmi felhő szolgáltatóknál oldanak meg. Ott tárolják a szekvenciákat, a metainformációkat, mindent. Amolyan jenki mentalitással, csapatmunkában dolgoznak. Tehát ha a Google-nek/Amazonnak van számítási kapacitása, akkor vesződjenek ők a biztonsági mentésekkel, rendszer frissítéssel és egyéb adminisztrációs gondokkal, a Broad meg vesződik a szekvenciák feldolgozásával.

Amiben más a koncepciójuk, az a szekvenciák tárolása. Egy BAM fájl ugyanis már nem elég az összes információ tárolására. Rákos mintáknál például vannak páciens adatok, a betegség előzményei de akár a könyvtár előkészítés, minta begyűjtés is fontos lehet. Ezeket a metainformációkat is tárolni kell. Egy hagyományos fájl alapú megközelítésben készítenek még egy fájlt (legrosszabb esetben egy Excel fájlt) ezek tárolására. De még ha adatbázist használnak is a tárolásra, akkor is a hagyományos fájl műveletek (másolás, átnevezés) könnyen inkonzisztens állapotba hozhatják azt.

Elsőre talán furcsa lehet ezt hallani, de ha belegondolunk, a mindennapi élet is hasonló kihívásokat tartogat. Vegyük csak a legegyszerűbb esetet, a fényképek tárolását. Normális esetben a fényképeket dátum és tematika szerint rendezzük egy könyvtár struktúrába. De bizonyára mindenki találkozott olyan esettel, hogy egy fénykép két csoportba is belefér. Mit csináljunk? Másoljuk két helyre? Hozzunk létre linket? Bármit választunk is, a hagyományos fájl alapú megközelítést használunk, annak minden korlátaival.

Szervertelen

Az Amazontól is tartott előadást egy ember. Ő az Amazon Lambda szolgáltatásról beszélt, ami - kicsit eltúlzott formában - a jó öreg callback függvények felhős megfelelője. Röviden arról van szó, hogy ahelyett, hogy szervereket konfigurálnánk, kódot írunk, majd beállítjuk, hogy a kód milyen esemény hatására fusson le (trigger).

Az előadáson a példa a fájl feltöltés volt. Adva van egy weboldal (természetesen ez is az Amazonon futott), ahol a felhasználó feltöltött egy képet. A kód elhelyezte azt az Amazon tárhelyére. Ezután a példát kiterjesztették egy mobilos alkalmazásra, ahol a kód változtatása nélkül, csupán egy új trigger hozzáadásával, immár mobilos felhő szolgáltatássá változtatták a programot.

Természetesen az egész csak az Amazon berkein belül, azok szolgáltatásaival összhangban működik, saját szerverrel mindez elképzelhetetlen. De hát ezért szervertelen, nem?

Az idén sem lehet ok a panaszra. Megjelent öt cikk, ami nem is rossz. Feldolgoztam 980 egysejtes RNA-seq-t (egy PhD hallgatónak kellett segíteni a dissszertációjához), 40 organoid RNA-seq-t, 1085 sejtvonalból származó RNA-seq-t (az RNA-seq pipeline teszteléséhez kellett), 196 egysejtes fitymabőr adatot (de csak 37 volt használható), 1224 exome szekvenálást (bár ez még folyamatban van), 4 humán teljes genomot, 9 egér teljes genomot és 6 microarray-t. Oh, az 550 GTEx adatról majdnem megfeledkeztem, bár az is folyamatban van.

Azokat az eseteket nem számoltam, ahol nem a teljes feldolgozást, csak egy részfeladatot csináltam, vagy segítettem másoknak kitalálni, miért akadtak el a programok az adatokon.

Összehasonításképp, a Sanger előtt egész életemben 67 szekvenálási adatot dolgoztam fel (gyarló memóriámnak tudjatok be +- 20-t). Szintén nem számolva azokat az eseteket, ahol csak részfeladatot hajtottam végre, vagy kapilláris szekvenálással kapcsolatos munka volt.

Demoscene fronton mind a megjelent posztok számában, mind az aktivitásban visszaesés volt tapasztalható. Csupán egyetlen produkció született, azt is eléggé lehúzta a közösség. Nem, mintha olyan sok pozitív kritika született volna más alkalommal, de legalább újabb partival bővült a gyűjtemény.

Nem csak a szekvenálási adatok száma növekedett, de a számítógépes technológiák száma is. Számomra kezd nehezen követhetővé válni, melyik mire is jó. Már két olyan előadáson is voltam, ahol csak kulcsszavakat jegyeztem le, mert az előadó úgy dobálózott velük, mintha ezek nyilvánvalóak lennének mindenki számára. Ez elég ijesztő volt számomra.

A másik döbbenetes dolog számomra az okos készülékek fejlődése. Én nem vagyok egy nagy mobil guru (a 8 éves Nokiám a példa rá), de az egyik Linux magazinban a személyi asszisztenseket hasonlították össze. Annak idején, mikor megjelent a Siri, az egyik Alma fan lehozta Amiga klubba. Az egész akkor egy viccre hasonlított. Alig értette, mit akartunk, lassú volt, béna. Most viszont az volt a kifogásuk a tesztelőknek, hogy egyes asszisztensek fárasztó vicceket mesélnek. Mi van? Komolyan ez egy probléma?

A trendektől kezdek eltávolodni, de azért remélem a tatter bloggal maradtok jövőre is.

Vége a hosszú, fáradságos munkának. Kész a kézirat is. A munkatársak, akik eddig együtt, vállvetve kutatták az ismeretlent, hirtelen ellenségekké vállnak és megkezdődik a harc, hogy milyen sorrendben legyenek a szerzők nevei feltűntetve.

De előtte még felmerül a kérdés: egyáltalán kit vegyünk be?

Az első kézenfekvő válasz: aki dolgozott az adott munkában. De mindenki tudja, hogy ez nem igaz. Az asszisztenseket nem vesszük bele, pedig ők is részt vettek a kísérletekben. A szerzősor végén pedig fel-feltűnnek olyan emberek, akik effektíve nem dolgoztak, de mégis nélkülözhetetlenek. Például mert adták a pénzt a világmegváltó vizsgálatokhoz vagy meghümmögték az eredményeket.

Akkor használjuk a következő definícót: akik nélkül nem jöhetett volna létre a cikk. Ez sem jó. A kutatók szülei sem szerepelnek a szerzőlistában (bár némelyik Nature cikknek annyi szerzője van, hogy a nagy számok törvénye alapján ez bekövetkezhet), sem a Microsoft Word fejlesztői. Valami pontosabb kéne.

Definiáljunk egy mérőszámot. Egy hányadost, ahol a számláló a hozzájárulásunk százaléka a munkához, a nevező pedig azoknak a cikkeknek a száma, amihez a fenti munka felhasználható. Ha ez a hányados átlép egy küszöbértéket, az illetőt bevesszük a cikkbe.

Tekintsünk el olyan apróságoktól, mint a "hozzájárulás százalékának" mérhetetlensége. Ez sem lenne teljesen objektív és a küszöbérték meghatározása sem olyan egyszerű. De tegyük félre ezeket és nézzük meg a mérőszám működését két szélsőséges esetben.

Tehát a Word fejlesztői mondjuk a munka 0.01%-t végezték, de ez a munka elméletileg az összes cikkhez hozzájárulás, hiszen majdnem mindegyik kéziratot Wordben írják. Tehát egy kicsi számot elosztunk egy nagy számmal: az eredmény a nullához közelít.

A PhD hallgató viszont a munka mondjuk 80%-t végzi, közben mással nem foglalkozik (hahaha). Ő átlépi a küszöböt.

Első körben a fenti mérőszám jónak tűnik, kellően objektív is. De ennek nincs jelentősége, mert a valóság az, hogy a levezető szerző teljesen önkényesen dönti el, kik legyenek bent.

Például nézzük ezt a cikket. Ebbe akár bele is vehettek volna. A munka alapjául szolgáló pályázatból fizettek, én csináltam a miRNS elemzést (ami valószínűleg újra megcsináltak) és később is rendszeresen küldtem a N. benthamiana homológ szekvenciákat nekik, mert a Blast keresés túl bonyolult volt. De eljöttem, elfelejtettek. Van egy olyan érzésem, nem ez lesz az egyetlen ilyen cikk.

Aztán itt van ez a cikk. Jóformán nem csináltam semmit, a munkát azután végezték, hogy eljöttem, csak egyfajta telefonos segítség voltam. Egyáltalán nem sértődtem volna meg, ha nem vesznek be. Mégis bevettek, mert úgy érezték a mérőszámom átlépi azt a bizonyos küszöbértéket.

Az egyik csoportvezetőnek YouTube nézés közben támadt az az ötlete, hogy karácsonyi dalokat algoritmizáljunk. Amolyan börtönös karatefilm szabályokat állított fel ("Az egyetlen szabály, hogy nincs szabály.") ésvárta a pályaműveket.

A folyamatábra nekem túl snassz volt, arra gondoltam, egy szándékosan túlbonyolított, nehezen olvasható kód működését kellene bemutatni, mintha egy hibakeresőben futna.

Elöször egy karácsonyi énekre volt szükségem, ami nagyon-nagyon redundáns. Meg is találtam. A második lépés egy Python kód megírása volt. Először megírtam szépen, hogy legyen egy referencia, majd fokozatosan elkezdtem a kódot tömöríteni és ezzel együtt nehezen értelmezhetővé tenni. De még így is inkrementálisan állította elő a szöveget. A végső verzióban viszont ezt is sikerült megváltoztatni, ami kicsit elrontotta a hibakereső-szerű bemutatót. (A refrénben ugyanis először elkészítem az ismétlődő részeket és utána szúrom be a nem ismétlődőket.)

A harmadik lépés demó-szerű lett, mert írtam még egy programot (C++), ami a zenére szikronizálva OpenGL-ben mutatja be a Python kód működését. Vizuálisan nem lett egy nagy eresztés, de nagyon jól szórakoztam az elkészítésével. Olyasmi volt, mintha Legóznék.

Végül videót is készítettem, mert a legtöbben almás PC-t használnak az intézetben, biztosan nem fogják lefordítani a kódot. Kipróbáltam az NVidia ShadowPlay felvevőjét, amit eredetileg a játékmenet rögzítésére találtak ki. Az elindítása nem volt teljesen triviális, mert nem tudtam, mire kell kattintani, ha aktiválni akarom a felvételi funkciót.

Legnagyobb meglepetésemre csak feketeséget vett fel, és csak akkor rögzítette a demót, amikor beállítottam a Desktop recording-ot. A végeredmény nem lett túl rossz, a tömörítési aránnyal is elégedett voltam. A textúrák kicsit kockásak lettek, talán további finomhangolással a minőséget is lehet javítani.

A nagy megmérettetés viszont elmaradt. A csoportvezető nem válaszolt az e-mailre és nem volt semmilyen visszajelzés. Egy hét elteltével is csak annyit válaszolt, hogy türelem. (Ellenálltam a késztetésnek, hogy azt írjam: türelmes vagyok, mint egy helpdesk kérés.) Többeket kérdeztem, egyáltalán tudnak-e valakiről, aki adott be valamit, de nemmel válaszoltak.

Remélem, azért Nektek tetszik:

Ha pedig valakit a kód érdekel, itt megtalálja.

Mikor elkezdtem a PhD-t, témavezetőm adott egy hatalmas füzetet és elmagyarázta, hogyan kell jegyzőkönyvet vezetni: minden nap új oldalt kezdeni, dátum a jobb felső sarokban, oldalszámozás a jobb alsóban és az előző munkafázis oldalszámát feltűntetni. Spirálfüzetet elfelejteni, mert abból csak kitépődnek a lapok.

A módszerrel a legnagyobb problémám a kereshetőség hiánya volt, a másik a limitált hely, amit kitölthetek. Legtöbb esetben elég volt egy oldal, de néha nagyon sokat kellett írni és ha volt még két-három gélkép (amit a megfelelő ragasztóval kellett beragasztani), akkor átlógtam a következő oldalra és borult a rendszer, vagy csak helyet pazaroltam.

Ezért kezdtem el számítógépes jegyzőkönyvet vezetni. Közönséges HTML oldalak voltak, a fenti logikát követték, oldalszámozás helyett linkekkel lehetett manőverezni. A formázási lehetőségek limitáltak voltak, eszközök alig akadtak (emlékszik még valaki a Dreamweaverre? Úgy értem, mikor még a Macromedia tulajdonában volt.).

Aztán ahogy átnyergeltem bioinformatikára, minden munkafolyamatot számítógéppel kezdtem végezni, és jött a LaTeX. A formázási problémák megoldódtak, a dokumentumok szép ábrákkal bővültek, de a kódok beillesztése rendszeres probléma volt. Ha csak simán beraktam, a LaTeX fordító értelmezni akarta. Ha nyers szövegként illesztettem be, túlnyúltak a sorok a lapszélen. Nehéz volt megosztani a munkarátsakkal, az e-mailes küldözgetés nem volt jó megoldás.

Később jött a Wikipedia/Confluence. Sokáig ez volt a legjobb módja a jegyzőkönyv készítésnek és bizonyos esetekben mind a mai napig. Könnyű formázni, könnyű megosztani. De nem védett a humán gyarlóságok ellen. Ugyanis nem vagyok egy jó jegyzőkönyv író. Ezt többször a fejemhez is vágták, mikor elhagytam az előző munkahelyemet. Miért nincs jegyzőkönyv arról, hogyan töltöttem fel a szekvenciákat az NCBI-ra? Miért nincsenek a szkriptek még jobban kommentelve? Hogyan csináltam ezt vagy azt? A probléma az, hogy ami számomra triviális, azt nem írtam le. Ha nagyon belemerültem a munkába és gyorsan követték egymást a lépések, a jegyzőkönyv és a valóság inkonzisztens állapotba került, ami csak később derült ki.

Milyen jó lenne egy olyan rendszer, ahol a munkavégzés egyben dokumentáció is lenne! A Bash saját history-ja végül is ilyen, de nem derül ki, mit miért csináltunk, ráadásul a fájlok kimenetét sem tartalmazza. A Screen-nek is van saját előzmény listája, ráadásul a kimenetet is képes tárolni a Ctrl + A és h billentyű kombinációval, ami igen viccesen néz ki egy less vagy vim kiadása után. Az R is tárol minden lépést, amit végrehajtottunk, tehát a hibás paraméterekkel végzett lépéseket, elgépeléseket, mindent. Ez sem tökéletes.

Jó volna egy olyan rendszer, ami egyesítené a Wikipédia formázási lehetőségeit az R/Bash előzménylistájával, de csak azokat a lépéseket tartalmazná, ami tényleg szükséges a munkához, nem a szerencsétlenkedéseket. Igen, van ilyen rendszer. Úgy hívják IPython.

A rendszert úgy kell elképzelni, mint egy Wikipédiát, amibe ha beleírjuk a parancsokat, azok a rendszer végrehajtja. A kimenetek is a dokumentumba kerülnek. A forráskódot is formázza, ami csak hab a tortán. Ez LaTeX alatt szinte lehetetlen lenne. De azért ez sem jelent mindenre megoldást. Először is csak Python alatt működik, én meg egy csomó programozási nyelvet használok. Illetve a napi munka során hármat: Bash-t, amiben a farmra küldöm ki a szekvencia feldolgozást. Azok kimenetéből Python szkriptekkel táblázatokat készítek, végül R-ben csinálom a végső elemzéseket.

A másik, sokkal nagyobb gond, hogy az IPython azon a gépen kell, hogy fusson, ahol végzem az elemzést. Ez legtöbbször egy távoli gép, amihez csak ssh kapcsolatom van. (Igen, lehet X-et alagutaztatni, de aki használt már IGV-t ilyen módon, az tudja, hogy ez minden, csak nem hatékony munkavégzés.)

Ugyan ez a probléma az RStudioval, knitr-el, csak épp R oldalon. De erre is van megoldás!

A Jupyter olyan, mint az IPython, de sokkal több nyelvet ismert. Igaz, egyszerre csak egyet lehet használni. Két részből áll. A kernel, ami a számítási feladatokat végzi. Ez lehet egy távoli gépen is, csak a megfelelő portok legyenek nyitva. A kliens pedig egy böngésző képes számítógépen fusson. Mindent a böngészőből irányíthatunk. A dokumentumunk különböző cellákra van osztva, alapvetően két típusuk van: leírás és kód. A leírásba szépen megfogalmazhatjuk, melyik lépést miért csináljuk, a kód pedig azonnal végrehajtja azt. Azonnal. Itt el is érkeztünk a rendszer egyik hátrányához.

Az egész az úgynevezett interaktív programozás keretében nyer értelmet. Tehát ha van egy lépés, ami három napig fut, akkor nem érdemes Jupyteren keresztül használni. A másik fontos dolog, hogy ez nem szoftver fejlesző környezet. A tesztelést és hibakeresést nem támogatja.

Maga a fájl egy JSON, ezért akár GitHub-on keresztül tárolhatjuk a verziókat. Ha elgépeltünk valamit, csak kijavítjuk azt és újra futtatjuk az adott cellát, nem az egész munkafolyamatot. Vizsgálatunkat exportálhatjuk HTML-be, PDF-be, de akár fel is tölthetjük egy szerverre, hogy megosszuk másokkal.

De a legjobb dologról még nem beszéltem. Ez egy olyan tulajdonság, amivel egyetlen hagyományos jegyzőkönyv sem rendelkezik. Mivel alapvetően böngészőben fut, lehetőség van interaktív elemeket beszúrni. A legjobb példa erre a plotty, Az ezzel elkészített grafikonok tetszőlegesen nagyíthatóak, mozgathatóak. A jegyzőkönyv olvasója kódolási ismeretek nélkül is tanulmányozhatja az eredményeket. Megszűnik a papír jelentette korlát és kód a dokumentációval tökéletes egységet alkot, megőrizve mindent örökké. Vagy legalábbis amíg van áram.

Már nagyon ritka, hogy egy munkafolyamathoz csak egyetlen programot használunk. Sok programot, sok különböző beállítást használunk, gyakran egymással összeegyeztethetetlen konfigurációkban.

Az is problémát jelenthet, ha nincs rendszergazdai hozzáférésünk az adott rendszerhez, de programokat kell rá telepítenünk. A megoldás a virtualizáció lehet. Ebben a posztban virtualizáció alatt értem mindazokat a technológiákat, amelyek lehetővé teszik, hogy több környezetet működtessünk, de nem feltétlenül párhuzamosan.

Természetesen csinálhatunk mindent kézzel. Linux alatt a PATH, LD_LIBRARY_PATH környezeti változókkal beállíthatjuk, hogy programjaink honnan olvassák be a fájlokat. Ha Python-t használunk, a modulok helyét a PYTHONPATH változó határozza meg. R esetében is van lehetőségünk újra definiálni a csomagok helyét: R_LIBS, R_LIBS_USER, R_LIBS_SITE. Ez a megoldás nem nyújt 100% védelmet a konfigurációk ütközése ellen, de limitált hozzáférésű rendszereknél egy lehetséges megoldás. Nekem például volt olyan problémám, hogy egy számítógép klaszteren a fő gépen más konfiguráció volt, mint a munkát végző csomópontokon. Ezekkel a trükkökkel viszont felül tudtam bírálni minden beállítást.

Természetesen az egész hardvert is emulálhatjuk. Ha nincs például PPC processzoros gépünk, a hardver emulációjával elérhetjük, hogy legyen. Újgenerációs Amiga-szerü rendszereken például futtathatunk régi, Motorola processzoros alkalmazásokat (Windows alatt a WinUAE teszi ugyan ezt). De említhetném a sokak által ismert VirtualBoxot is. Ez utóbbi segítségével az intézeti Windowsos laptopra telepítettem egy Linuxot. A géphez nincs rendszergazdai hozzáférésem, mert az IT nem engedi, de a virtuális Linuxhoz van.

Az egész hardver virtualizálása több erőforrást igényel, mintha csak egy operációs rendszert futtatnánk, de szerencsére van megoldás akkor is, ha csak a kernelt akarjuk virtualizálni. A Docker például kedvelt eszköz, pont ez a célja.

Könnyen összeállíthatjuk kedvenc szoftveres környezetünket, elküldhetjük másoknak és elfelejthetjük a telepítés során fellépő hibákat, mert a rendszer garantálja, hogy a csomagunkat használva mindenki ugyan azt kapja. Még webes szolgáltatások is vannak a csomagok tárolására.

Valamivel kevésbé flexibilis ugyan és nem is ez az eredeti célja, de a Conda is tartalmaz virtuális környezetek kezelésére parancsokat. A Conda igazi célja a csomagok egyszerű telepítése és frissítése, de a

conda create -n env
source activate env

létrehoz egy env nevű környezetet. Minden további parancs erre a környezetre fog vonatkozni. A Conda egyébként tartalmaz egy tisztán bioinformatikai csatornát biconda néven.

Elérkeztünk az RNA-seq legvitatottabb lépéséhez, az illesztéshez. Bárkivel beszéltem eddig RNA-seq-ről, mindenkinek megvolt a jól bejáratott véleménye arról, hogy melyik szoftvert kell használni és ebben a kérdésben nem ismert irgalmat. De én elárulom a nagy titkot. Csak nektek, ingyen. De aztán ne adjátok tovább, mert még mások is tudni fogják:

Tök mindegy. Először is, hiába használunk egy rendkívül kifinomult illesztőprogramot, ami a readek 100%-t illeszti, ha utána a kvantifikálás során a felét eldobjuk. Hiába akarunk egy elképesztően gyors alkalmazást futtatni, ha nincs elég memóriánk. Végezetül a legfontosabb indok: ha két gén expressziója jelentősen eltér, meg fogjuk találni azt, bármilyen programot is használunk.

Természetesen az illesztőprogramok között van különbség, nézzük meg, mik ezek. Még egy fontos dolgot meg kell említeni, ez pedig a felhasznált referencia típusa. Illeszthetünk transzkriptómra vagy genomra. Ha transzkriptómra illesztünk, nem kell foglalkoznunk azokkal a readekkel, melyek átérik az intronokat, tehát a DNS illesztésnél megismert programokat nyugodtan használhatjuk, például a BWA-t. Amire ügyelni kell, hogy a transzkriptóm egy gén több splice variánsát tartalmazhatja, ezért a paramétereket úgy válasszuk meg, hogy korrekt módon kezeljük a több helyre illeszkedő readeket. Erre most nem térek ki.

Ha genomra illesztünk, a többszörösen illeszkedő readek kevesebb gondot jelentenek (paralóg és pszeudó géneket leszámítva), cserébe meg kell küzdeni az intronokon átnyúló readek problémájával.

Bowtie

Ez a program elég egyszerű, még az indeleket sem kezeli. Pont ezért szeretik a mikroRNS elemzéseknél, ahol 20-24 nukleotid hosszú szekvenciákat kell a genomra pakolni, de azt gyorsan. Hátránya, hogy nagyon nagy genomokat nem képes kezelni. Embernél még nincs gond, de a búza genom gondot okoz neki.

TopHat

Kicsit öregecske program, már eljárt felette az idő. Lassú is szegény. Csak a történelmi hűség kedvéért említem meg. Annak idején a Tuxedo pipeline része volt, de a STAR megjelenésével sokan lecserélték. Ha az illesztésen kívül fúziós fehérjéket is akarunk keresni, a TopHat fusion a mi eszközünk.

GSNAP

Az EBI-ban az egy sejtes csoportok kedvenc illesztője. Változatos módon paraméterezhető. Sebességben a STAR mögött végez, de kicsivel több readet illeszt annál. Saját futtatásaim alapján viszont úgy láttam, hogy a génekre illeszkedő readek száma alacsonyabb, valamint több olyan readet is illeszt, ami három exonon is átnyúlik. Pontos számokat nem tudok mondani és az is elképzelhető, hogy paraméter tuningolással a fals illesztések száma is csökkenthető. Összességében teljesen korrekt program, memória szükséglete a STAR töredéke. Soha ne futtassuk alapértelmezett beállításokkal, mert használhatatlan eredményt produkál. A readeket annyi módon illeszti, ahány módon csak tudja, ami a sebességére is hátrányosan hat. Az EBI-os srácok a következő paraméterekre esküsznek:

gsnapl -t 4 -A sam -B 5 -n 1 -Q -N 1 -D referencia_dir -d referencia_nev read1.fastq read2.fastq

STAR

Ez az illesztő igazán univerzális. Hihetetlenül gyors, viszont a memória éhsége óriási, ami többszálú alkalmazásnál tovább növekszik. A Sangerben ez az RNA-seq pipeline része. A fejlesztője olyan opciókat is adott a programhoz, hogy a TopHathez hasonló kimenetet produkál, ezért egy TopHat alapú rendszerben fájdalommentesen beilleszthető. Rá épül egy fúziós fehérje kereső algoritmus, a STAR-Fusion. Ha van elég memóriánk érdemes ezt használni. Fontos, hogy alapértelmezett módon futtatva a nem illeszkedő readeket eldobja, ami a minőségi mutatók meghatározásánál gondot okozhat. Egy igazán sokoldalú paraméterezés a következő:

STAR --runThreadN 12 --outSAMstrandField intronMotif --outFileNamePrefix sample --outSAMattributes NH HI NM MD AS XS --twopassMode Basic --outSAMunmapped Within --chimSegmentMin 12 --chimJunctionOverhangMin 12 --alignSJDBoverhangMin 10 --alignMatesGapMax 200000 --alignIntronMax 200000 --chimSegmentReadGapMax parameter 3 --alignSJstitchMismatchNmax 5 -1 5 5 --limitBAMsortRAM 31532137230 --outSAMtype BAM SortedByCoordinate --readFilesCommand zcat --genomeDir reference --readFilesIn reads1.fastq.gz reads2.fastq.gz

Ezek a paraméterek STAR-Fusion és TopHat kompatibilisek.

HISat

Az új trónkövetelő. Erőforrás igénye alacsony, sebessége elképesztő. Rá épül az új Tuxedo pipeline, ami azért hagy még némi kívánni valót maga után. A program képes közvetlenül használni az SRA-t. Nem kell azokat Fastq-vá konvertálni.

Toby mindig meglepődött, ha a hülyeség váratlan helyen bukkant fel. Nem kellett volna neki. Negyvenkét évesen épp elég tapasztalata lehetett volna, hogy ne lepődjön meg olyan elemi dolgokon, mint, hogy vannak ostoba emberek.

Ez mégis csak a Broad Institute - gondolta magában. Ez egy hely, ahol a minőséget folyamatosan magasan tartják azzal, hogy mindenkinek be kell számolnia az eredményekről a ranglétra magasabb fokán álló embereknek. Ez pedig azt jelenti - okoskodott tovább Toby - hogy az arra érdemtelen embereket kiszelektálják. Egy molekuláris biológiával is foglalkozó tudományos intézetnél ez egyet kellene, hogy jelentsen az ostobaság eliminálásával.

Toby azért nem panaszkodhatott. Rengeteg brilliáns koponyával találkozott. Kollégái remek emberek voltak. Leszámítva Pault. Paul ugyanis munkája során mellőzött minden statisztikai eljárást és különböző Excel bűvészkedésekkel kereste a rákos mutációkat. Rendezgette az adatokat, kézzel eltávolított elemeket, melyek tapasztalata alapján "nem lehettek fontosak". Toby ezért igyekezett távol tartani magát tőle.

Most viszont Paul megtalálta és hozta az egyik haverját, Ulrichot is.

Toby elég sok expressziós adatot dolgozott fel élete során, ezért egy vizsgálathoz az ő segítségét kérték. Az elején minden rendben ment. Paul, Ulrich és az egyik PhD hallgatójuk, Grace meghívták egy megbeszélésre, ahol átbeszélték a kísérlet hátterét, a lehetséges feldolgozási lépéseket és felosztották a munkát. Ulrich határozott, kompetens személynek tűnt, aki értette a vizsgált betegség biológiai hátterét. Paul a megbeszélés alatt az Excellel játszott és büszkén mutogatta a nyers readszámokat, mint expressziós profilt, egészen addig, míg Toby fel nem világosította annak felesleges voltáról.

A megbeszélés után Grace rendszeresen felkereste Toby-t, hogy segítsen neki az analízis egyes lépéseinél. Gyorsan elsajátította az ismereteket, nem kellett neki kétszer elmondani semmit. Sőt, még Toby is ellesett egy-két R fogást. Grace néha panaszkodott Ulrichra, hogy ennyire képtelen egy témára koncentrálni, de Toby betudta ezt az ősidők óta fennálló hallgató-témavezető viszálynak és nem foglalkozott vele.

De a hülyeség nem egy masszív tömeg. A hülyeség alaktalan gáz. Nem állítják meg falak, gátak. A hülyeség átszivárog mindenen. Megtalálja a legkisebb rést is és könyörtelenül keresztül furakszik. Tobyhoz is elkezdett szivárogni, méghozzá furcsa e-mailek formájában. Az első egy metilációs heatmap volt PNG formátumban. Ulrich azt tudakolta, a minták hierarchikus elrendezése mennyire hasonlít az expressziós minták elrendezéséhez. Ha a kérés nem lett volna már így is elég furcsa, de a két kísérlet alig tartalmazott közös elemeket. Mintha valaki azt kérdezte volna: itt ez a szép piros paradicsom. Mennyire hasonlít erre a sárga paprikára? Termés mindkettő, tehát hasonlítanak. Egyik gömböly, másik nem. Tehát nem hasonlítanak.

Toby is hasonló dilemmával nézett szembe.

Egy másik alkalommal Ulrich kooperációs partnere további expressziós eredményeket adott, de CPM normalizálva, mert ők azt tarották a legjobbnak. Toby TPM normalizálást használt, ahogy korábban egy megbeszélés során eldöntötték. Talán csak az egy betű különbségnek köszönhetően, de Ulrich megkérdezte, össze lehet-e hasonlítani a kettőt. Az e-mail alján Toby látta a korábbi levelezést az említett partnerrel, ahol Ulrich már feltette ezt a kérdést és a választ is látta, ahol elmagyarázták neki, miben más a kettő. Ennek ellenére tőle ugyan azt megkérdezte. Toby nagy levegőt vett és a lehető legudvariasabban elmagyarázhat, hogy nem lehet őket összehasonlítani, de ő szívesen megcsinálja a CPM normalizálást, ha azt kérik tőle. Igyekezett más szavakat és kifejezéseket használni, mint amit a levél alján látott, mert abban bízott, így talán könnyebben megértik mondandóját.

De Ulrichet nem lehetett ilyen könnyen lerázni. Még egyszer visszakérdezett: A CPM ugyan az, mint a TPM? Toby keze megállt a billentyűzet felett. Azért nem szerette a hülyéket, mert a sok magyarázás és értetlenkedés után azt gondolta, ő maga a hülye, amiért nem képes érthetően kifejezni magát. Mint, amikor mindenki röhög egy poénon, csak egyvalaki néz kérdő szemmel, hogy most mi a csattanó?

Toby végül azzal nyugtatta magát, hogy különböző entitások különböző nevet szoktak kapni. Tehát ha más a neve, akkor más tulajdonságai vannak. Mint ahogy az sem mindegy, ki fekszik a nászéjszakán mellettünk: Józsi vagy Rózsi. Itt is csak egy betű a különbség.

Egyik ebédnél épp a fenti dilemmán töprengett, amikor mellé ült Katerina. Már korábban is beszélgettek és ahogy az lenni szokott, nagyon hamar munkára terelődött a szó. Tobynak meglepetésként hatott a hír, hogy Katerina Ulrich post-docja és Paul munkatársa. A nő elmondta neki, hogy főnöket képtelen egy dologra koncentrálni, folyamatosan egymással össze nem függő vizsgálatokat kér.

- Legtöbbször a könyvtárba menekülök előle - mondta Katerina. Ha meglát, rögtön eszébe jut valami hülyeség. Ha pedig haladni akarok a saját témámmal, otthonról dolgozom.

- Paul miért nem menekül el?

- Pault nagyon kedveli Ulrich, mert bármilyen ötlete is van Ulrichnak, Paul pillanatok alatt választ tud adni a táblázatok módosításával. Bármit össze tud hasonlítani, bármit kielemez.

- De Paul legtöbb válaszának semmi értelme.

- Igen, de a legtöbb kérdésnek sincs, amit Ulrich feltesz.

Toby mégsem szerette ilyen rövid ismerettség után kijelenteni emberekről, hogy diettánsok. Talán csak arról van szó - vélekedett. - hogy nem tudták feldolgozni azt a váltást, amit a számítógépek megjelenése okozott a biológiában. Talán próbálják megérteni ezt a számáukra idegen közeget, csak nem megy nekik. Hiszen egy világbajnok sprintertől sem lehet elvárni, hogy lefussa a maratont. De nem is állnak rajthoz.

Miután kiderült, hogy némi R szkripteléssel klikkelések ezreit tudom megspórolni, újabb búzás megbizatást kaptam. Maga a számítás nem volt bonyolult, hamar megvoltam vele. Utána viszont hosszú ideig nem történt semmi olyan, amihez a segítségem kellett volna.

A publikálás itt is nehezen ment. Talán három újságot is megjárt a kézirat, folyamatos átalakítást és újabb vizsgálatokat kellett végezni, de egyikhez sem kérték a segítségemet. Aztán a munkahely váltás még jobban elszigetelt a témától.

A helyzet pikantériája, hogy annyi új ismeretet tanultam itt, hogy ha most kellene megcsinálni a vizsgálatokat, sokkal többet tudnék a cikkhez hozzátenni. De ezen már nem lehet segíteni és a cikket sem lehet a végtelenségig csúsztatni, csak azért, mert van egy újabb ötletünk.

Mert ahhoz is érteni kell, hogy mikor hagyjunk abba egy munkát. Ha folyamatosan új kísérleteket csinálunk, új adatokat vonunk be, soha nem fogunk végezni és az a téves képzet alakul ki rólunk, hogy nem csinálunk semmit. Ez valahol igaz is, mert a hosszú munkával a beosztottak elunják magukat, az új projektek magasabb prioritást kapnak és végül minden ellaposodik, elhal. De akár a versenytársak is beelőzhetnek.

De ez a munka nem halt el, mert szerencsére nem hagyták elhalni. Nem úgy, mint más munkák, amikről azért nem írok, mert még reménykedem, hogy lesz belőlük valami. Vagy más kutatók, akik már 2014-ben beismerték egy intézet szintű összejövetelen, hogy már minden megvan a publikáláshoz, mégis további PhD hallgatókat emésztenek fel a kísérletek.

Egyszer volt, hol nem volt, valahol a CRT monitorok idején, de még a tapicskolós mobiltelefonok megjelenése előtt, a Napster fénykorában, volt egy zenelejátszó program, a WinAmp. Manapság, amikor gigabájtos tárak csücsülnek a zsebünkben, megmosolyogtató, hogy számítógépen tároltuk a zenéket, sőt még ott is játszottuk le őket (a prolik még CD-re is kiírták!), de nem volt más választásunk.

A WinAmp nagyszerű kis program volt, népszerűségben a Total Commanderrel (akkoriban Win Commander) vetekedett. Aztán történt valami rettenetes. Felütötte a fejét a túlprogramozás. A funkciók úgy potyogtak a kis szoftverbe, mint esőcseppek zápor idején. Kettőt pislogott a felhasználó és a program már filmeket is lejátszott, a külalak testre szabása pedig fontosabbá vált, mint maga a funkció. Az eszköz megszűnt az lenni, ami a célja volt és amiért szerették.

Ugyan ezt veszem észre az SRA Toolkit-el kapcsolatban is. A tárhely éhség csillapítására tömörítéseket vezettek be az NCBI-nál és elkészítették a kicsomagoló alkalmazást is. Kezdetben csak a Fastq konverzió volt a feladata, de mostanra egy felhő alapú, eloszott rendszer lett, ami a legegyszerűbb funkcióját is nyakatekert módon hajtja végre. A dokumentáció szerint a letöltött 180 MB nem müködöképes, szükséges hozzá a szerveren tárolt másik kódrészlet. Az, hogy milyen porton kommunikál, hogy esetleg megoldjuk a lehetséges hálózati problémákat, nem derül ki. A hibaüzenetek kriptikusak, nem szólnak semmit az aktuális problémáról, csak a virtuális adatbázis bejegyzésekről, amit valami hálózati alrendszeren keresztül akar végrehajtani. Ilyen szöveget manager meetingen mondanak, nem hibaüzenetekben.

A neten a legelterjedtebb megoldási forma: ne használd a programot, töltsd le a fastq-t az ENA-ról. Nekem nem volt szerencsém, a kérdéses szekvencia nem volt elérhető az említett formátumban. Kénytelen voltam megküzdeni az alkalmazással.

Habár elméletileg elég a szekvencia azonosítóját használni (az sra kiterjesztés nélkül megadott nevet azonosítónak tekinti a program és megpróbálja letölteni azt a szerverről), a gyakorlatban úgy tűnik a fastq-dump az egyetlen, amiből ez a funkció kimaradt. Miután letöltöttem mindent, folyamatosan időtúllépésre panaszkodott, mintha le akarna valamit tölteni. De ha már úgyis le akar valamit tölteni, letölthetné a szekvenciát is, nem? Mindegy, nem kell értelmet keresni egy túlprogramozott szoftverben.

Valami megmagyarázhatatlan csoda folytán a laptopomon működött minden. Csupán tárhely nem volt a szekvenciáknak. Ezért sshfs-el felcsatoltam a számítógép farm egyik gépét, majd oda mentettem mindent.

Később egy kolléga az NCBI embereit is bevonva kiderítette, hogy be kellett volna állítani a proxyt a vdb-config paranccsal. Csak azt nem értettem, máshol (git, wget, R, pip) miért nem kellett ilyesmit beállítani? Azok miért működnek nélküle? Bocsánat, megint értelmet keresek egy túlprogramozott szoftverben.

A szekvenálás minőségének meghatározása szinte teljesen megegyezik a DNS szekvenálásnál alkalmazott módszerekkel. Ez nem is meglepő, hiszen miután a minta előkészítés során az RNS-t átírtuk DNS-é, a gép ugyan úgy hajtja végre a szekvenálást. A FastQC ugyan úgy használható.

Az értékelés során viszont előfordulhat, hogy engedékenyebbek lehetünk. Ha kisRNS-eket szekvenálunk, az adaptor kontamináció nem hiba, a szekvenálás sajátosságából adódik. Egyszerűen a read hosszabb, mint a 20-24 nukleotidból álló molekula. Degradóm esetén úgyszintén.

De a felülreprezentált szekvenciák előfordulása sem feltétlenül hiba. Egy erősebben expresszálódó gén jelenlétét a program értelmezheti PCR hibának. Ha nem egy ismert adaptor a szekvencia, figyelmen kívül hagyhatjuk.

Attól, hogy a szekvencia maga jó, még nem jelenti azt, hogy értékelhető eredményt kapunk belőle. Illesztés után újabb minőségi ellenőrzést kell tartanunk. Különösen, ha egy sejtes szekvenálással van dolgunk. Erről lesz részletesebben is írás, most csak annyit mondok, ha a sejt feltárás során nem vagyunk elég óvatosak, az RNS bomlásnak indulhat és szekvenálás után csak annyit látunk, hogy kromoszómális génre alig, mitokondriális génre viszont töménytelen mennyiségű read illeszkedik. Ezen pedig a minta újraszekvenálása sem fog segíteni.

Érdemes tehát megnézni mennyi read esik arra a régióra, amit szekvenálni akarunk. KisRNS-eknél ezen kívül előszeretettel használják a méret eloszlásokat is. Egy sejtes szekvenálásnál több diagnosztikai ábrát is készítünk. Mivel nem ritka, hogy százas nagyságrendben szekvenálnak sejteket, ezért az X tengelyen az össz readszám szerepel, míg az Y tengelyen a százalékos mitokondriális génre illeszkedő readek száma vagy bármilyen génre illeszkedő readek aránya szerepel. Így viszonylag kevés ábrán sok mintát lehet áttekinteni. Ha az R-ben interaktív módon rajzoljuk ki a pontokat, klikkeléssel azonosíthatjuk a kiugró értékeket. Erre majd az egy sejtes posztban térek ki. Ha el nem felejtem.

Egy másik diagnosztikai eljárás a gén lefedettség (gene body coverage). Itt az összes génre illeszkedő lefedettséget ábrázolják, miután 100 nukleotid hosszúságúra normalizálják azokat. Ha bármelyik vég degradálásnak indul, ezen az ábrán könnyen észrevehető. Több, más statisztikával együtt az RSeqQC segítségével készíthetünk ilyen ábrákat.

A FastQC elég jó, de ha sok mintánk van, különböző, csak számunkra logikus könyvtár struktúrába rendezve, a FastQC használata elég macerás. Még akkor is, ha parancssorból is meghívható. Szerencsére a lusta elfoglalt embereknek fejlesztett multiQC megoldást nyújt erre.

Csupán egyetlen paramétert vár, a mintáinkat tartalmazó könyvtár nevét. A program rekurzívan bejárja az összes alkönyvtárat, a megtalált fájlok alapján megpróbálja kitalálni, milyen programokat futtattunk és HTML alapú összefoglalókat készít. Ennél egyszerűbb programot el sem lehet képzelni.

Már annyi RNA-seq adatot dolgoztam fel, hogy lassan én is elhiszem, hogy értek hozzá. Ez, és túláradó magabiztosságom együtt ezt a poszt sorozatot eredményezte, eredményezi.

Igyekszem valami rendszert felállítani, de annyi egymással összefüggő téma van, hogy akaratlanul lesznek átfedések. Azért remélem nem sűlyedünk a "Mindenféle dolog kusza összevisszasága" szintjére.

Az RNA-seq alapja az RNS szekvenálása. Először RNS-t izolálnak a mintából, ami lehet 1 sejt, néhány sejtes kolónia vagy még több sejtet tartalmazó minta. Ez "csak" olyan mértékben befolyásolja a bioinformatikai munkát, hogy mekkora amplifikációs hibával kell számolni. Mivel többféle RNS molekula létezik, az izolálás során megcélozhatunk bizonyos molekula típust és szekvenálhatunk mRNS-t, kisRNS-eket, riboszóma RNS-t, stb. A továbbiakban elsősorban az mRNS szekvenálásról lesz szó, mivel (szerintem) ez a legelterjedtebb RNA-seq alkalmazás.

Az mRNS-t a poliA vég segítségével viszonylag egyszerű kigyűjteni. A könyvtár készítés többi lépéséről nem írnék, mert ezt mások sokkal jobban tudják nálam, csupán egy részt emelnék ki, ami az analízist befolyásolhatja, ez pedig az irány specifikus (strand specific) kit alkalmazása. Ez röviden annyit jelent, hogy a szekvenálás során kapott read orientáció összefügghet annak az mRNS-nek az orientációjával, amiből származik. Az átfedő antiszensz transzkriptek esetén pontosabban meghatározható az expresszió. Itt a Sangerben az esetek 98%-ban a read orientáció megegyezik a transzkript orientációjával (first strand) és a veterán bioinformatikusok között van, aki már hallott olyanról, hogy valaki ellentétes irányú kitet használt (second strand).

Szemben a DNS szekvenálással, ahol a lefedettség jóval egyenletesebb, különböző transzkriptek lefedettsége eltérő lehet. Gyakran felmerülő kérdés, milyen mélységben szekvenáljuk mintáinkat? Természetesen attól függ, mire keressük a választ, de az Encode konzorcium ajánlása jó kiindulási alap. Ők egy standard emberi differenciál expressziós vizsgálathoz 30 millió readet elégnek tartanak, alacsonyan expresszálódó transzkriptek vizsgálata esetén 100-200 millió readet ajánlanak.

Ugyancsak gyakori vitakérdés az alkalmazni kívánt read hossza. Ugyan csak az analízistől függ. Ha csak gén szinten vagyunk kíváncsiak az expressziós különbségekre, akkor elég lehet a rövidebb méret is. Izoforma variánsok detektálása és referencia genommal nem rendelkező élőlények esetén a "hosszabb jobb" elv érvényesül.

Végezetül álljon itt néhány link, ahol nálam részletesebben mutatják be ezt a technológiát:

https://www.youtube.com/watch?v=0e1zWjgZODc

https://www.youtube.com/watch?v=Ji9nFCYl7Bk

http://rnaseq.uoregon.edu/

https://github.com/griffithlab/rnaseq_tutorial/wiki