Průvodce analýzou dat sekvenování DNA krok za krokem

Dr. Javier Quilez Oliete, zkušený bioinformatický konzultant na Kolabtree, poskytuje komplexního průvodce analýzou dat sekvenování DNA, včetně nástrojů a softwaru používaných ke čtení dat.

Úvod

Deoxyribonukleová kyselina (DNA) je molekula, která nese většinu genetické informace organismu. (U některých typů virů je nositelem genetické informace ribonukleová kyselina (RNA). Základními jednotkami molekul DNA jsou nukleotidy (obvykle označované písmeny A, C, G nebo T). Z koncepčního hlediska je sekvenování DNA proces čtení nukleotidů, které tvoří molekulu DNA (např. „GCAAACCAAT“ je řetězec 10 nukleotidů DNA). Současné technologie sekvenování produkují miliony takových čtení DNA v přiměřeném čase a za relativně nízkou cenu. Pro srovnání, náklady na sekvenování lidského genomu – genom je kompletní soubor molekul DNA v organismu – klesly pod hranici 100 dolarů a lze je provést během několika dnů. To kontrastuje s první iniciativou na sekvenování lidského genomu, která byla dokončena za deset let a jejíž náklady činily přibližně 2,7 miliardy dolarů.

Tato schopnost sekvenovat DNA s vysokou výkonností a nízkými náklady umožnila vývoj rostoucího počtu metod a aplikací založených na sekvenování. Například sekvenování celých genomů nebo jejich protein kódujících oblastí (dva přístupy známé jako sekvenování celého genomu, resp. exomu) u nemocných i zdravých jedinců může napovědět o změnách DNA způsobujících onemocnění. Také sekvenování RNA, která je přepisována z DNA – technika známá jako sekvenování RNA – se používá ke kvantifikaci genové aktivity a její změny za různých podmínek (např. neléčená versus léčená). Na druhé straně metody sekvenování zachycující konformaci chromozomů zjišťují interakce mezi blízkými molekulami DNA a pomáhají tak určit prostorové rozmístění chromozomů v buňce.

Společným znakem těchto a dalších aplikací sekvenování DNA je generování souborů dat v řádu gigabajtů a zahrnujících miliony přečtených sekvencí. Proto dávání smyslu experimentům s vysokokapacitním sekvenováním (HTS) vyžaduje značné možnosti analýzy dat. Naštěstí pro většinu typů dat HTS existují specializované výpočetní a statistické nástroje a relativně standardní pracovní postupy analýzy. Zatímco některé (počáteční) kroky analýzy jsou společné pro většinu typů sekvenačních dat, další navazující analýzy budou záviset na druhu dat a/nebo konečném cíli analýzy. Níže uvádím základní informace o základních krocích analýzy dat HTS a odkazuji na populární nástroje.

Některé z níže uvedených oddílů jsou zaměřeny na analýzu dat generovaných technologiemi sekvenování s krátkými čteními (většinou Illumina), protože ty historicky dominují trhu HTS. Rychle se však prosazují novější technologie, které generují delší čtení (např. Oxford Nanopore Technologies, PacBio). Vzhledem k tomu, že sekvenování s dlouhými čteními má některá specifika (např. vyšší chybovost), jsou pro analýzu tohoto druhu dat vyvíjeny specifické nástroje.

Kontrola kvality (QC) surových čtení

Dychtivý analytik začne analýzu ze souborů FASTQ; formát FASTQ je již dlouho standardem pro ukládání dat sekvenování s krátkými čteními. Soubory FASTQ v podstatě obsahují sekvenci nukleotidů a kvalitu volání na bázi pro miliony čtení. Ačkoli velikost souboru závisí na skutečném počtu čtení, soubory FASTQ jsou obvykle velké (v řádu megabajtů a gigabajtů) a komprimované. Je třeba poznamenat, že většina nástrojů, které používají soubory FASTQ jako vstup, je dokáže zpracovat v komprimovaném formátu, takže v zájmu úspory místa na disku se doporučuje je nerozbalovat. Konvenčně zde budu soubor FASTQ přirovnávat k sekvenačnímu vzorku.

FastQC je pravděpodobně nejoblíbenějším nástrojem pro provádění kontroly kvality surových čtení. Lze jej spustit prostřednictvím vizuálního rozhraní nebo programově. Zatímco první možnost může být pohodlnější pro uživatele, kteří se necítí dobře v prostředí příkazového řádku, druhá možnost nabízí nesrovnatelnou škálovatelnost a reprodukovatelnost (vzpomeňte si, jak zdlouhavé a náchylné k chybám může být ruční spuštění nástroje pro desítky souborů). Ať tak či onak, hlavním výstupem nástroje FastQC je soubor HTML uvádějící klíčové souhrnné statistiky o celkové kvalitě surových sekvenačních čtení z daného vzorku. Prohlížet desítky zpráv FastQC jednu po druhé je zdlouhavé a komplikuje to porovnávání jednotlivých vzorků. Proto můžete použít MultiQC, který shrnuje zprávy HTML z FastQC (a také z dalších nástrojů použitých následně, např. ořezávání adaptérů, zarovnávání) do jediné zprávy.

Informace o kvalitě vzorků jsou určeny k tomu, aby uživatel mohl posoudit, zda mají vzorky dobrou kvalitu a lze je proto použít pro další kroky, nebo je třeba je vyřadit. Bohužel neexistuje konsenzuální prahová hodnota založená na metrikách FastQC pro klasifikaci vzorků jako vzorků dobré nebo špatné kvality. Přístup, který používám, je následující. Očekávám, že všechny vzorky, které prošly stejným postupem (např. extrakce DNA, příprava knihovny), budou mít podobné statistiky kvality a většinu příznaků „vyhovuje“. Pokud mají některé vzorky nižší než průměrnou kvalitu, přesto je použiji v následné analýze s ohledem na tuto skutečnost. Na druhou stranu, pokud všechny vzorky v experimentu systematicky dostávají příznaky „varování“ nebo „neprošel“ ve více metrikách (viz tento příklad), mám podezření, že se v experimentu něco pokazilo (např. špatná kvalita DNA, příprava knihovny atd.), a doporučuji jej opakovat.

Ořezávání čtení

QC surových čtení pomáhá identifikovat problematické vzorky, ale nezlepšuje skutečnou kvalitu čtení. K tomu potřebujeme ořezat čtení, abychom odstranili technické sekvence a nekvalitní konce.

Technické sekvence jsou zbytky z experimentálního postupu (např. sekvenační adaptéry). Pokud takové sekvence sousedí se skutečnou sekvencí čtení, může zarovnání (viz níže) mapovat čtení na nesprávnou pozici v genomu nebo snížit důvěryhodnost daného zarovnání. Kromě technických sekvencí můžeme také chtít odstranit sekvence biologického původu, pokud jsou mezi čteními silně zastoupeny. Například neoptimální postupy přípravy DNA mohou ve vzorku zanechat vysoký podíl ribozomální RNA (rRNA) převedené na DNA. Pokud není tento typ nukleové kyseliny cílem sekvenačního experimentu, ponechání čtení odvozených od rRNA pouze zvýší výpočetní zátěž navazujících kroků a může zmást výsledky. Za zmínku stojí, že pokud jsou hladiny technických sekvencí, rRNA nebo jiných kontaminantů velmi vysoké, na což již pravděpodobně upozornila kontrola kvality, můžete chtít celý sekvenační vzorek vyřadit.

Při sekvenování s krátkými čteními se sekvence DNA určuje po jednom nukleotidu (technicky vzato, jeden nukleotid každý sekvenační cyklus). Jinými slovy, počet sekvenačních cyklů určuje délku čtení. Známým problémem metod sekvenování HTS je pokles přesnosti, s níž jsou určovány nukleotidy, s přibývajícími sekvenačními cykly. To se projevuje celkovým poklesem kvality volání na bázi, zejména ke konci čtení. Stejně jako se to děje u technických sekvencí, může snaha o zarovnání čtení, která obsahují nekvalitní konce, vést k nesprávnému umístění nebo špatné kvalitě mapování.

Pro odstranění technických/kontaminačních sekvencí a nekvalitních konců existují a jsou široce používány nástroje pro ořezávání čtení, jako jsou Trimmomatic a Cutadapt. Tyto nástroje v podstatě odstraní technické sekvence (interně dostupné a/nebo poskytnuté uživatelem) a oříznou čtení na základě kvality při maximalizaci délky čtení. Čtení, která po ořezání zůstanou příliš krátká, jsou vyřazena (příliš krátká čtení, např. <36 nukleotidů, komplikují krok zarovnání, protože budou pravděpodobně mapována na více míst v genomu). Můžete se podívat na procento čtení, která přežila ořezávání, protože vysoký podíl vyřazených čtení je pravděpodobně známkou špatné kvality dat.

Nakonec obvykle znovu spustím FastQC na ořezaných čteních, abych zkontroloval, zda byl tento krok účinný a systematicky zlepšil metriky QC.

Zarovnání

Až na výjimky (např. sestavení de novo) je zarovnání (označované také jako mapování) obvykle dalším krokem pro většinu typů dat a aplikací HTS. Zarovnání čtení spočívá v určení pozice v genomu, ze které pochází sekvence čtení (obvykle vyjádřeno jako chromozom:start-end). Proto v tomto kroku vyžadujeme použití referenční sekvence, na kterou zarovnáme/mapujeme čtení.

Výběr referenční sekvence bude určen více faktory. Jedním z nich je druh, z něhož sekvenovaná DNA pochází. Zatímco počet druhů, které mají k dispozici kvalitní referenční sekvenci, se zvyšuje, u některých méně prozkoumaných organismů tomu tak stále být nemusí. V těchto případech můžete chtít zarovnat čtení k evolučně blízkému druhu, pro který je k dispozici referenční genom. Protože například není k dispozici referenční sekvence pro genom kojota, můžeme pro zarovnání čtení použít sekvenci blízce příbuzného psa. Podobně můžeme ještě chtít zarovnat naše čtení k blízce příbuznému druhu, pro který existuje referenční sekvence vyšší kvality. Například byl sice publikován genom gibona, ten je však rozdělen na tisíce fragmentů, které plně nerekapitulují uspořádání tohoto genomu do desítek chromozomů; v takovém případě může být výhodné provést zarovnání pomocí referenční sekvence člověka.

Dalším faktorem, který je třeba vzít v úvahu, je verze sestavy referenční sekvence, protože s aktualizací a vylepšováním sekvence jsou vydávány nové verze. Důležité je, že souřadnice daného zarovnání se mohou mezi verzemi lišit. Například v prohlížeči genomu UCSC lze nalézt více verzí lidského genomu. U všech druhů rozhodně upřednostňuji přechod na nejnovější verzi sestavení, jakmile je plně uvolněna. Během přechodu to může způsobit určité nepříjemnosti, protože již existující výsledky budou vztaženy ke starším verzím, ale z dlouhodobého hlediska se to vyplatí.

Kromě toho záleží také na typu sekvenačních dat. Čtení generovaná z protokolů DNA-seq, ChIP-seq nebo Hi-C budou zarovnána na referenční sekvenci genomu. Na druhou stranu, protože RNA přepsaná z DNA je dále zpracována na mRNA (tj. odstraněny introny), mnoho čtení RNA-seq se nepodaří zarovnat na referenční sekvenci genomu. Místo toho je musíme buď zarovnat k referenčním sekvencím transkriptomu, nebo použít zarovnávače zohledňující rozdělení (viz níže), pokud jako referenční sekvenci použijeme sekvenci genomu. S tím souvisí i volba zdroje pro anotaci referenční sekvence, tedy databáze se souřadnicemi genů, transkriptů, centromer atd. Obvykle používám anotaci GENCODE, protože kombinuje komplexní anotaci genů a sekvencí transkriptů.

Byl vyvinut dlouhý seznam nástrojů pro zarovnávání sekvencí s krátkým čtením (viz část Zarovnávání sekvencí s krátkým čtením zde). Jejich přehled přesahuje rámec tohoto článku (podrobnosti o algoritmech těchto nástrojů naleznete zde). Podle mých zkušeností patří mezi nejoblíbenější Bowtie2, BWA, HISAT2, Minimap2, STAR a TopHat. Doporučuji, abyste si zarovnávač vybrali na základě zvážení klíčových faktorů, jako je typ HTS dat a aplikace, jakož i přijetí komunitou, kvalita dokumentace a počet uživatelů. Např. je třeba použít srovnávače jako STAR nebo Bowtie2, které si při mapování RNA-seq na genom uvědomují exon-exonové spoje.

Společná pro většinu mapovačů je potřeba indexovat sekvenci použitou jako referenční před vlastním zarovnáním. Tento krok může být časově náročný, ale pro každou referenční sekvenci jej stačí provést pouze jednou. Většina mapovačů ukládá zarovnání do souborů SAM/BAM, které se řídí formátem SAM/BAM (soubory BAM jsou binární verze souborů SAM). Zarovnání patří mezi výpočetně a časově nejnáročnější kroky při analýze sekvenačních dat a soubory SAM/BAM jsou těžké (řádově gigabajty). Proto je důležité se ujistit, že máte k dispozici potřebné zdroje (viz závěrečná část níže), abyste mohli zarovnání provést v rozumném čase a výsledky uložit. Stejně tak se vzhledem k velikosti a binárnímu formátu souborů BAM vyhněte jejich otevírání pomocí textových editorů; místo toho použijte unixové příkazy nebo specializované nástroje, jako je SAMtools.

Zarovnání

Řekl bych, že po zarovnání neexistuje jednoznačný společný krok, protože v tomto bodě se může každý typ dat HTS a aplikace lišit.

Běžnou následnou analýzou pro data DNA-seq je volání variant, tj. identifikace pozic v genomu, které se liší vzhledem k referenčnímu genomu a mezi jednotlivci. Oblíbeným analytickým rámcem pro tuto aplikaci je GATK pro jednonukleotidový polymorfismus (SNP) nebo malé inzerce/delece (indely) (obr. 2). Varianty obsahující větší kusy DNA (označované také jako strukturální varianty) vyžadují specializované metody volání (komplexní srovnání viz tento článek). Stejně jako u srovnávačů doporučuji vybrat správný nástroj s ohledem na klíčové faktory, jako je druh variant (SNP, indely nebo strukturální varianty), přijetí komunitou, kvalita dokumentace a počet uživatelů.

Pravděpodobně nejčastější aplikací RNA-seq je kvantifikace genové exprese. Historicky bylo třeba čtení zarovnat k referenční sekvenci a poté se počet čtení zarovnaných k danému genu nebo transkriptu používal jako zástupce pro kvantifikaci úrovně jeho exprese. Tento přístup zarovnání + kvantifikace se provádí pomocí nástrojů jako Cufflinks, RSEM nebo featureCounts. Přístup scuh je však stále více překonáván novějšími metodami implementovanými v softwaru jako Kallisto a Salmon. Z koncepčního hlediska není u těchto nástrojů nutné zarovnávat celou sekvenci čtení k referenční sekvenci. Místo toho stačí zarovnat pouze tolik nukleotidů, abychom měli jistotu, že čtení pochází z daného transkriptu. Zjednodušeně řečeno, přístup zarovnání + kvantifikace je omezen na jediný krok. Tento přístup je znám jako pseudomapování a výrazně zvyšuje rychlost kvantifikace genové exprese. Na druhou stranu mějte na paměti, že pseudomapování nebude vhodné pro aplikace, kde je potřeba úplné zarovnání (např. volání variant z dat RNA-seq).

Dalším příkladem rozdílů v krocích následné analýzy a potřebných nástrojích napříč aplikacemi založenými na sekvenování je ChIP-seq. Čtení generovaná touto technikou budou použita pro tzv. peak calling, který spočívá v detekci oblastí v genomu s výrazným přebytkem čtení, což naznačuje, kde je cílový protein vázán. Existuje několik programů pro volání vrcholů a tato publikace se jimi zabývá. Jako poslední příklad uvedu data Hi-C, v nichž se zarovnání používají jako vstup pro nástroje, které určují interakční matice a z nich 3D-charakteristiky genomu. Komentování všech testů založených na sekvenování přesahuje rámec tohoto článku (relativně úplný seznam viz tento článek).

Než začnete…

Zbývající část tohoto článku se dotýká aspektů, které nemusí být striktně považovány za kroky při analýze dat HTS a které jsou z velké části ignorovány. Naopak tvrdím, že je kapitálem, abyste se zamysleli nad otázkami položenými v tabulce 1, než začnete analyzovat údaje HTS (nebo skutečně jakýkoli druh údajů), a o těchto tématech jsem psal zde a zde.

Tabulka 1

Přemýšlejte o tom	Návrh postupu
Máte všechny informace o svém vzorku potřebné pro analýzu?	Shromažďujte systematicky metadata o experimentech
Budete schopni jednoznačně identifikovat svůj vzorek	Vytvořte systém, který každému vzorku přiřadí jedinečný identifikátor
Kde budou data a výsledky?	Strukturované a hierarchické uspořádání dat
Budete moci bezproblémově zpracovávat více vzorků	Škálovatelnost, paralelizace, automatická konfigurace a modularita kódu
Budete moci vy nebo kdokoli jiný reprodukovat výsledky?	Dokumentujte svůj kód a postupy!

Jak bylo uvedeno výše, surová data HTS a některé soubory generované při jejich analýze mají řádově gigabajty, takže není výjimečné, že projekt zahrnující desítky vzorků vyžaduje terabajty úložného prostoru. Kromě toho jsou některé kroky při analýze dat HTS výpočetně náročné (např. zarovnání). Úložná a výpočetní infrastruktura potřebná pro analýzu dat HTS je však důležitým faktorem, který se často přehlíží nebo se o něm nemluví. Jako příklad lze uvést, že v rámci nedávné analýzy jsme přezkoumali desítky publikovaných prací provádějících fenome-wide association analysis (PheWAS). Moderní PheWAS analyzují 100-1 000 genetických variant i fenotypů, což má za následek důležité ukládání dat a výpočetní výkon. A přesto se prakticky žádná z prací, které jsme přezkoumali, nevyjadřovala k infrastruktuře potřebné pro analýzu PheWAS. Není divu, že doporučuji, abyste předem naplánovali požadavky na úložiště a výpočetní techniku, s nimiž se budete potýkat, a podělili se o ně s komunitou.

Potřebujete pomoc s analýzou dat sekvenování DNA? Spojte se s bioinformatiky a odborníky na genomiku na volné noze na webu Kolabtree.