En trin-for-trin-guide til analyse af DNA-sekventeringsdata

Dr. Javier Quilez Oliete, en erfaren bioinformatik-konsulent på Kolabtree, giver en omfattende vejledning i analyse af DNA-sekventeringsdata, herunder værktøjer og software, der bruges til at læse data.

Indledning

Deoxyribonukleinsyre (DNA) er det molekyle, der bærer det meste af den genetiske information i en organisme. (I nogle typer virus bæres den genetiske information af ribonukleinsyre (RNA).) Nukleotider (konventionelt repræsenteret ved bogstaverne A, C, G eller T) er de grundlæggende enheder i DNA-molekyler. Begrebsmæssigt er DNA-sekventering en proces, hvor man aflæser de nukleotider, som et DNA-molekyle består af (f.eks. er “GCAAACCAAT” en DNA-streng på 10 nukleotider). De nuværende sekventeringsteknologier producerer millioner af sådanne DNA-reads på rimelig tid og til en relativt lav pris. Som reference kan nævnes, at prisen for sekventering af et menneskeligt genom – et genom er det komplette sæt af DNA-molekyler i en organisme – er faldet til under 100 USD, og det kan gøres på få dage. Dette står i kontrast til det første initiativ til sekventering af det menneskelige genom, som blev afsluttet på et årti og havde en omkostning på ca. 2,7 mia. dollars.

Denne evne til at sekventere DNA med høj gennemstrømning og lave omkostninger har gjort det muligt at udvikle et stigende antal sekventeringsbaserede metoder og applikationer. F.eks. kan sekventering af hele genomer eller deres proteinkodende regioner (to metoder, der er kendt som henholdsvis helgenom- og exomsekventering) hos syge og raske personer give et fingerpeg om sygdomsfremkaldende DNA-forandringer. Ligeledes anvendes sekventering af det RNA, der transskriberes fra DNA – en teknik, der er kendt som RNA-sekventering – til at kvantificere genaktivitet, og hvordan denne ændrer sig under forskellige forhold (f.eks. ubehandlet i forhold til behandling). På den anden side registrerer kromosomkonformitetsopsætningssekvenseringsmetoder interaktioner mellem nærliggende DNA-molekyler og bidrager således til at bestemme kromosomernes rumlige fordeling i cellen.

Fælles for disse og andre anvendelser af DNA-sekventering er genereringen af datasæt i størrelsesordenen gigabyte, der omfatter millioner af læste sekvenser. Derfor kræver det betydelige dataanalysekapaciteter at få mening ud af HTS-eksperimenter (high-throughput sequencing). Heldigvis findes der dedikerede beregningsmæssige og statistiske værktøjer og relativt standardiserede analysearbejdsgange for de fleste HTS-datatyper. Mens nogle af de (indledende) analysetrin er fælles for de fleste sekventeringsdatatyper, afhænger de mere efterfølgende analyser af datatypen og/eller det endelige mål med analysen. Nedenfor giver jeg en grundbog om de grundlæggende trin i analysen af HTS-data, og jeg henviser til populære værktøjer.

Nogle af nedenstående afsnit er fokuseret på analyse af data genereret fra short-read-sekventeringsteknologier (hovedsagelig Illumina), da disse historisk set har domineret HTS-markedet. Nyere teknologier, der genererer længere læsninger (f.eks. Oxford Nanopore Technologies, PacBio), vinder imidlertid hurtigt terræn. Da sekventering med lange læsninger har nogle særlige karakteristika (f.eks. højere fejlfrekvens), er der ved at blive udviklet specifikke værktøjer til analyse af denne type data.

Kvalitetskontrol (QC) af rå læsninger

Den ivrige analytiker vil starte analysen fra FASTQ-filer; FASTQ-formatet har længe været standard for lagring af short-read-sekventeringsdata. I det væsentlige indeholder FASTQ-filer nukleotidsekvensen og kvaliteten pr. base calling quality for millioner af læsninger. Selv om filstørrelsen afhænger af det faktiske antal læsninger, er FASTQ-filer typisk store (i størrelsesordenen megabytes og gigabytes) og komprimerede. Det skal bemærkes, at de fleste værktøjer, der bruger FASTQ-filer som input, kan håndtere dem i komprimeret format, så for at spare diskplads anbefales det ikke at udpakke dem for at spare diskplads. Som en konvention vil jeg her sidestille en FASTQ-fil med en sekventeringsprøve.

FastQC er sandsynligvis det mest populære værktøj til at udføre QC af de rå læsninger. Det kan køres via en visuel grænseflade eller programmatisk. Mens den første mulighed kan være mere bekvem for brugere, der ikke føler sig trygge ved kommandolinemiljøet, tilbyder sidstnævnte en uforlignelig skalerbarhed og reproducerbarhed (tænk på, hvor kedeligt og fejlbehæftet det kan være at køre værktøjet manuelt for titusindvis af filer). Uanset hvad er det vigtigste output af FastQC en HTML-fil, der rapporterer vigtige sammenfattende statistikker om den overordnede kvalitet af de rå sekventeringslæsninger fra en given prøve. Det er besværligt at inspicere snesevis af FastQC-rapporter én for én, og det komplicerer sammenligningen på tværs af prøverne. Derfor kan du måske bruge MultiQC, som samler HTML-rapporterne fra FastQC (samt fra andre værktøjer, der anvendes nedstrøms, f.eks. adaptertrimming, alignment) i en enkelt rapport.

QC-oplysninger er beregnet til at give brugeren mulighed for at vurdere, om prøverne har god kvalitet og derfor kan bruges til de efterfølgende trin, eller om de skal kasseres. Desværre findes der ikke en konsensusgrænse baseret på FastQC-metrikken til at klassificere prøver som værende af god eller dårlig kvalitet. Den fremgangsmåde, som jeg anvender, er følgende. Jeg forventer, at alle prøver, der har gennemgået den samme procedure (f.eks. DNA-ekstraktion, biblioteksforberedelse), har lignende kvalitetsstatistikker og et flertal af “bestået”-flag. Hvis nogle prøver har en lavere kvalitet end gennemsnittet, vil jeg stadig bruge dem i den efterfølgende analyse under hensyntagen til dette. På den anden side, hvis alle prøver i eksperimentet systematisk får “advarsel” eller “fail”-flag i flere metrikker (se dette eksempel), har jeg mistanke om, at noget gik galt i eksperimentet (f.eks. dårlig DNA-kvalitet, biblioteksforberedelse osv.), og jeg anbefaler, at det gentages.

Read trimming

QC af rå læsninger hjælper med at identificere problematiske prøver, men det forbedrer ikke læsningernes faktiske kvalitet. For at gøre det er vi nødt til at trimme læserne for at fjerne tekniske sekvenser og ender af lav kvalitet.

Tekniske sekvenser er rester fra den eksperimentelle procedure (f.eks. sekventeringsadaptere). Hvis sådanne sekvenser støder op til læsningens sande sekvens, kan en alignment (se nedenfor) kortlægge læsninger til den forkerte position i genomet eller mindske tilliden til en given alignment. Ud over tekniske sekvenser kan vi også ønske at fjerne sekvenser af biologisk oprindelse, hvis disse er meget udbredte blandt læsningerne. F.eks. kan suboptimale DNA-præparationsprocedurer efterlade en høj andel af DNA-konverteret ribosomalt RNA (rRNA) i prøven. Medmindre denne type nukleinsyre er målet for sekventeringseksperimentet, vil det at beholde læsninger, der stammer fra rRNA, blot øge den beregningsmæssige byrde for de efterfølgende trin og kan forvirre resultaterne. Det skal bemærkes, at hvis niveauet af tekniske sekvenser, rRNA eller andre forurenende stoffer er meget højt, hvilket sandsynligvis allerede vil være blevet fremhævet af QC, kan det være hensigtsmæssigt at kassere hele sekventeringsprøven.

I short-read-sekventering bestemmes DNA-sekvensen et nukleotid ad gangen (teknisk set et nukleotid for hver sekventeringscyklus). Med andre ord er det antallet af sekventeringscyklusser, der bestemmer læselængden. Et kendt problem ved HTS-sekventeringsmetoder er, at den nøjagtighed, hvormed nukleotiderne bestemmes, aftager, efterhånden som sekventeringscyklusserne akkumuleres. Dette afspejler sig i et generelt fald i kvaliteten af den kaldende kvalitet pr. base, især mod slutningen af læsningen. Som det sker med tekniske sekvenser, kan forsøg på at tilpasse læsninger, der indeholder ender af lav kvalitet, føre til fejlplacering eller dårlig kortlægningskvalitet.

For at fjerne tekniske/forurenende sekvenser og ender af lav kvalitet findes der værktøjer til trimning af læsninger som Trimmomatic og Cutadapt, som anvendes i vid udstrækning. I det væsentlige fjerner sådanne værktøjer tekniske sekvenser (internt tilgængelige og/eller leveret af brugeren) og trimmer læsninger baseret på kvalitet, samtidig med at læsningslængden maksimeres. Læsninger, der er for korte efter trimningen, kasseres (læsninger, der er for korte, f.eks. <36 nukleotider, komplicerer tilpasningsfasen, da de sandsynligvis vil være knyttet til flere steder i genomet). Du bør måske se på procentdelen af reads, der overlever trimningen, da en høj andel af kasserede reads sandsynligvis er et tegn på data af dårlig kvalitet.

Slutteligt kører jeg typisk FastQC igen på de trimmede reads for at kontrollere, at dette trin var effektivt og systematisk forbedrede QC-metrikken.

Alignment

Med undtagelser (f.eks. de novo assemblage) er alignment (også kaldet mapping) typisk det næste trin for de fleste HTS-datatyper og -applikationer. Læsejustering består i at bestemme den position i genomet, som sekvensen af læsningen stammer fra (typisk udtrykt som kromosom:start-end). På dette trin er det derfor nødvendigt at anvende en referencesekvens til at tilpasse/mappe read’erne på.

Valget af referencesekvensen vil blive bestemt af flere faktorer. For det første den art, som det sekventerede DNA stammer fra. Selv om antallet af arter med en tilgængelig referencesekvens af høj kvalitet er stigende, er dette måske stadig ikke tilfældet for nogle mindre undersøgte organismer. I disse tilfælde kan det være ønskeligt at tilpasse læsninger til en evolutivt nærliggende art, for hvilken der findes et referencegenom. Da der f.eks. ikke findes en referencesekvens for prærieulvens genom, kan vi bruge den nært beslægtede hunds genom til read-justeringen. På samme måde kan vi stadig ønske at tilpasse vores læsninger til en nært beslægtet art, for hvilken der findes en referencesekvens af højere kvalitet. For eksempel er gibbonens genom ganske vist blevet offentliggjort, men det er opdelt i tusindvis af fragmenter, som ikke fuldt ud gengiver organiseringen af dette genom i titusindvis af kromosomer; i det tilfælde kan det være en fordel at udføre tilpasningen ved hjælp af den menneskelige referencesekvens.

En anden faktor, der skal tages i betragtning, er versionen af referencesekvenssamlingen, da der udgives nye versioner, efterhånden som sekvensen opdateres og forbedres. Det er vigtigt, at koordinaterne for en given alignment kan variere mellem versioner. F.eks. kan der findes flere versioner af det menneskelige genom i UCSC Genome Browser. Uanset art går jeg stærkt ind for at migrere til den nyeste assemblageversion, når den er fuldt ud frigivet. Dette kan medføre nogle gener under overgangen, da allerede eksisterende resultater vil være relative til ældre versioner, men det betaler sig i det lange løb.

Dertil kommer, at typen af sekventeringsdata også har betydning. Læsninger, der genereres fra DNA-seq-, ChIP-seq- eller Hi-C-protokoller, vil blive afstemt med genomets referencesekvens. På den anden side vil mange RNA-seq-læsninger, da RNA, der er transskriberet fra DNA, viderebehandles til mRNA (dvs. introner fjernes), ikke blive afstemt med en genomreferencesekvens. I stedet skal vi enten tilpasse dem til transkriptomreferencesekvenser eller bruge split-aware aligners (se nedenfor), når vi bruger genomsekvensen som reference. Relateret til dette er valget af kilde til annotation af referencesekvensen, dvs. databasen med koordinaterne for generne, transkriptionerne, centromer osv. Jeg bruger typisk GENCODE-annotationen, da den kombinerer omfattende genannotation og transkriptsekvenser.

Der er blevet udviklet en lang række værktøjer til kortlæsningssekvensalignering (se afsnittet om kortlæsningssekvensalignering her). En gennemgang af dem ligger uden for rammerne af denne artikel (detaljer om algoritmerne bag disse værktøjer kan findes her). Min erfaring viser, at blandt de mest populære er Bowtie2, BWA, HISAT2, Minimap2, STAR og TopHat. Min anbefaling er, at du vælger din aligner på baggrund af nøglefaktorer som typen af HTS-data og anvendelse samt accept i samfundet, kvaliteten af dokumentationen og antallet af brugere. F.eks. har man brug for alignere som STAR eller Bowtie2, der er opmærksomme på exon-exon junctions, når man kortlægger RNA-seq til genomet.

Fælles for de fleste mappere er behovet for at indeksere den sekvens, der anvendes som reference, før den egentlige alignment finder sted. Dette trin kan være tidskrævende, men det behøver kun at blive udført én gang for hver referencesekvens. De fleste mappere gemmer alignments i SAM/BAM-filer, som følger SAM/BAM-formatet (BAM-filer er binære versioner af SAM-filer). Justeringen er et af de mest beregnings- og tidskrævende trin i analysen af sekventeringsdata, og SAM/BAM-filer er tunge (i størrelsesordenen gigabyte). Derfor er det vigtigt at sikre sig, at man har de nødvendige ressourcer (se det sidste afsnit nedenfor) til at køre tilpasningen på en rimelig tid og gemme resultaterne. På samme måde bør du på grund af BAM-filernes størrelse og binære format undgå at åbne dem med tekstredigeringsprogrammer; brug i stedet Unix-kommandoer eller dedikerede værktøjer som SAMtools.

Fra alignments

Jeg vil sige, at der ikke er et klart fælles trin efter alignmentet, da det er på dette punkt, hvor hver HTS-datatatype og applikation kan være forskellig.

En almindelig downstream-analyse for DNA-seq-data er variant calling, dvs. identifikation af positioner i genomet, der varierer i forhold til genomreferencen og mellem individer. En populær analyseramme til denne anvendelse er GATK til enkeltnukleotidpolymorfisme (SNP) eller små indsættelser/deletioner (indels) (figur 2). Varianter, der omfatter større dele af DNA (også kaldet strukturelle varianter), kræver dedikerede metoder til at kalde varianter (se denne artikel for en omfattende sammenligning). Som med alignere anbefaler jeg, at man vælger det rigtige værktøj under hensyntagen til nøglefaktorer som f.eks. slags varianter (SNP, indel eller strukturelle varianter), accept af fællesskabet, dokumentationens kvalitet og antallet af brugere.

Den mest hyppige anvendelse af RNA-seq er formentlig kvantificering af genekspression. Historisk set var det nødvendigt at tilpasse læsninger til referencesekvensen, og derefter blev antallet af læsninger tilpasset et givet gen eller transkript anvendt som en proxy til at kvantificere dets ekspressionsniveauer. Denne tilgang til alignment+kvantificering udføres af værktøjer som Cufflinks, RSEM eller featureCounts. Scuh-metoden er imidlertid i stigende grad blevet overgået af nyere metoder, der er implementeret i software som Kallisto og Salmon. Med sådanne værktøjer behøver den fulde sekvens af en læsning ikke at blive afstemt med referencesekvensen. I stedet behøver vi kun at tilpasse tilstrækkeligt mange nukleotider til at være sikre på, at en læsning stammer fra et givet transkript. Kort sagt kan man sige, at metoden med tilpasning+kvantificering er reduceret til et enkelt trin. Denne fremgangsmåde er kendt som pseudokortlægning og øger i høj grad hastigheden af genekspressionskvantificeringen. På den anden side skal man huske på, at pseudokortlægning ikke vil være egnet til anvendelser, hvor der er behov for en fuld tilpasning (f.eks. variant calling fra RNA-seq-data).

Et andet eksempel på forskellene i de efterfølgende analysetrin og de nødvendige værktøjer på tværs af sekventeringsbaserede anvendelser er ChIP-seq. Læsninger genereret med en sådan teknik vil blive brugt til peak calling, som består i at opdage regioner i genomet med et betydeligt overskud af læsninger, der indikerer, hvor målproteinet er bundet. Der findes adskillige peak callers, og i denne publikation gennemgås de forskellige. Som et sidste eksempel vil jeg nævne Hi-C-data, hvor alignments anvendes som input til værktøjer, der bestemmer interaktionsmatricerne og heraf genomets 3D-funktioner. Kommentarer til alle de sekventeringsbaserede assays ligger uden for rammerne af denne artikel (for en relativt komplet liste se denne artikel).

Hvor du begynder…

Den resterende del af denne artikel berører aspekter, som måske ikke strengt taget betragtes som trin i analysen af HTS-data, og som i vid udstrækning ignoreres. I modsætning hertil argumenterer jeg for, at det er vigtigt, at du tænker over de spørgsmål, der er stillet i tabel 1, før du begynder at analysere HTS-data (eller faktisk enhver form for data), og jeg har skrevet om disse emner her og her.

Tabel 1

Tænk over det	Forslag
Har du alle de oplysninger om din prøve, der er nødvendige for analysen?	Saml systematisk metadata om forsøgene
Kan du entydigt identificere din prøve?	Etabler et system til at tildele hver prøve en entydig identifikator
Hvor vil data og resultater blive opbevaret?	Struktureret og hierarkisk organisering af dataene
Vil du være i stand til at behandle flere prøver problemfrit?	Skalering, parallelisering, automatisk konfiguration og modularitet af koden
Vil du eller andre kunne reproducere resultaterne?	Dokumenter din kode og dine procedurer!

Som nævnt ovenfor er HTS-råddata og nogle af de filer, der genereres under analysen heraf, i størrelsesordenen gigabyte, så det er ikke usædvanligt, at et projekt, der omfatter titusindvis af prøver, kræver terabyte af lagerplads. Desuden er nogle trin i analysen af HTS-data beregningskrævende (f.eks. alignment). Den lager- og beregningsinfrastruktur, der kræves til analyse af HTS-data, er imidlertid en vigtig faktor, som ofte overses eller ikke diskuteres. Som et eksempel kan nævnes, at vi som led i en nylig analyse gennemgik snesevis af offentliggjorte artikler, der udførte phenome-wide association analysis (PheWAS). Moderne PheWAS analyserer 100-1.000-vis af både genetiske varianter og fænotyper, hvilket resulterer i betydelig datalagring og computerkraft. Alligevel kommenterede stort set ingen af de artikler, vi gennemgik, den infrastruktur, der er nødvendig for PheWAS-analysen. Det er ikke overraskende, at min anbefaling er, at I på forhånd planlægger de krav til lagring og databehandling, som I vil stå over for, og at I deler dem med fællesskabet.

Har I brug for hjælp til at analysere DNA-sekventeringsdata? Kom i kontakt med freelance bioinformatikere og genomeksperter på Kolabtree.