A cufflinks egy remek program, hogy a szekvenálásokból megállapítsuk, melyek lehetnek a kódoló régiók. Sajnos a program több, egymással összefüggő táblázatot generál, aminek a feldolgozása nem egyszerű feladat. Szerencsére az R programcsomag ismét a segítségünkre siet.
A cummeRbund csomag segítségével könnyen felderíthetjük a hatalmas adat mögött megbúvó összefüggéseket. Első lépésként töltsük be az adatokat:
library(cummeRbund)
adat = readCufflinks(dir="directory", genome="genom.fasta", gtfFile="annotacio", rebuild=T)
A parancs hatására adatainkat az adat változóval érhetjük el. A directory a cufflinks eredményeink kimeneti könyvtára. A genome.fasta az elemzéshez felhasznált genom, míg a gtfFile paraméterrel adhatjuk meg azt ismert transzkriptek annotációját. Amit nem látunk, az egy adatbázis generálás a háttérben. A szöveges állományok a gyorsabb kereshetőség érdekében egy SQLite adatbázisba lesznek elmentve. Ennek helye a directory-ban lesz.
Használata annyira egyszerű, hogy nem is kell érteni az R-hez, ha szép képeket és táblázatokat akarunk készíteni. Például ha normalizált expressziós értékek eloszlására vagyunk kíváncsiak, nem kell mást tenni, csak kiadni a
dens <- csDenstity(genes(adat))
dens
parancsot. A legtöbb parancs egyből rajzolja is az ábrákat. Annyira egyszerű használni, hogy csak a dokumentációt tudom ismételni. De hogy valami olyasmit is mutassak, ami nincs leírva a dokumentációba, nézzük, hogyan lehet az unalmas cufflinks azonosítókat gén szimbólumokra cserélni.
Az eredményül kapott szignifikáns gének neve ilyen furcsa karakterekkel kezdődik: XLOC_. Ha ezt egy biológus kezébe nyomjuk, az sikítani fog. Az annotáció segítségével viszont kicserélhetjük gén nevekre.
annot <- annotation(genes(adat))
sig.data <- as.data.frame(getSigTable(adat, alpha=0.05, level="genes"))
sig.data$gene_id <- row.names(sig.data)
sig.with.annot <- merge(annot, sig.data, all.y=T)
Mit csinál ez a kód? Először is létrehoz két data.frame típusú változót. Az egyik neve annot, a másiké sig.data. Az annot-nak van egy gene_id nevű oszlopa, ezt létrehozzuk a sig.data-ban is. Ezután a merge paranccsal összefűzzük őket. Az all.y segít, hogy ne veszítsünk adatot a sig.data-ból. Alapértelmezetten ugyanis a merge a két halmaz metszetét adja vissza.
Egy másik trükkel azt mutatom be, miként lehet GO analízist végezni. Ha olyan fajjal dolgozunk, amelynek megtalálható az annotációja a Bioconductorban, nincs nehéz dolgunk (ha nem található meg, az gáz):
go <- toTable(org.Hs.egGO) # get GO ids
go$term <- Term(go$go_id) # get GO terms
lipid <- go[grep("lipid", go$term),] # get GO ids by keyword
sym <- toTable(org.Hs.egSYMBOL) # get all the gene symbols
lipidsym <- merge(lipid, sym) # now the table contains
genesymbols, GO terms and GO IDs
x <- merge(lipidsym, annot, by.x = "symbol", by.y = "gene_short_name") #
put XLOC names into the mess
lipidgenes <- x$gene_id.y
lipidgenes <- getGenes(adat, lipidgenes) # select genes from
expression data
x <- expressionPlot(lipidgenes)
A kód működéséhez be kell tölteni a GO.db és az org.Hs.eg könyvtárakat. Ha más fajjal dolgozunk, akkor természetesen az adott fajnak megfelelő annotációs könyvtárra van szükségünk. Először táblázattá alakítjuk a GO annotációt, és hozzáfűzzük a táblázathoz a leírásokat. A grep parancs segítségével kiválasztjuk a minket érdeklő GO kategóriákat. Az annotációból kiszedjük a gén szimbólumokat. Erre a lépésre azért van szükség, mert ez az adat képez kapcsolatot a cummeRbund-os eredmények és a GO annotáció között. A két merge segítségével összefűzűnk minden adatot, kiválasztjuk a cummeRbund-os azonosítókat és máris van egy génlistánk GO alapján. Végül megnézzük az expressziós profilját a listánknak.
Remélem mindenki kedvet kapott a cummeRbundhoz.
