Domanda:
Come calcolare RPKM in R?
Iakov Davydov
2017-05-17 15:28:32 UTC
view on stackexchange narkive permalink

Ho i seguenti dati sui conteggi dei frammenti per ogni gene in 16 campioni:

  > str (espressione) 'data.frame': 42412 obs. di 16 variabili: $ 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 43877 3 1 1333 ...  

Come posso calcolare i valori RPKM da questi?

Cosa hai provato a risolvere questa domanda? Hai visitato Bioconductor?
Quattro risposte:
#1
+32
Konrad Rudolph
2017-05-17 15:46:42 UTC
view on stackexchange narkive permalink

Prima di tutto,

Non utilizzare RPKM .

Sono veramente deprecati perché una volta creano confusione si tratta di letture paired-end. Se non altro, usa FPKM , che matematicamente sono uguali ma usa un nome più corretto (contiamo le letture accoppiate separatamente? No, contiamo frammenti ).

Ancora meglio, utilizza TPM (= trascrizioni per milione) o un metodo di normalizzazione tra librerie appropriato. TMP è definito come:

$$ \ text {TPM} _ \ color {orchid} i = {\ color {dodgerblue} {\ frac {x_ \ color {orchidea} i} {{l_ \ text {eff}} _ \ color {orchidea} i}}} \ cdot \ frac {1} {\ sum_ \ color {pomodoro} j \ color {dodgerblue} {\ frac {x_ \ color {pomodoro} j} {{l_ \ text {eff}} _ \ color {tomato} j}}} \ cdot \ color {darkcyan} {10 ^ 6} $$

dove

  • $ \ color {orchid} i $: indice della trascrizione,
  • $ x_i $: conteggio grezzo della trascrizione,
  • $ \ color {pomodoro} j $ itera su tutte le trascrizioni (note),
  • $ \ color {dodgerblue} {\ frac {x_k} {{l_ \ text {eff}} _ k}} $: tasso di copertura dei frammenti per nucleobase ($ l_ \ text {eff} $ è il lunghezza effettiva),
  • $ \ color {darkcyan} {10 ^ 6} $: fattore di scala (= "per milioni").

Detto questo, FPKM può essere calcolato in R come segue. Tieni presente che la maggior parte del calcolo avviene nello spazio numerico trasformato nel log, per evitare l ' instabilità numerica:

  fpkm = function (counts , lunghezza_effettiva) {exp (log (conteggi) - log (lunghezze_effettiva) - log (somma (conteggi)) + log (1E9))}  

Qui, la lunghezza effettiva è la lunghezza della trascrizione meno la lunghezza media del frammento più 1; cioè tutte le possibili posizioni di un frammento medio all'interno della trascrizione, che è uguale al numero di tutti i frammenti distinti che possono essere campionati da una trascrizione.

Questa funzione gestisce una libreria Al tempo. Io ( e altri) sostengo che questo è il modo in cui le funzioni dovrebbero essere scritte. Se vuoi applicare il codice a più librerie, niente è più facile usando ‹dplyr›:

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

Tuttavia, i dati nella domanda non è in un formato dati ordinato, quindi dobbiamo prima trasformarlo di conseguenza utilizzando ‹tidyr›:

  tidy_expression = espressione% >% gather (Sample, Count) 

Questa equazione fallisce se tutti i tuoi conteggi sono zero; invece degli zeri otterrai un vettore di NaN. Potresti tenerne conto.


E ho detto che i TPM sono superiori, quindi ecco anche la loro funzione:

  tpm = function (counts, actual_lengths) {rate = log (counts) - log (actual_lengths) exp (rate - log (sum (exp (rate))) + log (1E6))}  
Posso fare una specie di meta domanda? Ho iniziato a vedere '%>%' nel codice R di recente e non l'avevo mai notato prima. Cosa fa esattamente?
@Greg [È una pipa] (https://github.com/tidyverse/magrittr). È in circolazione da un bel po 'di tempo (sebbene diverse librerie usassero simboli di operatori diversi; io stesso usavo `% |%`, simile allo shell pipe) ma solo di recente ha guadagnato la trazione principale, principalmente attraverso la libreria dplyr.
#2
+9
Iakov Davydov
2017-05-17 15:28:32 UTC
view on stackexchange narkive permalink

RPKM è definito come:

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

Come puoi vedere, devi hanno lunghezze di gene per ogni gene.

Supponiamo che geneLength sia un vettore che ha lo stesso numero di righe di data.frame e ogni valore del vettore corrisponde a un gene (riga) in expression.

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

Riguardo alla stabilità numerica

È suggerito in una delle risposte, che tutti i calcoli dovrebbero essere eseguiti in una scala trasformata in logaritmo. Secondo me non c'è assolutamente alcun motivo per farlo. IEEE binary64 memorizza un numero come numero binario 1.b_ {51} b {50} ... b_0 volte 2 ^ {e-1023}. La precisione relativa non dipende dal valore dell'esponente dato che un numero è compreso nell'intervallo [~ 10 ^ -308; 10 ^ 308].

In caso di RPKM possiamo uscire dall'intervallo solo se il numero totale di letture è intorno a 10 ^ 300, il che non è affatto realistico.

nel sito luminoso non c'è molto danno nemmeno nel fare calcoli nella scala logaritmica. Nel peggiore dei casi, perdi un po 'di precisione.

e poi scegliere la giusta lunghezza del gene stesso può essere un incubo a seconda dell'organismo ... Quindi forse dovresti anche sviluppare quelle che potrebbero essere le strategie per scegliere quello corretto.
@Mitra Penso che questo sia leggermente fuori portata. Ma se hai esperienza nel farlo, potresti aggiungere una risposta qui? Sarebbe davvero fantastico!
@Mitra ci ha pensato una seconda volta, questo merita una domanda a parte. Se ce ne sarà uno, lo collegherò da questa risposta.
Non credo sia fuori portata: hai bisogno della lunghezza del gene per calcolare l'RPKM. Altre risposte hanno aggirato questo problema utilizzando la lunghezza della trascrizione o dell'esone che possono essere facilmente recuperate dall'annotazione. Questo non è il caso della lunghezza del gene che stai usando nella tua formula.
#3
+1
arupgsh
2017-05-17 20:33:34 UTC
view on stackexchange narkive permalink

Se hai intenzione di eseguire un'analisi delle espressioni differenziali, probabilmente non avrai bisogno del calcolo RPKM.

RPK = Numero di letture mappate / lunghezza della trascrizione in kb (lunghezza della trascrizione / 1000 )

RPKM = RPK / numero totale di letture in milioni (numero totale di letture / 1000000)

L'intera formula insieme:

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

C = Numero di letture mappate a un gene

N = Letture mappate totali nel esperimento

L = lunghezza dell'esone in coppie di basi per un gene

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

Se stai cercando una soluzione più visiva (oltre alle altre risposte), NCI Genomic Data Commons (repository TCGA) offre una bella formula:

enter image description here

Dove:

RCg: numero di letture mappate al gene

RCpc: numero di letture mappate a tutti i geni codificanti proteine ​​

RCg75: il 75 ° percentile del valore del conteggio delle letture per i geni nel campione

L: lunghezza del gene in coppie di basi

Come altri hanno sottolineato, gli FPKM hanno alcuni problemi. GDC calcola anche i valori FPKM-UQ normalizzati per il quartile superiore. Quelli sono consigliati per il confronto tra campioni e l'analisi dell'espressione differenziale.

Il denominatore per FPKM-UQ è ancora specifico del campione, quindi questa non sarà una normalizzazione appropriata per l'analisi dell'espressione differenziale (è all'interno della normalizzazione del campione, non tra).
GDC sembra implicare il contrario: https://gdc.cancer.gov/about-data/data-harmonization-and-generation/genomic-data-harmonization/high-level-data-generation/rna-seq-quantification


Questa domanda e risposta è stata tradotta automaticamente dalla lingua inglese. Il contenuto originale è disponibile su stackexchange, che ringraziamo per la licenza cc by-sa 3.0 con cui è distribuito.
Loading...