Vraag:
Hoe RPKM in R te berekenen?
Iakov Davydov
2017-05-17 15:28:32 UTC
view on stackexchange narkive permalink

Ik heb de volgende gegevens van fragmentaantallen voor elk gen in 16 monsters:

  > str (uitdrukking) 'data.frame': 42412 obs. van 16 variabelen: $ sample1: int 4555 49122351 53 27 1 0 0 2513 ... $ sample2: int 2991 51 55 94 49 10 55 0 0978 ... $ sample3: int 3762 28136321 94 12 15 0 0 2181 ... $ sample4: int 4845 43193361 81 48 9 0 0 2883 ... $ sample5: int 2920 24104151 50 20 32 0 0 1743 ... $ sample6: int 4157 11135324 58 26 4 0 0 2364 ... $ sample7: int 3000 19155242 57 12 18 2 0 1946 ... $ sample8: int 5644 30227504 91 37 11 0 0 2988 ... $ sample9: int 2808 65247 93272 38 1108 1 0 1430 ... $ sample10: int 2458 37163 64150 29729 2 1 1049 ... $ sample11: int 2064 30123 51142 23637 0 0 1169 ... $ sample12: int 1945 63209 40171 41688 3 2749 ... $ sample13: int 2015 57432 82104 47948 4 0 1171 ... $ sample14: int 2550 54177 59201 36730 0 0 1474 ... $ sample15: int 2425 90279 73 358 34 1052 3 3 1027 ... $ sample16: int 2343 56365 67161 43 877 3 1 1333 ...  

Hoe bereken ik RPKM-waarden hieruit?

Wat heb je geprobeerd om deze vraag op te lossen? Heb je Bioconductor bezocht?
Vier antwoorden:
#1
+32
Konrad Rudolph
2017-05-17 15:46:42 UTC
view on stackexchange narkive permalink

Ten eerste:

Gebruik geen RPKM's .

Ze zijn echt verouderd omdat ze eenmaal verwarrend zijn het gaat om gepaarde eindes. Gebruik in ieder geval FPKM s, die wiskundig hetzelfde zijn, maar een correctere naam gebruiken (tellen we gepaarde lezingen afzonderlijk? Nee, we tellen fragmenten ).

Nog beter: gebruik TPM (= transcripties per miljoen) of een geschikte normalisatiemethode tussen bibliotheken. TMP wordt gedefinieerd als:

$$ \ text {TPM} _ \ color {orchid} i = {\ color {dodgerblue} {\ frac {x_ \ color {orchid} i} {{l_ \ text {eff}} _ \ kleur {orchidee} i}}} \ cdot \ frac {1} {\ sum_ \ kleur {tomaat} j \ kleur {dodgerblauw} {\ frac {x_ \ kleur {tomaat} j} {{l_ \ text {eff}} _ \ color {tomaat} j}}} \ cdot \ color {darkcyan} {10 ^ 6} $$

waar

  • $ \ color {orchid} i $: transcript index,
  • $ x_i $: transcript raw count,
  • $ \ color {tomaat} j $ herhaalt alle (bekende) transcripten,
  • $ \ color {dodgerblue} {\ frac {x_k} {{l_ \ text {eff}} _ k}} $: snelheid van fragmentdekking per nucleobase ($ l_ \ text {eff} $ is de effectieve lengte),
  • $ \ color {darkcyan} {10 ^ 6} $: schaalfactor (= "per miljoen").

Dat gezegd hebbende, kan FPKM als volgt worden berekend in R. Merk op dat het grootste deel van de berekening plaatsvindt in log-getransformeerde getallenruimte, om numerieke instabiliteit te voorkomen:

  fpkm = function (counts , effectieve_lengtes) {exp (log (counts) - log (effectieve_lengtes) - log (som (counts)) + log (1E9))}  

Hier, de effectieve lengte is de lengte van het transcript minus de gemiddelde fragmentlengte plus 1; dat wil zeggen, alle mogelijke posities van een gemiddeld fragment in het transcript, wat gelijk is aan het aantal verschillende fragmenten dat kan worden gesampled uit een transcript.

Deze functie behandelt één bibliotheek tegelijk. Ik ( en anderen) beweren dat dit de manier is waarop functies moeten worden geschreven. Als u de code op meerdere bibliotheken wilt toepassen, is niets eenvoudiger met ‹dplyr›:

  tidy_expression = tidy_expression% >% group_by (Sample)% >% mutate (FPKM = fpkm (Count, col_data $ Lengths))  

De gegevens in de vraag heeft geen netjes gegevensformaat, dus we moeten het eerst dienovereenkomstig transformeren met ‹tidyr›:

  tidy_expression = expression% >% verzamelen (Sample, Count) 

Deze vergelijking mislukt als al je tellingen nul zijn; in plaats van nullen krijg je een vector van NaNs. Misschien wil je daar rekening mee houden.


En ik zei dat TPM's superieur zijn, dus hier is ook hun functie:

  tpm = function (counts, effectieve_lengtes) {rate = log (counts) - log (effectieve_lengtes) exp (rate - log (som (exp (rate))) + log (1E6))}  
Kan ik een soort metavraag stellen? Ik ben onlangs '%>%' in R-code gaan zien en had het nog nooit eerder opgemerkt. Wat doet dat precies?
@Greg [Het is een pijp] (https://github.com/tidyverse/magrittr). Het bestaat al een behoorlijk lange tijd (hoewel verschillende bibliotheken verschillende operatorsymbolen gebruikten; ik gebruikte zelf `% |%`, vergelijkbaar met de shell pipe), maar kreeg pas sinds kort mainstream tractie, voornamelijk via de dplyr-bibliotheek.
#2
+9
Iakov Davydov
2017-05-17 15:28:32 UTC
view on stackexchange narkive permalink

RPKM wordt gedefinieerd als:

RPKM = numberOfReads / (geneLength / 1000 * totalNumReads / 1,000,000)

Zoals je kunt zien, moet je hebben genlengtes voor elk gen.

Laten we zeggen dat geneLength een vector is die hetzelfde aantal rijen heeft als uw data.frame , en elke waarde van de vector komt overeen met een gen (rij) in expressie.

  expression.rpkm <- data.frame (sapply (expressie, functie (kolom) 10 ^ 9 * column / geneLength / sum (column)))  

Met betrekking tot numerieke stabiliteit

Het wordt voorgesteld in een van de antwoorden, dat alle berekeningen moeten worden uitgevoerd in een log-getransformeerde schaal. Naar mijn mening is daar absoluut geen reden voor. IEEE binary64 slaat een getal op als binair getal 1.b_ {51} b {50} ... b_0 maal 2 ^ {e-1023}. De relatieve precisie is niet afhankelijk van de exponentwaarde, aangezien een getal in het bereik [~ 10 ^ -308; 10 ^ 308].

In het geval van RPKM kunnen we alleen buiten het bereik komen als het totale aantal leesbewerkingen ongeveer 10 ^ 300 is, wat helemaal niet realistisch is.

Aan op de heldere site is er ook niet veel kwaad in het uitvoeren van berekeningen in de log-schaal. In het ergste geval verlies je een beetje precisie.

en dan kan het kiezen van de juiste genlengte zelf een nachtmerrie zijn, afhankelijk van het organisme ... Dus misschien moet je ook de strategieën ontwikkelen om de juiste te kiezen.
@Mitra Ik denk dat dit enigszins buiten het bereik valt. Maar als u daar ervaring mee heeft, kunt u hier dan een antwoord toevoegen? Zou echt geweldig zijn!
@Mitra als ik er nog over nadenk, verdient dit een aparte vraag. Als er een zal zijn, zal ik het koppelen vanuit dit antwoord.
Ik denk niet dat het buiten de scope valt: je hebt de genlengte nodig om de RPKM te berekenen. Andere antwoorden omzeilden dit probleem door gebruik te maken van transcript of exon-lengte die gemakkelijk uit de annotatie kunnen worden opgehaald. Dat is niet het geval met de genlengte die u in uw formule gebruikt.
#3
+1
arupgsh
2017-05-17 20:33:34 UTC
view on stackexchange narkive permalink

Als u van plan bent om een ​​differentiële expressie-analyse uit te voeren, heeft u de RPKM-berekening waarschijnlijk niet nodig.

RPK = Aantal toegewezen leesbewerkingen / lengte van transcriptie in kb (lengte transcript / 1000 )

RPKM = RPK / totaal aantal lezingen in miljoen (totaal aantal lezingen / 1000000)

De hele formule samen:

RPKM = (10 ^ 9 * C) / (N * L) Waar,

C = Aantal leesbewerkingen toegewezen aan een gen

N = Totaal toegewezen leesbewerkingen in de experiment

L = lengte van exon in basenparen voor een gen

#4
  0
burger
2017-05-17 20:16:11 UTC
view on stackexchange narkive permalink

Als u op zoek bent naar een meer visuele oplossing (naast de andere antwoorden), biedt NCI Genomic Data Commons (TCGA-repository) een mooie formule:

enter image description here

Waar:

RCg: aantal leesbewerkingen toegewezen aan het gen

RCpc: aantal leesbewerkingen toegewezen aan alle eiwitcoderende genen

RCg75: de 75e percentiel afgelezen telwaarde voor genen in het monster

L: lengte van het gen in basenparen

Zoals anderen hebben opgemerkt, hebben FPKM's enkele problemen. GDC berekent ook FPKM-UQ-waarden die genormaliseerd zijn in het bovenste kwartiel. Die worden aanbevolen voor vergelijking tussen steekproeven en differentiële expressieanalyse.

De noemer voor FPKM-UQ is nog steeds steekproefspecifiek, dus dit zal geen geschikte normalisatie zijn voor differentiële expressie-analyse (deze valt binnen de steekproefnormalisatie, niet tussen).
GDC lijkt anders te suggereren: https://gdc.cancer.gov/about-data/data-harmonization-and-generation/genomic-data-harmonization/high-level-data-generation/rna-seq-quantification


Deze Q&A is automatisch vertaald vanuit de Engelse taal.De originele inhoud is beschikbaar op stackexchange, waarvoor we bedanken voor de cc by-sa 3.0-licentie waaronder het wordt gedistribueerd.
Loading...