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))}