Een stapsgewijze gids voor de analyse van DNA-sequencinggegevens

Dr. Javier Quilez Oliete, een ervaren bioinformaticaconsultant op Kolabtree, biedt een uitgebreide gids voor de analyse van DNA-sequencinggegevens, inclusief tools en software die worden gebruikt om gegevens te lezen.

Inleiding

Deoxyribonucleïnezuur (DNA) is de molecule die het grootste deel van de genetische informatie van een organisme draagt. (In sommige virustypen wordt de genetische informatie overgebracht door ribonucleïnezuur (RNA)). Nucleotiden (gewoonlijk aangeduid met de letters A, C, G of T) zijn de basiseenheden van DNA-moleculen. Conceptueel gezien is DNA-sequencing het proces van het lezen van de nucleotiden waaruit een DNA-molecule bestaat (bv. “GCAAACCAAT” is een DNA-streng van 10 nucleotiden). De huidige sequencing-technologieën produceren miljoenen van dergelijke DNA-lezingen in een redelijke tijd en tegen een relatief lage kostprijs. Ter vergelijking: de kosten voor het sequencen van een menselijk genoom – een genoom is de volledige set DNA-moleculen in een organisme – ligt nu onder de 100 dollar en het kan in een paar dagen gebeuren. Dit staat in contrast met het eerste initiatief om het menselijk genoom te sequencen, dat in tien jaar tijd werd voltooid en ongeveer 2,7 miljard dollar kostte.

Dit vermogen om DNA te sequencen tegen een hoge verwerkingscapaciteit en lage kosten heeft de ontwikkeling van een groeiend aantal op sequencing gebaseerde methoden en toepassingen mogelijk gemaakt. Zo kan het sequencen van volledige genomen of hun eiwitcoderende regio’s (twee benaderingen die respectievelijk bekend staan als “whole genome” en “exome sequencing”) bij ziekteveroorzakende en gezonde personen DNA-veranderingen aan het licht brengen die ziekten kunnen veroorzaken. Ook de sequentiebepaling van het RNA dat uit het DNA wordt getranscribeerd – een techniek die bekend staat als RNA-sequencing – wordt gebruikt om de genactiviteit te kwantificeren en na te gaan hoe deze verandert in verschillende omstandigheden (bv. onbehandeld versus behandeling). Anderzijds detecteren chromosoomconformatiesequencingmethoden interacties tussen DNA-moleculen in de buurt en helpen zo de ruimtelijke verdeling van chromosomen binnen de cel te bepalen.

Gemeenschappelijk aan deze en andere toepassingen van DNA-sequencing is het genereren van datasets in de orde van gigabytes en bestaande uit miljoenen gelezen sequenties. Daarom vereist het zinvol gebruik van high-throughput sequencing (HTS) experimenten aanzienlijke data-analyse mogelijkheden. Gelukkig bestaan er specifieke computationele en statistische tools en relatief standaard analyse workflows voor de meeste HTS datatypes. Sommige van de (initiële) analysestappen zijn gemeenschappelijk voor de meeste sequencing datatypes, maar meer downstream analyse zal afhangen van het soort data en/of het uiteindelijke doel van de analyse. Hieronder geef ik een primer op de fundamentele stappen in de analyse van HTS gegevens en verwijs ik naar populaire tools.

Enkele van de onderstaande secties zijn gericht op de analyse van gegevens die zijn gegenereerd met short-read sequencing-technologieën (meestal Illumina), omdat deze historisch gezien de HTS-markt hebben gedomineerd. Nieuwere technologieën die langere lezingen genereren (b.v. Oxford Nanopore Technologies, PacBio) winnen echter snel terrein. Aangezien long-read sequencing een aantal bijzonderheden heeft (b.v. hogere foutenpercentages), worden specifieke instrumenten ontwikkeld voor de analyse van dit soort gegevens.

Kwaliteitscontrole (QC) van raw reads

De gretige analist start de analyse vanuit FASTQ-bestanden; het FASTQ-formaat is lange tijd de standaard geweest om short-read sequencing-gegevens op te slaan. In essentie bevatten FASTQ bestanden de nucleotide sequentie en de per-base calling quality voor miljoenen reads. Hoewel de bestandsgrootte zal afhangen van het werkelijke aantal gelezen, FASTQ bestanden zijn meestal groot (in de orde van megabytes en gigabytes) en gecomprimeerd. De meeste programma’s die FASTQ bestanden als invoer gebruiken kunnen ze in gecomprimeerd formaat verwerken, dus, om schijfruimte te besparen, is het aanbevolen om ze niet te uncomprimeren. Als conventie, zal ik hier een FASTQ bestand gelijkstellen aan een sequencing sample.

FastQC is waarschijnlijk het meest populaire gereedschap om de QC van de ruwe leest uit te voeren. Het kan worden uitgevoerd via een visuele interface of programmatisch. Terwijl de eerste optie handiger kan zijn voor gebruikers die zich niet comfortabel voelen met de command-line omgeving, biedt de laatste onvergelijkbare schaalbaarheid en reproduceerbaarheid (denk aan hoe vervelend en foutgevoelig het kan zijn om het gereedschap handmatig uit te voeren voor tientallen bestanden). Hoe dan ook, de belangrijkste output van FastQC is een HTML-bestand met samenvattende statistieken over de algehele kwaliteit van de ruwe sequencing reads van een bepaald sample. Het één voor één inspecteren van tientallen FastQC rapporten is vervelend en bemoeilijkt de vergelijking tussen monsters. Daarom kunt u wellicht MultiQC gebruiken, dat de HTML-rapporten van FastQC (en van andere downstream gebruikte tools, bv. adapter trimmen, alignment) in één rapport samenvoegt.

QC-informatie is bedoeld om de gebruiker in staat te stellen te beoordelen of monsters van goede kwaliteit zijn en dus voor de volgende stappen kunnen worden gebruikt, of dat ze moeten worden weggegooid. Helaas is er geen consensusdrempel op basis van de FastQC-metriek om monsters als van goede of slechte kwaliteit te classificeren. De aanpak die ik gebruik is de volgende. Ik verwacht dat alle monsters die dezelfde procedure hebben doorlopen (bv. DNA-extractie, bibliotheekvoorbereiding) vergelijkbare kwaliteitsstatistieken en een meerderheid van “pass”-vlaggen hebben. Als sommige monsters van mindere dan gemiddelde kwaliteit zijn, gebruik ik ze toch in de downstream-analyse met dit in het achterhoofd. Aan de andere kant, als alle monsters in het experiment systematisch krijgen “waarschuwing” of “fail” vlaggen in meerdere metrieken (zie dit voorbeeld), vermoed ik dat er iets mis ging in het experiment (bijv. slechte DNA-kwaliteit, bibliotheek voorbereiding, enz.) en ik adviseer het te herhalen.

Read trimmen

QC van ruwe leest helpt om problematische monsters te identificeren, maar het verbetert niet de werkelijke kwaliteit van de leest. Om dit te doen, moeten we lezen trimmen om technische sequenties en lage-kwaliteit uiteinden te verwijderen.

Technische sequenties zijn overblijfselen van de experimentele procedure (bijv. sequencing adapters). Als dergelijke sequenties grenzen aan de ware sequentie van de afgelezen sequentie, kan de uitlijning (zie hieronder) de afgelezen sequenties op de verkeerde positie in het genoom plaatsen of het vertrouwen in een bepaalde uitlijning verminderen. Naast technische sequenties kunnen we ook sequenties van biologische oorsprong willen verwijderen als die in hoge mate aanwezig zijn onder de gelezen sequenties. Zo kunnen suboptimale DNA-preparatieprocedures een hoog aandeel van DNA-geconverteerd ribosomaal RNA (rRNA) in het monster achterlaten. Tenzij dit type nucleïnezuur is het doel van de sequencing experiment, het houden van leest afgeleid van rRNA zal gewoon verhoging van de computationele last van de downstream-stappen en kan verwarren de resultaten. Als het niveau van technische sequenties, rRNA of andere verontreinigingen erg hoog is, wat waarschijnlijk al door de QC is aangetoond, kunt u het hele sequencingmonster weggooien.

In short-read sequencing wordt de DNA-sequentie nucleotide voor nucleotide bepaald (technisch gesproken, nucleotide voor nucleotide per sequencingcyclus). Met andere woorden, het aantal sequencing-cycli bepaalt de leeslengte. Een bekend probleem bij HTS-sequencingmethoden is dat de nauwkeurigheid waarmee nucleotiden worden bepaald, afneemt naarmate het aantal sequencingcycli toeneemt. Dit komt tot uiting in een algemene daling van de aanroepkwaliteit per base, vooral tegen het einde van de lezing. Net als bij technische sequenties kan het proberen uit te lijnen van lezingen die uiteinden van lage kwaliteit bevatten, leiden tot misplaatsing of slechte mapping-kwaliteit.

Om technische/verontreinigende sequenties en uiteinden van lage kwaliteit te verwijderen, bestaan er lees-trimhulpmiddelen zoals Trimmomatic en Cutadapt, die op grote schaal worden gebruikt. In essentie verwijderen deze tools technische sequenties (intern beschikbaar en/of door de gebruiker aangeleverd) en trimmen ze de gelezen tekst op basis van kwaliteit, terwijl de leeslengte gemaximaliseerd wordt. Lezingen die na het trimmen te kort zijn gebleven, worden verwijderd (te korte lezingen, bv. <36 nucleotiden, bemoeilijken de alignmentstap, omdat ze waarschijnlijk op meerdere plaatsen in het genoom zullen voorkomen). U kunt kijken naar het percentage van de gelezen die overleven het trimmen, als een hoog percentage van weggegooide leest is waarschijnlijk een teken van slechte kwaliteit gegevens.

Ten slotte, ik typisch opnieuw FastQC uitvoeren op de getrimde leest om te controleren of deze stap was effectief en systematisch verbeterde de QC metrics.

Uitlijning

Op uitzonderingen na (bijv. de novo montage), uitlijning (ook wel aangeduid als mapping) is typisch de volgende stap voor de meeste HTS data types en toepassingen. Lees uitlijning bestaat in het bepalen van de positie in het genoom waarvan de volgorde van de gelezen afgeleid (meestal uitgedrukt als chromosoom:start-end). Vandaar dat in deze stap moeten we het gebruik van een referentie-sequentie uit te lijnen/map de gelezen op.

De keuze van de referentie-sequentie zal worden bepaald door meerdere factoren. Ten eerste de soort waarvan het gesequenteerde DNA afkomstig is. Hoewel het aantal soorten waarvan een referentiesequentie van hoge kwaliteit beschikbaar is, toeneemt, kan dit voor sommige minder bestudeerde organismen nog steeds niet het geval zijn. In die gevallen kan het wenselijk zijn de gelezen sequenties te aligneren met een soort die qua evolutie dichtbij staat en waarvoor een referentiegenoom beschikbaar is. Bijvoorbeeld, omdat er geen referentie sequentie is voor het genoom van de coyote, kunnen we die van de nauw verwante hond gebruiken voor de read alignment. Ook kan het zijn dat we onze gelezen willen uitlijnen op een nauw verwante soort waarvoor een referentiesequentie van hogere kwaliteit bestaat. Bijvoorbeeld, terwijl het genoom van de gibbon is gepubliceerd, dit is opgedeeld in duizenden fragmenten die niet volledig recapituleren de organisatie van dat genoom in tientallen chromosomen; in dat geval kan het uitvoeren van de uitlijning met behulp van de menselijke referentie-sequentie gunstig zijn.

Een andere factor om te overwegen is de versie van de referentie-sequentie assemblage, omdat nieuwe versies worden vrijgegeven als de sequentie wordt bijgewerkt en verbeterd. Belangrijk is dat de coördinaten van een bepaalde uitlijning tussen versies kunnen verschillen. Zo zijn er bijvoorbeeld meerdere versies van het menselijk genoom te vinden in de UCSC Genome Browser. In elk geval ben ik er een groot voorstander van om naar de nieuwste assemblageversie te migreren zodra die volledig is vrijgegeven. Dit kan wat overlast veroorzaken tijdens de overgang, omdat reeds bestaande resultaten zich zullen verhouden tot oudere versies, maar op de lange termijn loont het de moeite.

Bijkomend is ook het type sequencing-data van belang. Reads die zijn gegenereerd met DNA-seq-, ChIP-seq- of Hi-C-protocollen zullen worden afgestemd op de genoomreferentiesequentie. Aan de andere kant, als RNA getranscribeerd uit DNA verder wordt verwerkt tot mRNA (d.w.z. introns verwijderd), zullen veel RNA-seq gegevens niet uitlijnen met een genoom referentie sequentie. In plaats daarvan moeten we ofwel hen uit te lijnen met transcriptoom referentiesequenties of split-aware aligners (zie hieronder) te gebruiken bij gebruik van het genoom sequentie als referentie. Hiermee samenhangend is de keuze van de bron voor de annotatie van de referentiesequentie, d.w.z. de database met de coördinaten van de genen, transcripten, centromeren, enz. Ik gebruik meestal de GENCODE annotatie, omdat die een uitgebreide genannotatie en transcriptsequenties combineert.

Er is een lange lijst van short-read sequentie alignment tools ontwikkeld (zie de Short-read sequentie alignment sectie hier). Het bespreken van deze tools valt buiten het bestek van dit artikel (details over de algoritmen achter deze tools zijn hier te vinden). In mijn ervaring zijn Bowtie2, BWA, HISAT2, Minimap2, STAR en TopHat de meest populaire. Mijn aanbeveling is dat u uw aligner kiest op basis van belangrijke factoren zoals het type HTS data en toepassing, maar ook acceptatie door de gemeenschap, kwaliteit van de documentatie en het aantal gebruikers. Bijvoorbeeld, men heeft aligners zoals STAR of Bowtie2 nodig die zich bewust zijn van exon-exon junctions bij het in kaart brengen van RNA-seq naar het genoom.

Gemeenschappelijk voor de meeste mappers is de noodzaak om de sequentie die als referentie wordt gebruikt te indexeren voordat de eigenlijke uitlijning plaatsvindt. Deze stap kan tijdrovend zijn, maar het hoeft slechts eenmaal te worden gedaan voor elke referentiesequentie. De meeste mappers slaan de alignments op in SAM/BAM-bestanden, die het SAM/BAM-formaat volgen (BAM-bestanden zijn binaire versies van SAM-bestanden). De alignment is een van de meest rekenintensieve en tijdrovende stappen in de analyse van sequencing data en SAM/BAM-bestanden zijn zwaar (in de orde van gigabytes). Daarom is het belangrijk ervoor te zorgen dat u over de nodige middelen beschikt (zie de laatste paragraaf hieronder) om de uitlijning in een redelijke tijd uit te voeren en de resultaten op te slaan. Evenzo, vanwege de grootte en het binaire formaat van BAM-bestanden, vermijd ze te openen met tekstverwerkers; in plaats daarvan gebruik Unix commando’s of speciale tools zoals SAMtools.

Van de uitlijningen

Ik zou zeggen dat er niet een duidelijke gemeenschappelijke stap na de uitlijning, omdat op dit punt is waar elke HTS data type en toepassing kan verschillen.

Een veel voorkomende downstream-analyse voor DNA-seq-gegevens is variant calling, dat wil zeggen de identificatie van posities in het genoom die verschillen ten opzichte van het genoom referentie en tussen individuen. Een populair analysekader voor deze toepassing is GATK voor single nucleotide polymorfisme (SNP) of kleine invoegingen/deleties (indels) (figuur 2). Varianten die uit grotere stukken DNA bestaan (ook wel structurele varianten genoemd) vereisen specifieke calling methoden (zie dit artikel voor een uitgebreide vergelijking). Net als bij de aligners, adviseer ik bij het kiezen van het juiste gereedschap rekening te houden met belangrijke factoren zoals het soort varianten (SNP, indel of structurele varianten), acceptatie door de gemeenschap, kwaliteit van de documentatie en het aantal gebruikers.

Waarschijnlijk de meest frequente toepassing van RNA-seq is genexpressie kwantificering. Historisch gezien, gelezen moest worden uitgelijnd met de referentie-sequentie en vervolgens het aantal gelezen uitgelijnd met een bepaald gen of transcript werd gebruikt als een proxy om de expressieniveaus te kwantificeren. Deze aanpak van uitlijning en kwantificering wordt uitgevoerd door instrumenten zoals Cufflinks, RSEM of featureCounts. De scuh aanpak wordt echter steeds meer overtroffen door nieuwere methoden die zijn geïmplementeerd in software zoals Kallisto en Salmon. Met dergelijke hulpmiddelen hoeft de volledige sequentie van een lees niet te worden uitgelijnd met de referentiesequentie. In plaats daarvan hoeven alleen voldoende nucleotiden te worden uitgelijnd om er zeker van te zijn dat een afgelezen sequentie afkomstig is van een bepaald transcript. Eenvoudig gezegd: de aanpak uitlijning+kwantificering wordt tot één stap gereduceerd. Deze aanpak staat bekend als pseudo-mapping en verhoogt de snelheid van de genexpressie-kwantificering aanzienlijk. Aan de andere kant moet men bedenken dat pseudo-mapping niet geschikt is voor toepassingen waarbij de volledige uitlijning nodig is (b.v. variant calling van RNA-seq gegevens).

Een ander voorbeeld van de verschillen in de downstream analysestappen en de vereiste hulpmiddelen bij sequencing-gebaseerde toepassingen is ChIP-seq. Lezingen die met een dergelijke techniek worden gegenereerd, zullen worden gebruikt voor “peak calling”, waarbij regio’s in het genoom worden opgespoord met een significante overmaat aan lezingen die aangeeft waar het doeleiwit is gebonden. Er bestaan verscheidene “peak callers” en deze publicatie geeft een overzicht daarvan. Als laatste voorbeeld noem ik de Hi-C data, waarbij alignments worden gebruikt als input voor tools die de interactie matrices bepalen en daaruit de 3D-features van het genoom. Commentaar op alle op sequencing gebaseerde assays valt buiten het bestek van dit artikel (voor een betrekkelijk volledige lijst zie dit artikel).

Voordat u begint…

Het resterende deel van dit artikel raakt aan aspecten die wellicht niet strikt als stappen in de analyse van HTS-gegevens worden beschouwd en die grotendeels worden genegeerd. Ik stel daarentegen dat het van kapitaal belang is dat u nadenkt over de vragen die in tabel 1 worden gesteld voordat u begint met de analyse van HTS-gegevens (of van welk soort gegevens dan ook), en ik heb hier en hier over deze onderwerpen geschreven.

Tabel 1

Denk er eens over na	Voorgestelde actie
Krijgt u alle informatie van uw steekproef die nodig is voor de analyse?	Verzamel systematisch de metadata van de experimenten
Wilt u uw monster eenduidig kunnen identificeren?	Zet een systeem op om aan elk monster een unieke identificatiecode toe te kennen
Waar worden de gegevens en resultaten bewaard?	Gestructureerde en hiërarchische organisatie van de gegevens
Zult u in staat zijn om meerdere monsters naadloos te verwerken?	Kalibaarheid, parallellisatie, automatische configuratie en modulariteit van de code
Zult u of iemand anders in staat zijn om de resultaten te reproduceren?	Documenteer uw code en procedures!

Zoals hierboven vermeld, zijn de ruwe HTS-gegevens en sommige van de bestanden die bij de analyse daarvan worden gegenereerd in de orde van gigabytes, zodat het niet uitzonderlijk is dat een project met tientallen monsters terabytes aan opslagruimte vergt. Bovendien zijn sommige stappen in de analyse van HTS-gegevens rekenintensief (bv. uitlijning). De opslag- en computerinfrastructuur die nodig is voor de analyse van HTS-gegevens is echter een belangrijke overweging die vaak over het hoofd wordt gezien of niet wordt besproken. In het kader van een recente analyse hebben wij bijvoorbeeld tientallen gepubliceerde artikelen bestudeerd waarin fenomeenbrede associatieanalyses (PheWAS) werden uitgevoerd. Moderne PheWAS analyseren 100-1.000s van zowel genetische varianten als fenotypes, wat resulteert in belangrijke data-opslag en rekenkracht. En toch heeft vrijwel geen van de papers die we hebben bekeken commentaar gegeven op de infrastructuur die nodig is voor de PheWAS analyse. Het zal u dan ook niet verbazen dat ik u adviseer om van tevoren de opslag- en rekenvereisten te plannen waarmee u te maken krijgt en deze met de gemeenschap te delen.

Hulp nodig bij het analyseren van DNA-sequencinggegevens? Kom in contact met freelance bioinformatici en genomics-experts op Kolabtree.