Ein schrittweiser Leitfaden für die Analyse von DNA-Sequenzierungsdaten

Dr. Javier Quilez Oliete, ein erfahrener Bioinformatik-Berater auf Kolabtree, bietet einen umfassenden Leitfaden für die Analyse von DNA-Sequenzierungsdaten, einschließlich der zum Lesen der Daten verwendeten Werkzeuge und Software.

Einführung

Desoxyribonukleinsäure (DNA) ist das Molekül, das den größten Teil der genetischen Information eines Organismus trägt. (Bei einigen Virustypen wird die genetische Information durch Ribonukleinsäure (RNA) übertragen.) Nukleotide (üblicherweise durch die Buchstaben A, C, G oder T dargestellt) sind die Grundeinheiten von DNA-Molekülen. Bei der DNA-Sequenzierung werden die Nukleotide, aus denen ein DNA-Molekül besteht, ausgelesen (z. B. „GCAAACCAAT“ ist eine DNA-Kette mit 10 Nukleotiden). Mit den heutigen Sequenzierungstechnologien lassen sich Millionen solcher DNA-Lesevorgänge in angemessener Zeit und zu relativ geringen Kosten herstellen. Die Kosten für die Sequenzierung eines menschlichen Genoms – ein Genom ist der vollständige Satz von DNA-Molekülen in einem Organismus – liegen inzwischen unter der 100-Dollar-Grenze und können innerhalb weniger Tage durchgeführt werden. Dies steht im Gegensatz zur ersten Initiative zur Sequenzierung des menschlichen Genoms, die in einem Jahrzehnt abgeschlossen wurde und etwa 2,7 Milliarden Dollar kostete.

Diese Fähigkeit, DNA mit hohem Durchsatz und zu geringen Kosten zu sequenzieren, hat die Entwicklung einer wachsenden Zahl von auf Sequenzierung basierenden Methoden und Anwendungen ermöglicht. So kann beispielsweise die Sequenzierung ganzer Genome oder ihrer proteinkodierenden Regionen (zwei Ansätze, die als Ganzgenom- bzw. Exom-Sequenzierung bekannt sind) bei kranken und gesunden Personen auf krankheitsverursachende DNA-Veränderungen hinweisen. Auch die Sequenzierung der RNA, die von der DNA transkribiert wird – eine Technik, die als RNA-Sequenzierung bekannt ist – wird verwendet, um die Genaktivität zu quantifizieren und festzustellen, wie sich diese unter verschiedenen Bedingungen (z. B. unbehandelt und behandelt) verändert. Auf der anderen Seite werden mit Chromosomen-Konformations-Capture-Sequenzierungsmethoden Wechselwirkungen zwischen nahegelegenen DNA-Molekülen nachgewiesen und so die räumliche Verteilung von Chromosomen innerhalb der Zelle bestimmt.

Gemeinsam mit diesen und anderen Anwendungen der DNA-Sequenzierung ist die Erzeugung von Datensätzen in der Größenordnung von Gigabytes, die Millionen von Lesesequenzen umfassen. Daher erfordert die Auswertung von Hochdurchsatz-Sequenzierungsexperimenten (HTS) umfangreiche Datenanalysefähigkeiten. Glücklicherweise gibt es für die meisten HTS-Datentypen spezielle Berechnungs- und Statistik-Tools und relativ standardisierte Analyse-Workflows. Während einige der (anfänglichen) Analyseschritte für die meisten Sequenzierungsdatentypen gleich sind, hängt die weitere nachgeschaltete Analyse von der Art der Daten und/oder dem letztendlichen Ziel der Analyse ab. Im Folgenden gebe ich eine Einführung in die grundlegenden Schritte bei der Analyse von HTS-Daten und verweise auf gängige Tools.

Einige der folgenden Abschnitte konzentrieren sich auf die Analyse von Daten, die mit Short-Read-Sequenzierungstechnologien (hauptsächlich Illumina) erzeugt wurden, da diese in der Vergangenheit den HTS-Markt dominiert haben. Neuere Technologien, die längere Reads erzeugen (z. B. Oxford Nanopore Technologies, PacBio), sind jedoch auf dem Vormarsch. Da die Long-Read-Sequenzierung einige Besonderheiten aufweist (z. B. höhere Fehlerquoten), werden derzeit spezielle Tools für die Analyse dieser Art von Daten entwickelt.

Qualitätskontrolle (QC) von Rohdaten

Der eifrige Analytiker wird die Analyse mit FASTQ-Dateien beginnen; das FASTQ-Format ist seit langem der Standard für die Speicherung von Short-Read-Sequenzierungsdaten. Im Wesentlichen enthalten FASTQ-Dateien die Nukleotidsequenz und die Qualität des Callings pro Base für Millionen von Reads. Obwohl die Dateigröße von der tatsächlichen Anzahl der Reads abhängt, sind FASTQ-Dateien in der Regel groß (in der Größenordnung von Megabytes und Gigabytes) und komprimiert. Die meisten Tools, die FASTQ-Dateien als Eingabe verwenden, können sie in komprimiertem Format verarbeiten. Um Speicherplatz zu sparen, wird empfohlen, sie nicht zu dekomprimieren. Als Konvention werde ich hier eine FASTQ-Datei mit einer Sequenzierungsprobe gleichsetzen.

FastQC ist wahrscheinlich das beliebteste Tool zur Durchführung der Qualitätskontrolle von Rohdaten. Es kann über eine visuelle Schnittstelle oder programmatisch ausgeführt werden. Während die erste Option für Benutzer, die sich mit der Befehlszeilenumgebung nicht wohlfühlen, bequemer ist, bietet die zweite eine unvergleichliche Skalierbarkeit und Reproduzierbarkeit (man denke nur daran, wie mühsam und fehleranfällig es sein kann, das Tool für Dutzende von Dateien manuell auszuführen). In jedem Fall ist die Hauptausgabe von FastQC eine HTML-Datei mit den wichtigsten zusammenfassenden Statistiken über die Gesamtqualität der Rohsequenzierungs-Reads aus einer bestimmten Probe. Es ist mühsam, Dutzende von FastQC-Berichten einzeln zu prüfen, und es erschwert den Vergleich zwischen verschiedenen Proben. Daher sollten Sie MultiQC verwenden, das die HTML-Berichte von FastQC (sowie von anderen nachgeschalteten Tools, z. B. Adapter Trimming, Alignment) in einem einzigen Bericht zusammenfasst.

QC-Informationen sollen es dem Benutzer ermöglichen zu beurteilen, ob Proben eine gute Qualität aufweisen und daher für die nachfolgenden Schritte verwendet werden können oder ob sie verworfen werden müssen. Leider gibt es keinen einheitlichen Schwellenwert auf der Grundlage der FastQC-Metriken zur Einstufung von Proben in gute oder schlechte Qualität. Der von mir verwendete Ansatz ist der folgende. Ich erwarte, dass alle Proben, die dasselbe Verfahren durchlaufen haben (z. B. DNA-Extraktion, Bibliotheksvorbereitung), ähnliche Qualitätsstatistiken und eine Mehrheit von „Pass“-Flags aufweisen. Wenn einige Proben eine unterdurchschnittliche Qualität aufweisen, verwende ich sie dennoch in der nachgelagerten Analyse, wobei ich dies berücksichtige. Wenn andererseits alle Proben in einem Experiment systematisch „Warnung“ oder „Nicht bestanden“ in mehreren Metriken erhalten (siehe dieses Beispiel), vermute ich, dass bei dem Experiment etwas schief gelaufen ist (z. B. schlechte DNA-Qualität, Bibliotheksvorbereitung usw.), und ich empfehle, es zu wiederholen.

Read-Trimming

QC von Roh-Reads hilft, problematische Proben zu identifizieren, verbessert aber nicht die tatsächliche Qualität der Reads. Dazu müssen die Reads getrimmt werden, um technische Sequenzen und qualitativ schlechte Enden zu entfernen.

Technische Sequenzen sind Überbleibsel des experimentellen Verfahrens (z. B. Sequenzieradapter). Wenn solche Sequenzen an die tatsächliche Sequenz des Reads angrenzen, kann das Alignment (siehe unten) die Reads der falschen Position im Genom zuordnen oder das Vertrauen in ein bestimmtes Alignment verringern. Neben technischen Sequenzen kann es auch sinnvoll sein, Sequenzen biologischen Ursprungs zu entfernen, wenn diese in den Reads stark vertreten sind. So kann beispielsweise bei suboptimalen DNA-Präparationsverfahren ein hoher Anteil an DNA-umgewandelter ribosomaler RNA (rRNA) in der Probe zurückbleiben. Sofern diese Art von Nukleinsäure nicht das Ziel des Sequenzierungsexperiments ist, erhöht die Beibehaltung von Reads, die von rRNA stammen, nur den Rechenaufwand der nachgeschalteten Schritte und kann die Ergebnisse verfälschen. Wenn der Anteil an technischen Sequenzen, rRNA oder anderen Verunreinigungen sehr hoch ist, was wahrscheinlich bereits bei der Qualitätskontrolle festgestellt wurde, sollten Sie die gesamte Sequenzierprobe verwerfen.

Bei der Short-Read-Sequenzierung wird die DNA-Sequenz Nukleotid für Nukleotid bestimmt (technisch gesehen ein Nukleotid pro Sequenzierzyklus). Mit anderen Worten: Die Anzahl der Sequenzierzyklen bestimmt die Leselänge. Ein bekanntes Problem der HTS-Sequenzierungsmethoden ist die Abnahme der Genauigkeit, mit der die Nukleotide bestimmt werden, wenn die Sequenzierungszyklen zunehmen. Dies spiegelt sich in einer allgemeinen Abnahme der Qualität des Callings pro Base wider, insbesondere gegen Ende des Reads. Wie bei technischen Sequenzen kann der Versuch, Reads mit minderwertigen Enden zu alignieren, zu Fehlplatzierungen oder schlechter Mapping-Qualität führen.

Um technische/verunreinigende Sequenzen und minderwertige Enden zu entfernen, gibt es Read-Trimming-Tools wie Trimmomatic und Cutadapt, die weit verbreitet sind. Im Wesentlichen entfernen diese Tools technische Sequenzen (intern verfügbar und/oder vom Benutzer bereitgestellt) und trimmen Reads auf der Grundlage der Qualität bei gleichzeitiger Maximierung der Leselänge. Reads, die nach dem Trimming zu kurz sind, werden verworfen (zu kurze Reads, z. B. <36 Nukleotide, erschweren den Alignment-Schritt, da diese wahrscheinlich mehreren Stellen im Genom zugeordnet werden). Sie sollten sich den Prozentsatz der Reads ansehen, die das Trimming überleben, da eine hohe Rate an verworfenen Reads wahrscheinlich ein Zeichen für schlechte Datenqualität ist.

Schließlich führe ich in der Regel FastQC erneut auf den getrimmten Reads aus, um zu überprüfen, ob dieser Schritt effektiv war und die QC-Metriken systematisch verbessert hat.

Ausrichtung

Mit Ausnahmen (z. B. De-novo-Assembly) ist die Ausrichtung (auch als Mapping bezeichnet) bei den meisten HTS-Datentypen und -Anwendungen der nächste Schritt. Das Read-Alignment besteht darin, die Position im Genom zu bestimmen, von der die Sequenz des Reads stammt (normalerweise ausgedrückt als Chromosom:Start-Ende). Daher benötigen wir in diesem Schritt eine Referenzsequenz zum Alignment/Mapping der Reads.

Die Wahl der Referenzsequenz wird durch mehrere Faktoren bestimmt. Zum einen von der Spezies, von der die sequenzierte DNA stammt. Während die Zahl der Arten, für die eine qualitativ hochwertige Referenzsequenz zur Verfügung steht, zunimmt, kann dies bei einigen weniger erforschten Organismen noch nicht der Fall sein. In diesen Fällen sollten Sie die Reads an einer evolutiv nahen Art ausrichten, für die ein Referenzgenom verfügbar ist. Da es zum Beispiel keine Referenzsequenz für das Genom des Kojoten gibt, können wir für das Read-Alignment die Sequenz des nahe verwandten Hundes verwenden. Ebenso kann es sein, dass wir unsere Reads an einer eng verwandten Art ausrichten wollen, für die eine qualitativ hochwertigere Referenzsequenz existiert. Das Genom des Gibbons wurde zwar veröffentlicht, ist aber in Tausende von Fragmenten unterteilt, die die Organisation dieses Genoms in Dutzende von Chromosomen nicht vollständig wiedergeben. In diesem Fall kann es von Vorteil sein, das Alignment anhand der menschlichen Referenzsequenz durchzuführen.

Ein weiterer zu berücksichtigender Faktor ist die Version der Referenzsequenz, da neue Versionen veröffentlicht werden, wenn die Sequenz aktualisiert und verbessert wird. Wichtig ist, dass die Koordinaten eines bestimmten Alignments zwischen den Versionen variieren können. Im UCSC Genome Browser finden sich zum Beispiel mehrere Versionen des menschlichen Genoms. Bei allen Spezies empfehle ich dringend, auf die neueste Assembly-Version umzusteigen, sobald diese vollständig freigegeben ist. Dies kann während der Umstellung zu einigen Unannehmlichkeiten führen, da sich bereits vorhandene Ergebnisse auf ältere Versionen beziehen, aber langfristig zahlt es sich aus.

Außerdem spielt auch die Art der Sequenzierungsdaten eine Rolle. Reads, die aus DNA-seq-, ChIP-seq- oder Hi-C-Protokollen stammen, werden an der Genomreferenzsequenz ausgerichtet. Da jedoch die von der DNA transkribierte RNA zu mRNA weiterverarbeitet wird (d. h. die Introns werden entfernt), können viele RNA-seq-Reads nicht an einer Genomreferenzsequenz ausgerichtet werden. Stattdessen müssen wir sie entweder an Transkriptom-Referenzsequenzen ausrichten oder split-aware Aligner verwenden (siehe unten), wenn wir die Genomsequenz als Referenz verwenden. Damit verbunden ist die Wahl der Quelle für die Annotation der Referenzsequenz, d. h. der Datenbank mit den Koordinaten der Gene, Transkripte, Zentromere usw. In der Regel verwende ich die GENCODE-Annotation, da sie eine umfassende Genannotation und Transkriptsequenzen kombiniert.

Eine lange Liste von Short-Read-Sequenz-Alignment-Tools wurde entwickelt (siehe den Abschnitt Short-Read-Sequenz-Alignment hier). Sie zu besprechen, würde den Rahmen dieses Artikels sprengen (Details zu den Algorithmen hinter diesen Tools finden Sie hier). Nach meiner Erfahrung gehören zu den beliebtesten Tools Bowtie2, BWA, HISAT2, Minimap2, STAR und TopHat. Meine Empfehlung ist, dass Sie Ihren Aligner unter Berücksichtigung von Schlüsselfaktoren wie der Art der HTS-Daten und der Anwendung sowie der Akzeptanz in der Gemeinschaft, der Qualität der Dokumentation und der Anzahl der Nutzer auswählen. Man braucht z.B. Aligner wie STAR oder Bowtie2, die Exon-Exon-Verbindungen beim Mapping von RNA-seq auf das Genom berücksichtigen.

Gemeinsam mit den meisten Mappern ist die Notwendigkeit, die als Referenz verwendete Sequenz zu indizieren, bevor das eigentliche Alignment stattfindet. Dieser Schritt kann zeitaufwendig sein, muss aber nur einmal für jede Referenzsequenz durchgeführt werden. Die meisten Mapper speichern Alignments in SAM/BAM-Dateien, die dem SAM/BAM-Format entsprechen (BAM-Dateien sind binäre Versionen von SAM-Dateien). Das Alignment gehört zu den rechenintensivsten und zeitaufwändigsten Schritten bei der Analyse von Sequenzierungsdaten, und SAM/BAM-Dateien sind sehr groß (in der Größenordnung von Gigabytes). Daher müssen Sie sicherstellen, dass Sie über die erforderlichen Ressourcen verfügen (siehe den letzten Abschnitt unten), um das Alignment in einer angemessenen Zeit durchzuführen und die Ergebnisse zu speichern. Aufgrund der Größe und des Binärformats von BAM-Dateien sollten Sie es vermeiden, sie mit Texteditoren zu öffnen; verwenden Sie stattdessen Unix-Befehle oder spezielle Tools wie SAMtools.

Aus den Alignments

würde ich sagen, dass es keinen eindeutigen gemeinsamen Schritt nach dem Alignment gibt, da sich an diesem Punkt jeder HTS-Datentyp und jede Anwendung unterscheiden kann.

Eine gängige nachgelagerte Analyse für DNA-seq-Daten ist das Varianten-Calling, d. h. die Identifizierung von Positionen im Genom, die im Vergleich zur Genomreferenz und zwischen Individuen variieren. Ein beliebtes Analysesystem für diese Anwendung ist GATK für Einzelnukleotid-Polymorphismus (SNP) oder kleine Insertionen/Deletionen (Indels) (Abbildung 2). Varianten, die aus größeren DNA-Abschnitten bestehen (auch als Strukturvarianten bezeichnet), erfordern spezielle Calling-Methoden (siehe diesen Artikel für einen umfassenden Vergleich). Wie bei den Alignern empfehle ich, das richtige Tool unter Berücksichtigung von Schlüsselfaktoren wie der Art der Varianten (SNP, Indel oder Strukturvarianten), der Akzeptanz in der Gemeinschaft, der Qualität der Dokumentation und der Anzahl der Nutzer auszuwählen.

Die wahrscheinlich häufigste Anwendung von RNA-seq ist die Quantifizierung der Genexpression. In der Vergangenheit mussten die Reads an die Referenzsequenz angeglichen werden, und dann wurde die Anzahl der Reads, die an ein bestimmtes Gen oder Transkript angeglichen wurden, als Näherungswert für die Quantifizierung seiner Expressionswerte verwendet. Dieser Alignment+Quantifizierungsansatz wird von Tools wie Cufflinks, RSEM oder featureCounts durchgeführt. Der scuh-Ansatz wird jedoch zunehmend von neueren Methoden, die in Software wie Kallisto und Salmon implementiert sind, verdrängt. Bei diesen Tools muss nicht die gesamte Sequenz eines Read an die Referenzsequenz angeglichen werden. Stattdessen müssen wir nur genügend Nukleotide alignieren, um sicher zu sein, dass ein Read von einem bestimmten Transkript stammt. Vereinfacht ausgedrückt, wird der Ansatz von Alignment und Quantifizierung auf einen einzigen Schritt reduziert. Dieser Ansatz wird als Pseudo-Mapping bezeichnet und erhöht die Geschwindigkeit der Quantifizierung der Genexpression erheblich. Auf der anderen Seite ist zu beachten, dass Pseudo-Mapping nicht für Anwendungen geeignet ist, bei denen ein vollständiges Alignment erforderlich ist (z. B. Variantenaufruf aus RNA-seq-Daten).

Ein weiteres Beispiel für die Unterschiede bei den nachgeschalteten Analyseschritten und den erforderlichen Tools bei sequenzierungsbasierten Anwendungen ist ChIP-seq. Mit dieser Technik erzeugte Reads werden für das Peak-Calling verwendet, das darin besteht, Regionen im Genom mit einem signifikanten Überschuss an Reads aufzuspüren, die anzeigen, wo das Zielprotein gebunden ist. Es gibt mehrere Peak-Caller, die in dieser Veröffentlichung vorgestellt werden. Als letztes Beispiel erwähne ich Hi-C-Daten, bei denen Alignments als Input für Tools verwendet werden, die die Interaktionsmatrizen und daraus die 3D-Merkmale des Genoms bestimmen. Es würde den Rahmen dieses Artikels sprengen, auf alle sequenzbasierten Assays einzugehen (eine relativ vollständige Liste finden Sie in diesem Artikel).

Bevor Sie beginnen…

Der verbleibende Teil dieses Artikels berührt Aspekte, die vielleicht nicht unbedingt als Schritte bei der Analyse von HTS-Daten angesehen werden und die weitgehend ignoriert werden. Im Gegensatz dazu halte ich es für wichtig, dass Sie über die in Tabelle 1 gestellten Fragen nachdenken, bevor Sie mit der Analyse von HTS-Daten (oder überhaupt von Daten) beginnen, und ich habe über diese Themen hier und hier geschrieben.

Tabelle 1

Denken Sie darüber nach	Vorgeschlagene Maßnahmen
Haben Sie alle für die Analyse benötigten Informationen über Ihre Stichprobe?	Sammeln Sie systematisch die Metadaten der Experimente
Werden Sie in der Lage sein, Ihre Probe eindeutig zu identifizieren?	Einrichten Sie ein System, um jeder Probe eine eindeutige Kennung zuzuweisen
Wo werden Daten und Ergebnisse gespeichert?	Strukturierte und hierarchische Organisation der Daten
Werden Sie in der Lage sein, mehrere Proben nahtlos zu verarbeiten?	Skalierbarkeit, Parallelisierung, automatische Konfiguration und Modularität des Codes
Werden Sie oder andere Personen in der Lage sein, die Ergebnisse zu reproduzieren?	Dokumentieren Sie Ihren Code und Ihre Verfahren!

Wie bereits erwähnt, liegen die HTS-Rohdaten und einige der bei ihrer Analyse erzeugten Dateien in der Größenordnung von Gigabytes, so dass es nicht ungewöhnlich ist, dass ein Projekt mit Dutzenden von Proben Terabytes an Speicherplatz benötigt. Außerdem sind einige Schritte bei der Analyse von HTS-Daten rechenintensiv (z. B. das Alignment). Die für die Analyse von HTS-Daten erforderliche Speicher- und Datenverarbeitungsinfrastruktur ist jedoch ein wichtiger Aspekt, der oft übersehen oder nicht diskutiert wird. Als Beispiel haben wir im Rahmen einer kürzlich durchgeführten Analyse Dutzende von veröffentlichten Arbeiten durchgesehen, die phänomenweite Assoziationsanalysen (PheWAS) durchführen. Bei modernen PheWAS werden 100-1.000 genetische Varianten und Phänotypen analysiert, was zu einer erheblichen Datenspeicherung und Rechenleistung führt. Dennoch ging praktisch keine der von uns geprüften Arbeiten auf die für die PheWAS-Analyse erforderliche Infrastruktur ein. Es überrascht daher nicht, dass ich empfehle, die Speicher- und Rechenanforderungen im Voraus zu planen und sie mit der Gemeinschaft zu teilen.

Brauchen Sie Hilfe bei der Analyse von DNA-Sequenzierungsdaten? Nehmen Sie Kontakt mit freiberuflichen Bioinformatikern und Genomik-Experten auf Kolabtree auf.