inverkan av RNA-seq dataanalysalgoritmer på uppskattning av genuttryck och nedströms förutsägelse

FDA SEQC benchmark dataset

FDA SEQC-benchmark dataset (genuttryck omnibus anslutningsnummer GSE47792) inkluderar Parade RNA – seq-data som genereras med hjälp av Illumina HiSeq 2000-plattformen med läslängden på 100 nukleotider 7. Vi använde en delmängd av seqc-benchmark dataset sekvenserad på två platser-Beijing genomics Institute (BGI) och Mayo Clinic (maj). Vi använde fyra prover (dvs A, B, C och D), var och en med fyra replikatbibliotek förberedda på sekvenseringsplatserna. Prov A innehåller Universal Human Reference RNA( UHRR), prov B innehåller Human Brain Reference RNA (HBRR), prov C innehåller en blandning av A och B (75% A och 25% B) och prov D innehåller en blandning av A och B (25% A och 75% B). Vi använde data från två körfält i en enda flödescell för varje provreplikat. SEQC tillhandahöll också en kvantitativ PCR (qPCR) benchmark dataset som innehåller 20 801 gener analyserade med PrimePCR (Bio-Rad, Hercules, Kalifornien). Varje PrimePCR-gen analyserades en gång för vart och ett av de fyra proverna (dvs. A, B, C och D). FDA seqc benchmark dataset och prover sammanfattas i kompletterande tabeller S5 och S6.

neuroblastom och lung adenokarcinom dataset

vi använde en 176-prov neuroblastom dataset (en delmängd av en större 498-prov dataset; kallas SEQC-neuroblastom i detta manuskript) för att bedöma prestandan hos RNA-seq-rörledningar när det gäller genuttrycksbaserad förutsägelse av sjukdomsresultat. Dessa prover tillhandahölls av universitetets barnsjukhus i Köln och sekvenserades vid BGI med Illumina-plattformen48. Alla 176 prover togs från högriskpatienter som definierades som de antingen med neuroblastom i steg 4 och ålder > 18 månader eller med MYCN-förstärkta tumörer i något stadium eller ålder. SEQC-neuroblastoma dataset deponerades till genuttrycket Omnibus med anslutningsnummer GSE47792.

vi förutspådde två kliniska slutpunkter—händelsefri överlevnad (EFS), det vill säga förekomsten av händelser som framsteg, återfall eller död och total överlevnad (OS), det vill säga döden. För båda endpoints delades patienterna upp i två grupper (dvs. höga risker kontra låga risker). Högriskpatienter upplevde en händelse eller dog före en fördefinierad tröskelvärde för överlevnad, medan patienter med låg risk upplevde en händelse eller dog efter tröskeln, eller deras senaste uppföljning överskred tröskeln. Tröskelvärdena för överlevnadstid för EFS och OS var två respektive tre år. Tröskelvärdena valdes för att balansera antalet patienter med hög risk och låg risk. Uppgifter om SEQC-neuroblastom dataset finns i kompletterande tabell S9.

Vi använde också en 87-prov lung adenokarcinom RNA-seq dataset från Cancer Genome Atlas (TCGA) förvaret. Förutsägelsens slutpunkt var också överlevnad, och vi använde samma kriterier för att definiera högriskgrupper och lågriskgrupper med överlevnadstidsgränsen på två år. Den tvååriga tröskeln valdes för att balansera antalet patienter med hög risk och låg risk. Uppgifter om TCGA-lung-adenokarcinomdataset finns i kompletterande tabell S10.

filtrering av qPCR benchmark dataset för att producera en referensuppsättning gener

På grund av variabilitet i qPCR-mätningar och oenigheter mellan qPCR-plattformar7 filtrerade vi qPCR benchmark dataset för att behålla gener som uppvisade ”korrekt” beteende. Vi använde sedan dessa gener för att beräkna benchmark-mätvärdena (dvs. noggrannhet, precision, tillförlitlighet och reproducerbarhet). Sådan filtreringsprocessen sammanfattas i kompletterande Fig. S1.

Från och med den ursprungliga uppsättningen av 20,801-gener analyserade med PrimePCR filtrerade vi dessa gener för att behålla endast gener som kvantifierades som icke-noll (dvs detekterade) och med Ct-värden (cykeltröskel), 35 (35 indikerar detektering av endast en enda molekyl i ett prov). Filtrering av PrimePCR-data resulterade i 14 014 gener som också matchade AceView-transkriptomet som användes för att kartlägga SEQC-benchmark RNA-seq-dataset.

därefter filtrerade vi de 14 014 qPCR-generna för att behålla endast 12 610 gener som uppvisade rätt titreringsordning (till) och förväntade blandningsförhållanden (EMR). Detaljer om denna process finns i avsnittet” filtrering av qPCR-gener genom titreringsordning och förväntade blandningsförhållanden”.slutligen, eftersom vissa riktmärken som noggrannhet och precision är känsliga för noll – eller mycket låguttryckande gener, valde vi vidare gener som uttrycktes som icke-noll i alla replikater av alla prover av alla sekvenseringsställen och alla 278 RNA-seq-rörledningar. Den slutliga referensuppsättningen innehåller endast 10,222 qPCR-gener (kallad ”alla gener”) som användes för att beräkna alla tre benchmark-mätvärden för RNA-seq-rörledningar.

baserat på den tidigare studien är generna med lägre uttryck mer benägna att vara inkonsekventa bland rörledningar49. Således identifierade vi också en uppsättning låguttryckande gener i 10 222 gener baserat på det genomsnittliga qPCR-uttrycket av prover A, B, C och D. De lägsta 20% av de 10 222 generna (dvs. 2044 gener, kallade ”låguttryckande gener”) användes också för att beräkna samma uppsättning benchmark-mätvärden för RNA-seq-rörledningar. Denna design gjorde det möjligt för oss att undersöka förmågan hos RNA-seq-rörledningar vid uppskattning av låguttryckande genuttryck.

filtrering av qPCR-gener genom titreringsordning och förväntade blandningsförhållanden

SEQC-benchmark dataset (RNA-seq och qPCR) har unika egenskaper som möjliggör bedömning av kvantifieringskorrekthet. Efter att ha identifierat detekterbara (dvs. icke-noll-och Ct-kubi 35) och AceView-matchade qPCR-gener, använde vi två mätvärden (TO Och EMR) för att ytterligare filtrera riktmärket qPCR-dataset, vilket bara lämnar ”korrekta” qPCR-gener. Mätvärdena TO Och EMR fångar unika blandningsegenskaper för data, det vill säga

$$C= \frac{3}{4}A+\frac{1}{4}B\, \text{och }\,= \frac{1}{4}A+\frac{3}{4}B.$$

På grund av denna egenskap förväntas alla gener uttryckas i en av följande order, beroende på det relativa uttrycket av proverna A och B:

$$A\ge C\ge D\ge B\, \text{eller}\, a\le C\le D\le B.$$
$${\stackrel{-}{q}}_{s,\cdot, k}=\frac{1}{N}\sum_{n=1}^{n}{q}_{s ,n,k,}$$

uppsättningen qPCR-gener som följer rätt titreringsordning är

För en enda replikat qPCR-dataset (t. ex., PrimePCR-datasetet vi analyserade) kan den inneboende variationen i en enda qPCR-mätning resultera i vissa falska negativa gener som följer rätt till men misslyckas med att identifieras. Från litteraturen50, 51, variationskoefficienten för replikera qPCR-mätningar är i allmänhet 15% eller större, så vi använde detta nummer för att justera marginalen för att bestämma om en gen följer rätt till. Matematiskt beräknade vi intervallet plus och minus en standardavvikelse från varje qPCR-mätning och använde den som marginal. De reviderade ekvationerna för \({K}_{till}\) är följande:

$${K}_{TO}={K}_{TO,a\ge B}\cup {K}_{TO,a\le B,}$$

där \(a=1,15, b=0,85\)

förutom att, prover bör dessutom uppvisa ett specifikt blandningsförhållande. Med tanke på att förhållandet mellan proverna A och B är

$${R}_{A,B}=\frac{a}{b}$$

EMR mellan proverna C och D är

$$EM{R}_{C,D}=\frac{3z\cdot {R}_{A,B}+1}{z\cdot {R}_{A,B}+3}\cdot \frac{z+3}{3Z+1}$$
$${r}_{A,B}\in \left\equiv ,$$
$${R}_{C,D}\in \left\equiv \left,\text{ och}$$
$$em{r}_{c,d}\in \left\equiv ),$$

och bestämmer slutligen en uppsättning gener som uppfyller EMR-kriteriet enligt följande:

$${K}_{EMR}=\vänster\{k|\vänster({{R}_{C,D}^{nedre}\le {EMR}_{C,D}^{övre}|}_{{K, R}_{C,D}\ge EM{R}_{C,D}}\höger)\vee \vänster({{R}_{C,D}^{övre}\ge {EMR}_{C,D}^{lower}|}_{{k, r}_{c,d}\le em{r}_{c,d}}\right)\right\}$$

RNA-Seq dataanalys pipelines—kartläggning, kvantifiering och normalisering

vi undersökte 278 RNA-Seq-rörledningar som inkluderade tretton sekvenskartningsalgoritmer18,19,20,21,22,23,24,25,26,27,28,29, tre kategorier av uttryckskvantifieringsalgoritmer31, 32, 33 och sju uttrycksnormaliseringsmetoder. Tilläggstabeller S2-S4 sammanfattar alla alternativ som beaktas för varje rörledningskomponent (sekvenskartläggning, kvantifiering av uttryck och normalisering av uttryck). De tretton kartläggningsalgoritmerna som undersökts är Bowtie18, Bowtie219, BWA20, GSNAP21, Magic22 (en ny pipeline utvecklad av NCBI för SEQC-projektet: ftp://ftp.ncbi.nlm.nih.gov/repository/acedb/Software/Magic), MapSplice23, Novoalign (ett kommersialiserat paket utvecklat av Novocraft: https://www.novocraft.com/products/novoalign/), OSA24, RUM25, STAR26, Subread27, TopHat28 och WHAM29. Vissa använder oskärmad kartläggning av läsningar till transkriptomen, och några andra utför skarvad kartläggning till genomet. Magic använder både parallellt och jämför kvaliteten på varje anpassning för att hålla det bästa över flera mål. Kartläggningsalgoritmer kan bara rapportera unik kartläggning eller tillåta flera kartläggningsplatser per läsning. Kvantifieringsalgoritmer inkluderar enkla räknebaserade metoder (dvs. HTSeq31) och Poisson-distributionsbaserade probabilistiska metoder som tillämpas på antingen genomiska (dvs. Manschettlinks32) eller transkriptomiska kartdata (dvs. RSEM33). Magin, Rom och Subread (dvs., featureCounts52) rörledningar inkluderar inbäddade kvantifieringsmetoder som faller i kategorin enkla räknebaserade metoder. Normaliseringsmetoder inkluderar enkla skalningsmetoder (dvs. fragment per miljon mappade fragment , fragment per kilobas av genlängd per miljon mappade fragment , median och övre kvartil), robusta skalningsmetoder (dvs. relativt logguttryck och trimmat medelvärde av m-värden) och metoder inbäddade i specifika rörledningar (dvs. magiskt uttrycksindex).

Sequence mapping

Vi mappade sekvenser till varje referens i successiva steg med antingen oskärmade eller skarvade kartläggningsalgoritmer. Un-spliced mapping hänvisar till algoritmer som anpassar hela lässekvenser (t.ex. Bowtie2, BWA och Novoalign) medan spliced mapping hänvisar till algoritmer som delas läser in i segment för att rymma långa luckor eller introner i en läsning (t. ex. TopHat och MapSplice). I det första steget av oskärmad kartläggning försökte vi kartlägga alla parade sekvenser till ERCC/MT / rRNA-referensen (dvs., Extern RNA kontrollerar Konsortiesekvenser, mitokondriellt genom och ribosomala RNA-sekvenser). Alla omappade läspar mappades sedan till AceView-transkriptomen. Slutligen kartlades alla läspar som inte kartlades till antingen ERCC/MT/rRNA eller AceView-referenser till den mänskliga genomreferensen. Transkriptomiska kartläggningskoordinater översattes sedan till genomiska kartläggningskoordinater och slogs samman med de mänskliga genomkartningsresultaten för att producera de slutliga resultaten (kompletterande Fig. S21, vänster panel). Vi använde Bowtie2 som mapper för det första steget i alla skarvade kartläggningsledningar (kompletterande Fig. S21, höger panel). Skarvade kartläggningsalgoritmer antingen direkt mappade läsningar till det mänskliga genomet (t.ex. MapSplice och GSNAP) eller mappade hela icke-skarvade läsningar till transkriptomet och slog sedan samman dessa kartläggningsresultat med skarvade kartläggningsresultat av de återstående läsningarna till det mänskliga genomet (t. ex. TopHat och OSA). Kompletterande tabell S2 sammanfattar alla kartverktyg som undersökts i denna studie.

Bowtie2, Gsnap, Novoalign, TopHat och WHAM tillåter kontroll över antalet rapporterade mappningar per läspar. Som standard rapporterar dessa algoritmer vanligtvis en enda bästa kartläggningsplats per läspar. Vissa kvantifieringsalgoritmer kan dock använda information om flera tvetydiga kartläggningsplatser för att förbättra uppskattning av genuttryck. Förutom single-hit-rapportering genererade vi också kartläggningsresultat som rapporterade upp till 200 träffar per läsning (multi-hit). Vi inkluderade också Bowtie mapping pipeline med kartläggningsparametrar specifika för kvantifiering med RSEM, som beskrivs i följande avsnitt 33.

kommandoradsalternativ för alla sekvensinriktningsverktyg beskrivs i kompletterande anmärkning 1.

kvantifiering av genuttryck

kvantifieringssteget inkluderade tre kategorier av kvantifierare-räknebaserade kvantifierare (dvs. HTSeq och inbyggda kvantifierare för Magic, Rom och Subread pipelines), sannolikhetsmodellbaserade kvantifierare för genomisk kartläggning (dvs., Manschettknappar) och sannolikhetsmodellbaserade kvantifierare för transkriptomisk kartläggning (dvs. RSEM). De viktigaste egenskaperna hos dessa kvantifierare sammanfattas i kompletterande tabell S3. Manschettknappar är en Poisson-modellbaserad kvantifierare som uppskattar lästilldelningssannolikheter baserat på justeringsinformationen32. Det kan både montera transkript och kvantifiera gen-eller transkriptuttryck. I denna studie inaktiverade vi monteringsfunktionen och tillhandahöll GTF-filen genomanteckning som en kvantifieringsreferens. HTSeq är en na jacobve count-baserad kvantifierare som tilldelar mappade läser till genes31. HTSeq kan kvantifiera genuttryck, men inte transkriptuttryck. RSEM är också en Poisson modellbaserad kvantifierare som liknar konceptet Manschettlinks33. Information från multi-hit läser är viktigt för både manschettknappar och RSEM. Dessa algoritmer använder multi-hit läsa information för att mer exakt uppskatta gen eller transkript uttryck.

kartläggningsresultat från inriktningsledningar var inte alltid kompatibla med de tre kategorierna av kvantifierare. Manschettknappar kräver att justeringsresultat sorteras efter justeringskoordinater och multi-hit-läsningar markeras med ’NH’ – taggen i attributfältet i Sam-filen. HTSeq kräver att justeringsresultaten sorteras efter lästa namn och att’ NH ’ – taggen saknas i Sam-filen. RSEM kvantifierar endast transkriptomisk kartläggning, det vill säga läser mappade och rapporterade i transkriptomiska koordinater. Dessutom hanterar RSEM endast ospända anpassningar. Således krävs filtrering för att avlägsna gappade anpassningar. På grund av dessa krav förbehandlade vi alla anpassningsresultat före kvantifiering. Sammanfattningsvis var tjugo inriktningsrörledningar, inklusive skarvade, icke-skarvade, single-hit och multi-hit pipelines, lämpliga för räknebaserad kvantifiering. Sexton inriktningsledningar var lämpliga för manschettknappar, och endast tio var lämpliga för RSEM. RSEM är speciellt utformad för att fungera bra med Bowtie. Således inkluderade vi också denna inbäddade kartläggnings-och kvantifieringsrörledning.

kommandoradsalternativ för alla kvantifieringsverktyg beskrivs i kompletterande anmärkning 1.

genuttryck normalisering

RNA-seq data normalisering möjliggör jämförelse mellan prov. I allmänhet korrigerar normaliseringsmetoder bibliotekets storlek (dvs. det totala antalet läsningar i ett prov), vilket är den primära källan till varians mellan prov. Vi undersökte sju normaliseringsmetoder-fragment per miljon mappade fragment( FPM), fragment per kilobas av genlängd per miljon mappade fragment (FPKM), median (med.), övre kvartilen (UQ), relativt logguttryck (RLE), trimmat medelvärde för M-värden (TMM) och uttrycksindex (EIndex, som är specifikt för den magiska rörledningen) (se kompletterande tabell S4). Vi beskriver var och en av dessa normaliseringsmetoder baserat på följande matematiska beskrivning av SEQC-benchmark dataset.

$$\overline{x}_{S, \cdot ,k} = \frac{1}{n}\mathop \sum \limits_{n = 1}^{n} x_{s,n,k}$$

vi definierade uppsättningen nuvarande gener som

och den slutliga nuvarande genuppsättningen är

$$K_{p} = k_{p,BGI} \Cap K_{p,maj} .$$

vi använde samma uppsättning nuvarande gens för alla normaliseringsmetoder för en RNA-seq-pipeline.det totala antalet nuvarande gener för ett givet prov s och replikera n är

$$x_{s,n} = \mathop \sum \limits_{{k \In K_{p} }} x_{s,n,k}, $$

och det genomsnittliga totala antalet nuvarande gener för alla data från en enda plats är

$$\bar{x} = \frac{1}{4}\frac{1}{n}\mathop \sum \limits_{s} \mathop \sum \limits_{{n = 1}}^{n} x_{{s,n}}.$$

således definierade vi FPM-normaliserat uttryck för varje prov s, replikera n och gen k som

$$y_{s,n,k}^{FPM} = \frac{{x_{s,n,k} \cdot \overline{x}}}{{x_{s,n} }}.$$

Median – och övre kvartil-normaliserat uttryck för varje prov s, replikera n och gen k definieras sedan som

$$y_{s, n,k}^{med} = \frac{{x_{s,n,k} \cdot \tilde{x}}}{{\tilde{x}_{S,n} }}{\text{och }}y_{s,n,k}^{uq} = \frac{{x_{s,n,k} \cdot \hat{X}}}{{\hat{x}_{S,N} }}.{ } $$

För fpkm-normalisering definierade vi längden på en gen k som \(\ell_{k}\), vilket är längden på föreningen av alla exoner relaterade till genen enligt definitionen av AceView-transkriptomet. Den ursprungliga formuleringen av FPKM använde godtyckligt skalningsfaktorer på 1 kg 103 för genlängden och 1 kg 106 för det totala antalet mappade fragment. För att upprätthålla jämförbart dynamiskt intervall bland alla normaliseringsmetoder skalas vi istället av den genomsnittliga genlängden och det genomsnittliga totala antalet för alla nuvarande gener. Den genomsnittliga längden för alla närvarande gener är

$$\overline{\ell } = \frac{1}{{\left|{K_{p} } \right/}}\mathop \sum \limits_{{k \in K_{p} }} \ell_{k} .$$

således omskalas FPKM-normaliserat uttryck för varje prov s, replikera n och gen k är

$$y_{s, n, k}^{FPKM} = \frac{{x_{s,n,k} \cdot \overline{\ell } \cdot \overline{x}}}{{x_{S,N} \cdot \ell_{k} }}.$ $

TMM-och RLE-normaliseringsmetoderna liknar FPM-normaliseringen, men introducerar en ytterligare skalningsfaktor för att justera bibliotekets storlek. Vi använde edgerpaketet i R för att uppskatta en skalningsfaktor för varje provreplicate36,53. TMM-metoden väljer ett referensbibliotek från en pool av provreplikatbibliotek och beräknar sedan genvisa logguttrycksförhållanden (M-värden) och genvisa genomsnittliga logguttrycksvärden (a-värden) mellan målbiblioteket och referensbiblioteket. Extrema tal I M-värden och A-värden trimmas, och skalningsfaktorn för målbiblioteket är det vägda genomsnittet för återstående M-värden. RLE-metoden bestämmer en skalningsfaktor genom att först definiera medianbiblioteket som det genvisa geometriska medelvärdet över provrepliker35. Medianförhållandet för varje målbibliotek till medianbiblioteket tas som skalningsfaktor. TMM-och RLE-normaliserat uttryck för varje prov s, replikera n och gen k definieras sedan som:

där \(\hat {f} _ {s, n}^{TMM}\) och \(\hat{f}_{s,n}^{RLE}\) är skalningsfaktorn för prov s, replikera n.

RNA-seq pipeline performance metrics

Benchmark metrics för RNA-seq pipelines sammanfattas i kompletterande tabell S7.

noggrannhet mätt som avvikelse från qPCR-referenser

$${\stackrel {-} {y}}_{s,\cdot ,k}=\frac{1}{n}\sum_{n=1}^{n}{y}_{s,n,k}$$

givna prover A och B, den absoluta log-ratio avvikelsen för RNA-seq-baserat uttryck från qPCR-baserat uttryck för en gen k är

$$\delta_{\frac{a}{b},k} = \vänster / \log_2\vänster ( \frac{\bar{x}_{a,., k}} {\bar{x}_{B,., k}} \ höger)- \ log_2 \ vänster (\frac {\bar{q}_{A,., k}} {\bar{q}_{B,., k}} \right ) \right|, $$

och mätvärdet för slutlig noggrannhet definierades som medianen för alla \({\Delta }_{{\frac{A}{B},k}}\), \(k = 1 \ldots K\).

Precision mätt som variation av genuttryck över replikatbibliotek

Vi beräknade variationskoefficienten (CoV) för varje gen och varje prov över fyra replikatbibliotek enligt följande:

$$CoV_{s,k} = \frac{{sd\left( {x_{s, \cdot ,k} } \right)}}{{\overline{x}_{S, \cdot ,k} }},$$

tillförlitlighet mätt som intra-provkorrelation av genuttryck

tillförlitligheten hos ett mätsystem kan bedömas med intraclass korrelationskoefficient (ICC)54,55. ICC är tillämplig på mätningar som kan organiseras i grupper, och det beskriver hur liknande mätningar av samma grupp är varandra. Modern ICC-definition lånar ramverket för variansanalys (ANOVA), eller mer specifikt ANOVA med slumpmässiga effekter55. Typen av ANOVA beror på den experimentella designen och följer i allmänhet definitionen i Shrouts artikel publicerad 197955. ICC(1,1) och ICC (1,k) är baserade på envägs random effects-modellen och är tillämpliga på det fall att varje grupp bedöms av en annan uppsättning k-raters slumpmässigt utvalda från en större population av raters. ICC(2,1) och ICC(2,k) är baserade på tvåvägs random effects-modellen och är tillämpliga på det fall att ett slumpmässigt urval av k-bedömare väljs i förväg från en större population och varje bedömare bedömer varje grupp exakt en gång (dvs., varje rater bedömer n grupper helt och hållet). ICC (3,1) och ICC(3,k) är baserade på tvåvägsmodellen för blandade effekter och är tillämpliga på det fall att varje grupp bedöms av var och en av samma k-bedömare, som är de enda bedömarna i befolkningen. Den andra parametern i ICC(,) anger om ICC ska mäta tillförlitligheten hos en enda mätning eller genomsnittet av k-mätningar.

för seqc benchmark dataset med replikatbibliotek för varje prov,ICC(1,1) eller ICC (1, k) passade vårt mål eftersom, för en gen g, genuttryck av replikatbibliotek för olika prover (eller olika grupper i föregående sammanhang) inte bedömdes under exakt samma förhållanden (eller bedömdes av samma bedömare i föregående sammanhang). Vi valde att använda ICC (1, k) eftersom replikatbibliotek är tillgängliga för de flesta experiment. Matematiskt kan en enkelriktad slumpmässig effektmodell formuleras som

$$Y_{ij} = \mu + \alpha_{j} +{\varepsilon_{ij}}, $$
$$ICC\left( {1, k} \right) = \frac{BMS-WMS}{{BMS}},$$

Vi beräknade ICC för varje gen k, \(k = 1 \ldots k\) och använde sedan medianen för alla ICC: er som det slutliga måttet på tillförlitlighet.

Vi har också undersökt andra potentiella mätvärden, såsom Reproducerbarhet, som definieras som Spearman-korrelationen mellan två replikatbibliotek i samma prov (kompletterande anmärkning 2). Spearman-korrelationen varierade från 0.993 till 0,996 (kompletterande Fig. S8) använda AllGenes. Vi kasserade mätvärdet för Reproducerbarhet på grund av det relativt lilla dynamiska området.

utvärdera nyttan av benchmark metrics för RNA-Seq pipeline selection

vi rankade RNA-seq pipelines bas på den genomsnittliga rankningen av de tre benchmark metrics (dvs. noggrannhet, precision och tillförlitlighet). Vi utvärderade sedan nyttan av benchmark-mätvärdena genom att undersöka om bra och dåligt fungerande rörledningar identifierade baserat på benchmark-mätvärdena var informativa för att dra slutsatsen om prestanda för genuttrycksbaserad förutsägelse av sjukdomsresultat och statistisk signifikans av patientstratifiering för alla kliniska slutpunkter (dvs. seqc-neuroblastoma EFS och OS endpoints och TCGA-lung-adenokarcinomöverlevnad endpoint).

För det första, för de 278 representativa RNA-seq-rörledningarna som tillämpades på seqc-benchmark-datasetet, beräknade vi den genomsnittliga rankningen med hjälp av en delmängd av benchmark-mätvärdena som den slutliga prestandaindikatorn för varje rörledning. Totalt hade vi 6 metrics (3 benchmark metrics, 2 genuppsättningar), och vi undersökte 12 delmängder (4, 3) av de 6 metrics med hjälp av följande kriterier:

  1. (1)

    fyra kombinationer av de tre benchmark metrics med minst två i en delmängd—en kombination med alla tre benchmark metrics, tre kombinationer med två av de tre benchmark metrics.

  2. (2)

    tre delmängder bildade av mätvärden härledda från alla gener, de som härrör från låguttryckande gener eller en kombination av båda.

För det andra, för var och en av de 278 representativa RNA-seq-rörledningarna (156 för TCGA-lung-adenokarcinom survival endpoint), beräknade vi kapslade korsvalidering AUC och MCC, som beskrivs i avsnittet ”metod” ”neuroblastom och lung adenokarcinom prediktiv modellering”, vilket resulterade i 834 (468 för TCGA-lung-adenokarcinom survival endpoint) AUC och MCC-värden för varje klinisk endpoint (dvs., 278 rörledningar 3 klassificerare,eller 156 rörledningar 3 klassificerare) (Tilläggstabeller S11, S12). Vi modellerade också överlevnadsfunktioner med Kaplan-Meier-analys för varje rörledning, som beskrivs i avsnittet” metod ””Kaplan–Meier överlevnadsanalys”. För varje RNA-seq-pipeline sammanfattade vi prestanda för genuttrycksbaserad förutsägelse av sjukdomsresultat med både den genomsnittliga AUC och MCC över klassificerare och framgångsgraden för patientstratifiering (dvs., statistiskt signifikant separation av två Kaplan-Meier-kurvor) över alla iterationer och klassificerare i det kapslade korsvalideringsramen.

slutligen identifierade vi de bästa 10% goda rörledningarna och de nedre 10% dåliga rörledningarna baserat på den genomsnittliga rankningen för en delmängd av de tre riktmärkesmetoderna. Motsvarande prediktionsprestanda (dvs., AUC och MCC) för de godpresterande rörledningarna testades mot de Dålig utförande rörledningarna med hjälp av det ensidiga Wilcoxon rank-sum-testet med nollhypotesen att medianen för den tidigare gruppen inte var större än den senare gruppen.

neuroblastom och lungadenokarcinom prediktiv modellering

vi bedömde prestanda för 278 RNA-seq-rörledningar i termer av genuttrycksbaserat beslutsfattande med hjälp av seqc-neuroblastom dataset48. SEQC-neuroblastoma dataset och tillhörande kliniska endpoints sammanfattas i kompletterande tabell S9. RNA-seq-rörledningarna utvärderades i termer av att förutsäga neuroblastom-patientresultat för två kliniska slutpunkter med användning av kapslad korsvalidering (kompletterande Fig. S13) 56,57. Vi bedömde också på liknande sätt prestandan hos 156 RNA-seq-rörledningar applicerade på TCGA-lung-adenokarcinomdataset för att förutsäga sjukdomsresultat. Datasetet TCGA-lung-adenokarcinom och tillhörande klinisk slutpunkt sammanfattas i kompletterande tabell S10.

kapslad korsvalidering innebär utbildning och testning av en optimal prediktionsmodell. Detta uppnås med hjälp av trefaldig optimering eller inre korsvalidering, applicerad på träningsundermängden från den femfaldiga yttre korsvalideringen. När de slutliga optimala prediktionsmodellparametrarna (dvs. klassificeringshyperparametrarna och funktionsstorleken) identifieras tränas den slutliga modellen med hela träningsundermängden och testas sedan med den återstående vikten från den femfaldiga yttre korsvalideringen. Denna process upprepades för tio iterationer. Vi genomförde kapslad korsvalidering separat för var och en av de tre klassificerarna (dvs. logistisk regression och stödvektormaskiner) och använde minsta redundans, maximal relevans (mRMR) funktions urvalsmetod för att välja optimala funktionsstorlekar inom intervallet 5 till 40 med stegstorleken 558.

Kaplan-Meier survival analysis

för varje RNA-seq pipeline och klassificerare (dvs. 278 pipelines Microgaming 3 klassificerare för seqc-neuroblastoma endpoints och 156 pipelines Microgaming 3 klassificerare för TCGA-lung-adenokarcinom survival endpoint) modellerade vi Kaplan–Meier survival functions baserat på de förutspådda etiketterna för varje prov. Vi använde sedan det två-tailed log-rank-testet för att avgöra om uppskattade överlevnadskurvor för varje förutsagd patientgrupp var statistiskt olika.

analys av varians och beräkning av bidraget från varje RNA-seq-rörlednings faktor till den totala rörledningens varians

Vi använde analys av varians (ANOVA) för att avgöra om varje RNA-seq-rörlednings faktor väsentligt bidrar till variansen för var och en av de tre riktmärkena (dvs. noggrannhet, precision och tillförlitlighet) samt till variansen för prediktionsprestanda (dvs. AUC och MCC). För var och en av de tre benchmark-mätvärdena använde vi en linjär modell (r-funktionen ”lm”) för att passa data från alla 278 rörledningar med metriken som den beroende variabeln och RNA-seq-pipelinefaktorerna som oberoende kategoriska variabler. Vi betraktade följande faktorer som oberoende kategoriska variabler-kartläggningsalgoritm, kartläggningsstrategi (dvs. spliced vs un-spliced), kartläggningsrapportering (dvs. single-hit vs multi-hit), kvantifieringsalgoritm och normaliseringsalgoritm. Vi inkluderade alla faktorer och deras tvåvägsinteraktioner i den linjära modellen. För var och en av förutsägelsens slutpunkter tillämpade vi samma teknik för att passa data från alla 278 rörledningar med genomsnittlig AUC eller MCC som den beroende variabeln och samma uppsättning RNA-seq-rörlednings faktorer som oberoende kategoriska variabler. Vi genomförde sedan ANOVA på den linjära modellen (r-funktionen ”anova”). ANOVA beräknar en ”summan av kvadrater” (dvs varians) som tillskrivs varje faktor eller interaktion och använder ett F-test för att avgöra om variansen är statistiskt signifikant. Vi beräknade procenten att varje faktor eller interaktion bidrar till den totala variansen genom att beräkna förhållandet ”summan av kvadrater” för varje faktor till den totala summan av kvadrater.

regressionsanalys

vi undersökte förhållandet mellan inriktningsprofiler eller genuttrycksfördelningsegenskaper och benchmark-mätvärden. Inriktningsprofilerna inkluderade det totala antalet mappade fragment, det totala antalet läsningar som spänner över den introniska regionen, det totala antalet läsningar med Infogningar eller raderingar, det totala antalet perfekt matchade läsningar, det totala antalet läsningar med högst en felaktighet och antalet felaktigheter per mappad läsning. Varje justeringsalgoritm representerades av den genomsnittliga statistiken över 2 sekvenseringsplatser, 4 prover, 4 replikera bibliotek och 2 banor. Med hjälp av” MASS ” -paketet i R antog vi m-estimation with Huber weighting-metoden för att passa robusta linjära regressionsmodeller mellan en beroende variabel (benchmark metric performance) och en förklarande variabel (en justeringsprofil). M-estimering med Huber weighting approach är en regressionsmetod som är robust i närvaro av avvikare. Genuttrycksfördelningsegenskaperna inkluderade lägre kvartil, median, övre kvartil, maximal, interkvartilområde, standardavvikelse, skevhet, kurtos och entropi av en genuttrycksfördelning. Vi använde samma m-uppskattning med Huber viktningsmetod för att passa robusta linjära regressionsmodeller och rapporterade sedan det återstående standardfelet för varje modell.

Disclaimer

de synpunkter som presenteras i denna artikel återspeglar inte nödvändigtvis nuvarande eller framtida åsikt eller politik US Food and Drug Administration. Varje omnämnande av kommersiella produkter är för förtydligande och inte avsedd som ett godkännande.

Lämna ett svar

Din e-postadress kommer inte publiceras. Obligatoriska fält är märkta *