Przewodnik krok po kroku po analizie danych z sekwencjonowania DNA

Dr Javier Quilez Oliete, doświadczony konsultant ds. bioinformatyki na Kolabtree, przedstawia kompleksowy przewodnik po analizie danych z sekwencjonowania DNA, w tym narzędzia i oprogramowanie używane do odczytu danych.

Wprowadzenie

Kwas dezoksyrybonukleinowy (DNA) jest cząsteczką, która przenosi większość informacji genetycznej organizmu. (W niektórych typach wirusów, informacja genetyczna jest przenoszona przez kwas rybonukleinowy (RNA)). Nukleotydy (umownie reprezentowane przez litery A, C, G lub T) są podstawowymi jednostkami cząsteczek DNA. Koncepcyjnie, sekwencjonowanie DNA jest procesem odczytywania nukleotydów składających się na cząsteczkę DNA (np. „GCAAACCAAT” jest 10-nukleotydowym ciągiem DNA). Obecne technologie sekwencjonowania wytwarzają miliony takich odczytów DNA w rozsądnym czasie i przy stosunkowo niskich kosztach. Dla porównania, koszt sekwencjonowania ludzkiego genomu – genom to kompletny zestaw cząsteczek DNA w organizmie – spadł poniżej 100 dolarów i można to zrobić w ciągu kilku dni. To kontrastuje z pierwszą inicjatywą sekwencjonowania ludzkiego genomu, która została zakończona w ciągu dekady i miała koszt około 2,7 miliardów dolarów.

Ta zdolność do sekwencjonowania DNA przy wysokiej przepustowości i niskich kosztach umożliwiła rozwój rosnącej liczby metod i zastosowań opartych na sekwencjonowaniu. Na przykład, sekwencjonowanie całych genomów lub ich regionów kodujących białka (dwa podejścia znane odpowiednio jako sekwencjonowanie całego genomu i sekwencjonowanie egzomu) u osób chorych i zdrowych może podpowiedzieć chorobotwórcze zmiany DNA. Również sekwencjonowanie RNA, które jest przepisywane z DNA – technika znana jako sekwencjonowanie RNA – jest używana do ilościowego określenia aktywności genów i jej zmian w różnych warunkach (np. nieleczone kontra leczone). Z drugiej strony, metody sekwencjonowania wychwytu konformacji chromosomów wykrywają interakcje pomiędzy pobliskimi cząsteczkami DNA i w ten sposób pomagają określić przestrzenne rozmieszczenie chromosomów w komórce.

Wspólną cechą tych i innych zastosowań sekwencjonowania DNA jest generowanie zbiorów danych rzędu gigabajtów i zawierających miliony odczytanych sekwencji. Dlatego też, nadanie sensu eksperymentom sekwencjonowania o wysokiej przepustowości (HTS) wymaga znacznych możliwości analizy danych. Na szczęście, dla większości typów danych HTS istnieją dedykowane narzędzia obliczeniowe i statystyczne oraz stosunkowo standardowe procesy analizy. Podczas gdy niektóre z (początkowych) kroków analizy są wspólne dla większości typów danych sekwencjonowania, dalsze etapy analizy będą zależały od rodzaju danych i/lub ostatecznego celu analizy. Poniżej przedstawiam elementarz podstawowych kroków analizy danych HTS i odsyłam do popularnych narzędzi.

Niektóre z poniższych sekcji koncentrują się na analizie danych wygenerowanych za pomocą technologii sekwencjonowania krótkich odczytów (głównie Illumina), ponieważ to one historycznie zdominowały rynek HTS. Jednakże, nowsze technologie generujące dłuższe odczyty (np. Oxford Nanopore Technologies, PacBio) szybko zyskują na popularności. Ponieważ sekwencjonowanie długich odczytów ma pewne cechy szczególne (np. wyższy wskaźnik błędów), opracowywane są specjalne narzędzia do analizy tego rodzaju danych.

Kontrola jakości (QC) surowych odczytów

Niecierpliwy analityk rozpocznie analizę od plików FASTQ; format FASTQ jest od dawna standardem przechowywania danych sekwencjonowania krótkich odczytów. W istocie, pliki FASTQ zawierają sekwencję nukleotydów i jakość wywoływania per-baza dla milionów odczytów. Chociaż rozmiar pliku zależy od rzeczywistej liczby odczytów, pliki FASTQ są zazwyczaj duże (rzędu megabajtów i gigabajtów) i skompresowane. Warto zauważyć, że większość narzędzi, które używają plików FASTQ jako danych wejściowych, może obsługiwać je w skompresowanym formacie, więc w celu zaoszczędzenia miejsca na dysku zaleca się ich nie dekompresować. Jako konwencję będę tu utożsamiał plik FASTQ z próbką sekwencjonowania.

FastQC jest prawdopodobnie najpopularniejszym narzędziem do przeprowadzania QC surowych odczytów. Można go uruchomić za pomocą interfejsu wizualnego lub programowo. Podczas gdy pierwsza opcja może być wygodniejsza dla użytkowników, którzy nie czują się komfortowo w środowisku wiersza poleceń, druga oferuje nieporównywalną skalowalność i odtwarzalność (pomyślmy, jak żmudne i podatne na błędy może być ręczne uruchamianie narzędzia dla dziesiątek plików). Tak czy inaczej, głównym wynikiem działania FastQC jest plik HTML, który przedstawia kluczowe statystyki podsumowujące ogólną jakość surowych odczytów sekwencjonowania z danej próbki. Przeglądanie dziesiątek raportów FastQC jeden po drugim jest żmudne i komplikuje porównywanie próbek. Dlatego warto skorzystać z MultiQC, który agreguje raporty HTML z FastQC (jak również z innych narzędzi używanych w dalszej kolejności, np. przycinanie adapterów, wyrównywanie) w jeden raport.

Informacje o jakości QC mają pozwolić użytkownikowi ocenić, czy próbki mają dobrą jakość i mogą być w związku z tym użyte do dalszych kroków, czy też należy je odrzucić. Niestety, nie istnieje próg konsensusu oparty na metryce FastQC, który pozwoliłby zaklasyfikować próbki jako dobrej lub złej jakości. Podejście, którego używam jest następujące. Oczekuję, że wszystkie próbki, które przeszły przez tę samą procedurę (np. ekstrakcja DNA, przygotowanie biblioteki) będą miały podobną statystykę jakości i większość flag „pass”. Jeżeli niektóre próbki mają jakość niższą niż przeciętna, to pamiętając o tym, nadal będę ich używał w dalszych analizach. Z drugiej strony, jeśli wszystkie próbki w eksperymencie systematycznie otrzymują flagi „warning” lub „fail” w wielu metrykach (zobacz ten przykład), podejrzewam, że coś poszło nie tak w eksperymencie (np. zła jakość DNA, przygotowanie biblioteki, itp.) i zalecam jego powtórzenie.

Przycinanie odczytów

QC surowych odczytów pomaga zidentyfikować problematyczne próbki, ale nie poprawia rzeczywistej jakości odczytów. W tym celu musimy przyciąć odczyty, aby usunąć sekwencje techniczne i końce niskiej jakości.

Sekwencje techniczne są pozostałościami po procedurze eksperymentalnej (np. adaptery sekwencjonowania). Jeśli takie sekwencje sąsiadują z prawdziwą sekwencją odczytu, wyrównanie (patrz niżej) może zmapować odczyty do niewłaściwej pozycji w genomie lub zmniejszyć zaufanie do danego wyrównania. Poza sekwencjami technicznymi, możemy również chcieć usunąć sekwencje pochodzenia biologicznego, jeśli są one silnie obecne wśród odczytów. Na przykład, nieoptymalne procedury przygotowania DNA mogą pozostawić w próbce wysoki odsetek rybosomalnego RNA (rRNA) przekształconego w DNA. O ile ten rodzaj kwasu nukleinowego nie jest celem eksperymentu sekwencjonowania, zachowanie odczytów pochodzących z rRNA tylko zwiększy obciążenie obliczeniowe dalszych etapów i może zafałszować wyniki. Należy zauważyć, że jeśli poziomy sekwencji technicznych, rRNA lub innych zanieczyszczeń są bardzo wysokie, co prawdopodobnie zostało już podkreślone przez QC, można odrzucić całą próbkę sekwencjonowania.

W sekwencjonowaniu krótkich odczytów sekwencja DNA jest określana po jednym nukleotydzie na raz (technicznie, po jednym nukleotydzie w każdym cyklu sekwencjonowania). Innymi słowy, liczba cykli sekwencjonowania określa długość odczytu. Znanym problemem metod sekwencjonowania HTS jest spadek dokładności, z jaką oznaczane są nukleotydy, w miarę narastania liczby cykli sekwencjonowania. Znajduje to odzwierciedlenie w ogólnym spadku jakości wywołania per-baza, szczególnie pod koniec odczytu. Podobnie jak w przypadku sekwencji technicznych, próba wyrównania odczytów zawierających końce niskiej jakości może prowadzić do błędnego rozmieszczenia lub niskiej jakości mapowania.

Aby usunąć sekwencje techniczne/zanieczyszczające i końce niskiej jakości, istnieją i są szeroko stosowane narzędzia do przycinania odczytów, takie jak Trimmomatic i Cutadapt. W istocie, narzędzia te usuwają sekwencje techniczne (dostępne wewnętrznie i/lub dostarczone przez użytkownika) i przycinają odczyty na podstawie jakości, maksymalizując ich długość. Odczyty, które są zbyt krótkie po przycięciu są odrzucane (odczyty zbyt krótkie, np. <36 nukleotydów, komplikują etap wyrównania, ponieważ będą one prawdopodobnie mapowane do wielu miejsc w genomie). Warto przyjrzeć się procentowi odczytów, które przetrwały przycinanie, ponieważ wysoki odsetek odrzuconych odczytów jest prawdopodobnie oznaką złej jakości danych.

Na koniec, zwykle ponownie uruchamiam FastQC na przyciętych odczytach, aby sprawdzić, czy ten krok był skuteczny i systematycznie poprawiał metryki QC.

Wyrównanie

Z wyjątkiem wyjątków (np. de novo assembly), wyrównanie (określane również jako mapowanie) jest zwykle następnym krokiem dla większości typów danych HTS i aplikacji. Wyrównanie odczytów polega na określeniu pozycji w genomie, z której pochodzi sekwencja odczytu (zwykle wyrażana jako chromosom:początek-koniec). Stąd, na tym etapie wymagamy użycia sekwencji referencyjnej do wyrównania/mapowania odczytów.

Wybór sekwencji referencyjnej będzie zdeterminowany przez wiele czynników. Jednym z nich jest gatunek, z którego pochodzi sekwencjonowane DNA. Podczas gdy liczba gatunków z dostępną sekwencją referencyjną wysokiej jakości wzrasta, może to nadal nie dotyczyć niektórych mniej zbadanych organizmów. W takich przypadkach można chcieć wyrównać odczyty do ewolucyjnie bliskiego gatunku, dla którego dostępny jest genom referencyjny. Na przykład, ponieważ nie ma sekwencji referencyjnej dla genomu kojota, do wyrównania odczytów możemy użyć genomu blisko spokrewnionego psa. Podobnie, możemy chcieć wyrównać nasze odczyty do blisko spokrewnionego gatunku, dla którego istnieje sekwencja referencyjna o wyższej jakości. Na przykład, chociaż genom gibona został opublikowany, jest on rozbity na tysiące fragmentów, które nie odtwarzają w pełni organizacji tego genomu na dziesiątki chromosomów; w tym przypadku, przeprowadzenie wyrównania przy użyciu ludzkiej sekwencji referencyjnej może być korzystne.

Innym czynnikiem do rozważenia jest wersja sekwencji referencyjnej, ponieważ nowe wersje są wydawane w miarę uaktualniania i ulepszania sekwencji. Co ważne, współrzędne danego wyrównania mogą się różnić w zależności od wersji. Na przykład, wiele wersji ludzkiego genomu można znaleźć w UCSC Genome Browser. W przypadku każdego gatunku, zdecydowanie opowiadam się za migracją do najnowszej wersji złożenia, gdy tylko zostanie ona w pełni wydana. Może to powodować pewne niedogodności podczas przejścia, ponieważ już istniejące wyniki będą odnosiły się do starszych wersji, ale opłaca się to na dłuższą metę.

Poza tym, typ danych sekwencjonowania również ma znaczenie. Odczyty wygenerowane z protokołów DNA-seq, ChIP-seq lub Hi-C będą wyrównane do sekwencji referencyjnej genomu. Z drugiej strony, ponieważ RNA transkrybowane z DNA jest dalej przetwarzane na mRNA (tj. usuwane są introny), wiele odczytów z RNA-seq nie będzie wyrównanych do sekwencji referencyjnej genomu. Zamiast tego, musimy albo wyrównać je do sekwencji referencyjnych transkryptomu, albo użyć alignerów świadomych podziału (patrz poniżej), gdy jako odniesienie używamy sekwencji genomu. Wiąże się to z wyborem źródła anotacji sekwencji referencyjnej, czyli bazy danych zawierającej współrzędne genów, transkryptów, centromerów itp. Zazwyczaj korzystam z anotacji GENCODE, ponieważ łączy ona kompleksową anotację genów i sekwencji transkryptów.

Opracowano długą listę narzędzi do wyrównywania sekwencji short-read (patrz sekcja Short-read sequence alignment tutaj). Ich przegląd wykracza poza ramy tego artykułu (szczegóły na temat algorytmów stojących za tymi narzędziami można znaleźć tutaj). Z mojego doświadczenia wynika, że do najbardziej popularnych należą Bowtie2, BWA, HISAT2, Minimap2, STAR i TopHat. Zalecam, aby przy wyborze alignera kierować się kluczowymi czynnikami, takimi jak rodzaj danych HTS i ich zastosowanie, a także akceptacja przez społeczność, jakość dokumentacji i liczba użytkowników. Np. potrzebne są alignery takie jak STAR lub Bowtie2, które są świadome połączeń ekson-ekson podczas mapowania RNA-seq do genomu.

Wspólną cechą większości maperów jest potrzeba indeksowania sekwencji używanej jako referencja przed właściwym wyrównaniem. Ten krok może być czasochłonny, ale musi być wykonany tylko raz dla każdej sekwencji referencyjnej. Większość mapperów przechowuje wyrównania w plikach SAM/BAM, które są zgodne z formatem SAM/BAM (pliki BAM są binarną wersją plików SAM). Wyrównanie jest jednym z najbardziej obliczeniowych i czasochłonnych kroków w analizie danych sekwencjonowania, a pliki SAM/BAM są ciężkie (rzędu gigabajtów). Dlatego ważne jest, aby upewnić się, że masz wymagane zasoby (patrz ostatnia sekcja poniżej), aby uruchomić wyrównanie w rozsądnym czasie i przechowywać wyniki. Podobnie, ze względu na rozmiar i binarny format plików BAM, unikaj otwierania ich za pomocą edytorów tekstowych; zamiast tego użyj komend uniksowych lub dedykowanych narzędzi, takich jak SAMtools.

Od wyrównania

Powiedziałbym, że nie ma wyraźnego wspólnego kroku po wyrównaniu, ponieważ w tym punkcie każdy typ danych HTS i aplikacja mogą się różnić.

Powszechną analizą downstream dla danych DNA-seq jest wywoływanie wariantów, czyli identyfikacja pozycji w genomie, które różnią się w stosunku do odniesienia do genomu i pomiędzy osobami. Popularnym narzędziem analitycznym dla tego zastosowania jest GATK dla polimorfizmu pojedynczych nukleotydów (SNP) lub małych insercji/delecji (indeli) (Rysunek 2). Warianty obejmujące większe fragmenty DNA (określane również jako warianty strukturalne) wymagają dedykowanych metod wywoływania (kompleksowe porównanie znajduje się w tym artykule). Podobnie jak w przypadku alignerów, zalecam wybór odpowiedniego narzędzia, biorąc pod uwagę kluczowe czynniki, takie jak rodzaj wariantów (SNP, indel lub warianty strukturalne), akceptację przez społeczność, jakość dokumentacji i liczbę użytkowników.

Prawdopodobnie najczęstszym zastosowaniem RNA-seq jest kwantyfikacja ekspresji genów. Historycznie, odczyty musiały być wyrównane do sekwencji referencyjnej, a następnie liczba odczytów wyrównanych do danego genu lub transkryptu była używana jako przybliżenie do ilościowego określenia poziomu jego ekspresji. Takie podejście alignment+quantification jest wykonywane przez takie narzędzia jak Cufflinks, RSEM czy featureCounts. Podejście to jest jednak coraz częściej wypierane przez nowsze metody zaimplementowane w takich programach jak Kallisto czy Salmon. Koncepcyjnie, w tych narzędziach pełna sekwencja odczytu nie musi być wyrównywana do sekwencji referencyjnej. Zamiast tego wystarczy dopasować do siebie tyle nukleotydów, by mieć pewność, że dany odczyt pochodzi z danego transkryptu. Mówiąc prościej, podejście wyrównanie+kwantyfikacja jest zredukowane do jednego kroku. Podejście to znane jest jako pseudo-mapping i znacznie zwiększa szybkość kwantyfikacji ekspresji genów. Z drugiej strony, należy pamiętać, że pseudo-mapowanie nie będzie odpowiednie dla aplikacji, gdzie potrzebne jest pełne wyrównanie (np. wywoływanie wariantów z danych RNA-seq).

Innym przykładem różnic w dalszych etapach analizy i wymaganych narzędzi w aplikacjach opartych na sekwencjonowaniu jest ChIP-seq. Odczyty wygenerowane taką techniką posłużą do wywoływania pików, które polega na wykryciu w genomie regionów o znacznej nadwyżce odczytów, wskazujących na miejsce związania docelowego białka. Istnieje kilka peak callerów, które zostały omówione w niniejszej publikacji. Jako ostatni przykład podam dane Hi-C, w których alignacje są wykorzystywane jako dane wejściowe dla narzędzi wyznaczających macierze interakcji, a na ich podstawie cechy 3D genomu. Komentarz do wszystkich testów opartych na sekwencjonowaniu wykracza poza zakres tego artykułu (względnie kompletna lista znajduje się w tym artykule).

Zanim zaczniesz…

Pozostała część tego artykułu porusza aspekty, które mogą nie być ściśle rozważane jako etapy analizy danych HTS i które są w dużej mierze ignorowane. W przeciwieństwie do tego twierdzę, że kapitalne jest przemyślenie pytań postawionych w Tabeli 1 przed rozpoczęciem analizy danych HTS (lub jakiegokolwiek rodzaju danych w istocie), i pisałem na te tematy tutaj i tutaj.

Tabela 1

Pomyśl o tym	Proponowane działanie
Czy masz wszystkie informacje o swojej próbie potrzebne do analizy?	Zbieraj systematycznie metadane eksperymentów
Czy będziesz w stanie jednoznacznie zidentyfikować swoją próbkę?	Wprowadź system przypisywania każdej próbce unikalnego identyfikatora
Gdzie będą znajdować się dane i wyniki?	Strukturalna i hierarchiczna organizacja danych
Czy będziesz w stanie płynnie przetwarzać wiele próbek?	Skalowalność, równoległość, automatyczna konfiguracja i modularność kodu
Czy Ty lub ktokolwiek inny będzie w stanie odtworzyć wyniki?	Dokumentuj swój kod i procedury!

Jak wspomniano powyżej, surowe dane HTS i niektóre pliki generowane podczas ich analizy są rzędu gigabajtów, więc nie jest niczym wyjątkowym, że projekt obejmujący dziesiątki próbek wymaga terabajtów pamięci. Ponadto, niektóre etapy analizy danych HTS są intensywne obliczeniowo (np. wyrównywanie). Jednakże, pamięć masowa i infrastruktura obliczeniowa wymagana do analizy danych HTS jest ważnym czynnikiem, który jest często pomijany lub nie jest omawiany. Jako przykład, w ramach ostatniej analizy, przejrzeliśmy dziesiątki opublikowanych prac, w których przeprowadzono analizę asocjacyjną całego fenomu (PheWAS). Nowoczesne PheWAS analizują 100-1,000 zarówno wariantów genetycznych jak i fenotypów, co wymaga dużej ilości danych i mocy obliczeniowej. A jednak praktycznie żadna z przejrzanych przez nas prac nie zawierała komentarza na temat infrastruktury potrzebnej do analizy PheWAS. Nie jest więc zaskoczeniem, że zalecam, abyście z góry zaplanowali wymagania dotyczące przechowywania danych i mocy obliczeniowej, z którymi będziecie musieli się zmierzyć i podzielili się nimi ze społecznością.

Potrzebujesz pomocy w analizie danych z sekwencjonowania DNA? Skontaktuj się z niezależnymi bioinformatykami i ekspertami w dziedzinie genomiki na Kolabtree.