Informatie

Wat is de meest geschikte manier om genexpressiegegevens te normaliseren?


Deze vraag komt omdat het lezen van een paper over normalisatie van genexpressiegegevens niet duidelijk is of de methode voor het normaliseren van de gegevens alleen voor RNA-Seq-gegevens is of ook voor microarrays kan worden toegepast.

Voor RNA-Seq-gegevens zijn er normalisatiemethoden die corrigeren voor het GC-inhoudseffect of andere effecten op genniveau. Is het logisch om deze effecten te overwegen bij de normalisatie van microarray-genexpressiegegevens?


Normalisatie van expressiegegevens is een groot onderwerp en er worden regelmatig nieuwe methoden gepubliceerd. Als je zoiets als dit benadert, wil je over het algemeen mensen bekijken die soortgelijke dingen hebben gedaan als jij, en als je eenmaal begrijpt waarom ze deden wat ze deden, kun je vragen wat je moet doen om je vragen te beantwoorden. Houd altijd je biologische vraag in gedachten. Als u bijvoorbeeld QTL's meet, moet u veel voorzichtiger zijn dan wanneer u alleen op zoek bent naar genen die zijn beïnvloed door een knock-out-mutatie.

Over het algemeen wil je heel verschillende methoden gebruiken voor RNAseq- en Microarray-gegevens. De twee gegevenstypen volgen totaal verschillende distributies (RNAseq geeft u telgegevens, microarraygegevens geven u continue signalen) en hebben verschillende soorten technische ruis die ze beïnvloeden (GC-inhoud heeft invloed op beide, maar op een andere manier). Sommige methoden kunnen op beide worden gebruikt, maar omvatten meestal het dwingen van de gegevens in een andere vorm (bijv. mapping telt tot een normale verdeling). Het limmapakket voor R kan beide aan, met verschillende distributies, en is een goed begin. Er bestaan ​​nieuwere, zogenaamd betere methoden voor RNAseq, die ik niet persoonlijk heb gebruikt.


Over het algemeen wilt u voor RNA-seq-gegevens niet corrigeren voor GC-inhoud of andere effecten op genniveau (bijv. lengte) omdat u expressiewaarden vergelijkt tussen aandoeningen BINNEN een gen. Om deze reden wordt aanbevolen om onbewerkte tellingen te gebruiken en niet genormaliseerde waarden zoals FPKM. Zie paragraaf 2.7 van de edgeR-gebruikershandleiding.

Deze recente benchmark die RNA-seq-kwantificatiemethoden vergelijkt, is misschien de moeite van het bekijken waard.


Trainingsmodules

De volgende stap in de RNA-seq-workflow is de differentiële expressieanalyse. Het doel van differentiële expressietesten is om te bepalen welke genen op verschillende niveaus tussen omstandigheden tot expressie worden gebracht. Deze genen kunnen biologisch inzicht bieden in de processen die worden beïnvloed door de aandoening(en) van belang.

De stappen die in het grijze vak hieronder worden beschreven, hebben we al besproken, en we zullen nu doorgaan met het beschrijven van de stappen in een end-to-end workflow voor RNA-seq differentiële expressie op genniveau.

Dus wat vertegenwoordigen de telgegevens eigenlijk? De telgegevens die worden gebruikt voor differentiële expressieanalyse vertegenwoordigen het aantal sequentielezingen dat afkomstig is van een bepaald gen. Hoe hoger het aantal tellingen, hoe meer reads geassocieerd zijn met dat gen, en de veronderstelling dat er een hoger expressieniveau van dat gen in het monster was.

De analysestappen voor differentiële expressie worden in het onderstaande stroomschema groen weergegeven. Ten eerste moeten de telgegevens worden genormaliseerd om rekening te houden met verschillen in bibliotheekgroottes en RNA-samenstelling tussen monsters. Vervolgens zullen we de genormaliseerde tellingen gebruiken om enkele plots voor QC op gen- en monsterniveau te maken. Ten slotte wordt de differentiële expressieanalyse uitgevoerd met behulp van uw interessetool.


Normalisatiemethoden voor miRNA-kwantificering – Ask TaqMan #40

Wat betekenen endogene controles, exogene controles en gemiddelde expressiewaarde voor u?

Je hebt het geraden! Vandaag beantwoorden we uw vragen over normalisatiemethoden voor miRNA-kwantificering.

MicroRNA's of miRNA's zijn klein

Niet-coderende RNA's van 22 nucleotiden die genexpressie reguleren. microRNA's kunnen worden gekwantificeerd door real-time PCR met behulp van TaqMan-assays. In een microRNA-expressie-experiment kunnen variatie in de hoeveelheid uitgangsmateriaal, monsterverzameling, RNA-voorbereiding en -kwaliteit en efficiëntie van omgekeerde transcriptie bijdragen aan kwantificeringsfouten. Om deze redenen is het belangrijk om de juiste normalisatiecontroles te gebruiken bij het kwantificeren van miRNA's.

Er zijn drie soorten normalisatiemethoden die gewoonlijk worden gebruikt voor miRNA-analyse door qPCR - endogene controles, exogene controles en gemiddelde expressiewaarde-normalisatie, of "globale gemiddelde normalisatie"

Normalisatie met behulp van endogene controlegenen is momenteel de meest nauwkeurige methode om mogelijke verschillen in RNA-invoer of RT-efficiëntiebias te corrigeren. Exogene controles of "spike-ins" worden meestal gebruikt om de extractie-efficiëntie of de hoeveelheid monsterinvoer te bewaken voor moeilijke monsters zoals plasma/serum of andere biovloeistoffen. Grootschalige miRNA-expressieprofileringsstudies kunnen gebruik maken van globale gemiddelde normalisatie, waarbij het berekende gemiddelde van alle miRNA's in een bepaald monster als normalisatiemiddel wordt gebruikt.

Historisch gezien werden niet-coderende RNA's zoals snRNA's en snoRNA's gebruikt als endogene normalisatoren voor miRNA-kwantificering.

Maar meer recentelijk zijn belangrijke opinieleiders in de miRNA-gemeenschap om de volgende redenen afgestapt van het gebruik van snoRNA's/snRNA's als endogene controles:

  • Ze zijn groter dan miRNA's
  • Ze 'spiegelen' niet de fysiochemische eigenschappen van miRNA
  • Ze hebben een andere cellulaire verwerking en verschillende functies dan miRNA's
  • De expressieniveaus van snoRNA en snRNA zijn recentelijk geassocieerd met kanker en prognose.

De miRNA-gemeenschap heeft ook gesuggereerd dat de ideale endogene controle genexpressie heeft die relatief constant en matig overvloedig is in een verscheidenheid aan weefsels en celtypen en behandeling.

miRNA's die uniform tot expressie worden gebracht, kunnen worden gebruikt als een endogene controle. Er zijn verschillende miRNA's waarvan is aangetoond dat ze in de literatuur en in experimentele onderzoeken tot expressie worden gebracht op relatief constante niveaus in veel verschillende weefseltypes (zie tabel). Deze kunnen werken als goede endogene controles voor uw experimentele conditie.

Trouwens, het wordt aanbevolen om 2 of meer van deze miRNA's te valideren als endogene controles voor de doelcel, het weefsel of de behandeling die u gebruikt, omdat geen enkele controle kan fungeren als een universele endogene controle voor alle experimentele omstandigheden.

Bovendien kunnen synthetische miRNA-moleculen worden gebruikt als spike-in-controles en zijn ze uiterst nuttig als exogene controles in moeilijke monsters zoals serum/plasma.

Zoals de naam al aangeeft, zijn spike-ins of exogene controles synthetische RNA-oligonucleotiden die aan het monster worden toegevoegd.

Een spike-in-controle moet een doelsequentie zijn die niet in uw monster aanwezig is. Ath-miR-159a (uitspraakbegeleiding nodig) is bijvoorbeeld niet aanwezig bij mensen, dus het is een goede exogene controle voor mensen.

Er zijn 3 verschillende normalisatiemethoden, exogene controles, endogene controles en gemiddelde expressiewaardenormalisatie en met deze normalisatiemethoden kunt u bepaalde aspecten van uw experimentele proces controleren bij het analyseren van miRNA's door qPCR.


Resultaten

Quantro: test voor globale verschillen in verdelingen tussen groepen

Overweeg een set onbewerkte gegevens met hoge doorvoer x ik vertegenwoordigen l ∈ (1, …, N k) monsters in elk van de k ∈ (1, …, K) groepen (N t totale monsters) van een genexpressie- of DNA-methyleringsexperiment. Wij nemen aan x ik heeft een gemeenschappelijke distributie (x ikk) waar k is de theoretische verdeling voor de k e groep. We definiëren ( _^ <-1>) als de waargenomen kwantielverdeling voor de l e monster in de k e groep. Als eerste stap gebruiken we een ANOVA om te testen of het gemiddelde van de medianen van de verdelingen tussen groepen verschilt en de mediaan om de steekproeven dienovereenkomstig te normaliseren. Laat ( >_<.k>^<-1>=frac<1>_^_^<-1>> ) de kwantielverdeling zijn, gemiddeld over alle steekproeven in de k e groep en laat ( >_<..>^<-1>=frac<1>frac<1>_^K_^_^<-1>>> ) de kwantielverdeling zijn, gemiddeld over alle steekproeven en groepen.

Om de verschillen tussen twee verdelingen te kwantificeren, gebruiken we de afstand van Mallow [28], die wordt gedefinieerd als de afstand tussen twee kansverdelingen over een regio (vergelijking S1 in aanvullend bestand 1). We definiëren de totale variantie van de verdelingen als de som van gekwadrateerde verschillen tussen F ik − 1 en ( >_<..>^ <-1>) met behulp van Mallow's afstand (in het geval waarin P = 2) als:

De totale variantie kan worden ontleed (Vgl. S2–7 in aanvullend bestand 1) in de variantie tussen groepen (SS tussen) en de variantie binnen groepen (SS binnenin):

We stellen voor om een ​​gegevensgestuurde teststatistiek te gebruiken, ook wel: F quantro, om te testen op globale verschillen in de verdelingen tussen de K groepen. De nulhypothese is dat er geen globale verschillen zijn in de verdelingen tussen de groepen en de alternatieve hypothese is dat ten minste één groep verschilt van de rest.

Als er geen globale verschillen zijn in de verdelingen tussen de groepen (vanwege technische of biologische variatie), kunnen we een globale aanpassingsmethode toepassen, zoals kwantielnormalisatie, om ongewenste technische variatie te verwijderen. Als er globale verschillen zijn in de verdelingen tussen de groepen, is kwantielnormalisatie mogelijk geen geschikte normalisatietechniek, afhankelijk van de bron van variatie (technische of biologische variatie).

De F quantro teststatistiek (Vgl. S8 in aanvullend bestand 1) is een verhouding van de gemiddelde kwadratische fout tussen groepen (MEVROUW tussen) naar de gemiddelde kwadratische fout binnen groepen (MEVROUW binnenin):

We gebruiken permutatietests om de statistische significantie van F quantro en verwerp de nulhypothese als de P waarde (Vgl. S9 in aanvullend bestand 1) van de permutatietest is minder dan een α significantieniveau.

Gerichte en globale veranderingen in genexpressie

wij hebben gesolliciteerd quantro naar verschillende openbaar beschikbare gegevenssets voor genexpressie op basis van zowel microarray- als RNA-Seq-platforms (tabel S1 in aanvullend bestand 1) om te onderzoeken gericht en globaal verschillen in verdeling over groepen. We gebruikten een α = 0,05 significantieniveau als de drempel om te testen op globaal veranderingen in de verdeling over groepen. Voorbeelden van gerichte veranderingen in verdelingen over groepen zijn de genexpressie van monsters uit de Yoruba (YRI)-populatie, gestratificeerd naar genotype op basis van een expression kwantitatieve trait loci (eQTL) (P = 0,917 Fig. 2a Figuur S1 in aanvullend bestand 1), monsters van twee inteeltmuisstammen (P = 0,245 Figuur S2 in aanvullend bestand 1), monsters van alveolaire macrofagen van niet-rokers, rokers en patiënten met astma (P = 0,562 Fig. 2b Figuur S3 in aanvullend bestand 1), monsters van bronchiale poetsbeurten van personen met en zonder chronische obstructieve longziekte (P = 0,218 Figuur S4 in aanvullend bestand 1) en monsters uit twee hersengebieden bij patiënten met de ziekte van Parkinson (P = 0.264 Afbeelding S5 in aanvullend bestand 1). In alle bovenstaande voorbeelden wordt kwantielnormalisatie geschikt geacht omdat er geen globale verschillen in de verdelingen tussen groepen werden gedetecteerd op het significantieniveau α = 0,05.

Bij het vergelijken van de genexpressie van twee weefsels, vonden we opvallende globale verschillen in de verdelingen tussen hersen- en leverweefsels (P = 0,004 Afb. 2c Afbeelding S6 in Aanvullend bestand 1). We hebben meerdere onderzoeken van de Gene Expression Omnibus (GEO) overwogen om elk weefsel weer te geven om te voorkomen dat batcheffecten [29] van verschillende onderzoeken van GEO worden verward met verschillen in weefsels. We vergeleken ook de genexpressie van normale en tumormonsters. We verkregen meerdere onderzoeken van GEO en vonden globale verschillen in de verdelingen tussen de normale en tumormonsters van long (P < 0,001 Afb. 2d), borst (P < 0,001), prostaat (P < 0,001), schildklier (P < 0,001), maag (P < 0,001) en leverweefsels (P = 0,044) (figuren S7-12 in aanvullend bestand 1). We vonden ook globale veranderingen in de verdeling van leverweefsel tussen vier groepen patiënten (controle, gezonde obesitas, steatosis en nash-monsters) uit een onderzoek naar de genexpressie van niet-alcoholische leververvetting (P = 0,004 Afbeelding S13 in Aanvullend bestand 1).

Gerichte en globale veranderingen in DNA-methylatie

Naast genexpressie hebben we drie openbaar beschikbare DNA-methyleringsgegevenssets overwogen. We ontdekten geen globale verschillen in de verdeling van vetweefsel van patiënten voor en na zes maanden inspanning (P = 0.132 Fig. 3a Figuur S14 in Aanvullend bestand 1) en pancreasweefsels van niet-diabetische en type 2 diabetes (P = 0,069 Afbeelding S15 in aanvullend bestand 1). In tegenstelling tot, quantro ontdekte globale verschillen in de verdelingen over zes gezuiverde celtypes uit volbloed (P < 0.001 Fig. 3b Figuur S16 in aanvullend bestand 1), wat relevant kan zijn voor de studies die de celsamenstelling van volbloed schatten met behulp van DNA-methylatie [30, 31].

Biologische variatie in distributies van onbewerkte DNA-methylatie-microarrays. een Voorbeeld van gerichte veranderingen in distributies: onbewerkte bètawaarden van N = 46 arrays waarin vetweefselmonsters van gezonde mannen voor en na 6 maanden inspanning worden vergeleken. B Voorbeeld van globale veranderingen in distributies: onbewerkte bètawaarden van N = 35 arrays waarin zes gezuiverde celtypen uit volbloed worden vergeleken: CD14+-monocyten (Mono), CD19+ B-cellen (Bcell), CD4+ T-cellen (CD4T), CD56+ natuurlijke killercellen (NK), CD8+ T-cellen (CD8T), en granulocyten (oma)

Quantro verbetert de nauwkeurigheid van het detecteren van differentieel gemethyleerde CpG's

Hier evalueren we de prestaties van globale normalisatiemethoden in de context van: gericht en globaal veranderingen in distributies met als doel het detecteren van differentieel gemethyleerde CpG's. We hebben een Monte Carlo-simulatiestudie uitgevoerd om te illustreren hoe het gebruik van globale normalisatiemethoden, zoals kwantielnormalisatie, niet altijd geschikt is en de F quantro teststatistieken kunnen de keuze voor normalisatie leiden. Voor de simulatiestudie simuleren we DNA-methylatie-arrays met als doel differentieel gemethyleerde CpG's te detecteren, maar merk op dat deze resultaten ook vertalen voor differentiële genexpressie. We vergelijken naïef het gebruik van kwantielnormalisatie met het gebruik van quantro om de beslissing om ofwel kwantielnormalisatie ofwel geen normalisatie te gebruiken te begeleiden om de kosten te beoordelen van het gebruik van globale normalisatiemethoden in de context van distributies met globale verschillen.

Als er slechts een minderheid van differentieel gemethyleerde CpG's is, vermindert kwantielnormalisatie de vertekening en de gemiddelde kwadratische fout (MSE) bij het detecteren van echte verschillen tussen groepen monsters omdat het ongewenste technische variatie verwijdert (figuren S21 en S22 in aanvullend bestand 1). Naarmate het aantal differentieel gemethyleerde CpG's toeneemt, zal kwantielnormalisatie zowel de ongewenste technische als interessante biologische variatie verwijderen, wat resulteert in een hogere bias en MSE bij het detecteren van differentiële methylatie. Daarentegen is het gebruik van quantro detecteert deze globale verschillen en vermindert daarom de vertekening en MSE in vergelijking met het gebruik van kwantielnormalisatie (figuren S21 en S22 in aanvullend bestand 1). Evenzo wordt het aantal valse ontdekkingen verminderd bij het gebruik van quantro om de normalisatiekeuze te begeleiden in het geval dat er globale verschillen zijn tussen groepen. Als u bijvoorbeeld een 450K DNA-methylatie-array overweegt als er slechts een klein aantal differentieel gemethyleerde CpG's is (1% van CpG's of 4500 CpG's), quantro en kwantielnormalisatie vergelijkbaar zijn in het aantal valse ontdekkingen (respectievelijk 873 en 873), maar als er globale verschillen zijn in de verdelingen tussen groepen (10% van CpG's of 45.000 CpG's), quantro is in staat om die globale verschillen te detecteren en het aantal valse ontdekkingen te verminderen in vergelijking met kwantielnormalisatie (respectievelijk 4887 en 6583) (Figuur S23 in aanvullend bestand 1). Gebruik makend van quantro geeft onderzoekers een datagestuurd hulpmiddel om te testen of globale normalisatiemethoden geschikt zijn, zoals kwantielnormalisatie, wat kan resulteren in grotere vertekening, MSE en meer valse ontdekkingen bij het detecteren van differentieel gemethyleerde CpG's in de context van globale verschillen in distributies.

Bovendien hebben we rekening gehouden met de werkelijke positieve en valse positieve snelheid van het gebruik van kwantielnormalisatie en het gebruik van quantro om de keuze van normalisatie te begeleiden, terwijl de drempel van het aantal geselecteerde top differentieel gemethyleerde CpG's wordt gevarieerd. Als er slechts een klein aantal differentieel gemethyleerde CpG's zijn, kwantielnormalisatie en quantro zijn vergelijkbaar in prestaties, maar wanneer het aandeel differentieel gemethyleerde CpG's toeneemt, kan kwantielnormalisatie geen globale verschillen tussen de groepen detecteren, wat resulteert in een lagere gevoeligheid en specificiteit (Figuur S24 in aanvullend bestand 1). Gebruik makend van quantro als een hulpmiddel om te bepalen welk type normalisatiebenadering moet worden gebruikt, resulteert in een hogere gevoeligheid en specificiteit bij het detecteren van echte differentieel gemethyleerde CpG's in vergelijking met naïef gebruik van kwantielnormalisatie.


De impact van normalisatiemethoden op RNA-Seq-gegevensanalyse

High-throughput sequencing-technologieën, zoals de Illumina Hi-seq, zijn krachtige nieuwe hulpmiddelen voor het onderzoeken van een breed scala aan biologische en medische problemen. Enorme en complexe datasets die door de sequencers worden geproduceerd, creëren een behoefte aan de ontwikkeling van statistische en computationele methoden die de analyse en het beheer van gegevens kunnen aanpakken. De gegevensnormalisatie is een van de meest cruciale stappen van gegevensverwerking en dit proces moet zorgvuldig worden overwogen omdat het een diepgaand effect heeft op de resultaten van de analyse. In dit werk richten we ons op een uitgebreide vergelijking van vijf normalisatiemethoden met betrekking tot sequencing-diepte, veel gebruikt voor transcriptome sequencing (RNA-seq) gegevens, en hun impact op de resultaten van genexpressie-analyse. Op basis van deze studie stellen we een universele workflow voor die kan worden toegepast voor de selectie van de optimale normalisatieprocedure voor een bepaalde dataset. De beschreven workflow omvat de berekening van de bias- en variantiewaarden voor de controlegenen, de gevoeligheid en specificiteit van de methoden en classificatiefouten, evenals het genereren van de diagnostische plots. Het combineren van bovenstaande informatie vergemakkelijkt de selectie van de meest geschikte normalisatiemethode voor de bestudeerde datasets en bepaalt welke methoden onderling uitwisselbaar kunnen worden gebruikt.

1. Inleiding

In de afgelopen jaren is RNA-sequencing (RNA-seq)-technologie een steunpilaar geworden van biomedisch onderzoek en een aantrekkelijk alternatief voor microarrays. RNA-seq-technologieën hebben verschillende voordelen ten opzichte van microarrays, waaronder minder ruis (de technische vooroordelen die inherent zijn aan microarray-technologie zijn niet aanwezig in RNA-seq-experimenten) [1], de mogelijkheid om alternatieve splicing-isovormen te detecteren [2, 3], en de kracht om nieuwe genen, genpromotors, isovormen en allelspecifieke expressie te detecteren [4]. De dalende kosten van next-generation sequencing (NGS) is een extra argument voor de keuze voor RNA-seq in plaats van op microarray gebaseerde genexpressie-analyse. Ondanks de lagere bias in het RNA-seq-experiment, zijn er nog enkele bronnen van systematische variatie die moeten worden geëlimineerd uit RNA-seq-gegevens vóór de differentiële expressie (DE) -analyse. Deze variaties omvatten met name verschillen tussen monsters, zoals bibliotheekgrootte (sequentiëringsdiepte) [5] of verschillen binnen monsters, bijvoorbeeld in genlengte [6], guanine-cytosine (GC)-gehalte [7, 8], of ongewenste variatie geïntroduceerd door batch-effect [9]. Ervaring met microarray-gegevens heeft herhaaldelijk aangetoond dat normalisatie ervoor moet zorgen dat expressieschattingen beter vergelijkbaar zijn tussen kenmerken (genen) en monsters. Er zijn echter nog veel vragen over de impact van de normalisatiemethode op de resultaten van RNA-seq data-analyse. Het belang van normalisatie van RNA-seq-gegevens werd aangetoond in [10]. Hun belangrijkste bevinding was dat de keuze van de normalisatieprocedure de resultaten van DE-analyse beïnvloedt: de gevoeligheid varieert meer tussen normalisatieprocedures dan tussen teststatistieken. Auteurs van [10, 11] toonden aan dat de normalisatiestap vragen oproept, en dat er nog steeds behoefte is aan het geven van bruikbare praktische tips en het opstellen van duidelijke richtlijnen voor onderzoekers die niet zeker weten hoe ze een normalisatiemethode moeten kiezen. Hiertoe stellen we voor om een ​​combinatie van grafische en statistische methoden toe te passen om de impact van de specifieke normalisatie op de resultaten van DE-analyse te vergelijken. Ons onderzoek heeft betrekking op vijf normalisatiemethoden die veel worden gebruikt voor de normalisatie van RNA-seq-gegevens: Trimmed Mean of

-waarden, bovenste kwartiel, mediaan, kwantiel en PoissonSeq-normalisatie geïmplementeerd in respectievelijk R-pakketten edgeR (v3.2.4), DESeq (v1.12.1), EBSeq (v1.3.1) en PoissonSeq 3 (v1.1.2). De vergelijking was gebaseerd op de analyse van drie datasets. Twee daarvan zijn openbaar beschikbare datasets (Bodymap- en Cheung-gegevens) en één is de dataset (AML-gegevens) afkomstig van een van onze projecten (nog niet gepubliceerd). In dit artikel schetsen we een eenvoudige en effectieve methode voor het vergelijken van verschillende normalisatiebenaderingen en laten we zien hoe onderzoekers de resultaten van differentiële expressieanalyse kunnen verbeteren door in de normalisatiestap verschillende aspecten op te nemen op basis van de biologie of informatica.

2. materialen en methoden

In deze sectie beschrijven we de te vergelijken normalisatiemethoden en de datasets die in ons onderzoek zijn gebruikt. Vervolgens stellen we de criteria voor voor het vergelijken van de impact van de normalisatiemethoden op de resultaten van DE-analyse.

2.1. Normalisatiemethoden

Sinds de opkomst van de RNA-seq-technologie zijn er een aantal normalisatiemethoden ontwikkeld. In ons werk richten we ons voornamelijk op een vergelijking van vijf van de meest populaire normalisatiemethoden die worden gebruikt voor DE-analyse van RNA-seq-gegevens, geïmplementeerd in vier Bioconductor-pakketten: Trimmed Mean of -values ​​(TMM) [11] en Upper Quartile (UQ) [10], beide geïmplementeerd in het edgeR Bioconductor-pakket [12], Median (DES) geïmplementeerd in het DESeq Bioconductor-pakket [13], Quantile (EBS) [10] geïmplementeerd in het EBSeq Bioconductor-pakket [14] en PoissonSeq (PS ) normalisatie geïmplementeerd in het PoissonSeq-pakket [15]. Alle pakketten zijn verkrijgbaar bij CRAN (http://cran.r-project.org/web/packages) en Bioconductor (http://www.bioconductor.org/packages/release/bioc).

Omdat de basisbron van variaties tussen monsters het verschil in bibliotheekgrootte is (RNA-monsters kunnen naar verschillende diepten worden gesequenced), wordt aan elke bibliotheek een normalisatiefactor toegewezen. Er zijn verschillende manieren om een ​​normalisatiefactor te berekenen. We beschouwen vijf verschillende methoden die hieronder worden beschreven.

Laten we aannemen dat we hebben

duiden leestellingen aan voor gen

is de schaalfactor voor het e monster.

De eerste gepresenteerde methode berekent een getrimd gemiddelde van -waarden tussen elk paar monsters. Deze methode werd voorgesteld door Robinson en Oshlack [11] en is gebaseerd op de hypothese dat de meeste genen niet differentieel tot expressie worden gebracht. De auteurs definieerden een normalisatiefactor voor een bestudeerd monster met een referentiemonster als volgt:

. zijn, respectievelijk, het totale aantal uitlezingen voor monster th en referentiemonster , en

vertegenwoordigt de set genen met niet bijgesneden en

(absolute expressieniveaus) waarden (volgens [11] is het getrimde gewogen gemiddelde het gemiddelde na het verwijderen van het bovenste en onderste percentage van de gegevenswaarden worden met 30% getrimd en waarden worden met 5% getrimd). Volgens de veronderstelling dat de meeste genen geen DE zijn,

De volgende methode, uitgewerkt in [10], is Upper Quartile (UQ) geïmplementeerd in het edgeR-pakket. Hier wordt de schaalfactor berekend vanaf het 75e percentiel van de tellingen voor elke bibliotheek na het verwijderen van transcripten, die in alle bibliotheken nul zijn. Het heeft de volgende vorm:

is het bovenste kwartiel van monster , dat is het monster met genormaliseerde tellingen en .

Anders en Huber [13] suggereerden de Median-normalisatiemethode die is geïmplementeerd in het DESeq Bioconductor-pakket. Deze methode gaat uit van dezelfde aanname als de TMM-methode (de meeste genen zijn geen DE). Een schaalfactor voor een gegeven monster neemt de mediaan van de verhoudingen van de waargenomen tellingen van het monster tot het geometrische gemiddelde over monsters (d.w.z. een pseudoreferentiemonster):

Het is ook mogelijk om kwantitatieve normalisatie uit te voeren over monsters, zoals vaak wordt gedaan in het geval van microarray-gegevens [10]. Hier hebben we Quantile-normalisatie gebruikt die is geïmplementeerd in het EBSeq Bioconductor-pakket [14]. Deze normalisatiemethode schat de sequentiediepte van een experiment met een bovenste kwartiel van zijn tellingen en heeft de volgende vorm:

zijn de bovenste kwartielen van respectievelijk het e monster en het e monster. Deze methode lijnt de verdeling van alle monsters uit.

Ten slotte hebben we ook een normalisatiemethode (PS) getest die is voorgesteld in [15], geïmplementeerd in het PoissonSeq-pakket. Het heeft de volgende formule:

waar is de set genen gevonden door goodness-of-fit-statistieken van de vorm:

is de Total Count (TC)-schaalfactor. We kozen genen voor de set waarvoor de waarden van

Verder hebben we de prestaties van alle hierboven genoemde normalisatiemethoden vergeleken met niet-genormaliseerde gegevens, aangeduid met "onbewerkte gegevens" (RD).

Daarnaast werd ook gekeken naar een andere bron van variatie, gerelateerd aan GC-gehalte en batch-effect. Verkregen resultaten onthulden dat opname van GC-inhoudsanalyse niet bijdroeg aan de voorgestelde vergelijkingsstrategie en geen invloed had op de samenvattingsrangschikking (zie figuur S5 en tabel S4 in aanvullend materiaal online beschikbaar op //dx.doi.org/10.1155/ 2015/621690). Daarom is deze normalisatiestrategie niet meegenomen in verdere analyse.

De ongewenste variatie die voortkomt uit batcheffecten, zoals bemonsteringstijd, andere technologie, kan worden aangepast met behulp van methoden die zijn geïmplementeerd in sva-pakket [9]. We hadden alleen informatie over de batches in de AML-gegevens die op twee bemonsteringsdatums konden worden geïntroduceerd. Daarom hebben we een dergelijke aanpak alleen voor AML-gegevens overwogen. Alvorens aanvullende normalisatie toe te voegen, hebben we gecontroleerd of de aanwezigheid van batcheffecten in AML-datasets bestond. Het detecteren van de aanwezigheid van batcheffecten in AML-gegevens werd op twee manieren bereikt. We hebben een hiërarchische clustering gegeven (zie figuur S8) samen met hoofdcomponentenanalyse (zie figuur S9). Onze resultaten bevestigden niet het bestaan ​​van een batcheffect in verband met de bemonsteringsdatum in het geval van AML-gegevens. Daarom werd deze normalisatiestrategie in verdere analyse verwaarloosd.

2.2. Data bronnen

De vijf normalisatiemethoden werden vergeleken op basis van de analyse van drie echte RNA-seq-gegevenssets. Twee datasets die in deze studie zijn opgenomen, zijn verkregen uit openbaar beschikbare bronnen (Bodymap- en Cheung-gegevens) en één is verkregen uit het project dat is uitgevoerd aan het Institute of Bioorganic Chemistry in Poznan (AML-gegevens). De Bodymap-gegevensset werd gepubliceerd in [16], waar het transcriptoom van niet-getransformeerde menselijke borstepitheelcellen met verwijzing naar Illumina Bodymap-gegevens verzameld uit normale weefsels werd bestudeerd. De Cheung-dataset [17] komt uit de studie van genexpressie van menselijke B-cellen van individuen die tot grote families behoren. Het doel van de studie was de identificatie van polymorfe transregulatoren. Beide datasets zijn gedeponeerd in de Recount-database op http://bowtie-bio.sourceforge.net/recount/ [18].

De AML-dataset is afkomstig van ons nog niet gepubliceerde RNA-seq-experiment. We bestudeerden genexpressieprofielen in 30 perifere bloed (PB) en beenmerg (BM) monsters verkregen van 25 volwassen AML (acute myeloïde leukemie) patiënten die genezen waren aan de afdeling Hematologie en Beenmergtransplantatie, Poznan University of Medical Sciences. Twee monsters per baan werden gesequenced bij het Institute of Bioorganic Chemistry in Poznan, met behulp van single-read flow-cell (SR-FC) en Genome Analyzer IIx (GAIIx, Illumina). Als controle werden één BM-monster en een pool van 12 PB-monsters verkregen van gezonde vrijwilligers gesequenced op één extra baan. De bibliotheken werden voorbereid tot 4 μg totaal RNA met behulp van de TruSeq RNA Sample Preparation Kit v2 (Illumina), volgens de instructies van de fabrikant, en gevalideerd met een DNA 1000-chip (Agilent). Elke bibliotheek genereerde ongeveer 20 miljoen 50-nt-lange reads, verwerkt in CASAVA, FastQC en de NGS QC Toolkit. De uitlezingen werden toegewezen aan het menselijke referentiegenoom UCSC hg19 met TopHat-run (v2.0.6).

Vóór alle berekeningen hebben we de datasets gefilterd om een ​​vergelijkbare hoeveelheid genen in elke dataset te verkrijgen. Uit de Cheung- en Bodymap-gegevenssets kozen we alleen die genen waarvoor het gemiddelde van de tellingen voor alle monsters groter was dan 0, terwijl we in de AML-gegevensset 50 kozen als de grenswaarde voor het gemiddelde van de tellingen in alle monsters. De samengevatte informatie over elke gegevensset wordt weergegeven in tabel 1, terwijl de details te vinden zijn in aanvullende tabel S1 en figuur S1. Zoals we kunnen zien, varieerden de datasets met de steekproefomvang en genaantallen, evenals de genexpressieniveaus. In de Cheung-dataset overheersen genen met lage expressieniveaus. In de Bodymap-dataset wordt de grootste groep gevormd door de genen met een hoog expressieniveau, terwijl in de AML-dataset de grootste groep wordt gevormd door de genen met een gemiddeld expressieniveau. Een dergelijke verscheidenheid aan datasets maakt het mogelijk om verschillen in prestaties van normalisatiemethoden te onthullen in het geval van gegevens met een verschillende structuur van genen.

2.3. Analytische criteria voor vergelijking van normalisatiemethoden
2.3.1. Selectie van de Housekeeping-genen (HG)

Het idee om HG te gebruiken is voortgekomen uit het eerdere onderzoek naar het microarray-experiment. In dat experiment werden dye-swapped microarrays gebruikt. Gekozen kleurstoffen speelden de rol van huishoudgenen. In weloverwogen RNA-seq-experimenten hebben we niet direct de lijsten met huishoudgenen. Daarom hebben we besloten om ons onderzoek naar de impact van normalisatiemethoden uit te voeren op basis van de analytische versie van HG-lijsten, genen die op dezelfde manier tot expressie worden gebracht in monsters.

De vierkantswortel van de gemiddelde vierkantsfout werd als maat gebruikt. De huishoudgenen werden geselecteerd op basis van de onbewerkte gegevens met behulp van de volgende formule:

waar is het gemiddelde van de tellingen voor gen. We besloten om voor alle datasets huishoudgenen te kiezen op basis van een bepaalde afkapwaarde en pasten een lineaire transformatie van de resultaten toe op het interval

met behulp van min-max normalisatie:

Vervolgens selecteerden we 1% van alle genen met de laagste waarden als huishoudgenen. Hoewel het aantal huishoudgenen in elke dataset hetzelfde was, werden voor elke set verschillende genen opgenomen. De tabellen en staafdiagrammen met betrekking tot de selectie van huishoudgenen zijn te vinden in aanvullende materialen (tabel S2 en figuur S2). Deze duiden op dezelfde feiten als voor het geval waarbij alle genen in aanmerking worden genomen. Er is enige variabiliteit in datasets met betrekking tot het aantal HG-genen voor elk niveau van uitlezingen. In alle datasets kunnen we zien dat er het hoogste aantal HG-genen is met een aantal reads in genen onder 500 (relatief gemiddeld abundantieniveau).

2.3.2. Bias en variantie

Om de normalisatiemethoden te evalueren die worden gebruikt voor de verwerking van RNA-seq-gegevens, hebben we het bias- en variantiecriterium toegepast dat is voorgesteld door Argyropoulos et al. [19] voor de analyse van dubbelkanaals microarray-gegevens. We hebben deze methode eerder aangepast voor eenkanaals microarray-gegevens [20]. In dit artikel transformeren we de methode om geschikt te zijn voor de RNA-seq-gegevens. We kunnen de volgende formules van vooringenomenheid en variantie gebruiken als respectievelijke proxy's voor nauwkeurigheid en precisie:

waar staat voor lezen telt voor

het controlegen van het e monster is de gemiddelde telling voor het e controlegen en is de gemiddelde waarde van voor het e controlegen, en True Log Ratio is de werkelijke waarde van voor het e controlegen uit de definitie van HG. De nauwkeurigheid van een meting is de mate waarin de metingen van een grootheid de werkelijke waarde benaderen. De nauwkeurigheid van een meting is de mate waarin herhaalde metingen onder ongewijzigde omstandigheden dezelfde resultaten laten zien. Op basis van de definitie moet elk gen, dat als HG wordt beschouwd, hetzelfde aantal tellingen in elk monster hebben, evenals de gemiddelde waarde van het aantal tellingen voor dit gen. De True Log Ratio is dus gelijk aan 0 en de formule (10) voor de bias reduceert tot de root mean square error (RMSE): de ratio's (bias en variantie) voor elke normalisatiemethode zijn het gemiddelde van de bias- en variantiewaarden berekend voor alle controle (huishoudelijke) genen. De normalisatiemethode "A" zou de voorkeur hebben boven de methode "B" als deze gepaard gaat met de kleinste vertekening en variantie [19].

2.3.3. Differentiële expressieanalyse

Ons doel was om erachter te komen hoe normalisatiemethoden de resultaten van differentiële expressie beïnvloeden. Daarom werd na toepassing van elke normalisatiemethode differentiële expressieanalyse uitgevoerd met behulp van de edgeR-methode uit het edgeR Bioconductor-pakket [12]. Deze methode is gekozen omdat het betrouwbare prestaties liet zien in een breed scala aan experimenten en het gemakkelijk is om schaalfactoren in de statistische test op te nemen. De edgeR-methode is speciaal ontwikkeld om de spreiding van tellinggegevens te modelleren en is ontworpen voor oververspreide RNA-seq-gegevens. In het kort wordt aangenomen dat tellingen een negatieve binomiale verdeling zijn, volgens

waar is het expressieniveau van het gen in het monster, is de normalisatiefactor in het monster en is de dispersie van het gen. De methode van differentiële expressie-analyse, geïmplementeerd in het edgeR-pakket, breidt de exacte test van Fisher uit.

2.3.4. Voorspellingsfouten

Discriminerende analyse werd toegepast om de effectiviteit van de monsterclassificatie te bepalen op basis van gevonden DEG's die in elke dataset na elke normalisatieprocedure werden geïdentificeerd. Verschillende classificatiemethoden kunnen leiden tot verschillende voorspellingsfouten. We schatten de classificatiefouten op basis van de vijf classificaties: naïeve Bayes, neuraal netwerk,

-nearest buur, en ondersteunen vectormachines en willekeurig bos die worden weergegeven in Tabel 2. Bij het analyseren van elke dataset, hebben we leave-one-out kruisvalidatie (LOOCV) gebruikt om de schatting te verkrijgen van de voorspellerfouten van de testset die resulteert uit verschillende classificaties gebruiken. Voor een dataset met

monsters, deze methode omvat afzonderlijke runs. Voor elke run wordt een aantal steekproeven, minus één gegevenspunt, gebruikt om het model te trainen en vervolgens wordt een voorspelling uitgevoerd op het resterende gegevenspunt. De totale voorspellingsfout is de som van de fouten voor alle runs [21]. Als input voor de classifiers in de discriminantanalyse kozen we deze informatieve genen, met een hoog onderscheidingsvermogen die differentieel tot expressie worden gebracht. De set genen werd verkregen door een genselectieproces uit een teststatistiek na elke normalisatieprocedure.

Alle berekeningen en diagnostische plots werden uitgevoerd met R 3.0.2 [22].

3. Resultaten

Het doel van onze studie was om verschillende normalisatiebenaderingen te vergelijken en de workflow te schetsen die zal helpen bij het selecteren van de normalisatiemethode die geschikt is voor een bepaalde dataset. We hebben besloten om vijf van de meest gebruikte methoden voor RNA-seq-gegevensnormalisatie te testen: Trimmed Mean of -values ​​(TMM), Upper Quartile (UQ), Mediaan (DES), Quantile (EBS) en PoissonSeq (PS), beschreven in detail in de sectie Methoden. Alle methoden zijn getest op drie verschillende datasets: Bodymap-, Cheung- en AML-gegevens. De details met betrekking tot datasets worden beschreven in hoofdstuk 2.

3.1. Differentiële expressie (DE) analyse

Het belangrijkste doel van veel RNA-seq-experimenten is het identificeren van genen die differentieel tot expressie worden gebracht tussen de vergeleken omstandigheden. Door drie bovengenoemde datasets te analyseren, hebben we de directe impact van de hierboven beschreven normalisatiemethoden op de resultaten van DE-analyse gecontroleerd. Na elke normalisatie hebben we de lijsten met differentieel tot expressie gebrachte genen (DEG's) bepaald met behulp van de statistische test uit het edgeR-pakket. In het geval van elke dataset vergeleken we genexpressieniveaus tussen twee soorten biologische monsters en rangschikten we de genen volgens aangepast

waarden. Genen met aangepaste waarden < 0,05 werden geselecteerd als differentieel uitgedrukt.

De resultaten van DE-analyse kunnen worden vergeleken op basis van het aantal en de inhoud van DEG's. Visualisatie van de DEG's wordt weergegeven in figuur 1. De staafdiagrammen tonen de bijdrage van genen met bepaalde expressieniveaus in DEG-lijsten voor alle datasets, geselecteerd uit gegevens die zijn ingediend bij vijf geteste methoden van normalisatie en onbewerkte gegevens (RD). Zoals we kunnen zien, verschilt de invloed van normalisatiemethoden tussen datasets. In het geval van de Cheung-dataset resulteren alle methoden in een vergelijkbaar aantal DEG's en hun significante deel bestaat uit genen met lage expressieniveaus. Een meer evenwichtige bijdrage van DEG's werd waargenomen in de Bodymap-gegevensset. Het aantal DEG's met een gemiddeld expressieniveau is iets hoger dan het aantal zeer sterk of zwak tot expressie gebrachte genen. In de AML-dataset overheersen genen met een gemiddeld expressieniveau duidelijk op de lijst van DEG's. Hier zijn er ook de meest significante verschillen in het aantal DEG's tussen verschillende genormaliseerde gegevens. De meest beperkende methode leek TMM te zijn, terwijl de hoogste aantallen DEG's werden verkregen met EBS- en PS-methoden. De bijdrage van alle groepen genen is hetzelfde voor elke normalisatiemethode.

-as worden de methoden van normalisatie weergegeven, terwijl de

-as vertegenwoordigt het aantal DEG's bepaald na elke normalisatieprocedure. De balkkleuren vertegenwoordigen de groepen genen met een bepaald expressieniveau.

MA-plots, beschikbaar in aanvullende materialen (Figuur S3), werden gegenereerd voor extra visualisatie van DEG's. Ze geven de relatie weer tussen het basisgemiddelde en log2FC van de tellingen. De resultaten laten zien dat DEG's in elke MA-plot na normalisatie iets meer verspreid zijn in vergelijking met de onbewerkte gegevens. De locatie van DEG verschilt echter afhankelijk van de normalisatiemethode voor AML- en Bodymap-gegevenssets. In het geval van AML-gegevens is het ook de moeite waard om erop te wijzen dat we meer DEG's kunnen waarnemen die overexpressie zijn dan onderexpressie.

Omdat het moeilijk is om te beoordelen welke normalisatiemethode geschikter is voor een dataset op basis van MA-plotanalyse, is een nauwkeurigere verificatie noodzakelijk. Hiervoor hebben we bias- en variantiewaarden berekend.

3.2. Bias- en variantieberekening

Op basis van details die in de vorige sectie zijn beschreven, hebben we huishoudgenen geselecteerd waarvoor we bias- en variantiewaarden hebben berekend. In navolging van het idee beschreven in de studie [19], gingen we ervan uit dat de meest geschikte normalisatiemethode die is die de laagste bias- en variantiewaarden voor de controlegenen genereert. De bias- en variantiewaarden werden berekend volgens de formules (11) en (12) voor alle huishoudgenen die afzonderlijk voor elke dataset waren geselecteerd, zoals beschreven in Paragraaf 2. Vervolgens, voor elke dataset, het gemiddelde van de bias- en variantiewaarden voor elke normalisatiemethode en voor alle controlegenen werd berekend. Tabellen 3 en 4 geven de rangorde weer van normalisatiemethoden op basis van deze bias- en variantiewaarden. De methode met de laagste bias of variantie werd gerangschikt als 1. De hoogste bias- en variantiewaarden, althans in de Cheung- en Bodymap-gegevenssets, werden waargenomen voor de niet-genormaliseerde onbewerkte gegevens (RD), ter vergelijking in de tabellen opgenomen. Rekening houdend met al deze datasets is de conclusie dat de beste resultaten werden verkregen met behulp van DES-, EBS- en PS-normalisatiemethoden, en de methoden die de hoogste bias- en variantiewaarden genereerden, waren TMM en UQ. In het geval van AML-gegevens leidde de toepassing van de TMM-methode tot een toename van de vertekening en de variantie die aanwezig is in RD. Het is vermeldenswaard dat de verschillen tussen hen in alle datasets klein zijn.

3.3. Gevoeligheid en specificiteit

De gevoeligheid en de specificiteit van de normalisatiemethoden werden onderzocht met behulp van de AML-dataset, op basis van onze eerdere ervaring met de analyse van microarray-gegevens, beschreven in [20], evenals bewijs uit de literatuur [23, 24]. Eerst werden de sets genen geselecteerd als positieve en negatieve controles. De positieve controles bestonden uit de genen die sterke kandidaten waren voor DEG's. Voor de set positieve controles selecteerden we genen die werden gevalideerd door real-time polymerasekettingreactie (PCR) -analyse of beschreven in de literatuur als overexpressie (of, minder vaak, onderexpressie) in AML of onrijpe hematopoëtische cellen [23]. De negatieve controles omvatten de genen die geen DEG's zijn (hun expressieniveaus mogen niet verschillen tussen deze twee soorten monsters) [24]. In totaal werden 44 genen gekozen als positieve controles en 44 genen werden gekozen als negatieve controles. De lijsten met genen zijn te vinden in Tabel S3. De gevoeligheid en specificiteit van de normalisatiemethoden werden respectievelijk berekend als een percentage positieve controles die aanwezig waren en het percentage negatieve controles dat afwezig waren in elke lijst van differentieel tot expressie gebrachte genen. De normalisatiemethoden met de hoogste waarden van specificiteit tonen beter niet-differentieel tot expressie gebrachte genen, terwijl aan de andere kant de methoden met de hoogste gevoeligheidswaarden wijzen op een grote kans om DEG's te vinden die echt DEG's zijn. De resultaten van deze analyse zijn weergegeven in tabel 5.

Tabel 5 laat zien dat de specificiteitswaarden voor EBS- en PS-methoden aanzienlijk lager zijn dan voor de overige. De gevoeligheidswaarden verschilden minder tussen de methoden en waren over het algemeen laag, tussen 10 en 31%. Bovendien kregen we voor TMM- en EBS-methoden de meest uiteenlopende resultaten: de hoogste specificiteit maar de laagste sensitiviteit voor TMM, en het tegenovergestelde voor EBS. In het geval van specificiteit kunnen we zien dat de meeste methoden specificiteitswaarden van meer dan 80% opleveren.

3.4. Voorspellingsfouten

Tabel 6 geeft een vergelijking van de vijf normalisatiemethoden waarbij voor elke dataset de verschillende aantallen informatieve genen werden gebruikt op basis van vijf classifiers en LOOCV. In elk geval selecteerden we het aantal differentieel tot expressie gebrachte genen als 75% van het aantal monsters. Daarom hebben we voor de Cheung-, Bodymap- en AML-gegevenssets respectievelijk 30 (

) differentieel tot expressie gebrachte genen. De resultaten in de tabel suggereren dat PS, DES en EBS in alle datasets beter presteren dan TMM en UQ.

3.5. Diagnostische plots

Naast de analytische methoden voor het vergelijken van normalisatiemethoden, stellen we voor om aanvullende determinanten te gebruiken die nuttig kunnen zijn voor de afwijzing van de meest gebruikte methoden die duidelijk falen. Het is mogelijk dat de normalisatiemethoden verschillende resultaten opleveren voor verschillende datasets. Daarom stellen we voor om de volgende workflow toe te passen op basis van diagnostische plots om te bepalen welke normalisatiemethode optimaal is voor een specifieke dataset. Hier richten we ons alleen op de AML-dataset. In figuur 2 (a) kunnen we de verschillen observeren die zijn geïntroduceerd in normalisatiefactoren die door elke normalisatiemethode zijn verkregen. Uit deze figuur kunnen we concluderen dat normalisatiecoëfficiënten bepaald met TMM- en UQ-methoden zich groeperen en afwijken van de rest van normalisatiemethoden.

De resultaten in Tabel 6 kunnen worden samengevat rekening houdend met de gemiddelde fouten uitgedrukt in cijfers. Voor de AML-dataset kunnen de prestaties van de vijf normalisatiemethoden worden gerangschikt. De percentages classificatiefouten verkregen door normalisatiemethoden zouden kunnen worden gebruikt om een ​​95% betrouwbaarheidsinterval te verkrijgen van het gemiddelde van het aandeel classificatiefouten voor de normalisatie. De overeenkomstige staafdiagram die betrouwbaarheidsintervallen uitzet wordt gegeven in figuur 2(b). Deze grafiek geeft aan dat de TMM-, UQ- en DES-methoden beter presteren dan de twee andere methoden met betrekking tot classificatieprestaties.

3.6. Gemeenschappelijke DEG's

Om het aantal DE-genen en het aantal algemene DE-genen te vergelijken dat wordt gevonden bij de normalisatiemethoden die voor een bepaalde dataset zijn uitgevoerd, genereren we ballonplots (Figuur 2(c) en Figuur S5) en Venn-diagrammen (zie Figuur S4). Ballonplots vertegenwoordigen percentages van het aantal algemeen gedetecteerde differentieel tot expressie gebrachte genen tussen de e en e methoden. Eerst berekenen we in de e cel een deel van de gemeenschappelijke detectie met betrekking tot de e methode:

, is het aantal differentieel tot expressie gebrachte genen dat gewoonlijk wordt gedetecteerd met de e methode en de e methode, en is het aantal differentieel tot expressie gebrachte genen dat wordt gedetecteerd met de e methode. Vervolgens nemen we de gemiddelde procentuele waarde van gemeenschappelijke genen tussen elk paar methoden. De meest geprefereerde methode is die met het hoogste aantal gemeenschappelijke DEG's.

Voor de AML-gegevensset (Figuur 2(c)) wordt het laagste aantal gemeenschappelijke DEG's met andere methoden geproduceerd door de TMM-normalisatiemethode en het beste door de EBS-methode. Het laagste aantal gemeenschappelijke DEG's werd verkregen voor de EBS-TMM- en PS-TMM-vergelijkingen. Uit deze analyse concludeerden we dat de set genen geïdentificeerd als DEG's niet stabiel is en afhankelijk is van zowel de methode van datanormalisatie als de dataset zelf. In het geval van een AML-dataset is er echter een set van 227 genen die gemeenschappelijk zijn voor alle geteste methoden (aanvullende figuur S4) en deze genen kunnen worden beschouwd als de sterkste kandidaten voor DEG's (figuur S4). In het geval van de Cheung-dataset merkten we op dat de meeste methoden op dezelfde manier lijken te presteren. Voor de Bodymap-dataset leverden alle normalisatiemethoden iets lagere percentages van algemeen gedetecteerde DEG's op dan in het geval van de Cheung-dataset. Voor de onbewerkte gegevens (RD) was het aantal gemeenschappelijke genen voor elk paar methoden minder dan 50% (zie figuur S5). De Venn-diagrammen voor de Cheung- en Bodymap-gegevenssets bevestigden deze tendens ook (Figuur S4).

3.7. Groepering van normalisatiemethoden

Clustering is een andere benadering om de prestaties van de normalisatiemethoden te vergelijken. Op basis van de overeenkomsten tussen de DE-genlijsten, hebben we dendrogrammen gegenereerd om gemakkelijk te zien welke methoden zich bij elkaar voegen. De exacte mate van overeenkomst was de DE-genrangschikking. In het geval van elke gegevensset werden voor vijf varianten van genormaliseerde gegevens, evenals voor onbewerkte gegevens (RD), bepaalde sets van differentieel tot expressie gebrachte genen verkregen. Eerst kozen we genen die alle zes DEG-lijsten gemeen hebben. Omdat we 20 genen gemeen hebben tussen zes methoden, hebben we de normalisatiemethoden geclusterd op basis van deze gemeenschappelijke DEG's. Vervolgens hebben we voor alle methoden deze genen gerangschikt, waardoor we ranglijsten van genen hebben verkregen. Op basis van deze lijsten hebben we de afstandsmatrix berekend met behulp van Euclidische afstand en plotdendrogrammen (Figuur 2 (d)). De dendrogrammen werden geconstrueerd op basis van hiërarchische clustering met behulp van de methode van Ward. Dit criterium is een andere benadering die de DEG-lijsten vergelijkt. Het vorige criterium (gemeenschappelijke DEG's) geeft ons de informatie over het percentage van DEG dat gemeenschappelijk is per paar methoden, terwijl dit criterium de informatie geeft over de overeenkomst tussen de methoden op basis van de volgorde van gemeenschappelijke DEG's in elke lijst.

3.8. Overzicht van rangen

Door alle criteria te combineren die in de vorige secties zijn beschreven, willen we bepalen welke van de vijf geteste normalisatiemethoden geschikt zou zijn voor de AML-gegevensset. Tabel 7 geeft een overzicht van de rangen die zijn verkregen op basis van de bias- en variantiewaarden, de voorspellingsfouten, gevoeligheid, specificiteit en het aantal veelvoorkomende DEG's voor de AML-gegevensset. In het geval dat alle criteria in tabel 7 even belangrijk zijn, kan de onderzoeker de beslissing baseren op de uiteindelijke rangorde die wordt berekend als het gemiddelde van de rangschikkingen die afzonderlijk voor elk criterium zijn vastgesteld met behulp van gekozen normalisatiemethoden.

4. Discussie

In de praktijk blijft normalisatie van high-throughput data nog steeds een belangrijke vraag en heeft het veel aandacht gekregen in de literatuur. Het toenemende aantal normalisatiemethoden maakt het voor wetenschappers moeilijk om te beslissen welke methode voor welke bepaalde dataset moet worden gebruikt. Op basis van de resultaten die in dit artikel en in [25] worden gepresenteerd, kunnen we concluderen dat normalisatie de analyse van differentiële expressie beïnvloedt, daarom is een belangrijk aspect hoe de meest verstandige methode voor de gegevens te kiezen. In dit artikel hebben we laten zien dat, afhankelijk van de datastructuur, de invloed van normalisatie verschilt (zie figuren 1 en 2(d)). In ons werk hebben we aangetoond dat op basis van sommige van deze criteria de keuze van de normalisatiemethode geschikter en robuuster kan zijn en meer automatisch kan worden gemaakt. Coëfficiënten zoals bias en variantie kunnen worden beschouwd als criteria voor een vergelijking van normalisatiemethoden. Het is de moeite waard erop te wijzen dat in ons onderzoek in de meeste gevallen inclusief de normalisatie de bias- en variantiewaarden vermindert in vergelijking met onbewerkte gegevens, wat de noodzaak van normalisatie bevestigt. Wanneer de verschillen tussen de bias- en variantiewaarden significant zijn, geeft het gebruik van rangen van deze waarden de werkelijke verschillen nauwkeuriger weer. In ons onderzoek zijn de verschillen tussen de verkregen bias- en variantiewaarden voor alle methoden in alle datasets echter klein. In dat geval weerspiegelt het gebruik van rangen niet de werkelijke verschillen tussen methoden en zijn aanvullende criteria nodig. De diagnostische plots zouden als dergelijke aanvullende determinanten kunnen dienen en kunnen nuttig zijn voor de afwijzing van de meest toonaangevende methoden die duidelijk falen. Uit ons onderzoek blijkt dat het gebruik van de TMM-methode in de meeste gevallen slecht wordt weergegeven. Deze conclusie komt niet overeen met de evaluatie van [26, 27]. Uit hun onderzoek bleek dat het gebruik van de TMM-methode leidde tot goede prestaties op gesimuleerde datasets. Een reden voor de onenigheid van deze resultaten kan te wijten zijn aan verschillende benaderingen van vergelijking en het gebruik van criteria die niet door andere auteurs worden overwogen. Een andere reden voor verschillen in de conclusies kan voortkomen uit het aantal biologische replica's dat in onze studie is gebruikt of kan verband houden met de specifieke datasets die in deze artikelen worden gebruikt.

Verder vinden we dat andere conclusies beschreven in [26] consistent zijn met onze eigen resultaten. Deze resultaten bevestigen de bevredigende resultaten van de DESeq-methode. In tabel 7 hebben we de samenvatting van de rangen voor de AML-dataset gepresenteerd. Het is mogelijk dat de normalisatie van andere datasets andere resultaten kan opleveren dan die verkregen met de AML-dataset (in dit artikel hebben we laten zien dat, afhankelijk van de datastructuur, de invloed van normalisatie verschilt). Over het algemeen richt elk criterium dat in het document wordt voorgesteld zich op verschillende aspecten van vergelijking. Afhankelijk van de hoofddoelstellingen van het onderzoek kunnen sommige criteria nuttiger zijn, bijvoorbeeld als de impact ligt in een goede voorspelling op basis van gekozen genen, zal het belangrijkste aspect voorspellingsfouten zijn. In sommige gevallen zijn de resultaten van criteria echter niet doorslaggevend. In een dergelijke situatie raden we de toepassing van de volgende workflow aan om te bepalen welke normalisatiemethode optimaal is voor een specifieke dataset: (i) normaliseer de gegevens met behulp van weloverwogen methoden, (ii) bereken de "bias" en "variantie" en rangschik de methoden op basis van deze waarden, (iii) na elke normalisatie differentiële analyse uitvoeren en DEG-lijsten bepalen die door elke normalisatiemethode zijn gevonden, (iv) een subset van genen selecteren die kunnen dienen als positieve en negatieve controles om de gevoeligheid en specificiteit van normalisatiemethoden te onderzoeken en rangschik de methoden op basis van deze criteria, (v) bereken het percentage van het gemiddelde van de voorspellingsfouten die zijn verkregen met behulp van gekozen classificaties voor DEG's die door elke normalisatiemethode worden gevonden en rangschik ze, (vi) teken Venn-diagrammen of ballongrafieken op basis van het aantal differentieel tot expressie gebrachte genen en rangschik de methoden op basis van het aantal gemeenschappelijke DEG-waarden, en (vii) kies op basis van de samenvatting van de rangen de meest geschikte normalisatie me van de onderzochte dataset. We merken hier dat de normalisatiemethode de expressieresultaten kan beïnvloeden, wat leidt tot foutieve DE-analyse, daarom is het erg belangrijk om in dit stadium van de analyse moeite te doen. Ten slotte wilden we de aandacht vestigen op het feit dat ons artikel niet duidelijk aangeeft welke normalisatiemethode de beste is, maar het geeft een nieuwe kijk op hoe de normalisatie voor RNA-Seq-gegevensanalyse moet worden gekozen om foutieve DE-analyse te voorkomen. Het kan niet alleen worden toegepast op methoden met betrekking tot sequencing-diepte, maar het voorgestelde algoritme is ook geschikt om normalisaties te vergelijken die rekening houden met andere bronnen van ongewenste variatie.

5. Conclusies

In de studie werden nieuwe coëfficiënten zoals bias en variantie voorgesteld als objectieve criteria voor een vergelijking van normalisatiemethoden. Concluderend suggereren onze resultaten dat, afhankelijk van de RNA-seq-gegevensstructuur en de toegepaste methode, de invloed van normalisatie verschilt. De gepresenteerde criteria, met name bias en variantie, kunnen echter de keuze van de optimale normalisatiemethode voor een specifieke dataset ondersteunen.

Belangenverstrengeling

De auteurs verklaren dat er geen belangenconflict is met betrekking tot de publicatie van dit artikel.

Dankbetuigingen

De auteurs willen de recensent bedanken voor de grondige evaluatie van het artikel en de waardevolle opmerkingen. Dit werk werd ondersteund door het KNOW-programma van het ministerie van Wetenschap en Hoger Onderwijs van de Republiek Polen.

Aanvullende materialen

In aanvullende materialen bieden we aanvullende details over datasets die in het gepresenteerde onderzoek zijn gebruikt. Verder bieden we hier aanvullende resultaten die niet zijn opgenomen in de hoofdtekst voor vergelijking van normalisatiemethoden, die andere bronnen van variatie bevatten, zoals GC-inhoud en batchcorrectie.

Referenties

  1. S. Zhao, W.-P. Fung-Leung, A. Bittner, K. Ngo en X. Liu, "Vergelijking van RNA-Seq en microarray in transcriptoomprofilering van geactiveerde T-cellen", PLoS ONE, vol. 9, nee. 1, artikel-ID e78644, 2014. Bekijk op: Publisher-site | Google geleerde
  2. E. T. Wang, R. Sandberg, S. Luo et al., "Alternatieve isovormregulatie in transcriptomen van menselijk weefsel", Natuur, vol. 456, nee. 7221, blz. 470-476, 2008. Bekijk op: Publisher Site | Google geleerde
  3. Q. Pan, O. Shai, L.J. Lee, B.J. Frey en B.J. Blencowe, "Diep onderzoek van alternatieve splitsingscomplexiteit in het menselijke transcriptoom door high-throughput sequencing," Natuurgenetica, vol. 40, nee. 12, pp. 1413-1415, 2008. Bekijk op: Publisher Site | Google geleerde
  4. W. M. Landau en P. Liu, "Dispersieschatting en het effect ervan op testprestaties in RNA-seq-gegevensanalyse: een op simulatie gebaseerde vergelijking van methoden," PLoS ONE, vol. 8, nee. 12, Artikel-ID e81415, 2013. Bekijk op: Publisher-site | Google geleerde
  5. A. Mortazavi, B.A. Williams, K. McCue, L. Schaeffer en B. Wold, "In kaart brengen en kwantificeren van zoogdiertranscriptomen door RNA-Seq," Natuurmethoden, vol. 5, nee. 7, pp. 621-1628, 2008. Bekijk op: Publisher Site | Google geleerde
  6. A. Oshlack en M. J. Wakefield, "Vertekening van transcriptielengte in RNA-seq-gegevens verwart systeembiologie," Biologie Direct, vol. 4, artikel 14, 2009. Bekijk op: Publisher Site | Google geleerde
  7. J.K. Pickrell, J.C. Marioni, A.A. Pai et al., "Mechanismen begrijpen die ten grondslag liggen aan menselijke genexpressievariatie met RNA-sequencing," Natuur, vol. 464, nee. 7289, blz. 768–772, 2010. Bekijk op: Publisher Site | Google geleerde
  8. D. Risso, K. Schwartz, G. Sherlock et al., "GC-inhoudnormalisatie voor RNA-seq", Tech. Rep. 291, Afdeling Biostatistiek, Universiteit van Californië, Berkeley, 2011. Bekijk op: Google Scholar
  9. J. T. Leek, R. B. Scharpf, H. C. Bravo et al., "De wijdverbreide en kritieke impact van batcheffecten in gegevens met hoge doorvoer aanpakken", Natuur beoordelingen Genetica, vol. 11, nee. 10, blz. 733–739, 2010. Bekijk op: Publisher Site | Google geleerde
  10. J.H. Bullard, E. Purdom, K.D. Hansen en S. Dudoit, "Evaluatie van statistische methoden voor normalisatie en differentiële expressie in mRNA-Seq-experimenten," BMC Bio-informatica, vol. 11, artikel 94, 2010. Bekijk op: Publisher Site | Google geleerde
  11. M. D. Robinson en A. Oshlack, "Een schaalnormalisatiemethode voor differentiële expressie-analyse van RNA-seq-gegevens," Genoombiologie, vol. 11, nee. 3, artikel r25, 2010. Bekijk op: Publisher Site | Google geleerde
  12. M. D. Robinson, D. J. McCarthy en G. K. Smyth, "EdgeR: een biogeleiderpakket voor differentiële expressie-analyse van digitale genexpressiegegevens," Bio-informatica, vol. 26, nee. 1, blz. 139-140, 2010. Bekijk op: Publisher Site | Google geleerde
  13. S. Anders en W. Huber, "Differentiële expressieanalyse voor sequentietellingsgegevens," Genoombiologie, vol. 11, nee. 10, artikel R106, 2010. Bekijk op: Publisher Site | Google geleerde
  14. N. Leng, J. Dawson, J. Thomson et al., "EBSeq: een empirisch Bayes-hiërarchisch model voor gevolgtrekking in RNA-seq-experimenten", Tech. Rep. 226, Universiteit van Wisconsin, 2012. Bekijk op: Google Scholar
  15. J. Li, D. M. Witten, I. M. Johnstone en R. Tibshirani, "Normalisatie, testen en schatting van valse ontdekkingssnelheden voor RNA-sequencinggegevens", Biostatistieken, vol. 13, nee. 3, pp. 523–538, 2012. Bekijk op: Publisher Site | Google geleerde
  16. YW Asmann, B.M. Necela, K.R. Kalari et al., "Detectie van overtollige fusietranscripten als biomarkers of ziektespecifieke therapeutische doelen bij borstkanker," Kankeronderzoek, vol. 72, nee. 8, blz. 1921-1928, 2012. Bekijk op: Publisher Site | Google geleerde
  17. V.G. Cheung, R.R. Nayak, I.X. Wang et al., "Polymorfe cis- en trans-regulatie van menselijke genexpressie", PLoS Biologie, vol. 8, nee. 9, artikel-ID e1000480, 2010. Bekijk op: Publisher-site | Google geleerde
  18. A.C. Frazee, B. Langmead en J.T. Leek, "Recount: een multi-experimentbron van analyseklare RNA-seq-gentellingsdatasets," BMC Bio-informatica, vol. 12, artikel 449, 2011. Bekijk op: Publisher Site | Google geleerde
  19. C. Argyropoulos, A. A. Chatziioannou, G. Nikiforidis, A. Moustakas, G. Kollias en V. Aidinis, "Operationele criteria voor het selecteren van een cDNA-microarray-gegevensnormalisatie-algoritme," Oncologische rapporten, vol. 15, nee. 4, blz. 983-996, 2006. Bekijken op: Google Scholar
  20. B. Uszczynska, J. Zyprych-Walczak, L. Handschuh et al., "Analyse van boutique-arrays: een universele methode voor de selectie van de optimale procedure voor gegevensnormalisatie," International Journal of Molecular Medicine, vol. 32, nee. 3, pp. 668–684, 2013. Bekijk op: Publisher Site | Google geleerde
  21. D. Chen, Z. Liu, X. Ma en D. Hua, "Genen selecteren op basis van teststatistieken," Journal of Biomedicine and Biotechnology, vol. 2005, nee. 2, pp. 132-138, 2005. Bekijk op: Publisher Site | Google geleerde
  22. R Ontwikkeling kernteam, R: Een taal en omgeving voor statistische berekeningen, R Foundation for Statistical Computing, Wenen, Oostenrijk, 2011.
  23. M. Goswami, N. Hensel, B.D. Smith et al., "Expressie van vermeende doelen van immunotherapie bij acute myeloïde leukemie en gezonde weefsels", Leukemie, vol. 28, nee. 5, pp. 1167-1170, 2014. Bekijk op: Uitgeverssite | Google geleerde
  24. E. Eisenberg en E. Y. Levanon, "Human house-genen, herzien," Trends in genetica, vol. 29, nee. 10, pp. 569–574, 2013. Bekijk op: Uitgeverssite | Google geleerde
  25. F. Seyednasrollah, A. Laiho en L.L. Elo, "Vergelijking van softwarepakketten voor het detecteren van differentiële expressie in RNA-seq-onderzoeken," Briefings in bio-informatica, vol. 16, nee. 1, pp. 59-70, 2015. Bekijk op: Publisher Site | Google geleerde
  26. M.-A. Dillies, A. Rau, J. Aubert et al., "Een uitgebreide evaluatie van normalisatiemethoden voor Illumina high-throughput RNA-sequencing data-analyse," Briefings in bio-informatica, vol. 14, nee. 6, pp. 671–683, 2013. Bekijk op: Publisher Site | Google geleerde
  27. E. Maza, P. Frasse, P. Senin, M. Bouzayen en M. Zouine, "Vergelijking van normalisatiemethoden voor differentiële genexpressie-analyse in RNA-Seq-experimenten: een kwestie van relatieve grootte van bestudeerde transcriptomen," Communicatieve & Integratieve Biologie, vol. 6, nee. 6, artikel-ID e25849, 2013. Bekijk op: Publisher-site | Google geleerde

Auteursrechten

Copyright © 2015 J. Zyprych-Walczak et al. Dit is een open access-artikel dat wordt gedistribueerd onder de Creative Commons Attribution-licentie, die onbeperkt gebruik, distributie en reproductie in elk medium toestaat, op voorwaarde dat het originele werk correct wordt geciteerd.


Analysemethoden

CB Johansson, . K. Roeser, in uitgebreide biomaterialen, 2011

3.313.7.8.2 Genexpressie

Genexpressieanalyses naderen het onderzoeksveld van biomaterialen. De eerste rapporten onthulden de mogelijkheid om functionele hechting en de biologische mechanismen gerelateerd aan de biomechanische gegevens te bestuderen door informatie te combineren die is verzameld uit genexpressieresultaten van eiwitten die verband houden met botremodellering in teruggewonnen weefsel. 65 Andere genexpressiestudies waarbij zowel het implantaat als het weefsel rondom het implantaat zijn onderzocht, onthullen interessante gegevens die belangrijke kennis zullen toevoegen aan het gebied van osseo-integratie die voorheen niet mogelijk was. 66 Nog een ander recent uitgevoerd onderzoek betrof het verwijderen van koppelgegevens gekoppeld aan de hoeveelheid verschillende botvorming en botresorptie, evenals ontstekingsgenen, zowel op het teruggewonnen implantaat als in het omringende weefsel. 41 In een onderzoek waarin machinaal bewerkte (controle) en gezandstraalde (test) implantaten in het dijbeen van muizen werden vergeleken, vonden we opmerkelijke histologische verschillen in botvorming, en de overeenkomstige resultaten van de Northern-blot-analyse toonden een significant hogere expressie van osteocalcine-mRNA voor de testgroep ( Afbeelding 12 ) (niet-gepubliceerde gegevens). Dit onderzoeksgebied krijgt nu meer aandacht en zal naar verwachting nog verder groeien. Nogmaals, om zoveel mogelijk tests uit te kunnen voeren op hetzelfde monster, moet een zorgvuldig onderzoeksontwerp worden geschetst voordat het onderzoek begint.

Afbeelding 12 . (a) Northern-blot-analyse van machinaal bewerkte en gezandstraalde implantaten die in het dijbeen van de muis zijn geplaatst. De osteocalcine-mRNA-expressie werd gevolgd van dag 1 tot dag 7. (b) Overeenkomstige HE (hematoxyline-eosine) gekleurde paraffinesectie op dag 3 en 1 week (dag 7) voor beide groepen. De implantaatdiameter is in millimeters.


Veelgestelde vragen over het WGCNA-pakket

Op deze pagina vindt u een lijst met veelgestelde vragen en onze veelgestelde antwoorden. Lees deze voordat u ons een e-mail stuurt over een probleem. Deze FAQ is voor het laatst bijgewerkt op 24 december 2017.

Vragen over gegevensanalyse

We raden af ​​om WGCNA te proberen op een dataset die uit minder dan 15 samples bestaat. In een typische high-throughput-omgeving zullen correlaties op minder dan 15 monsters gewoon te veel ruis zijn om het netwerk biologisch zinvol te laten zijn. Indien mogelijk, zou men minstens 20 monsters moeten hebben, zoals bij alle analysemethoden, meer monsters leiden meestal tot robuustere en verfijndere resultaten.

Probesets of genen kunnen worden gefilterd op gemiddelde expressie of variantie (of hun robuuste analogen zoals mediane en mediane absolute afwijking, MAD), aangezien genen met een lage of niet-variabele gewoonlijk ruis vertegenwoordigen. Of het beter is om te filteren op gemiddelde expressie of variantie is een kwestie van debat, beide hebben voor- en nadelen, maar wat nog belangrijker is, ze hebben de neiging om vergelijkbare sets genen eruit te filteren, aangezien gemiddelde en variantie meestal gerelateerd zijn.

We Niet doen raden aan om genen te filteren op differentiële expressie. WGCNA is ontworpen als een analysemethode zonder toezicht die genen clustert op basis van hun expressieprofielen. Het filteren van genen op differentiële expressie zal leiden tot een reeks gecorreleerde genen die in wezen een enkele (of enkele sterk gecorreleerde) modules zullen vormen. Het maakt ook de aanname van de schaalvrije topologie volledig ongeldig, dus het kiezen van zacht drempelvermogen door schaalvrije topologie-fit zal mislukken.

  • Ondertekende netwerken. De keuze van ondertekende versus niet-ondertekende netwerken is complex, maar over het algemeen geven we de voorkeur aan ondertekende (of "ondertekende hybride") netwerken boven niet-ondertekende netwerken. Om ondertekende netwerken te construeren, gebruikt u argument type = "signed" of type = "signed hybrid" in functies zoals nauwkeurigheidMeasures, adjacency, chooseOneHubInEachModule, chooseTopHubInEachModule, nearNeighborConnectivity, nearNeighborConnectivityMS, orderBranchesUsingHubGenes, softConnectivity help en mogelijk andere (zie het helpbestand voor elke functie in geval van twijfel). Sommige functies gebruiken het argument netwerktype opmerkelijke voorbeelden van netwerktypes zijn blockwiseModules, blockwiseConsensusModules, blockwiseIndividualTOMs, consensusTOM, intramodularConnectivity, modulePreservation, pickSoftThreshold, TOMsimilarityFromExpr, vectorTOM maar er zijn er ook andere. Nogmaals, lees bij twijfel het helpbestand.
  • Robuuste correlatie. De standaard correlatiemethode in alle functies in WGCNA is de standaard Pearson-correlatie. Over het algemeen raden we (en gebruiken we het zelf) de biweight-middencorrelatie aan als een robuust alternatief, tenzij er een goede reden is om aan te nemen dat er geen uitbijtermetingen zijn. Dit is geïmplementeerd in de WGCNA-functie bicor. Veel WGCNA-functies gebruiken het argument corFnc waarmee men een alternatieve correlatiefunctie kan specificeren voor de standaard cor en bicor is een optie. Aanvullende argumenten voor de correlatiefunctie kunnen worden opgegeven met het argument corOptions (afhankelijk van de functie kan voor dit argument een van de twee alternatieve vormen nodig zijn, zie de help voor elke functie voor details). In bepaalde functies, met name die van de blockwise-familie, kan de correlatiefunctie niet direct als een functie worden gespecificeerd, maar moet men het argument corType gebruiken om ofwel Pearson of biweight mid-correlatie te specificeren.

  • Beperking van het aantal uitgesloten uitbijters: argument maxPOutliers. De standaardversie van de dubbelgewicht middencorrelatie, beschreven in Langfelder en Horvath (2011) (link naar artikel), kan ongewenste resultaten opleveren wanneer de gegevens een bimodale verdeling hebben (bijv. wanneer een genexpressie sterk afhankelijk is van een binaire variabele zoals ziektestatus of genotype) of wanneer een van de variabelen die deel uitmaken van de correlatie zelf binair (of ordinaal) is. Om deze reden raden we ten zeerste aan om het argument maxPOutliers = 0,05 of 0,10 te gebruiken wanneer de biweight-middencorrelatie wordt gebruikt. Dit argument dwingt bicor er in wezen toe om nooit meer dan het gespecificeerde aantal steekproeven als uitbijters te beschouwen.
  • Omgaan met binaire gegevens. Bij het relateren van high-throughput data x aan binaire variabele y zoals steekproefkenmerken, kan men argument robustY = FALSE gebruiken om de robuuste behandeling voor het y argment van bicor uit te schakelen. Dit resulteert in een hybride robuuste Pearson-correlatie zoals beschreven in Langfelder en Horvath (2011). De hybride correlatie kan ook worden gebruikt wanneer een van de ingangen numeriek is maar waarvan bekend is dat deze geen uitbijters heeft.

Ja. Wat WGCNA betreft, is het werken met (goed genormaliseerde) RNA-seq-gegevens niet echt anders dan het werken met (goed genormaliseerde) microarray-gegevens.

We raden aan om objecten te verwijderen waarvan het aantal constant laag is (bijvoorbeeld het verwijderen van alle objecten met een telling van minder dan bijvoorbeeld 10 in meer dan 90% van de steekproeven), omdat dergelijke laag uitgedrukte kenmerken de neiging hebben om ruis en correlaties weer te geven op basis van tellingen die meestal nul zijn, zijn niet echt zinvol. De werkelijke drempels moeten gebaseerd zijn op het experimentele ontwerp, de sequentiediepte en het aantal monsters.

We raden dan een variantie-stabiliserende transformatie aan. Pakket DESeq2 implementeert bijvoorbeeld de functie variantanceStabilizingTransformation die we nuttig hebben gevonden, maar men zou ook kunnen beginnen met genormaliseerde tellingen (of RPKM/FPKM-gegevens) en deze log-transformeren met log2(x+1) . Voor sterk uitgedrukte kenmerken zijn de verschillen tussen volledige variantiestabilisatie en een eenvoudige log-transformatie klein.

Of je nu RPKM, FPKM of gewoon genormaliseerde tellingen gebruikt, maakt niet veel uit voor WGCNA-analyse, zolang alle monsters zijn verwerkt dezelfde manier. Deze normalisatiemethoden maken een groot verschil als men expressie van gen A wil vergelijken met expressie van gen B, maar WGCNA berekent correlaties waarvoor gen-gewijze schaalfactoren geen verschil maken. (Samplegewijze schaalfactoren doen dat natuurlijk, dus monsters moeten wel worden genormaliseerd.)

Als gegevens uit verschillende batches komen, raden we aan om te controleren op batcheffecten en deze zo nodig aan te passen. We gebruiken ComBat voor het verwijderen van batch-effecten, maar andere methoden zouden ook moeten werken.

Ten slotte controleren we gewoonlijk kwantielspreidingsdiagrammen om er zeker van te zijn dat er geen systematische verschuivingen tussen steekproeven zijn als steekproefkwantielen correlaties vertonen (wat ze gewoonlijk doen), kwantielnormalisatie kan worden gebruikt om dit effect te verwijderen.

Heterogeniteit van gegevens kan van invloed zijn op elke statistische analyse, en nog meer op een niet-gecontroleerde analyse zoals WGCNA. Welke eventuele wijzigingen aan de analyse moeten worden aangebracht, hangt in belangrijke mate af van of de heterogeniteit (of de onderliggende driver) als "interessant" wordt beschouwd voor de vraag die de analist probeert te beantwoorden, of niet. Als je geluk hebt, is de belangrijkste oorzaak van steekproefverschillen de behandeling/conditie één-onderzoeken, in welk geval WGCNA op de gegevens kan worden toegepast zoals ze zijn. Helaas zijn de drijfveren voor heterogeniteit vaak oninteressant en moet hiervoor worden gecorrigeerd. Dergelijke factoren kunnen technisch zijn (batcheffecten, technische variabelen zoals postmorteminterval enz.) of biologisch (bijvoorbeeld verschillen in geslacht, weefsel of soort).

Als men een categorische bron van variatie heeft (bijv. geslacht of weefselverschillen) en het aantal monsters in elke categorie groot genoeg is (minstens 30 bijvoorbeeld) om een ​​netwerk in elke categorie afzonderlijk te construeren, kan het de moeite waard zijn om uit te voeren een consensusmodule-analyse (Tutorial II, zie WGCNA Tutorials). Omdat deze analyse in elke categorie afzonderlijk een netwerk construeert, heeft de variatie tussen categorieën geen invloed op de analyse.

Als het gewenst is om voor alle steekproeven één netwerk te bouwen, moet worden gecorrigeerd voor de ongewenste of oninteressante bronnen van grote variatie in de data. Voor categorische (ordinale) factoren raden we aan om de functie ComBat te gebruiken (uit het pakket sva ). Gebruikers die nog nooit ComBat hebben gebruikt, moeten het helpbestand voor ComBat lezen en het sva-vignet doorlopen (typ vignette("sva") bij de R-prompt) om er zeker van te zijn dat ze ComBat correct gebruiken.

Voor continue bronnen van variatie (bijvoorbeeld postmortem interval), kan men eenvoudige lineaire regressie gebruiken om de gegevens aan te passen. Er zijn mogelijk meer geavanceerde methoden die ook het gebruik van covariaten mogelijk maken en beschermen tegen overcorrectie.

Welke methode ook wordt gebruikt, we waarschuwen de gebruiker dat het verwijderen van ongewenste variatiebronnen nooit perfect is en in sommige gevallen kan leiden tot het verwijderen van echt interessant signaal, en in zeldzame gevallen kan het een vals associatiesignaal introduceren. Daarom moeten alleen bronnen met relatief grote variatie worden verwijderd.

Ten eerste moet de gebruiker ervoor zorgen dat variabelen (probesets, genen enz.) niet zijn gefilterd op differentiële expressie met betrekking tot een monsterkenmerk. Zie punt 2 hierboven voor details over gunstige en schadelijke filtergenen of probesets.

Als de schaalvrije topologie-fit-index geen waarden boven 0,8 bereikt voor redelijke bevoegdheden (minder dan 15 voor niet-ondertekende of ondertekende hybride netwerken en minder dan 30 voor ondertekende netwerken) en de gemiddelde connectiviteit relatief hoog blijft (in de honderden of hoger) , is de kans groot dat de gegevens een sterke driver vertonen die een subset van de steekproeven globaal anders maakt dan de rest. Het verschil veroorzaakt een hoge correlatie tussen grote groepen genen, wat de aanname van de schaalvrije topologiebenadering ongeldig maakt.

Gebrek aan schaalvrije topologie-fit op zichzelf maakt de gegevens niet ongeldig, maar moet zorgvuldig worden onderzocht. Het helpt altijd om de monsterclusterboom en alle technische of biologische monsterinformatie eronder uit te zetten, zoals in figuur 2 van Tutorial I, sectie 1 sterke clusters in de clusteringboom geven globaal verschillende groepen monsters aan. Het kan het resultaat zijn van een technisch effect zoals een batcheffect, biologische heterogeniteit (bijvoorbeeld een dataset bestaande uit monsters van 2 verschillende weefsels), of sterke veranderingen tussen omstandigheden (bijvoorbeeld in een tijdreeks). Er moet zorgvuldig worden onderzocht of er steekproefheterogeniteit is, wat de heterogeniteit drijft en of de gegevens moeten worden aangepast (zie vorig punt).

Als het ontbreken van schaalvrije topologie-fit blijkt te worden veroorzaakt door een interessante biologische variabele die men niet wil verwijderen (dwz de gegevens aanpassen voor), kan de juiste soft-thresholding power worden gekozen op basis van het aantal monsters zoals in de onderstaande tabel. Deze tabel is in december 2017 bijgewerkt om de resulterende netwerken conservatief te maken.

Aantal monsters Niet-ondertekende en ondertekende hybride netwerken Ondertekende netwerken
Minder dan 20 9 18
20-30 8 16
30-40 7 14
meer dan 40 6 12

Veel van de WGCNA-functies gebruiken meerdere argumenten die verschillende subtiliteiten in netwerkconstructie en module-identificatie beheersen. Over het algemeen proberen we standaardinstellingen te bieden die redelijk goed werken in de meest voorkomende situaties. In sommige gevallen merken we echter na verloop van tijd dat een andere instelling geschikter is. In de meeste gevallen behouden we de oude standaard voor reproduceerbaarheid.

Fouten bij het uitvoeren van zelfstudievoorbeelden

Deze fout treedt bijna altijd op omdat R het WGCNA-pakket niet kon laden. Typ in R

Als u een fout krijgt in de functie pickSoftThreshold , raadpleeg dan item 1 in Runtime-fouten. Stuur anders Peter Langfelder een e-mail. Het is zeker mogelijk dat de tutorials nog steeds onontdekte bugs bevatten, of dat onze wijzigingen aan het pakket of wijzigingen in R de beledigende tutorial hebben verbroken.

De meest waarschijnlijke boosdoener is de grootte van uw dataset. Met name secties 2.a van Tutorials I en II gaan ervan uit dat je minder dan 5000 probes in je dataset hebt. Als u meer dan dat heeft, kijk dan naar de bijbehorende paragraaf 2.c (Omgaan met grote datasets). Pas het argument maxBlockSize aan aan de mogelijkheden van uw computer. De details worden beschreven in Paragraaf 2.c van de bijbehorende zelfstudie.

Runtime-fouten

Hoewel de meeste moderne processors meerdere kernen hebben, maken sommige omgevingen, met name clusters (zoals de Sun Grid Engine-clusters), slechts één processorkern beschikbaar voor elk proces. Pogingen om meerdere threads te starten, leiden vaak tot foutmeldingen vergelijkbaar met dit voorbeeld: thread 0 kon niet succesvol worden gestart. Foutcode: 11.

Als dit gebeurt, schakel dan threading uit (bijvoorbeeld door de functie

Verschillende Mac-gebruikers hebben malloc-fouten gemeld, zoals:

R(9073,0xa013dfa0) malloc: *** mmap(size=95006720) mislukt (foutcode=12)
*** fout: kan regio niet toewijzen
*** stel een breekpunt in malloc_error_break in om te debuggen

Voor zover we weten, is dit een onechte en ongevaarlijke boodschap en kan deze worden genegeerd. Voor alle C- en gcc-savvy Mac-gebruikers die er zijn, laat het ons weten als je een oorzaak en oplossing vindt.

  • De gegevens bevatten te veel ontbrekende gegevens. We proberen de code bestand te maken tegen ontbrekende gegevens, maar soms missen we iets. Probeer als test de ontbrekende gegevens toe te rekenen en voer een testrun uit met geïmputeerde gegevens.
  • De gegevens kunnen te weinig monsters of sondes bevatten.
  • De invoergegevens zijn niet numeriek. Dit is vaak een subtiel en moeilijk te vatten probleem. Bij het inlezen van gegevens uit teksttabellen (door tabs gescheiden .txt- of komma-gescheiden .csv-bestanden), kan R een dataframe converteren naar een lijst en/of sommige kolommen converteren naar een teken of een factor. Dit kan bijvoorbeeld gebeuren omdat ontbrekende gegevens in uw gegevensbestand worden gecodeerd als NULL of N/A , terwijl R NA verwacht en al het andere wordt behandeld als een tekenreeks.Het is het beste om dergelijke problemen bij de bron op te lossen, maar het is misschien niet altijd mogelijk en de gebruiker moet een versie van
  • Last but not least, uw gegevens kunnen prima in orde zijn, maar leiden tot geen modules, of misschien slechts één module, en dergelijke gevallen worden mogelijk niet correct afgehandeld door de code (d.w.z. de code bevat fouten). Stuur Peter Langfelder een e-mail als u een dergelijke situatie tegenkomt.

Sommige fouten kunnen worden veroorzaakt door conversie tussen matrices en dataframes. Voor de meeste doeleinden zijn een data.frame en een matrix equivalent en zullen veel functies met beide typen argumenten worden uitgevoerd. Een grote uitzondering is de verwerking van kolomnamen. Hoewel kolomnamen van matrices vrij willekeurige kolomnamen zijn, moeten gegevensframes beginnen met een letter, onderstrepingsteken of een punt. Dit kan een probleem zijn als bijvoorbeeld de expressiegegevens zijn opgeslagen in een matrix en de probeset-ID's beginnen met een nummer, bijvoorbeeld "1552612_at" . Wanneer geconverteerd naar een dataframe, zal R een "X" toevoegen aan elke ongeldige kolomnaam, waardoor het voorbeeld "X1552612_at" wordt. Dergelijke wijzigingen kunnen fouten veroorzaken, bijvoorbeeld in de functie plotNetworkHeatmap en andere.

We hebben verspreide meldingen ontvangen van crashes op Mac OSX 10.6.x (10.5.x-systemen vertonen deze fout niet). De symptomen zijn verschillende harde crashes (bevriezen, segmentatiefouten en dergelijke) die optreden wanneer u bijna alles probeert uit te voeren na het laden van het WGCNA-pakketversie 1.00 en lager. De boosdoener bleek het bij het verkeerde eind te hebben en/of Tcl/Tk verkeerd te hebben geïnstalleerd. Om deze reden hebben we de Tcl/Tk-afhankelijkheid verwijderd vanaf WGCNA versie 1.10. Installeer de nieuwe versie. Als uw problemen aanhouden, neem dan contact op met Peter Langfelder.

Installatie problemen

Voordat u tijd besteedt aan het oplossen van een installatieprobleem met het gedownloade pakket, kunt u overwegen het pakket vanuit CRAN te installeren. We begrijpen dat het upgraden van R een beetje gedoe kan lijken, maar uiteindelijk is het de moeite waard.

Vanaf R-versie 2.14.0 is de pakketimpute ingetrokken bij CRAN en is deze nu exclusief verkrijgbaar bij Bioconductor. Typ de volgende regels in een R-sessie om het te installeren:

Dit zou het pakket moeten installeren, maar let op eventuele fouten.

Een veelvoorkomende oorzaak van deze fout is dat wanneer de gebruiker het bestand opslaat, het besturingssysteem het zal decomprimeren of uitpakken. Meestal betekent dit dat de .zip- of .tar.gz-bundel wordt gedecomprimeerd en uitgepakt, waardoor het bestand onbruikbaar wordt voor R. Mac OS X lijkt bijvoorbeeld het gzip-bestand automatisch te decomprimeren. De oplossing is om het bestand op schijf op te slaan zoals het is, zonder dat een programma zoals WinZip het kan aanraken. R zal het pakket zelf decomprimeren en uitpakken. Op een Mac moet je misschien een terminal openen, naar de map gaan waar je het bestand hebt opgeslagen en typ

De beste oplossing is om je R bij te werken naar de nieuwste versie, dan gewoon R uit te voeren en het commando install.packages("WGCNA") te gebruiken. Als je om wat voor reden dan ook je R niet kunt of wilt updaten, bekijk dan de installatie-instructies en zorg ervoor dat je de vereiste XCode-tools hebt geïnstalleerd.

Sommige gebruikers hebben gemeld dat een pakket met de naam gfortran.pkg ook nodig is. Dit is mogelijk een nieuwe functie van R vanaf versie 2.9.0.

Update: Vanaf versie 1.10 vereist WGCNA geen qvalue. Installeer de nieuwste versie van WGCNA. Dit zou alle installatieproblemen van Tcl/Tk en qvalue moeten oplossen.

Compatibiliteitsproblemen met versies

Algemene vragen

Bij het construeren van een netwerk op basis van een gegevensset van een typische genomische grootte (dwz tussen 10 000 en 30 000 genen of andere variabelen), is de meest tijdrovende stap de berekening van de Topologische Overlap Matrix, waarbij matrices worden vermenigvuldigd met tienduizenden rijen en kolommen. Met een standaard R-distributie kan dit meerdere uren duren, zelfs op een modern werkstation, aangezien matrixvermenigvuldiging in standaard R geen voordeel haalt uit multithreading (parallelle uitvoering). Het is mogelijk om dit proces met een factor 10-100 te versnellen door een voor snelheid geoptimaliseerde Basic Linear Algebra Subprograms (BLAS)-bibliotheek te installeren en R ertegen te compileren. Het proces van het compileren van R tegen een verbeterde BLAS-bibliotheek wordt beschreven in de installatie- en beheerhandleiding van R. Het compileren van R op Linux- en Unix-smaken is meestal relatief eenvoudig en duidelijk. Op Mac OSX en (meer nog) op Windows vereist het de installatie van extra tools en pakketten. Hoewel het handig is om beheerdersrechten te hebben om R te compileren en te installeren, is dit meestal niet nodig. Zie de R installatie- en beheerhandleiding voor volledige details.

Een deel van de WGCNA-code is geschreven om te profiteren van parallelle uitvoering om de berekeningen te versnellen. Er zijn twee belangrijke parallelle berekeningsmechanismen die door WGCNA worden gebruikt: de gecompileerde functies cor en bicor gebruiken POSIX-threading voor een deel, maar niet voor de hele berekening. Deze parallelle code is alleen beschikbaar op platforms met POSIX-threads (verschillende Linux- en Unix-varianten en Mac OS). Het is niet beschikbaar op Windows. Sommige functies (zoals pickSoftThreshold) kunnen POSIX-threads gebruiken waar ze beschikbaar zijn en SNOW-clusters gebruiken waar multi-threading niet beschikbaar is. De gebruikers moeten zich ervan bewust zijn dat POSIX en parallelle uitvoering van clusters heel verschillend zijn en niet uitwisselbaar zijn, in feite wijzen veel clusteromgevingen slechts een enkele kern toe aan elk proces ("taak") en pogingen om extra threads te starten leidt tot fouten (zie Runtime-fouten onderstaand).

    Roep binnen R de functie . aan

als u bijvoorbeeld 2 cores heeft (of 2 cores wilt gebruiken),

Houd er rekening mee dat deze instelling geen invloed heeft op de multi-threading-status van de onderliggende BLAS-bibliotheek.


Achtergrond

Kwantitatieve reverse-transcriptase-polymerase-kettingreactie (qRT-PCR) is een voorkeursmethode geworden voor genexpressieonderzoeken in klinische monsters, met name voor doelen met een lage kopie die van belang zijn en voor monsters van beperkte grootte [1-3]. In vergelijking met microarrays [4], profiteert qRT-PCR van een breed dynamisch bereik, gevoeligheid en maakt nauwkeurige kwantificering mogelijk [5, 6].

Om veranderingen in het expressieniveau van doelgenen door qRT-PCR nauwkeurig te kwantificeren, moet men echter normalisatie toepassen voor heterogeniteit in klinische monsters en ook voor variabiliteit die wordt geïntroduceerd tijdens RNA-extractie en cDNA-synthese [1, 7]. Naast normalisatie naar steekproefomvang en totaal RNA, vertegenwoordigt normalisatie met behulp van endogene referentiegenen een relevante benadering [3]. Referentiegenen zouden idealiter constitutief tot expressie moeten worden gebracht door alle celtypen en mogen niet worden beïnvloed door ziekte en experimentele procedures. Tot op heden is er nog geen universeel referentiegen geïdentificeerd. Huishoudgenen (HKG's) zijn de meest gebruikte referentiegenen [1]. Hoewel HKG's door elke cel tot expressie worden gebracht, varieert hun expressie tussen verschillende celtypen/organen [8, 9]. Het gebruik van HKG's als referentiegenen voor een bepaald monstertype moet daarom worden gevalideerd.

Tot nu toe zijn slechts enkele referentiegenen gevalideerd voor cellen uit het ademhalingscompartiment, specifiek GNB2L1 werd gevalideerd voor bronchoalveolaire macrofagen bij patiënten met chronische obstructieve longziekte (COPD) [10] en GAPDH (glyceraldehyde-3-fosfaatdehydrogenase) voor niet-kleincellige longkanker [11]. De meerderheid van de studies die zijn gepubliceerd over qRT-PCR in de longomgeving, gebruiken een algemene benadering van normalisatie tegen GAPDH of ACTB (bèta-actine) [12-16]. Deze "traditionele" referentiegenen zijn echter al ongeschikt bevonden voor het normaliseren van mRNA-niveaus in astmatische luchtwegen [17, 18] en ook voor expressiestudies met bronchoalveolaire macrofagen [10].

Om geschikte referentiegenen te identificeren voor qRT-PCR-normalisatie in de setting van het bronchoalveolaire compartiment, was ons doel daarom om HKG's te identificeren met de meest stabiele mRNA-expressie in bronchoalveolaire (BAL) cellen. Onze keuze voor kandidaat-HKG's was gebaseerd op 1) hun algemeen gebruik in eerdere qRT-PCR-experimenten (ACTB, GAPDH, G6PD), 2) stabiele expressie in verschillende menselijke weefsels in microarray-experimenten (ARF1, CANX, GPS1, PSMB2, PSMD2) [ 8, 9] en 3) stabiele expressie in bronchoalveolaire macrofagen en perifere neutrofielen (GNB2L1, RPL32) [10, 19]. Om rekening te houden met variaties van het BAL-cellulaire profiel bij verschillende luchtwegaandoeningen, hebben we de stabiliteit van HKG's mRNA-expressie bestudeerd bij eenenzeventig proefpersonen over een spectrum van longpathologieën. Naast BAL cellulair profiel en type longpathologie, werden vier variabelen onderzocht op hun mogelijke invloed op mRNA-expressie van bestudeerde HKG's: roken, geslacht, behandeling en leeftijd. Verder werd de stabiliteit van mRNA-expressie van alle tien HKG's gevalideerd in het tweede, onafhankelijke BAL-cohort bestaande uit zeventien controlepersonen en drieënzestig sarcoïdosepatiënten met speciale nadruk op patiëntensubgroepen. Ten slotte hebben we door onderzoek naar mRNA-expressie van twee cytokinen waarvan bekend is dat ze geassocieerd zijn met sarcoïdose, INFG (interferon-gamma) en CCL2/MCP-1, praktisch bewijs geleverd dat normalisatie met gevalideerde referentiegenen in klinische monsters een absolute voorwaarde is voor het verkrijgen van klinisch onbevooroordeelde geldige informatie van qRT-PCR.


Kader 1: Vergelijkingen van microarrays en sequencing voor analyse van genexpressie

Er zijn nu verschillende vergelijkingen gemaakt van RNA-seq- en microarray-gegevens. Deze omvatten proof-of-principle demonstraties van het sequencing-platform [2, 31, 32], toegewijde vergelijkingsstudies [34, 75-77] en de ontwikkeling van analysemethodologie [10]. De resultaten zijn unaniem: sequencing heeft een hogere gevoeligheid en dynamisch bereik, gekoppeld aan een lagere technische variatie. Bovendien hebben vergelijkingen een sterke overeenstemming aangetoond tussen microarrays en sequencing in metingen van zowel absolute als differentiële expressie. Desalniettemin zijn microarrays zeer succesvol geweest en zijn ze nog steeds zeer succesvol in het ondervragen van het transcriptoom in veel biologische omgevingen. Voorbeelden zijn het definiëren van de cel van oorsprong voor subtypes van borstkanker [78] en het onderzoeken van het effect van evolutie op genexpressie in Drosophila [79].

Microarrays en sequencing hebben elk hun eigen specifieke vooroordelen die van invloed kunnen zijn op het vermogen van een platform om DE te meten. Het is algemeen bekend dat kruishybridisatie van microarray-probes de expressiemetingen op een niet-uniforme manier beïnvloedt [80, 81] en dat de sequentie-inhoud de gemeten probe-intensiteiten beïnvloedt [82]. Ondertussen hebben verschillende onderzoeken een GC-bias waargenomen in RNA-seq-gegevens [45] en kan RNA-seq lijden aan ambiguïteit in de afbeelding van paraloge sequenties. Bovendien is er een groter statistisch vermogen om veranderingen bij hogere tellingen te detecteren (een tweevoudig verschil van 200 reads tot 100 reads is bijvoorbeeld meer statistisch significant dan 20 reads tot 10, onder de nulhypothese van geen verschil). RNA-seq als associatie tussen DE en genlengte, een effect dat niet aanwezig is in microarraygegevens [66, 68]. Andere onderzoeken geven aan dat specifieke sequentieprotocollen vooroordelen produceren in de gegenereerde uitlezingen, die gerelateerd kunnen zijn aan de sequentiesamenstelling en afstand langs het transcript [49, 50, 83, 84]. Het is bijvoorbeeld gevonden dat bibliotheekvoorbereiding voor kleine RNA's de reeks waargenomen sequenties sterk beïnvloedt [85]. Bovendien zijn benaderingen van transcriptoomassemblage noodzakelijkerwijs bevooroordeeld door het expressieniveau omdat er minder informatie beschikbaar is voor genen die op een laag niveau tot expressie worden gebracht [11, 14]. Veel van deze vooroordelen worden nog onderzocht en slimme statistische methoden die deze kennis benutten, kunnen mogelijk verbeteringen aan bestaande methoden opleveren.

Naast het grotere dynamische bereik en de gevoeligheid van RNA-seq hebben verschillende aanvullende factoren bijgedragen aan de snelle opname van sequencing voor differentiële expressie-analyse. Ten eerste zijn microarrays eenvoudigweg niet beschikbaar voor veel niet-modelorganismen (Affymetrix biedt bijvoorbeeld microarrays voor ongeveer 30 soorten [86]). Daarentegen zijn genomen en sequentie-informatie direct beschikbaar voor duizenden soorten [87]. Bovendien, zelfs wanneer genomen niet beschikbaar zijn, kan RNA-seq nog steeds worden uitgevoerd en kan het transcriptoom nog steeds worden ondervraagd (een recent onderzoek gebruikte bijvoorbeeld RNA-seq om de celoorsprong van de gezichtstumor van Tasmaanse duivel [88]) te onderzoeken. Ten tweede geeft sequencing ongekende details over transcriptionele kenmerken die arrays niet kunnen, zoals nieuwe getranscribeerde regio's, allelspecifieke expressie, RNA-bewerking en een uitgebreide mogelijkheid om alternatieve splicing vast te leggen. Een recent RNA-seq-onderzoek [11] kon bijvoorbeeld verschillende voorbeelden van isovormomschakeling tijdens celdifferentiatie laten zien, en RNA-seq werd gebruikt om de ouder-van-oorsprong-expressie in muizenhersenen aan te tonen [5].

Sequencing is natuurlijk niet zonder uitdagingen. De kosten van het platform kunnen voor sommige onderzoeken beperkend zijn. Met de uitbreiding van de totale sequentiecapaciteit en de mogelijkheid om te multiplexen, zullen de kosten per monster om voldoende sequentiediepte te genereren echter binnenkort vergelijkbaar zijn met die van microarrays. De kosten van informatica voor het onderbrengen, verwerken en analyseren van de gegevens zijn echter aanzienlijk [89]. Onderzoekers met beperkte toegang tot computerpersoneel en middelen kunnen ervoor kiezen om microarrays te gebruiken omdat de procedures voor gegevensanalyse relatief volwassen zijn. Ten slotte is het duidelijk dat data-analysemethodologieën voor het rangschikken van gegevens nog geruime tijd zullen blijven evolueren.


Validatie van Tuba1a als geschikte interne controle voor normalisatie van genexpressie-analyse tijdens de ontwikkeling van de muislong

De expressieverhouding tussen het geanalyseerde gen en een intern controlegen is de meest gebruikte normalisatiemethode voor kwantitatieve RT-PCR (qRT-PCR) expressie-analyse. Het ideale referentiegen voor een specifiek experiment is het gen waarvan de expressie niet wordt beïnvloed door de verschillende geteste experimentele omstandigheden. In deze studie valideren we de toepasbaarheid van vijf veelgebruikte referentiegenen tijdens verschillende stadia van de ontwikkeling van de muislong. De stabiliteit van expressie van vijf verschillende referentiegenen (Tuba1a, Actb Gapdh, Rn18S en Hist4h4) werd berekend binnen vijf experimentele groepen met behulp van het statistische algoritme van geNorm-software. Over het algemeen vertoonde Tuba1a de minste variabiliteit in expressie tussen de verschillende stadia van longontwikkeling, terwijl Hist4h4 en Rn18S de maximale variabiliteit in hun expressie vertoonden. Expressie-analyse van twee longspecifieke markers, oppervlakteactief proteïne C (SftpC) en Clara celspecifiek 10 kDA-eiwit (Scgb1a1), genormaliseerd naar elk van de vijf hier geteste referentiegenen, bevestigde onze resultaten en toonde aan dat een onjuiste keuze van het referentiegen kan leiden tot artefacten. Bovendien zal een combinatie van twee interne controles voor normalisatie van expressieanalyse tijdens longontwikkeling de nauwkeurigheid en betrouwbaarheid van de resultaten vergroten.

Figuren

Genstructuur van markers vaak ...

Genstructuur van markers die vaak worden gebruikt als interne controles voor expressieanalyse tijdens ...

Karakterisering van primerparen ontworpen...

Karakterisering van primerparen ontworpen voor expressieanalyse van de interne controle tijdens...

Sequentie-identiteit van de PCR...

Sequentie-identiteit van de PCR-producten voor de vijf referentiegenen geëvalueerd in ...

Expressieanalyse van de interne…

Expressie-analyse van de genen voor interne controle die in deze studie zijn geëvalueerd tijdens verschillende ...

Berekening van de gemiddelde expressiestabiliteit…

Berekening van de gemiddelde expressiestabiliteit en paarsgewijze variatiecoëfficiënt voor elk van de...

Expressie-analyse van twee cel...

Expressie-analyse van twee cellijnmarkers tijdens de ontwikkeling van de longen van de muis met behulp van verschillende ...

Expressie-analyse van twee cel...

Expressie-analyse van twee cellijnmarkers tijdens de ontwikkeling van de longen van de muis met behulp van een ...


Methoden:

Voorbereiding van plantenmonsters

Een gemeenschapsstandaard diploïde inteeltlijn van Brachypodium distachyon, Bd21, werd gebruikt voor alle experimentele behandelingen. Het palea en het lemma werden voorzichtig met een fijne pincet van de rijpe zaden afgepeld. De gestripte zaden werden gesteriliseerd door ze gedurende 10 minuten te weken in een oplossing van 10% huishoudbleekmiddel (5,25% NaOCl) aangevuld met 0,1% Tween 20 met af en toe schudden. De gesteriliseerde zaden werden driemaal grondig gespoeld met steriel dubbel gedestilleerd water en op Murashige en Skoog (MS)-agarplaten geplaatst (hierna MS-agarplaten) (4,3 g/l MS-zouten met vitamines, 3% sucrose, pH 5,8, en 0,4% fytagel). De zaden werden gedurende 2 dagen koud behandeld bij 4°C om de ontkieming te synchroniseren en men liet ze ontkiemen in een kweekkamer die werd gecontroleerd op 25°C met een dagelijkse fotoperiodieke cyclus van 16 uur licht en 8 uur donker.

Vier weken oude Brachypodium-planten werden gebruikt voor behandelingen met groeihormonen en abiotische stress. Geoogste plantenmonsters werden ingevroren in vloeibare stikstof en bewaard bij -80°C tot RNA-extractie. Voor genexpressiestudies in verschillende plantenweefsels werden geschikte plantenweefsels geoogst van 5 weken oude, volgroeide planten. Voor genexpressiestudies in verschillende ontwikkelingsstadia werden planten geoogst op 7 dagen na ontkieming (DAG) (vroege vegetatieve fase), 12 (late vegetatieve fase), 20 (overgangsfase) en 30 (reproductieve fase) (zie aanvullend bestand 2 ).

Groeihormonen en abiotische stress

Voor groeihormoonbehandelingen werden 4 weken oude Brachypodium-planten overgebracht naar de MS-vloeistofculturen aangevuld met IAA (50 M), brassinolide (BL, 50 M), zeatine (50 M), ABA (100 M), GA (50 M) M), 1-aminocyclopropaan-1-carbonzuur (ACC, 50 M), SA (100 M), of met MeJA (100 M) en gedurende 5 uur onder zacht schudden geïncubeerd. De nepzaailingen werden op dezelfde manier geïncubeerd in een MS-vloeistof maar zonder toevoeging van groeihormonen.

Voor zout- en droogtestressbehandelingen werden 4 weken oude Brachypodium-planten overgebracht naar de MS-vloeistofculturen die waren aangevuld met respectievelijk 300 mM NaCl of 400 mM mannitol en zachtjes geschud gedurende 5 uur. Voor koude- en warmtebehandelingen werden de zaailingen respectievelijk 5 uur bij 4°C of 2 uur bij 42°C geïncubeerd.

Selectie van potentiële referentiegenen in Brachypodium

Om potentiële Brachypodium-homologen van de Arabidopsis- of rijstgenen te identificeren die gewoonlijk worden gebruikt als interne controles voor genexpressiestudies, hebben we versie 0.52 van de HarvEST:Brachypodium-software http://harvest-web.org opgevraagd, die 6 verschillende bibliotheken van Brachypodium weergeeft. Alle gensequenties werden verkregen uit de GenBank dbEST-database. De HarvEST:Brachypodium-software bevat de beste BLASTX-hits van de UniProt en de rijst- en Arabidopsis-genomen (respectievelijk TIGR-versie 5, februari 2007 en TAIR-versie 7, april 2007).

De HarvEST is in de eerste plaats een EST-software voor het bekijken van databases die de nadruk legt op de genfunctie en is gericht op vergelijkende genomica en het ontwerp van oligonucleotiden, met als doel diverse onderzoeksactiviteiten te ondersteunen, zoals het ontwerp van microarray-inhoud, functionele annotatie en fysieke en genetische mapping http: //oogst.ucr.edu. De HarvEST:Brachypodium is de meest recente en gestandaardiseerde EST-database voor Brachypodium en kan dus worden gebruikt om sequentie-uitlijningen te onderzoeken en te bepalen waar individuele sequenties betrouwbaar afwijken van een consensussequentie.

Geselecteerde Brachypodium EST's van potentiële referentiegenen werden gebruikt om primers te ontwerpen. Een set qRT-PCR-primers met hoge efficiëntie voor negen individuele referentiegenen werd ontworpen met behulp van de Primer3-software (versie 0.4.0) [24]. De primers zijn ontworpen om smelttemperaturen te hebben in een bereik van 50 – 60°C, afhankelijk van individuele genen. Hun sequenties zijn samengevat in aanvullend bestand 1.

Totale RNA-extractie en primaire cDNA-synthese

Totaal RNA werd geëxtraheerd uit geschikte plantenmonsters met behulp van de RNeasy Plant-minikit (Qiagen, Valencia, CA) volgens de procedure van de fabrikant. De kwaliteit en integriteit van de RNA-monsters werden geëvalueerd door absorptiemetingen en door elektroforetische analyse met behulp van het Labwork Image Acquisition and Analysis Program (Media Cybernetics, San Diego, CA). Alle RNA-monsters die in de qRT-PCR-reacties werden gebruikt, vertoonden een absorptieverhouding van 260/280 nm van 1,8 – 2,2. RNA-monsters met een verhouding van ≈ 2 zijn over het algemeen geschikt voor daaropvolgende enzymatische reacties. Voorafgaand aan RT-PCR en qRT-PCR werden totale RNA-monsters voorbehandeld met een RNase-vrij DNase I om verontreinigend genomisch DNA te elimineren. Het primaire cDNA werd gesynthetiseerd uit ongeveer 3 μg totaal RNA met behulp van het MMLV first-strand synthesesysteem (Promega, Madison, WI) en de oilgo-dT en willekeurige primers in een reactievolume van 40 l volgens de procedure van de fabrikant.

Om elke contaminatie van genomisch DNA in de RNA-preparaten uit te sluiten, werden de RNA-monsters en het genomische DNA parallel onderworpen aan PCR-amplificaties van de ARR4 en GAPDH gensequenties (30 cycli) en de PCR-producten werden vergeleken. Er werden geen zichtbare amplificaties van genomische DNA's gedetecteerd uit de RNA-monsters (zie aanvullend bestand 15).

RT-PCR en qRT-PCR

Eén l van het primaire cDNA-synthesereactiemengsel werd genomen voor daaropvolgende PCR-amplificatie door RT-PCR of qRT-PCR. RT-PCR-runs werden routinematig uitgevoerd gedurende 20 tot 35 cycli, afhankelijk van het lineaire bereik van PCR-amplificatie voor elk gen. Elke PCR-cyclus omvatte incubaties bij 94°C gedurende 30 seconden, bij 55°C gedurende 1 minuut en bij 72°C gedurende 5 minuten. Na de laatste cyclus werd nog een cyclus van 10 min bij 72 °C uitgevoerd om het trimmen van onvolledige polymerisaties mogelijk te maken. Positieve en negatieve controlegenen werden in de reactiesets opgenomen om de haalbaarheid van de testomstandigheden te verzekeren. De gebruikte RT-PCR-primers staan ​​vermeld in aanvullend bestand 1.

qRT-PCR-reacties werden uitgevoerd in blokken met 96 putjes met een Applied Biosystems 7500 Real-Time PCR-systeem met behulp van de SYBR Green I-mastermix in een reactievolume van 25 l, dat 1 μl van het primaire cDNA-reactiemengsel bevat, 2X SYBR Groene PCR Master Mix (Applied Biosystems, Foster City, CA), en een primerpaar. De gebruikte primers worden vermeld in aanvullend bestand 1. Het gebruikte tweestaps thermische cyclusprofiel was 15 s bij 95°C en 1 minuut bij 60°C. Alle qRT-PCR-reacties werden uitgevoerd in biologische duplicaten, die elk werden gebruikt voor RNA-extractie gevolgd door qRT-PCR in drievoud. De laatste drempelcyclus (Ct) waarden waren het gemiddelde van zes waarden (biologische duplicaten, elk met drievoud). De vergelijkende ΔΔCt methode werd gebruikt om de relatieve hoeveelheden van elk geamplificeerd product in de monsters te evalueren. de Ct werd automatisch bepaald voor elke reactie door het Applied Biosystems 7500 Real-Time PCR-systeem ingesteld met standaardparameters. De specificiteit van de qRT-PCR-reacties werd bepaald door smeltcurveanalyse van de geamplificeerde producten met behulp van de standaardmethode die in het systeem is geïnstalleerd.

Elke qRT-PCR-reactieset bevatte een negatieve controle met water in plaats van cDNA. Dubbele metingen werden gemiddeld en de gemiddelde waarden werden gebruikt voor verdere berekeningen.

Bepaling van expressiestabiliteit van referentiegenen

De expressiestabiliteit van elk referentiegen werd geanalyseerd met behulp van de softwarepakketten geNorm (versie 3.5) en NormFinder (versie 0.953), die ook zijn geïntegreerd in de GenEx-software (versie 4.3.5, http://www.multid.se). De geNorm-software berekent de genexpressiestabiliteit (M) voor een referentiegen als de gemiddelde paarsgewijze variatie V voor het gen met alle andere geteste referentiegenen. Door stapsgewijze uitsluiting van het gen met de hoogste M-waarde kunnen de geteste genen worden gerangschikt volgens de stabiliteit van hun expressiepatronen.

De NormFinder-software is een algoritme voor het identificeren van het optimale normalisatiegen onder een reeks kandidaatgenen. Het rangschikt de set kandidaat-normalisatiegenen volgens de stabiliteit van hun expressiepatronen in een bepaalde steekproefset onder een bepaald experimenteel ontwerp. Daarom kan het expressiegegevens analyseren die zijn verkregen via elke kwantitatieve methode, zoals qRT-PCR en op microarray gebaseerde expressieprofilering. De laagste stabiliteitswaarde vertegenwoordigt de meest stabiele genexpressie binnen de onderzochte genenset.


Bekijk de video: DNA translatie, transcriptie en eiwitsynthese (December 2021).