FDA SEQC-Benchmark-Datensätze
Der FDA SEQC-Benchmark-Datensatz (Gene Expression Omnibus Accession Number GSE47792) enthält gepaarte RNA-seq-Daten, die mit der Illumina HiSeq 2000-Plattform mit der Leselänge von 100 Nukleotiden7 generiert wurden. Wir verwendeten eine Teilmenge des SEQC-Benchmark-Datensatzes, der an zwei Standorten sequenziert wurde – Beijing Genomics Institute (BGI) und Mayo Clinic (MAYO). Wir verwendeten vier Proben (d. H. A, B, C und D) mit jeweils vier Replikatbibliotheken, die an den Sequenzierungsstellen vorbereitet wurden. Probe A enthält die Universal Human Reference RNA (UHRR), Probe B enthält die Human Brain Reference RNA (HBRR), Probe C enthält eine Mischung aus A und B (75% A und 25% B) und Probe D enthält eine Mischung aus A und B (25% A und 75% B). Wir verwendeten Daten von zwei Spuren einer einzelnen Durchflusszelle für jedes Probenreplikat. Der SEQC lieferte auch einen quantitativen PCR-Benchmark-Datensatz (qPCR), der 20.801 Gene enthält, die mit PrimePCR (Bio-Rad, Hercules, Kalifornien) getestet wurden. Jedes PrimePCR-Gen wurde einmal für jede der vier Proben (d. H. A, B, C und D) getestet. Die FDA SEQC Benchmark-Datensätze und Proben sind in ergänzenden Tabellen S5 und S6 zusammengefasst.
Neuroblastom- und Lungenadenokarzinom-Datensätze
Wir verwendeten einen Neuroblastom-Datensatz mit 176 Stichproben (eine Teilmenge eines größeren Datensatzes mit 498 Stichproben; in diesem Manuskript als SEQC-Neuroblastom bezeichnet), um die Leistung von RNA-seqc in Bezug auf die genexpressionsbasierte Vorhersage des Krankheitsausganges zu bewerten. Diese Proben wurden von der Universitätskinderklinik Köln zur Verfügung gestellt und am BGI mit der Illumina-Plattform sequenziert48. Alle 176 Proben wurden von Hochrisikopatienten entnommen, die entweder mit Neuroblastom im Stadium 4 und Alter > 18 Monate oder mit MYCN-amplifizierten Tumoren in jedem Stadium oder Alter definiert waren. Der SEQC-Neuroblastom-Datensatz wurde mit der Zugangsnummer GSE47792 in den Gene Expression Omnibus hinterlegt.
Wir haben zwei klinische Endpunkte vorhergesagt — das ereignisfreie Überleben (EFS), dh das Auftreten von Ereignissen wie Fortschritt, Rückfall oder Tod, und das Gesamtüberleben (OS), dh den Tod. Für beide Endpunkte wurden die Patienten in zwei Gruppen eingeteilt (d. H. Hohe Risiken versus niedrige Risiken). Hochrisikopatienten erlebten ein Ereignis oder starben vor einer vordefinierten Überlebenszeitschwelle, während Patienten mit niedrigem Risiko ein Ereignis erlebten oder nach der Schwelle starben oder ihre letzte Nachsorge die Schwelle überschritt. Die Überlebenszeitschwellen für EFS und OS lagen bei zwei bzw. drei Jahren. Die Schwellenwerte wurden gewählt, um die Anzahl der Patienten mit hohem und niedrigem Risiko auszugleichen. Einzelheiten des SEQC-Neuroblastom-Datensatzes sind in der ergänzenden Tabelle S9 angegeben.
Wir verwendeten auch einen 87-Proben-Lungenadenokarzinom-RNA-seq-Datensatz aus dem Cancer Genome Atlas (TCGA) Repository. Der Vorhersageendpunkt war auch das Überleben, und wir verwendeten die gleichen Kriterien, um Hochrisiko- und Niedrigrisikogruppen mit der Überlebenszeitschwelle von zwei Jahren zu definieren. Die Zweijahresschwelle wurde gewählt, um die Anzahl der Patienten mit hohem und niedrigem Risiko auszugleichen. Einzelheiten des TCGA-Lungenadenokarzinom-Datensatzes sind in der ergänzenden Tabelle S10 angegeben.
Filtern des qPCR-Benchmark-Datensatzes, um einen Referenzsatz von Genen zu erzeugen
Aufgrund der Variabilität der qPCR-Messungen und der Meinungsverschiedenheiten zwischen den qPCR-Plattformen7 haben wir den qPCR-Benchmark-Datensatz gefiltert, um Gene beizubehalten, die „korrektes“ Verhalten zeigten. Wir haben dann diese Gene verwendet, um die Benchmark-Metriken (d. H. Genauigkeit, Präzision, Zuverlässigkeit und Reproduzierbarkeit) zu berechnen. Ein solches Filterverfahren ist in der ergänzenden Fig. S1.
Beginnend mit dem anfänglichen Satz von 20.801 Genen, die mit PrimePCR getestet wurden, filterten wir diese Gene, um nur Gene zu behalten, die als ungleich Null quantifiziert wurden (d. h. detektiert wurden) und mit Ct-Werten (Zyklusschwelle) ≤ 35 (35 zeigt den Nachweis nur eines einzelnen Moleküls in einer Probe an). Das Filtern von PrimePCR-Daten führte zu 14.014 Genen, die auch mit dem AceView-Transkriptom übereinstimmten, das für die Kartierung des SEQC-Benchmark-RNA-seq-Datensatzes verwendet wurde.Anschließend filterten wir die 14.014 qPCR-Gene, um nur 12.610 Gene zu erhalten, die die richtige Titrationsreihenfolge (TO) und die erwarteten Mischungsverhältnisse (EMR) aufwiesen. Einzelheiten zu diesem Prozess finden Sie im Abschnitt „Filtern von qPCR-Genen nach Titrationsreihenfolge und erwarteten Mischungsverhältnissen“.Schließlich, da einige Benchmark-Metriken wie Genauigkeit und Präzision empfindlich auf Null- oder sehr niedrig exprimierende Gene reagieren, haben wir weiter Gene ausgewählt, die in allen Replikaten aller Proben aller Sequenzierungsstellen und aller 278 RNA-seq-Pipelines als Nicht-Null exprimiert wurden. Der endgültige Referenzsatz enthält nur 10.222 qPCR-Gene (als „alle Gene“ bezeichnet), die zur Berechnung aller drei Benchmark-Metriken für RNA-seq-Pipelines verwendet wurden.
Basierend auf der vorherigen Studie sind die Gene mit geringerer Expression eher inkonsistent unter den Pipelinen49. So identifizierten wir auch eine Reihe von niedrig exprimierenden Genen in den 10.222 Genen basierend auf der durchschnittlichen qPCR-Expression der Proben A, B, C und D. Die niedrigsten 20% der 10.222 Gene (d. H. 2044 Gene, bezeichnet als „niedrig exprimierende Gene“) wurden auch verwendet, um den gleichen Satz von Benchmark-Metriken für RNA-seq-Pipelines zu berechnen. Dieses Design ermöglichte es uns, die Fähigkeit von RNA-seq-Pipelines zur Schätzung der niedrig exprimierenden Genexpression zu untersuchen.
Filterung von qPCR-Genen nach Titrationsreihenfolge und erwarteten Mischungsverhältnissen
Die SEQC-Benchmark-Datensätze (RNA-seq und qPCR) haben einzigartige Eigenschaften, die die Beurteilung der Quantifizierungskorrektheit ermöglichen. Nach der Identifizierung nachweisbarer (d. H. Nicht Null und Ct ≤ 35) und AceView-übereinstimmender qPCR-Gene verwendeten wir zwei Metriken (TO und EMR), um den Benchmark-qPCR-Datensatz weiter zu filtern, wobei nur „korrekte“ qPCR-Gene übrig blieben. Die TO- und EMR-Metriken erfassen eindeutige Mischeigenschaften der Daten, d. h.
Aufgrund dieser Eigenschaft wird erwartet, dass alle Gene in einer der folgenden Ordnungen exprimiert werden, abhängig von der relativen Expression der Proben A und B:
Die Menge der qPCR-Gene, die der korrekten Titrationsreihenfolge folgen, ist
Für einen einzelnen replizierten qPCR-Datensatz (z., dem PrimePCR-Datensatz, den wir analysiert haben), kann die inhärente Variabilität einer einzelnen qPCR-Messung zu einigen falsch negativen Genen führen, die dem korrekten TO folgen, aber nicht identifiziert werden können. Aus der Literatur50,51 ist der Variationskoeffizient für Replikat-qPCR-Messungen im Allgemeinen 15% oder größer, so dass wir diese Zahl verwendet haben, um den Spielraum für die Bestimmung anzupassen, ob ein Gen der richtigen RICHTUNG folgt. Mathematisch haben wir den Bereich von plus und minus einer Standardabweichung aus jeder qPCR-Messung berechnet und als Marge verwendet. Die überarbeiteten Gleichungen für \({K}_{TO}\) lauten wie folgt:
wobei \(a=1,15, b=0,85\)
Neben TO sollten Proben zusätzlich ein bestimmtes Mischungsverhältnis aufweisen. Da das Verhältnis zwischen den Proben A und B
ist, ist die EMR zwischen den Proben C und D
und bestimmt schließlich einen Satz von Genen, die das EMR-Kriterium wie folgt erfüllen:
RNA-seq—Datenanalyse-Pipelines – Kartierung, Quantifizierung und Normalisierung
Wir untersuchten 278 RNA-seq-Pipelines, die dreizehn Sequenzabbildungsalgorithmen enthielten 18,19,20,21,22,23,24,25,26,27,28,29, drei Kategorien von Ausdrucksquantifizierungsalgorithmen31,32,33 und sieben Ausdrucksnormalisierungsmethoden. Ergänzende Tabellen S2–S4 fassen alle Optionen zusammen, die für jede Pipelinekomponente in Betracht gezogen werden (Sequenzzuordnung, Expressionsquantifizierung und Expressionsnormalisierung). Die dreizehn untersuchten Mapping-Algorithmen sind Bowtie18, Bowtie219, BWA20, GSNAP21, Magic22 (eine neue Pipeline, die von NCBI für das SEQC-Projekt entwickelt wurde: ftp://ftp.ncbi.nlm.nih.gov/repository/acedb/Software/Magic ), MapSplice23, Novoalign (ein kommerzialisiertes Paket, das von Novocraft entwickelt wurde: https://www.novocraft.com/products/novoalign/), OSA24, RUM25, STAR26, Subread27, TopHat28 und WHAM29. Einige verwenden eine nicht gespleißte Zuordnung von Lesevorgängen zum Transkriptom, und einige andere führen eine gespleißte Zuordnung zum Genom durch. Magic verwendet beides parallel und vergleicht die Qualität jeder Ausrichtung, um die besten Ergebnisse für mehrere Ziele zu erzielen. Zuordnungsalgorithmen können nur eindeutige Zuordnungen melden oder mehrere Zuordnungsorte pro Lesevorgang zulassen. Quantifizierungsalgorithmen umfassen einfache zählbasierte Methoden (dh HTSeq31) und Poisson-verteilungsbasierte probabilistische Methoden, die entweder auf genomische (dh Cufflinks32) oder transkriptomische Kartierungsdaten (dh RSEM33) angewendet werden. Die Magie, RUM und Subread (dh., featureCounts52) Pipelines umfassen eingebettete Quantifizierungsmethoden, die in die Kategorie der einfachen zählbasierten Methoden fallen. Normalisierungsmethoden umfassen einfache Skalierungsmethoden (d. H. Fragmente pro Million abgebildeter Fragmente, Fragmente pro Kilobase Genlänge pro Million abgebildeter Fragmente, Median und oberes Quartil), robuste Skalierungsmethoden (d. H. Relative Log-Expression und getrimmter Mittelwert von m-Werten) und Methoden, die in spezifische Pipelines eingebettet sind (d. H. Magischer Expressionsindex).
Sequenzzuordnung
Wir haben Sequenzen in aufeinanderfolgenden Schritten mit un-Spliced- oder Spliced-Mapping-Algorithmen auf jede Referenz abgebildet. Un-spliced Mapping bezieht sich auf Algorithmen, die ganze Lesesequenzen ausrichten (z. B. Bowtie2, BWA und Novoalign), während sich Spliced Mapping auf Algorithmen bezieht, die Lesevorgänge in Segmente aufteilen, um lange Lücken oder Introns in einem Lesevorgang aufzunehmen (z. B. TopHat und MapSplice). Im ersten Schritt des un-gespleißten Mappings haben wir versucht, alle gepaarten Endsequenzen der ERCC / MT / rRNA-Referenz zuzuordnen (d. h., Externe RNA steuert Dna-Sequenzen, das mitochondriale Genom und ribosomale RNA-Sequenzen). Alle nicht zugeordneten Lesepaare wurden dann dem AceView-Transkriptom zugeordnet. Schließlich wurden alle Lesepaare, die weder den ERCC / MT / rRNA- noch den AceView-Referenzen zugeordnet waren, der menschlichen Genomreferenz zugeordnet. Transkriptomische Kartierungskoordinaten wurden dann in genomische Kartierungskoordinaten übersetzt und mit den Ergebnissen der humanen Genomkartierung zusammengeführt, um die Endergebnisse zu erhalten (Ergänzende Abb. S21, linke Seite). Wir haben Bowtie2 als Mapper für den ersten Schritt aller gespleißten Mapping-Pipelines verwendet (Ergänzende Abb. S21, rechte Seite). Gespleißte Abbildungsalgorithmen haben Lesevorgänge entweder direkt auf das menschliche Genom abgebildet (z. B. MapSplice und GSNAP) oder ganze nicht gespleißte Lesevorgänge auf das Transkriptom abgebildet und diese Abbildungsergebnisse dann mit gespleißten Abbildungsergebnissen der verbleibenden Lesevorgänge auf das menschliche Genom (z. B. TopHat und OSA) zusammengeführt. Ergänzende Tabelle S2 fasst alle in dieser Studie untersuchten Mapping-Tools zusammen.
Bowtie2, GSNAP, Novoalign, TopHat und WHAM ermöglichen die Kontrolle über die Anzahl der gemeldeten Zuordnungen pro Lesepaar. Standardmäßig melden diese Algorithmen normalerweise einen einzelnen besten Zuordnungsort pro Lesepaar. Einige Quantifizierungsalgorithmen können jedoch Informationen über mehrere mehrdeutige Abbildungsorte verwenden, um die Schätzung der Genexpression zu verbessern. So haben wir neben dem Single-Hit-Reporting auch Mapping-Ergebnisse generiert, die bis zu 200 Treffer pro Lesevorgang (Multi-Hit) auswiesen. Wir haben auch die Bowtie-Mapping-Pipeline mit Mapping-Parametern aufgenommen, die für die Quantifizierung mit RSEM spezifisch sind, wie im folgenden Abschnitt beschrieben33.
Die Befehlszeilenoptionen für alle Sequenzausrichtungswerkzeuge sind im ergänzenden Hinweis 1 beschrieben.
Quantifizierung der Genexpression
Die Quantifizierungsstufe umfasste drei Kategorien von Quantifizierern — zählbasierte Quantifizierer (d. h. HTSeq und integrierte Quantifizierer für die Pipelines Magic, RUM und Subread), wahrscheinlichkeitsmodellbasierte Quantifizierer für die genomische Kartierung (d. h., Manschettenknöpfe) und wahrscheinlichkeitsmodellbasierte Quantifizierer für transkriptomisches Mapping (d. h. RSEM). Die wichtigsten Merkmale dieser Quantifizierer sind in der ergänzenden Tabelle S3 zusammengefasst. Cufflinks ist ein auf dem Poisson-Modell basierender Quantifizierer, der die Wahrscheinlichkeiten der Lesezuweisung basierend auf den Ausrichtungsinformationen schätzt32. Es ist in der Lage, sowohl Transkripte zu assemblieren als auch Gen- oder Transkriptexpressionen zu quantifizieren. In dieser Studie haben wir die Assemblierungsfunktion deaktiviert und die GTF-Datei für Genomanmerkungen als Quantifizierungsreferenz bereitgestellt. HTSeq ist ein naïve Count-basierter Quantifizierer, der zugeordnete Lesevorgänge genes31 zuweist. HTSeq ist in der Lage, die Genexpression zu quantifizieren, aber nicht die Transkriptexpression. RSEM ist auch ein Poisson-modellbasierter Quantifizierer, dessen Konzept Cufflinks33 ähnelt. Informationen aus Multi-Hit-Lesevorgängen sind sowohl für Manschettenknöpfe als auch für RSEM wichtig. Diese Algorithmen verwenden Multi-Hit-Leseinformationen, um die Gen- oder Transkriptexpression genauer abzuschätzen.
Mapping-Ergebnisse von Alignment-Pipelines waren nicht immer mit den drei Kategorien von Quantifizierern kompatibel. Manschettenknöpfe erfordern, dass die Ausrichtungsergebnisse nach Ausrichtungskoordinaten sortiert sind und Multi-Hit-Lesevorgänge mit dem Tag ‚NH‘ im Attributfeld der SAM-Datei markiert sind. HTSeq erfordert, dass die Ausrichtungsergebnisse nach Lesenamen sortiert sind und dass das Tag ‚NH‘ in der SAM-Datei fehlt. RSEM quantifiziert nur transkriptomisches Mapping, das heißt, liest abgebildet und berichtet in transkriptomischen Koordinaten. Darüber hinaus behandelt RSEM nur un-gapped Alignments. Daher ist eine Filterung erforderlich, um gapped Alignments zu entfernen. Aufgrund dieser Anforderungen haben wir alle Ausrichtungsergebnisse vor der Quantifizierung vorverarbeitet. Zusammenfassend waren zwanzig Ausrichtungspipelines, einschließlich gespleißter, nicht gespleißter, Single-Hit- und Multi-Hit-Pipelines, für die zählbasierte Quantifizierung geeignet. Sechzehn Ausrichtungsrohre waren für Manschettenknöpfe geeignet, und nur zehn waren für RSEM geeignet. RSEM wurde speziell entwickelt, um gut mit Bowtie zu arbeiten. Daher haben wir auch diese eingebettete Mapping- und Quantifizierungspipeline einbezogen.
Die Befehlszeilenoptionen für alle Quantifizierungswerkzeuge sind im ergänzenden Hinweis 1 beschrieben.
Normalisierung der Genexpression
Die Normalisierung der RNA-seq-Daten ermöglicht einen Vergleich zwischen den Proben. Im Allgemeinen korrigieren Normalisierungsmethoden die Bibliotheksgröße (d. H. Die Gesamtzahl der Lesevorgänge in einer Stichprobe), die die Hauptquelle für die Varianz zwischen den Stichproben ist. Wir untersuchten sieben Normalisierungsmethoden — Fragmente pro Million kartierter Fragmente (FPM), Fragmente pro Kilobase Genlänge pro Million kartierter Fragmente (FPKM), Median (Med.), oberes Quartil (UQ), relativer logarithmischer Ausdruck (RLE), getrimmter Mittelwert der M-Werte (TMM) und Ausdrucksindex (EIndex, der für die Magic-Pipeline spezifisch ist) (siehe ergänzende Tabelle S4). Wir beschreiben jede dieser Normalisierungsmethoden basierend auf der folgenden mathematischen Beschreibung des SEQC-Benchmark-Datensatzes.
Wir definierten die Menge der gegenwärtigen Gene als
und die endgültige gegenwärtige Genmenge ist
Wir haben den gleichen Satz von present gens für alle Normalisierungsmethoden für eine RNA-seq-Pipeline verwendet.
Die Gesamtzahl der vorhandenen Gene für eine gegebene Probe s und replizieren n ist
und die durchschnittliche Gesamtzahl der vorhandenen Gene für alle Daten von einer einzelnen Site ist
Daher haben wir die FPM-normalisierte Expression für jede Probe s, jedes Replikat n und jedes Gen k als
Die Median- und oberquartilnormalisierte Expression für jede Probe s, jedes Replikat n und jedes Gen k wird dann definiert als
Für die FPKM-Normalisierung haben wir die Länge eines Gens k als \(\ell_{k}\) definiert, was die Länge der Vereinigung aller Exons ist, die mit dem Gen verwandt sind, wie durch das AceView-Transkriptom definiert. Die ursprüngliche Formulierung von FPKM verwendete willkürlich Skalierungsfaktoren von 1 × 103 für die Genlänge und 1 × 106 für die Gesamtzahl der abgebildeten Fragmente. Um einen vergleichbaren Dynamikumfang unter allen Normalisierungsmethoden zu erhalten, skalierten wir stattdessen mit der durchschnittlichen Genlänge und der durchschnittlichen Gesamtzahl für alle vorhandenen Gene. Die durchschnittliche Länge aller vorhandenen Gene ist
Somit ist die neu skalierte FPKM-normalisierte Expression für jede Probe s, replizieren n und Gen k
Die TMM- und RLE-Normalisierungsmethoden ähneln der FPM-Normalisierung, führen jedoch einen zusätzlichen Skalierungsfaktor ein, um die Bibliotheksgröße anzupassen. Wir haben das edgeR-Paket in R verwendet, um einen Skalierungsfaktor für jedes Sample replicate36,53 zu schätzen. Die TMM-Methode wählt eine Referenzbibliothek aus einem Pool von Probenreplikatbibliotheken aus und berechnet dann genweise Log-Expressionsverhältnisse (M-Werte) und genweise durchschnittliche Log-Expressionswerte (A-Werte) zwischen der Zielbibliothek und der Referenzbibliothek. Extreme Zahlen in M-Werten und A-Werten werden getrimmt, und der Skalierungsfaktor für die Zielbibliothek ist der gewichtete Durchschnitt der verbleibenden M-Werte. Die RLE-Methode bestimmt einen Skalierungsfaktor, indem zuerst die Medianbibliothek als genweises geometrisches Mittel über Probenreplikate definiert35. Als Skalierungsfaktor wird das Medianverhältnis jeder Zielbibliothek zur Medianbibliothek herangezogen. TMM- und RLE-normalisierte Expression für jede Probe s, replizieren n und Gen k werden dann definiert als:
wobei \(\hat{f}_{s,n}^{TMM}\) und \(\hat{f}_{s,n}^{RLE}\) der Skalierungsfaktor für Probe s, replizieren n sind.
RNA-seq-Pipeline-Leistungsmetriken
Benchmark-Metriken für RNA-seq-Pipelines sind in der ergänzenden Tabelle S7 zusammengefasst.
Genauigkeit gemessen als Abweichung von qPCR-Referenzen
Bei den Proben A und B ist die absolute Log-Ratio-Abweichung der RNA-seq-basierten Expression von der gen k ist
und die endgültige Genauigkeitsmetrik wurde als Median aller \({\Delta }_{{\frac{A}{B},k}}\), \(k = 1 \ldots K\) definiert.
Präzision gemessen als Variation der Genexpression über Replikatbibliotheken hinweg
Wir berechneten den Variationskoeffizienten (CoV) für jedes Gen und jede Probe über vier Replikatbibliotheken wie folgt:
Zuverlässigkeit gemessen als Intra-Sample-Korrelation der Genexpression
Die Zuverlässigkeit eines Messsystems kann anhand des intraclass correlation coefficient (ICC)54,55 beurteilt werden. ICC ist auf Messungen anwendbar, die in Gruppen organisiert werden können, und beschreibt, wie ähnlich Messungen derselben Gruppe zueinander sind. Die moderne ICC-Definition entlehnt den Rahmen der Varianzanalyse (ANOVA), genauer gesagt ANOVA mit zufälligen Effekten55. Die Art der ANOVA hängt vom experimentellen Design ab und folgt im Allgemeinen der Definition in Shrouts Artikel aus dem Jahr 197955. ICC (1,1) und ICC (1,k) basieren auf dem Einweg-Zufallseffektmodell und sind auf den Fall anwendbar, dass jede Gruppe von einem anderen Satz von k-Ratern bewertet wird, die zufällig aus einer größeren Population von Ratern ausgewählt wurden. ICC (2,1) und ICC (2,k) basieren auf dem Zwei-Wege-Zufallseffektmodell und sind auf den Fall anwendbar, dass eine Zufallsstichprobe von k Bewertern aus einer größeren Population vorausgewählt wird und jeder Bewerter jede Gruppe genau einmal bewertet (d. h., jeder Rater bewertet insgesamt n Gruppen). ICC (3,1) und ICC (3,k) basieren auf dem Zwei-Wege-Mischeffektmodell und sind auf den Fall anwendbar, dass jede Gruppe von jeweils denselben k-Bewertern bewertet wird, die die einzigen Bewerter in der Bevölkerung sind. Der zweite Parameter in ICC (,) gibt an, ob der ICC die Zuverlässigkeit einer einzelnen Messung oder den Durchschnitt von k Messungen messen soll.
Für den SEQC-Benchmark-Datensatz mit Replikatbibliotheken für jede Probe entsprach ICC (1,1) oder ICC(1,k) unserem Ziel, da für ein Gen g die Genexpression von Replikatbibliotheken für verschiedene Proben (oder verschiedene Gruppen im vorherigen Kontext) nicht unter genau den gleichen Bedingungen bewertet wurden (oder von denselben Bewertern im vorherigen Kontext bewertet wurden). Wir haben uns für ICC (1,k) entschieden, da für die meisten Experimente Replikatbibliotheken verfügbar sind. Mathematisch kann ein Einweg-Zufallseffektmodell wie folgt formuliert werden:
Wir haben ICC für jedes Gen k berechnet, \ (k = 1 \ ldots K\), und verwendete dann den Median aller ICCs als letztes Maß für die Zuverlässigkeit.
Wir haben auch andere potenzielle Metriken untersucht, wie z. B. die Reproduzierbarkeit, die als Spearman-Korrelation zwischen zwei Replikatbibliotheken derselben Probe definiert ist (Ergänzende Anmerkung 2). Die Spearman-Korrelation lag im Bereich von 0.993 bis 0,996 (Ergänzende Fig. S8) unter Verwendung von AllGenen. Wir haben die Reproduzierbarkeitsmetrik wegen des relativ kleinen Dynamikbereichs verworfen.
Bewertung des Nutzens der Benchmark-Metriken für die Auswahl der RNA-Seq-Pipeline
Wir haben die RNA-seq-Pipelines basierend auf dem durchschnittlichen Rang der drei Benchmark-Metriken (d. h. Genauigkeit, Präzision und Zuverlässigkeit) eingestuft. Anschließend bewerteten wir den Nutzen der Benchmark-Metriken, indem wir untersuchten, ob die auf der Grundlage der Benchmark-Metriken identifizierten Pipelines mit guter und schlechter Leistung aussagekräftig waren, um die Leistung der auf Genexpression basierenden Vorhersage des Krankheitsausganges und der statistischen Signifikanz der Patienten abzuleiten Schichtung für alle klinischen Endpunkte (d. H. die SEQC-Neuroblastom-EFS- und OS-Endpunkte und den TCGA-Lungenadenokarzinom-Überlebensendpunkt).Zunächst berechneten wir für die 278 repräsentativen RNA-seq-Pipelines, die auf den SEQC-Benchmark-Datensatz angewendet wurden, den durchschnittlichen Rang unter Verwendung einer Teilmenge der Benchmark-Metriken als endgültigen Leistungsindikator für jede Pipeline. Insgesamt hatten wir 6 Metriken (3 Benchmark—Metriken × 2 Gensätze) und untersuchten 12 Teilmengen (4 × 3) der 6 Metriken anhand der folgenden Kriterien:
- (1)
Vier Kombinationen der drei Benchmark-Metriken mit mindestens zwei in einer Teilmenge – eine Kombination mit allen drei Benchmark-Metriken, drei Kombinationen mit zwei der drei Benchmark-Metriken.
- (2)
Drei Teilmengen, die durch Metriken gebildet werden, die von allen Genen abgeleitet sind, die von niedrig exprimierenden Genen abgeleitet sind, oder eine Kombination aus beiden.
Zweitens berechneten wir für jede der 278 repräsentativen RNA-seq-Pipelines (156 für den TCGA-Lungenadenokarzinom-Überlebensendpunkt) verschachtelte Kreuzvalidierungs-AUC und MCC, wie im Abschnitt „Methode“ „Vorhersagemodellierung von Neuroblastom und Lungenadenokarzinom“ beschrieben, was zu 834 (468 für den TCGA-Lungenadenokarzinom-Überlebensendpunkt) AUC- und MCC-Werten für jeden klinischen Endpunkt (d. h., 278 Rohrleitungen × 3 Sichter oder 156 Rohrleitungen × 3 Sichter) (Ergänzende Tabellen S11, S12). Wir modellierten auch Überlebensfunktionen mithilfe der Kaplan-Meier-Analyse für jede Pipeline, wie im Abschnitt „Methode“ „Kaplan–Meier-Überlebensanalyse“ beschrieben. Für jede RNA-seq-Pipeline haben wir die Leistung der genexpressionsbasierten Vorhersage des Krankheitsausganges unter Verwendung der durchschnittlichen AUC und MCC über Klassifikatoren und der Erfolgsrate der Patientenstratifizierung (d. h., statistisch signifikante Trennung zweier Kaplan-Meier-Kurven) über alle Iterationen und Klassifikatoren im verschachtelten Kreuzvalidierungsframework.
Schließlich haben wir die oberen 10% der Pipelines mit guter Leistung und die unteren 10% der Pipelines mit schlechter Leistung basierend auf dem durchschnittlichen Rang einer Teilmenge der drei Benchmark-Metriken identifiziert. Die entsprechende Vorhersageleistung (d.h., AUC und MCC) der Pipelines mit guter Leistung wurde mit dem einseitigen Wilcoxon-Rangsummentest gegen den der Pipelines mit schlechter Leistung getestet mit der Nullhypothese, dass der Median der ersteren Gruppe nicht größer war als der der letzteren Gruppe.
Prädiktive Modellierung von Neuroblastom und Lungenadenokarzinom
Wir haben die Leistung von 278 RNA-seqc in Bezug auf genexpressionsbasierte Entscheidungen anhand des SEQC-Neuroblastom-Datensets48 bewertet. Der SEQC-Neuroblastom-Datensatz und die damit verbundenen klinischen Endpunkte sind in der ergänzenden Tabelle S9 zusammengefasst. Die RNA-seq-Pipelines wurden hinsichtlich der Vorhersage der Neuroblastom-Patientenergebnisse für zwei klinische Endpunkte unter Verwendung einer verschachtelten Kreuzvalidierung bewertet (Ergänzende Abb. S13)56,57. In ähnlicher Weise bewerteten wir auch die Leistung von 156 RNA-seq-Pipelines, die auf den TCGA-Lungenadenokarzinom-Datensatz angewendet wurden, um das Krankheitsergebnis vorherzusagen. Der TCGA-Lungenadenokarzinom-Datensatz und der zugehörige klinische Endpunkt sind in der ergänzenden Tabelle S10 zusammengefasst.
Die verschachtelte Kreuzvalidierung beinhaltet das Training und Testen eines optimalen Vorhersagemodells. Dies wird unter Verwendung der dreifachen Optimierung oder inneren Kreuzvalidierung erreicht, die von der fünffachen äußeren Kreuzvalidierung auf die Trainingsteilmenge angewendet wird. Sobald die endgültigen optimalen Vorhersagemodellparameter (d. h. Die Klassifikator-Hyperparameter und die Merkmalsgröße) identifiziert sind, wird das endgültige Modell unter Verwendung der gesamten Trainingsteilmenge trainiert und dann unter Verwendung der verbleibenden Falte aus der fünffachen äußeren Kreuzvalidierung getestet. Dieser Vorgang wurde für zehn Iterationen wiederholt. Wir haben die verschachtelte Kreuzvalidierung für jeden der drei Klassifikatoren separat durchgeführt (d. H., adaptive Boosting, logistische Regression und Support Vector Machines) und verwendete die mRMR-Merkmalsauswahlmethode (Minimum Redundancy, Maximum Relevance), um optimale Merkmalsgrößen aus dem Bereich von 5 bis 40 mit der Schrittgröße von 558 auszuwählen.
Kaplan–Meier-Überlebensanalyse
Für jede RNA-seq-Pipeline und jeden Klassifikator (d. h. 278 Pipelines × 3-Klassifikatoren für die SEQC–Neuroblastom-Endpunkte und 156 Pipelines × 3-Klassifikatoren für den TCGA-Lungenadenokarzinom-Überlebensendpunkt) modellierten wir Kaplan-Meier-Überlebensfunktionen basierend auf den vorhergesagten Markierungen jeder Probe. Wir verwendeten dann den Two-Tailed Log-Rank-Test, um festzustellen, ob die geschätzten Überlebenskurven für jede vorhergesagte Patientengruppe statistisch unterschiedlich waren.
Analyse der Varianz und Berechnung des Beitrags jedes RNA-seq-Pipeline-Faktors zur Gesamt-Pipeline-Varianz
Wir haben die Varianzanalyse (ANOVA) verwendet, um festzustellen, ob jeder RNA-seq-Pipeline-Faktor signifikant zur Varianz jeder der drei Benchmark-Metriken (d. h. Genauigkeit, Präzision und Zuverlässigkeit) sowie zur Varianz der Vorhersageleistung (d. h. AUC und MCC) beiträgt. Für jede der drei Benchmark-Metriken verwendeten wir ein lineares Modell (R-Funktion „lm“), um die Daten aus allen 278 Pipelines unter Verwendung der Metrik als abhängige Variable und der RNA-seq-Pipelinefaktoren als unabhängige kategoriale Variablen anzupassen. Wir haben die folgenden Faktoren als unabhängige kategoriale Variablen betrachtet – Mapping-Algorithmus, Mapping-Strategie (d. H. gespleißt vs. nicht gespleißt), Mapping-Reporting (d. H. Single-Hit vs. Multi-Hit), Quantifizierungsalgorithmus und Normalisierungsalgorithmus. Wir haben alle Faktoren und ihre wechselseitigen Wechselwirkungen in das lineare Modell einbezogen. Für jeden der Vorhersageendpunkte haben wir die gleiche Technik angewendet, um die Daten aus allen 278 Pipelines unter Verwendung der durchschnittlichen AUC oder MCC als abhängige Variable und des gleichen Satzes von RNA-seq-Pipeline-Faktoren als unabhängige kategoriale Variablen anzupassen. Wir haben dann die ANOVA auf dem linearen Modell durchgeführt (R-Funktion „Anova“). ANOVA berechnet eine „Summe der Quadrate“ (d. H. Varianz), die jedem Faktor oder jeder Interaktion zugeordnet wird, und verwendet einen F-Test, um festzustellen, ob die Varianz statistisch signifikant ist. Wir haben den Prozentsatz berechnet, den jeder Faktor oder jede Interaktion zur Gesamtvarianz beiträgt, indem wir das Verhältnis von „Summe der Quadrate“ für jeden Faktor zur Gesamtsumme der Quadrate berechnet haben.
Regressionsanalyse
Wir untersuchten die Beziehung zwischen Alignment-Profilen oder Genexpressionsverteilungsmerkmalen und Benchmark-Metriken. Die Ausrichtungsprofile umfassten die Gesamtzahl der abgebildeten Fragmente, die Gesamtzahl der Lesevorgänge, die den intronischen Bereich überspannen, die Gesamtzahl der Lesevorgänge mit Einfügungen oder Löschungen, die Gesamtzahl der perfekt übereinstimmenden Lesevorgänge, die Gesamtzahl der Lesevorgänge mit höchstens einer Fehlanpassung und die Anzahl der Fehlanpassungen pro abgebildetem Lesevorgang. Jeder Ausrichtungsalgorithmus wurde durch die durchschnittlichen Statistiken über 2 Sequenzierungsstellen, 4 Stichproben, 4 Replikatbibliotheken und 2 Spuren dargestellt. Unter Verwendung des „MASSEN“ -Pakets in R haben wir die M-Schätzung mit einem Gewichtungsansatz übernommen, um robuste lineare Regressionsmodelle zwischen einer abhängigen Variablen (Benchmark-Metrikleistung) und einer erklärenden Variablen (einem Ausrichtungsprofil) anzupassen. Die M-Schätzung mit Huber-Gewichtungsansatz ist eine Regressionsmethode, die in Gegenwart von Ausreißern robust ist. Die Merkmale der Genexpressionsverteilung umfassten das untere Quartil, den Median, das obere Quartil, das Maximum, den Interquartilbereich, die Standardabweichung, die Schiefe, die Kurtose und die Entropie einer Genexpressionsverteilung. Wir verwendeten dieselbe M-Schätzung mit einem Gewichtungsansatz, um robuste lineare Regressionsmodelle anzupassen, und berichteten dann über den verbleibenden Standardfehler für jedes Modell.
Haftungsausschluss
Die in diesem Artikel dargestellten Ansichten spiegeln nicht unbedingt die aktuelle oder zukünftige Meinung oder Politik der US-amerikanischen Food and Drug Administration wider. Jede Erwähnung kommerzieller Produkte dient der Klarstellung und ist nicht als Billigung gedacht.