Az immunoprecipitáción alapuló szelvenálási módszerek sokkal körültekintőbb ellenőrzést igényelnek, mint a "hagyományos" DNS vagy RNS szekvenálás. Bár a címben ChIP-seq szerepel, az itt leírtak szinte minden immunoprecipitáción alapuló módszerre (CUT&RUN, ChIP-exo, Dip-seq) alkalmazhatóak.
Az első lépés természetesen a kapott FASTQ fájlok ellenőrzése például a FastQC programmal. Itt ugyan azok az elvek érvényesek, mint bármilyen más szekvenálás ellenőrzésénél, bár a szekvencia duplikáció magasabb lehet, mivel a teljes genomnak csak egy kis részét szekvenáljuk. Ha a FASTQ fájlok megfelelőek, akkor más módszereknél megnyugodhatunk, hogy az elemzéssel nem várható komoly gond (hacsak nincs kontamináció), viszont ChIP-seq akkor is lehet használhatatlan, ha a FASTQ fájlok tökéletesek.
A szekvenciák illesztése ugyan úgy zajlik, mint más esetben. Amennyiben adaptor szennyezést tapasztaltunk, természetesen azokat levághatjuk, de a BWA úgysem illeszti az idegen részeket. Az illesztett readek aránya ugyancsak hasznos információ, 90% alatt nem szabadna lennie.
Utána jöhet a peak keresés. Ezt legtöbbször a macs3-al végezzük, de meg kell jegyezni, hogy bár ez a legelterjedtebb peak kereső, vannak más alternatívák, és lehetnek olyan kísérleti elrendezések, ahol a macs3 nem a legjobb választás. Ökölszabályként elmondhatjuk, hogy 10 ezer alatti peakszám valami problémát jelent.
Az immunoprecipitáció nem egy egyszerű lépés, ezért minden mintát kétszer szekvenálnak meg. Van egy úgynevezett "input", ahol a szekvenálás immunprecipitáció nélkül történik. A párja természetesen az "IP", ahol a korábban említett lépést nem hagyják ki. Ez a két minta alkot egy egységet. Bár vannak bátrabb bioinformatikusok, akik input nélkül elemeznek, és tőlem is többször kérték, hogy "azért nézzük meg input nélkül is...", tudni kell, hogy a sokkal jobb jel-zaj arányt kapunk, kevesebb lesz a fals pozitív, ha inputot is használunk.
A peak fájlok BED formátumúak, a fő információ a kromoszóma koordináta, a legtöbb program azt használja, még akkor is, ha a macs3 további értékeket tárol.
ChIPQC
Ez egy nagyon jó R csomag lenne, ha nem tartalmazna 11 éves elavult függőséget, ami miatt a legújabb Bioconductor nem tartalmazza. Nekem is csak azért működik, mert kézzel forgattam a csomagot és a kérdéses függőséget. Használata nagyon egyszerű:
sample <- ChIPQCsample("ip.bam", peaks="NA_peaks.narrowPeak", annotation="hg19")
A két minta közül az IP-t kell megadni illesztésként, a peak fájlt, valamint a genom annotációját. Táblázatos formában megkapjuk az összes hasznos minőségi mutatót a következő paranccsal:
QCmetrics(sample)
Lássunk egy kimenetet:
Reads Map% Filt% Dup% ReadL FragL RelCC SSD RiP%
3320201.000 100.000 37.700 0.000 51.000 174.000 0.683 0.237 2.210
Ez egy tipikusan rossz kísérletből származik. Habár a readek 100%-a illeszkedik a genomra (ez a duplikáció kiszűrése utáni BAM fájl, de eredetileg is 96% volt. És ezért 0% a Dup oszlop), a peakekre csak a readek 2,21% esik (RiP, read in peaks), ami nagyon rossz. Az SSD a normalizált lefedettségből készített hisztogram szórása, amire a dokumentáció azt írta, "minél nagyobb, annál jobb". Nos, jó kísérleteknél ez 20-40 is lehet, de tipikusan 2-3. A 0.237 nem elfogadható. Végül a RelCC egy olyan mérőszám, ami a forward és reverz readek lefedettsége között számol korrelációt.
DiffBind
Ez a csomag különbségeket keres a kezelt mintában, viszont van pár nagyon hasznos funkciója, amivel a minták minőségét vizsgálhatjuk. Ha vannak ismétléseink, és nem csak egy mintával akarjuk megváltani a világot, akkor fontos tudni, hogy van-e batch effekt, nem történt-e minta összecserélés, vagy bármi olyan, amit a laborosok kerek perec képtelenségnek tartanak (és amiről utólag kiderül, hogy nem is olyan képtelenség).
A program használatához készíteni kell egy DataFrame-et, ami tartalmazza a mintákat, illesztett BAM-ok, peak-ek elérési útját, kezelést, ismétlések számát, . Ezen kívül bármilyen más információt felvihetünk, de ez a minimum, ami szükséges.
sampleSheet <- data.frame(SampleID=c("c1","c2","em1", "em2"),
Tissue=c("unknown", "unknown","unknown","unknown"),
Treatment=c("control","control","treatment","treatment"), Replicate=c(1,2,1,2),
bamReads=c("bams/c1ip.bam", "bams/c2ip.bam", "bams/em1ip.bam", "bams/em2ip.bam"),
bamControl=c("bams/c1input.bam", "bams/c2input.bam", "bams/em1input.bam","bams/em2input.bam"),
Peaks=c("peaks/C1/c1_peaks.narrowPeak", "peaks/C2/c2_peaks.narrowPeak", "peaks/EM1/em1_peaks.narrowPeak"),
PeakCaller=c("narrow","narrow","narrow","narrow"))
Ezután létrehozunk egy adatbázist a fenti információkkal:
sample <- dba(sampleSheet = sampleSheet)
Ha megvagyunk az adatbázis elkészítésével, jöhetnek az érdekes plotok. Először is egy PCA:
dba.plotPCA(samples, label=DBA_ID)
Ha a minták nem a kísérleti elrendezés szerint csoportosulnak, akkor valami gond van.
Homer
A Homer egy mindent az egyben csomag. Nagyon könnyű használni, a nehézségét az adja, hogy annyi minden van beleépítve, hogy elsőre a felhasználó azt sem tudja, mit futtasson. Dedikált QC elemzés nincs benne, nem ír ki számokat, de az elemzés során ad szöveges megjegyzéseket a minőséggel kapcsolatban. Például a peakCall parancs adhat ilyen eredményeket:
Guessing sample is ChIP-Seq - uneven enrichment between same strand and
different strands - may have problems such as clonal amplification.
IGV
Bár az IGV nem számol semmit, segít, hogy egy általános képünk legyen az illesztésről. Ha szemmel nem látunk feldúsuló readeket, akkor a programok sem fognak peakeket találni.