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.
