En steg-för-steg-guide till analys av DNA-sekvensdata

Dr. Javier Quilez Oliete, en erfaren bioinformatikkonsult på Kolabtree, ger en omfattande guide till analys av DNA-sekvensdata, inklusive verktyg och programvara som används för att läsa data.

Introduktion

Deoxyribonukleinsyra (DNA) är den molekyl som bär den mesta av den genetiska informationen i en organism. (I vissa typer av virus bärs den genetiska informationen av ribonukleinsyra (RNA).) Nukleotider (konventionellt representerade av bokstäverna A, C, G eller T) är de grundläggande enheterna i DNA-molekyler. Begreppsmässigt är DNA-sekvensering processen att läsa av de nukleotider som ingår i en DNA-molekyl (t.ex. ”GCAAACCAAT” är en DNA-sträng med 10 nukleotider). Med nuvarande sekvenseringsteknik produceras miljontals sådana DNA-avläsningar på en rimlig tid och till en relativt låg kostnad. Som referens kan nämnas att kostnaden för att sekvensera ett mänskligt genom – ett genom är den fullständiga uppsättningen av DNA-molekyler i en organism – har sjunkit till under 100 dollar och kan göras på några dagar. Detta står i kontrast till det första initiativet för att sekvensera det mänskliga genomet, som avslutades på ett decennium och hade en kostnad på cirka 2,7 miljarder dollar.

Denna förmåga att sekvensera DNA med hög genomströmning och låg kostnad har möjliggjort utvecklingen av ett växande antal sekvenseringsbaserade metoder och tillämpningar. Genom att sekvensera hela genomer eller deras proteinkodande regioner (två metoder som kallas helgenomsekvensering respektive exomsekvensering) hos sjuka och friska individer kan man till exempel hitta tips på sjukdomsorsakande DNA-förändringar. Sekvensering av det RNA som transkriberas från DNA – en teknik som kallas RNA-sekvensering – används också för att kvantifiera genaktivitet och hur denna förändras under olika förhållanden (t.ex. obehandlad kontra behandling). Å andra sidan upptäcker metoder för sekvensering av kromosomkonformationer interaktioner mellan närliggande DNA-molekyler och hjälper därmed till att bestämma kromosomernas rumsliga fördelning i cellen.

Gemensamt för dessa och andra tillämpningar av DNA-sekvensering är genereringen av datamängder i storleksordningen gigabyte och som består av miljontals avlästa sekvenser. För att kunna förstå experiment med sekvensering med högt genomflöde (HTS) krävs därför omfattande kapacitet för dataanalys. Lyckligtvis finns det särskilda beräknings- och statistikverktyg och relativt standardiserade analysarbetsflöden för de flesta HTS-datatyper. Vissa av de (inledande) analysstegen är gemensamma för de flesta typer av sekvenseringsdata, men mer efterföljande analyser beror på typen av data och/eller det slutliga målet med analysen. Nedan ger jag en grundkurs i de grundläggande stegen i analysen av HTS-data och hänvisar till populära verktyg.

En del av avsnitten nedan är inriktade på analysen av data som genererats från sekvenseringstekniker med korta läsningar (främst Illumina), eftersom dessa historiskt sett har dominerat HTS-marknaden. Nyare tekniker som genererar längre läsningar (t.ex. Oxford Nanopore Technologies, PacBio) vinner dock snabbt mark. Eftersom sekvensering med långa läsningar har vissa särdrag (t.ex. högre felprocent) utvecklas särskilda verktyg för analys av denna typ av data.

Kvalitetskontroll (QC) av råa läsningar

Den ivriga analytikern kommer att börja analysen från FASTQ-filer; FASTQ-formatet har länge varit standard för att lagra sekvenseringsdata med korta läsningar. I huvudsak innehåller FASTQ-filer nukleotidsekvensen och kvaliteten per bas för miljontals läsningar. Även om filstorleken beror på det faktiska antalet läsningar är FASTQ-filer vanligtvis stora (i storleksordningen megabyte och gigabyte) och komprimerade. Det bör noteras att de flesta verktyg som använder FASTQ-filer som indata kan hantera dem i komprimerat format, så för att spara diskutrymme rekommenderas det att inte packa upp dem. Som en konvention kommer jag här att likställa en FASTQ-fil med ett sekvenseringsprov.

FastQC är sannolikt det mest populära verktyget för att utföra QC av de råa läsningarna. Det kan köras via ett visuellt gränssnitt eller programmatiskt. Det första alternativet kan vara bekvämare för användare som inte känner sig bekväma med kommandoradsmiljön, men det senare erbjuder ojämförlig skalbarhet och reproducerbarhet (tänk på hur tråkigt och felbenäget det kan vara att manuellt köra verktyget för tiotals filer). Hur som helst är det viktigaste resultatet av FastQC en HTML-fil som rapporterar viktig sammanfattande statistik om den övergripande kvaliteten på de obearbetade sekvenseringsavläsningarna från ett visst prov. Att granska tiotals FastQC-rapporter en och en är tråkigt och försvårar jämförelsen mellan olika prover. Därför kanske du vill använda MultiQC, som aggregerar HTML-rapporterna från FastQC (liksom från andra verktyg som används nedströms, t.ex. adaptertrimmning, justering) till en enda rapport.

QC-informationen är tänkt att göra det möjligt för användaren att bedöma om proverna har god kvalitet och därför kan användas för de efterföljande stegen eller om de måste kasseras. Tyvärr finns det inte något tröskelvärde som bygger på FastQC-metrikerna för att klassificera proverna som av god eller dålig kvalitet. Det tillvägagångssätt som jag använder är följande. Jag förväntar mig att alla prover som har genomgått samma förfarande (t.ex. DNA-extraktion, framställning av bibliotek) ska ha liknande kvalitetsstatistik och en majoritet av ”godkänd” flaggor. Om vissa prover har lägre kvalitet än genomsnittet kommer jag ändå att använda dem i analysen i efterföljande led med detta i åtanke. Å andra sidan, om alla prover i experimentet systematiskt får ”varningsflaggor” eller ”fail”-flaggor i flera mätvärden (se det här exemplet), misstänker jag att något gick fel i experimentet (t.ex. dålig DNA-kvalitet, bibliotekspreparering osv.) och rekommenderar att det upprepas.

Read trimming

QC av råa reads hjälper till att identifiera problematiska prover, men det förbättrar inte den faktiska kvaliteten på reads. För att göra det måste vi trimma reads för att ta bort tekniska sekvenser och ändar av låg kvalitet.

Tekniska sekvenser är rester från det experimentella förfarandet (t.ex. sekvenseringsadaptrar). Om sådana sekvenser gränsar till avläsningens sanna sekvens kan anpassningen (se nedan) kartlägga avläsningar till fel position i genomet eller minska förtroendet för en viss anpassning. Förutom tekniska sekvenser kan vi också vilja ta bort sekvenser av biologiskt ursprung om dessa är mycket närvarande bland läsningarna. Till exempel kan suboptimala förfaranden för DNA-preparering lämna en hög andel DNA-konverterat ribosomalt RNA (rRNA) i provet. Om inte denna typ av nukleinsyra är målet för sekvenseringsexperimentet, kommer att behålla reads som härrör från rRNA bara att öka den beräkningsmässiga bördan för stegen nedströms och kan förvirra resultaten. Om nivåerna av tekniska sekvenser, rRNA eller andra föroreningar är mycket höga, vilket förmodligen redan har framhävts av QC, kan det vara lämpligt att kassera hela sekvenseringsprovet.

I short-read-sekvensering bestäms DNA-sekvensen en nukleotid i taget (tekniskt sett en nukleotid varje sekvenseringscykel). Med andra ord är det antalet sekvenseringscykler som bestämmer avläsningslängden. Ett känt problem med HTS-sekvenseringsmetoderna är att noggrannheten med vilken nukleotiderna bestäms minskar när sekvenseringscyklerna ackumuleras. Detta återspeglas i en övergripande minskning av kvaliteten på anropet per bas, särskilt mot slutet av avläsningen. Precis som med tekniska sekvenser kan försök att anpassa läsningar som innehåller ändar av låg kvalitet leda till felplacering eller dålig kartläggningskvalitet.

För att avlägsna tekniska/kontaminerande sekvenser och ändar av låg kvalitet finns det verktyg för trimning av läsningar som Trimmomatic och Cutadapt, som används i stor utsträckning. I huvudsak tar sådana verktyg bort tekniska sekvenser (internt tillgängliga och/eller tillhandahållna av användaren) och trimmar avläsningar baserat på kvalitet samtidigt som de maximerar avläsningslängden. Läsningar som är för korta efter trimningen kasseras (alltför korta läsningar, t.ex. <36 nukleotider, komplicerar anpassningssteget eftersom de sannolikt kommer att kartlägga flera platser i genomet). Du kanske vill titta på andelen reads som överlever trimningen, eftersom en hög andel kasserade reads sannolikt är ett tecken på data av dålig kvalitet.

Slutligt kör jag vanligtvis FastQC på nytt på de trimmade läsningarna för att kontrollera att detta steg var effektivt och systematiskt förbättrade QC-metrikerna.

Alignment

Med undantag (t.ex. de novo assemblering) är alignment (även kallad mappning) typiskt sett nästa steg för de flesta HTS-datatyper och applikationer. Utläsningsanpassning består i att bestämma den position i genomet från vilken sekvensen för utläsningen härstammar (vanligtvis uttryckt som kromosom: start- och slutpunkt). I detta steg kräver vi därför att vi använder en referenssekvens för att anpassa/mappa avläsningarna till.

Valet av referenssekvens bestäms av flera faktorer. För det första den art som det sekvenserade DNA:t kommer från. Även om antalet arter med en tillgänglig referenssekvens av hög kvalitet ökar, kan detta fortfarande inte vara fallet för vissa mindre studerade organismer. I dessa fall kanske du vill anpassa läsningar till en evolutionärt nära art för vilken ett referensgenom finns tillgängligt. Eftersom det till exempel inte finns någon referenssekvens för prärievargens genom kan vi använda den närbesläktade hundens genom för läsarjusteringen. På samma sätt kan vi fortfarande vilja anpassa våra läsningar till en närbesläktad art för vilken det finns en referenssekvens av högre kvalitet. Även om till exempel gibbonens genom har publicerats är detta uppdelat i tusentals fragment som inte fullt ut återger organisationen av detta genom i tiotals kromosomer; i det fallet kan det vara fördelaktigt att utföra anpassningen med hjälp av den mänskliga referenssekvensen.

En annan faktor att ta hänsyn till är versionen av referenssekvenssamlingen, eftersom nya versioner släpps i takt med att sekvensen uppdateras och förbättras. Det är viktigt att koordinaterna för en viss justering kan variera mellan olika versioner. Till exempel finns flera versioner av det mänskliga genomet i UCSC Genome Browser. I alla arter är jag starkt för att migrera till den senaste versionen av sammansättningen när den är helt tillgänglig. Detta kan orsaka vissa olägenheter under övergången, eftersom redan existerande resultat kommer att vara relativa till äldre versioner, men det lönar sig i det långa loppet.

För övrigt spelar typen av sekvenseringsdata också roll. Avläsningar som genereras från DNA-seq-, ChIP-seq- eller Hi-C-protokoll kommer att anpassas till genomets referenssekvens. Eftersom RNA som transkriberats från DNA bearbetas vidare till mRNA (dvs. introner avlägsnas) kommer många RNA-seq-avläsningar å andra sidan att misslyckas med att anpassa sig till en referenssekvens för genomet. I stället måste vi antingen anpassa dem till transkriptomreferenssekvenser eller använda split-aware aligners (se nedan) när vi använder genomsekvensen som referens. Relaterat till detta är valet av källa för annotering av referenssekvensen, dvs. databasen med koordinaterna för gener, transkript, centromerer osv. Jag använder vanligtvis GENCODE-annotationen eftersom den kombinerar omfattande genannotationer och transkriptsekvenser.

En lång rad verktyg för kortläsningssekvensanpassning har utvecklats (se avsnittet Kortläsningssekvensanpassning här). Att granska dem ligger utanför den här artikelns räckvidd (detaljer om algoritmerna bakom dessa verktyg finns här). Enligt min erfarenhet är Bowtie2, BWA, HISAT2, Minimap2, STAR och TopHat bland de mest populära. Min rekommendation är att du väljer din alignerare med hänsyn till nyckelfaktorer som typ av HTS-data och tillämpning samt acceptans i samhället, dokumentationens kvalitet och antalet användare. Man behöver t.ex. aligners som STAR eller Bowtie2 som är medvetna om exon-exon-övergångar när man mappar RNA-seq till genomet.

För de flesta kartläggare är det nödvändigt att indexera den sekvens som används som referens innan den faktiska anpassningen äger rum. Detta steg kan vara tidskrävande, men det behöver bara göras en gång för varje referenssekvens. De flesta kartläggare lagrar anpassningar i SAM/BAM-filer som följer SAM/BAM-formatet (BAM-filer är binära versioner av SAM-filer). Anpassningen är ett av de mest beräknings- och tidskrävande stegen i analysen av sekvenseringsdata och SAM/BAM-filerna är tunga (i storleksordningen gigabyte). Det är därför viktigt att se till att du har de resurser som krävs (se det sista avsnittet nedan) för att köra anpassningen på en rimlig tid och lagra resultaten. På samma sätt bör man på grund av BAM-filernas storlek och binära format undvika att öppna dem med textredigerare, utan istället använda Unix-kommandon eller särskilda verktyg som SAMtools.

Från anpassningarna

Jag skulle vilja säga att det inte finns något tydligt gemensamt steg efter anpassningen, eftersom det är vid den här punkten som varje HTS-datatyp och tillämpning kan skilja sig åt.

En vanlig nedströmsanalys för DNA-seq-data är variant calling, dvs. identifiering av positioner i genomet som varierar i förhållande till genomreferensen och mellan individer. Ett populärt analysramverk för denna tillämpning är GATK för polymorfism av en enda nukleotid (SNP) eller små instick/avskiljningar (indels) (figur 2). Varianter som består av större delar av DNA (även kallade strukturella varianter) kräver särskilda metoder för att kalla in varianter (se den här artikeln för en omfattande jämförelse). Liksom när det gäller aligners rekommenderar jag att man väljer rätt verktyg med hänsyn till nyckelfaktorer som typ av varianter (SNP, indel eller strukturella varianter), acceptans av samhället, dokumentationens kvalitet och antalet användare.

Det vanligaste användningsområdet för RNA-seq är troligen kvantifiering av genuttryck. Historiskt sett behövde läsningar anpassas till referenssekvensen och sedan användes antalet läsningar som anpassats till en viss gen eller transkript som en proxy för att kvantifiera dess uttrycksnivåer. Detta tillvägagångssätt för anpassning och kvantifiering utförs av verktyg som Cufflinks, RSEM eller featureCounts. Scuh-metoden har dock i allt högre grad överträffats av nyare metoder som genomförs i programvaror som Kallisto och Salmon. Med sådana verktyg behöver inte hela sekvensen av en läsning anpassas till referenssekvensen. I stället behöver vi bara anpassa tillräckligt många nukleotider för att vara säkra på att en avläsning kommer från ett visst transkript. Enkelt uttryckt kan man säga att metoden för anpassning och kvantifiering reduceras till ett enda steg. Detta tillvägagångssätt är känt som pseudokartläggning och ökar kraftigt hastigheten på kvantifieringen av genuttryck. Å andra sidan bör man komma ihåg att pseudomappning inte lämpar sig för tillämpningar där en fullständig anpassning behövs (t.ex. variantbestämning från RNA-seq-data).

Ett annat exempel på skillnaderna i de efterföljande analysstegen och de nödvändiga verktygen i olika sekvenseringsbaserade tillämpningar är ChIP-seq. Läser som genereras med denna teknik kommer att användas för peak calling, vilket innebär att man upptäcker regioner i genomet med ett betydande överskott av läsningar som indikerar var målproteinet är bundet. Det finns flera peak callers och i denna publikation ges en översikt över dem. Som ett sista exempel nämner jag Hi-C-data, där anpassningar används som indata för verktyg som bestämmer interaktionsmatriserna och, utifrån dessa, genomets 3D-funktioner. Att kommentera alla sekvenseringsbaserade analyser ligger utanför den här artikelns räckvidd (för en relativt fullständig lista se den här artikeln).

För att börja…

Den återstående delen av den här artikeln berör aspekter som kanske inte strikt betraktas som steg i analysen av HTS-data och som till stor del ignoreras. Däremot hävdar jag att det är viktigt att du tänker på de frågor som ställs i tabell 1 innan du börjar analysera HTS-data (eller vilken typ av data som helst förvisso), och jag har skrivit om dessa ämnen här och här.

Tabell 1

Tänk på saken	Föreslagna åtgärder
Har du all den information om ditt prov som behövs för analysen?	Samla systematiskt in metadata om försöken
Kommer du att entydigt kunna identifiera ditt prov?	Inrätta ett system för att tilldela varje prov en unik identitetsbeteckning
Var kommer data och resultat att finnas?	Strukturerad och hierarkisk organisering av data
Kommer du att kunna bearbeta flera prover sömlöst?	Skalerbarhet, parallellisering, automatisk konfiguration och moduläritet hos koden
Kommer du eller någon annan att kunna reproducera resultaten?	Dokumentera din kod och dina förfaranden!

Som nämnts ovan är HTS-rådata och vissa av de filer som genereras under analysen av dem i storleksordningen gigabyte, så det är inte ovanligt att ett projekt som omfattar tiotals prover kräver terabyte lagring. Dessutom är vissa steg i analysen av HTS-data beräkningsintensiva (t.ex. anpassning). Den lagrings- och beräkningsinfrastruktur som krävs för att analysera HTS-data är dock ett viktigt övervägande som ofta förbises eller inte diskuteras. Som ett exempel kan nämnas att vi som en del av en nyligen genomförd analys granskade tiotals publicerade artiklar som utförde fenomövergripande associationsanalyser (PheWAS). Moderna PheWAS analyserar 100-1 000-talet både genetiska varianter och fenotyper, vilket resulterar i betydande datalagring och datorkraft. Ändå kommenterade praktiskt taget ingen av de artiklar vi granskade den infrastruktur som behövs för PheWAS-analysen. Inte överraskande är min rekommendation att ni i förväg planerar de lagrings- och datakrav som ni kommer att ställas inför och att ni delar med er av dem till samhället.

Behövs det hjälp med att analysera DNA-sekvenseringsdata? Kom i kontakt med frilansande bioinformatiker och genomik-experter på Kolabtree.