Lépésről lépésre útmutató a DNS-szekvenálási adatok elemzéséhez

Dr. Javier Quilez Oliete, a Kolabtree tapasztalt bioinformatikai tanácsadója átfogó útmutatót nyújt a DNS-szekvenálási adatok elemzéséhez, beleértve az adatok olvasásához használt eszközöket és szoftvereket.

Bevezetés

A dezoxiribonukleinsav (DNS) az a molekula, amely egy szervezet genetikai információinak nagy részét hordozza. (Egyes vírustípusokban a genetikai információt a ribonukleinsav (RNS) hordozza). A DNS-molekulák alapegységei a nukleotidok (amelyeket hagyományosan A, C, G vagy T betűkkel jelölnek). Fogalmilag a DNS-szekvenálás a DNS-molekulát alkotó nukleotidok leolvasásának folyamata (pl. a “GCAAACCAAT” egy 10 nukleotidból álló DNS-szál). A jelenlegi szekvenálási technológiák több millió ilyen DNS-olvasatot állítanak elő ésszerű idő alatt és viszonylag alacsony költséggel. Referenciaként egy emberi genom szekvenálásának költsége – a genom a DNS-molekulák teljes halmaza egy szervezetben – leesett a 100 dolláros határról, és ez néhány nap alatt elvégezhető. Ez ellentétben áll az emberi genom szekvenálására irányuló első kezdeményezéssel, amely egy évtized alatt készült el, és körülbelül 2,7 milliárd dollárba került.

Ez a nagy áteresztőképességű és alacsony költségű DNS-szekvenálási képesség lehetővé tette egyre több szekvenáláson alapuló módszer és alkalmazás kifejlesztését. Például a teljes genomok vagy azok fehérjéket kódoló régióinak szekvenálása (két megközelítés, amelyet teljes genom-, illetve exom-szekvenálásnak neveznek) beteg és egészséges egyéneknél betegséget okozó DNS-változásokra utalhat. A DNS-ből átírt RNS szekvenálása – ez az RNS-szekvenálásnak nevezett technika – arra is alkalmas, hogy számszerűsítsük a génaktivitást és azt, hogy ez hogyan változik különböző körülmények között (pl. kezeletlen versus kezelés). Másrészt a kromoszóma-konformáció rögzítő szekvenálási módszerek a közeli DNS-molekulák közötti kölcsönhatásokat detektálják, és így segítenek meghatározni a kromoszómák sejten belüli térbeli eloszlását.

A DNS-szekvenálás ezen és más alkalmazásainak közös jellemzője a gigabájt nagyságrendű, több millió leolvasott szekvenciát tartalmazó adathalmazok létrehozása. Ezért a nagy áteresztőképességű szekvenálási (HTS) kísérletek értelmezéséhez jelentős adatelemzési képességekre van szükség. Szerencsére a legtöbb HTS-adattípushoz léteznek dedikált számítási és statisztikai eszközök és viszonylag szabványos elemzési munkafolyamatok. Míg a (kezdeti) elemzési lépések egy része a legtöbb szekvenálási adattípus esetében közös, a további elemzések az adatok típusától és/vagy az elemzés végső céljától függenek. Az alábbiakban bemutatom a HTS-adatok elemzésének alapvető lépéseit, és hivatkozom a népszerű eszközökre.

Az alábbi szakaszok egy része a rövid leolvasású szekvenálási technológiák (főként az Illumina) által generált adatok elemzésére összpontosít, mivel történelmileg ezek domináltak a HTS-piacon. A hosszabb leolvasásokat generáló újabb technológiák (pl. Oxford Nanopore Technologies, PacBio) azonban gyorsan teret hódítanak. Mivel a hosszú leolvasású szekvenálásnak vannak sajátosságai (pl. magasabb hibaarány), az ilyen típusú adatok elemzésére speciális eszközöket fejlesztenek ki.

A nyers leolvasások minőségellenőrzése (QC)

A lelkes elemző a FASTQ fájlokból kezdi az elemzést; a FASTQ formátum már régóta a rövid leolvasású szekvenálási adatok tárolásának szabványa. A FASTQ fájlok lényegében több millió leolvasás nukleotidszekvenciáját és bázisonkénti hívásminőségét tartalmazzák. Bár a fájl mérete a leolvasások tényleges számától függ, a FASTQ fájlok jellemzően nagyok (megabájt és gigabájt nagyságrendűek) és tömörítettek. Megjegyzendő, hogy a FASTQ fájlokat bemenetként használó legtöbb eszköz tömörített formátumban is képes kezelni azokat, így a lemezterület megtakarítása érdekében ajánlott nem kitömöríteni őket. Konvencióként itt egy FASTQ-fájlt egy szekvenálási mintának fogok megfeleltetni.

A FastQC valószínűleg a legnépszerűbb eszköz a nyers leolvasások QC-jének elvégzésére. Futtatható vizuális felületen keresztül vagy programozottan. Míg az első lehetőség kényelmesebb lehet azoknak a felhasználóknak, akik nem érzik jól magukat a parancssori környezetben, az utóbbi összehasonlíthatatlan skálázhatóságot és reprodukálhatóságot kínál (gondoljunk csak arra, milyen fárasztó és hibakényes lehet kézzel futtatni az eszközt több tíz fájl esetében). Bárhogy is legyen, a FastQC fő kimenete egy HTML-fájl, amely egy adott mintából származó nyers szekvenálási leolvasások általános minőségére vonatkozó legfontosabb összefoglaló statisztikákat közöl. A FastQC-jelentések tucatjainak egyenkénti vizsgálata fárasztó, és megnehezíti a minták közötti összehasonlítást. Ezért érdemes lehet használni a MultiQC-t, amely a FastQC-ből (valamint más, a későbbi lépésekben használt eszközökből, pl. adapter-trimmelésből, igazításból) származó HTML-jelentéseket egyetlen jelentésben összesíti.

AQC-információk célja, hogy a felhasználó meg tudja ítélni, hogy a minták jó minőségűek-e, és ezért felhasználhatók-e a további lépésekhez, vagy el kell őket vetni. Sajnos a FastQC metrikák alapján nem létezik konszenzusos küszöbérték a minták jó vagy rossz minőségűvé minősítésére. Az általam alkalmazott megközelítés a következő. Elvárom, hogy minden olyan minta, amely ugyanazon az eljáráson ment keresztül (pl. DNS extrakció, könyvtárkészítés), hasonló minőségi statisztikával rendelkezzen, és többségében “megfelelt” jelzővel legyen ellátva. Ha néhány minta az átlagosnál gyengébb minőségű, akkor is ezt szem előtt tartva használom fel őket a későbbi elemzésben. Másrészt, ha a kísérletben az összes minta szisztematikusan több metrikában is “warning” vagy “fail” jelzést kap (lásd ezt a példát), akkor arra gyanakszom, hogy valami rosszul sikerült a kísérletben (pl. rossz DNS-minőség, könyvtár-előkészítés stb.), és javaslom a kísérlet megismétlését.

Read trimming

A nyers leolvasások QC-je segít a problémás minták azonosításában, de nem javítja a leolvasások tényleges minőségét. Ehhez a leolvasásokat trimmelnünk kell, hogy eltávolítsuk a technikai szekvenciákat és a gyenge minőségű végeket.

A technikai szekvenciák a kísérleti eljárásból visszamaradt részek (pl. szekvenálási adapterek). Ha az ilyen szekvenciák a leolvasás valódi szekvenciájával szomszédosak, az igazítás (lásd alább) a leolvasásokat a genom rossz pozíciójához rendelheti, vagy csökkentheti az adott igazítás megbízhatóságát. A technikai szekvenciák mellett a biológiai eredetű szekvenciákat is eltávolíthatjuk, ha ezek nagymértékben jelen vannak a leolvasások között. Például a nem optimális DNS-előkészítési eljárások nagy arányban hagyhatnak a mintában DNS-konvertált riboszomális RNS-t (rRNS). Hacsak nem ez a típusú nukleinsav a szekvenálási kísérlet célpontja, az rRNS-ből származó leolvasások megtartása csak növeli a későbbi lépések számítási terhét, és megzavarhatja az eredményeket. Megjegyzendő, hogy ha a technikai szekvenciák, az rRNS vagy más szennyeződések szintje nagyon magas, amit valószínűleg a minőségellenőrzés már kiemelt, érdemes lehet a teljes szekvenálási mintát eldobni.

A rövid leolvasású szekvenálás során a DNS-szekvencia meghatározása nukleotidonként történik (technikailag minden szekvenálási ciklusban egy nukleotidot). Más szóval, a szekvenálási ciklusok száma határozza meg a leolvasás hosszát. A HTS-szekvenálási módszerek ismert problémája, hogy a nukleotidok meghatározásának pontossága a szekvenálási ciklusok halmozódásával csökken. Ez a bázisonkénti hívásminőség általános csökkenésében mutatkozik meg, különösen a leolvasás vége felé. Ahogy a technikai szekvenciák esetében is történik, a gyenge minőségű végeket tartalmazó leolvasások összehangolásának kísérlete rossz elhelyezéshez vagy gyenge leképezési minőséghez vezethet.

A technikai/szennyező szekvenciák és a gyenge minőségű végek eltávolítására léteznek és széles körben használatosak az olyan leolvasás-trimmelő eszközök, mint a Trimmomatic és a Cutadapt. Ezek az eszközök lényegében eltávolítják a (belsőleg rendelkezésre álló és/vagy a felhasználó által biztosított) technikai szekvenciákat, és a minőség alapján trimmelik a leolvasásokat, miközben maximalizálják a leolvasás hosszát. A trimmelés után túl rövidnek bizonyuló olvasatokat elvetik (a túl rövid olvasatok, pl. <36 nukleotid, megnehezítik az illesztési lépést, mivel ezek valószínűleg a genom több helyére is leképeződnek). Érdemes megnézni a trimmelést túlélő olvasatok százalékos arányát, mivel az elvetett olvasatok magas aránya valószínűleg rossz minőségű adatokra utal.

Végül általában újra lefuttatom a FastQC-t a levágott olvasatokon, hogy ellenőrizzem, hogy ez a lépés hatékony volt-e és szisztematikusan javította-e a QC-metrikákat.

Alignment

A kivételektől eltekintve (pl. de novo assembly), a legtöbb HTS-adattípus és alkalmazás esetében jellemzően az alignment (más néven mapping) a következő lépés. A leolvasás összehangolása abból áll, hogy meghatározzuk a genomban azt a pozíciót, ahonnan a leolvasás szekvenciája származik (jellemzően kromoszóma:start-end-ként kifejezve). Ennél a lépésnél tehát szükség van egy referenciaszekvencia használatára a leolvasások illesztéséhez/térképezéséhez.

A referenciaszekvencia kiválasztását több tényező határozza meg. Egyrészt a faj, amelyből a szekvenált DNS származik. Bár egyre több olyan fajról áll rendelkezésre jó minőségű referenciaszekvencia, néhány kevésbé vizsgált szervezet esetében ez még mindig nem biztos, hogy így van. Ezekben az esetekben érdemes a leolvasásokat olyan evolúciósan közeli fajhoz igazítani, amelyhez referencia genom áll rendelkezésre. Mivel például a prérifarkas genomjához nem áll rendelkezésre referenciaszekvencia, a leolvasások összehangolásához a közeli rokon kutya genomját használhatjuk. Hasonlóképpen előfordulhat, hogy olvasatainkat egy olyan közeli rokon fajhoz szeretnénk igazítani, amelyre létezik jobb minőségű referenciaszekvencia. Például, bár a gibbon genomját közzétették, ez több ezer töredékre van bontva, amelyek nem teljesen reprodukálják az adott genom több tíz kromoszómára való szerveződését; ebben az esetben előnyös lehet, ha az igazítást az emberi referenciaszekvenciával végezzük el.

Egy másik figyelembe veendő tényező a referenciaszekvencia-összeállítás verziója, mivel a szekvencia frissítése és javítása során új verziók jelennek meg. Fontos, hogy egy adott illesztés koordinátái változhatnak a különböző verziók között. Például a humán genom több változata is megtalálható a UCSC Genome Browserben. Bármely faj esetében erősen javaslom, hogy a legújabb összeállítási verzióra térjünk át, amint az teljes egészében megjelent. Ez okozhat némi kellemetlenséget az átmenet során, mivel a már meglévő eredmények a régebbi verziókhoz lesznek viszonyítva, de hosszú távon kifizetődő.

Ez mellett a szekvenálási adatok típusa is számít. A DNS-seq, ChIP-seq vagy Hi-C protokollokból generált leolvasásokat a genom referenciaszekvenciájához igazítják. Másfelől, mivel a DNS-ből átírt RNS-t mRNS-é dolgozzák tovább (azaz az intronokat eltávolítják), sok RNS-seq leolvasás nem fog igazodni a genom referenciaszekvenciájához. Ehelyett vagy transzkriptom-referenciaszekvenciákhoz kell igazítanunk őket, vagy a genomszekvenciát referenciaként használva split-érzékeny igazítókat kell használnunk (lásd alább). Ehhez kapcsolódik a referenciaszekvencia annotációjának forrása, azaz a gének, transzkriptek, centromerek stb. koordinátáit tartalmazó adatbázis kiválasztása. Én jellemzően a GENCODE annotációt használom, mivel az átfogó génannotációt és a transzkriptek szekvenciáit egyesíti.

A rövid olvasatú szekvenciaillesztési eszközök hosszú listáját fejlesztették ki (lásd a Rövid olvasatú szekvenciaillesztés című részt itt). Ezek áttekintése meghaladja e cikk kereteit (az ezen eszközök mögött álló algoritmusok részletei itt találhatók). Tapasztalataim szerint a legnépszerűbbek közé tartozik a Bowtie2, a BWA, a HISAT2, a Minimap2, a STAR és a TopHat. Azt javaslom, hogy az igazítót olyan kulcsfontosságú tényezők figyelembevételével válassza ki, mint a HTS-adatok és az alkalmazás típusa, valamint a közösség általi elfogadottság, a dokumentáció minősége és a felhasználók száma. Például olyan illesztőprogramokra van szükség, mint a STAR vagy a Bowtie2, amelyek tisztában vannak az exon-exon csomópontokkal, amikor az RNS-seq-et a genomra illesztik.

A legtöbb mappernek közös jellemzője, hogy a tényleges illesztés előtt indexelni kell a referenciaként használt szekvenciát. Ez a lépés időigényes lehet, de csak egyszer kell elvégezni minden egyes referenciaszekvencia esetében. A legtöbb leképező SAM/BAM fájlokban tárolja az illesztéseket, amelyek a SAM/BAM formátumot követik (a BAM fájlok a SAM fájlok bináris változatai). Az illesztés a szekvenálási adatok elemzésének legszámításigényesebb és legidőigényesebb lépései közé tartozik, és a SAM/BAM fájlok nagy tömegűek (gigabájt nagyságrendűek). Ezért fontos, hogy biztosítsuk a szükséges erőforrásokat (lásd az alábbi utolsó szakaszt) az igazítás ésszerű időn belüli lefuttatásához és az eredmények tárolásához. Hasonlóképpen, a BAM-fájlok mérete és bináris formátuma miatt kerülje a szövegszerkesztőkkel való megnyitást; ehelyett használjon Unix-parancsokat vagy dedikált eszközöket, például a SAMtools-t.

Az igazításokból

Az igazítás után azt mondanám, hogy nincs egyértelmű közös lépés, mivel ezen a ponton az egyes HTS-adattípusok és alkalmazások különbözhetnek.

A DNS-seq-adatok közös utólagos elemzése a variánshívás, azaz a genomban azon pozíciók azonosítása, amelyek a genom-referenciához képest és az egyének között eltérnek. Egy népszerű elemzési keretrendszer erre az alkalmazásra a GATK az egynukleotid-polimorfizmus (SNP) vagy a kis inszerciók/deléciók (indelek) esetében (2. ábra). A nagyobb DNS-tömböket tartalmazó variánsok (más néven strukturális variánsok) külön erre a célra kifejlesztett hívó módszereket igényelnek (átfogó összehasonlításért lásd ezt a cikket). Az igazítókhoz hasonlóan a megfelelő eszköz kiválasztásához is azt tanácsolom, hogy olyan kulcsfontosságú tényezőket vegyünk figyelembe, mint a változatok fajtája (SNP, indel vagy strukturális változatok), a közösség általi elfogadottság, a dokumentáció minősége és a felhasználók száma.

Az RNA-seq talán leggyakoribb alkalmazása a génexpresszió számszerűsítése. Történelmileg a leolvasásokat a referenciaszekvenciához kellett igazítani, majd az adott génhez vagy transzkripthez igazított leolvasások számát helyettesítő értékként használták az expressziós szintek számszerűsítésére. Ezt az összehangolás+kvantitatív megközelítést olyan eszközökkel végzik, mint a Cufflinks, az RSEM vagy a featureCounts. A scuh megközelítést azonban egyre inkább felülmúlják az olyan szoftverekben megvalósított újabb módszerek, mint a Kallisto és a Salmon. Koncepcionálisan az ilyen eszközökkel a leolvasás teljes szekvenciáját nem kell a referenciaszekvenciához igazítani. Ehelyett csak annyi nukleotidot kell összehangolnunk, hogy biztosak lehessünk abban, hogy a leolvasás egy adott transzkriptumból származik. Egyszerűen fogalmazva, az összehangolás+kvantitatív megközelítés egyetlen lépésre redukálódik. Ezt a megközelítést pszeudo-térképezésnek nevezik, és nagymértékben növeli a génexpresszió mennyiségi meghatározásának sebességét. Másfelől ne feledjük, hogy a pszeudo-térképezés nem alkalmas olyan alkalmazásokhoz, ahol teljes igazításra van szükség (pl. variánsmeghatározás RNS-seq adatokból).

Egy másik példa a downstream elemzési lépések és a szükséges eszközök közötti különbségekre a szekvenáláson alapuló alkalmazások között a ChIP-seq. Az ilyen technikával generált leolvasásokat csúcsmeghatározásra használjuk, amely a genom azon régióinak felderítéséből áll, amelyekben jelentős többlet leolvasás van, ami azt jelzi, hogy a célfehérje hol kötődik. Számos peak caller létezik, és ez a kiadvány áttekinti őket. Utolsó példaként megemlítem a Hi-C adatokat, amelyekben az összehangolásokat olyan eszközök bemeneteként használják, amelyek meghatározzák a kölcsönhatási mátrixokat, és ezekből a genom 3D-jellemzőit. Az összes szekvenáláson alapuló vizsgálat kommentálása meghaladja e cikk kereteit (a viszonylag teljes listát lásd ebben a cikkben).

Előtte…

A cikk hátralévő része olyan szempontokat érint, amelyeket nem feltétlenül tekintünk szigorúan a HTS-adatok elemzésének lépéseinek, és amelyeket nagyrészt figyelmen kívül hagyunk. Ezzel szemben én amellett érvelek, hogy az 1. táblázatban feltett kérdéseken el kell gondolkodni, mielőtt a HTS-adatok (vagy valójában bármilyen adat) elemzésébe kezdünk, és ezekről a témákról itt és itt írtam.

1. táblázat

Gondolkodj el rajta	A javasolt intézkedés
Megvan a mintádról az elemzéshez szükséges összes információ?	Szisztematikusan gyűjtse össze a kísérletek metaadatait
Egyértelműen azonosítani tudja majd a mintáját?	Szervezzen egy olyan rendszert, amely minden mintához egyedi azonosítót rendel
Hol lesznek az adatok és az eredmények?	Az adatok strukturált és hierarchikus szervezése
A több minta zökkenőmentes feldolgozására lesz képes?	A kód skálázhatósága, párhuzamosítása, automatikus konfigurálása és modularitása
Az eredményeket Ön vagy bárki más képes lesz reprodukálni?	Dokumentálja a kódját és az eljárásait!

Amint fentebb említettük, a HTS nyers adatok és az elemzésük során keletkező fájlok egy része gigabájtos nagyságrendű, így nem kivételes, hogy egy több tíz mintát tartalmazó projekt terabájtnyi tárhelyet igényel. Emellett a HTS-adatok elemzésének egyes lépései számításigényesek (pl. az igazítás). A HTS-adatok elemzéséhez szükséges tárolási és számítási infrastruktúra azonban fontos szempont, és gyakran figyelmen kívül hagyják vagy nem tárgyalják. Egy nemrégiben végzett elemzés részeként például több tíz publikált, fenoméra kiterjedő asszociációs elemzést (PheWAS) végző cikket tekintettünk át. A modern PheWAS-ok 100-1000 genetikai variánst és fenotípust elemeznek, ami jelentős adattárolást és számítási teljesítményt eredményez. Ennek ellenére az általunk áttekintett cikkek közül gyakorlatilag egyik sem foglalkozott a PheWAS-elemzéshez szükséges infrastruktúrával. Nem meglepő módon az a javaslatom, hogy előre tervezze meg a felmerülő tárolási és számítástechnikai követelményeket, és ossza meg azokat a közösséggel.

Segítségre van szüksége a DNS-szekvenálási adatok elemzéséhez? Vegye fel a kapcsolatot szabadúszó bioinformatikusokkal és genomikai szakértőkkel a Kolabtree oldalon.