HTML

Az élet kódjai

Csináld maga. Senki nem csinálja meg helyetted.

Friss topikok

Majdnem mindent az RNA-seq-ről (4.rész)

2017.01.16. 00:58 Travis.CG

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.

Szólj hozzá!

Címkék: bioinformatika

A bejegyzés trackback címe:

https://cybernetic.blog.hu/api/trackback/id/tr5212069567

Kommentek:

A hozzászólások a vonatkozó jogszabályok  értelmében felhasználói tartalomnak minősülnek, értük a szolgáltatás technikai  üzemeltetője semmilyen felelősséget nem vállal, azokat nem ellenőrzi. Kifogás esetén forduljon a blog szerkesztőjéhez. Részletek a  Felhasználási feltételekben és az adatvédelmi tájékoztatóban.

Nincsenek hozzászólások.
süti beállítások módosítása