Frage:
Wie berechnet man RPKM in R?
Iakov Davydov
2017-05-17 15:28:32 UTC
view on stackexchange narkive permalink

Ich habe die folgenden Daten der Fragmentanzahl für jedes Gen in 16 Proben:

  > str (Ausdruck) 'data.frame': 42412 obs. von 16 Variablen: $ sample1: int 4555 49 122 351 53 27 1 0 0 2513 ... $ sample2: int 2991 51 55 94 49 10 55 0 0 978 ... $ sample3: int 3762 28 136 321 94 12 15 0 0 2181 ... $ sample4: int 4845 43 193 361 81 48 9 0 0 2883 ... $ sample5: int 2920 24 104 151 50 20 32 0 0 1743 ... $ sample6: int 4157 11 135 324 58 26 4 0 0 2364 ... $ sample7: int 3000 19 155 242 57 12 18 2 0 1946 ... $ sample8: int 5644 30 227 504 91 37 11 0 0 2988 ... $ sample9: int 2808 65 247 93 272 38 1108 1 0 1430 ... $ sample10: int 2458 37 163 64 150 29 729 2 1 1049 ... $ sample11: int 2064 30 123 51 142 23 637 0 0 1169 ... $ sample12: int 1945 63 209 40 171 41 688 3 2 749 ... $ sample13: int 2015 57 432 82 104 47 948 4 0 1171 ... $ sample14: int 2550 54 177 59 201 36 730 0 0 1474 ... $ sample15: int 2425 90 279 73 358 34 1052 3 3 1027 ... $ sample16: int 2343 56 365 67 161 43 877 3 1 1333 ...  

Wie berechne ich RPKM-Werte aus diesen?

Was haben Sie versucht, diese Frage zu lösen? Haben Sie Bioconductor besucht?
Vier antworten:
#1
+32
Konrad Rudolph
2017-05-17 15:46:42 UTC
view on stackexchange narkive permalink

Zunächst einmal

Verwenden Sie keine RPKMs .

Sie sind wirklich veraltet, weil sie einmal verwirrend sind es kommt zu Paired-End-Reads. Wenn überhaupt, verwenden Sie FPKMs , die mathematisch identisch sind, aber einen korrekteren Namen verwenden (zählen wir gepaarte Lesevorgänge separat? Nein, wir zählen Fragmente ). P. >

Noch besser, verwenden Sie TPM (= Transkripte pro Million) oder eine geeignete bibliotheksübergreifende Normalisierungsmethode. TMP ist definiert als:

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

wobei

  • $ \ color {orchid} i $: Transkriptionsindex,
  • $ x_i $: Transkript-Rohzahl,
  • $ \ color {tomate} j $ iteriert über alle (bekannten) Transkripte,
  • $ \ color {dodgerblue} {\ frac {x_k} {{l_ \ text {eff}} _ k}} $: Rate der Fragmentabdeckung pro Nukleobase ($ l_ \ text {eff} $ ist die effektive Länge),
  • $ \ color {darkcyan} {10 ^ 6} $: Skalierungsfaktor (= "pro Million").

Das heißt, FPKM kann in R wie folgt berechnet werden. Beachten Sie, dass der größte Teil der Berechnung im logarithmisch transformierten Zahlenraum erfolgt, um numerische Instabilität zu vermeiden:

  fpkm = Funktion (Anzahl , effektive_Längen) {exp (log (Anzahl) - log (effektive_Längen) - log (Summe (Anzahl)) + log (1E9))}  

Hier die effektive Länge ist die Transkriptlänge minus der mittleren Fragmentlänge plus 1; Das heißt, alle möglichen Positionen eines durchschnittlichen Fragments innerhalb des Transkripts, was der Anzahl aller unterschiedlichen Fragmente entspricht, die aus einem Transkript entnommen werden können.

Diese Funktion verarbeitet eine Bibliothek zu einer Zeit. Ich ( und andere) argumentiere, dass Funktionen so geschrieben werden sollten. Wenn Sie den Code auf mehrere Bibliotheken anwenden möchten, ist mit ‹dplyr›:

nichts einfacher
  tidy_expression = tidy_expression% >% group_by (Beispiel)% >% mutate (FPKM = fpkm (Anzahl, col_data $ Lengths))  

Die Daten in der Frage ist nicht im ordentlichen Datenformat, daher müssen wir es zuerst mit ‹tidyr› entsprechend transformieren:

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

Diese Gleichung schlägt fehl, wenn alle Ihre Zählungen Null sind. Anstelle von Nullen erhalten Sie einen Vektor von NaNs. Vielleicht möchten Sie das erklären.


Und ich erwähnte, dass TPMs überlegen sind, daher hier auch ihre Funktion:

  tpm = Funktion (Anzahl, effektive_Längen) {rate = log (Anzahl) - log (effektive_Längen) exp (rate - log (Summe (exp (Rate))) + log (1E6))}   vor>
Kann ich eine Art Meta-Frage stellen? Ich habe kürzlich angefangen, '%>%' im R-Code zu sehen und hatte es noch nie bemerkt. Was macht das genau?
@Greg [Es ist eine Pipe] (https://github.com/tidyverse/magrittr). Es gibt es schon seit geraumer Zeit (obwohl verschiedene Bibliotheken unterschiedliche Operatorsymbole verwendeten; ich selbst habe "% |%" verwendet, ähnlich wie bei der Shell-Pipe), aber erst vor kurzem hat sich der Mainstream durchgesetzt, hauptsächlich durch die dplyr-Bibliothek.
#2
+9
Iakov Davydov
2017-05-17 15:28:32 UTC
view on stackexchange narkive permalink

RPKM ist definiert als:

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

Wie Sie sehen können, müssen Sie dies tun haben Genlängen für jedes Gen.

Angenommen, geneLength ist ein Vektor mit der gleichen Anzahl von Zeilen wie Ihr data.frame und jedem Wert des Vektors entspricht einem Gen (Zeile) in Ausdruck .

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

In Bezug auf die numerische Stabilität

Dies wird in einer der Antworten vorgeschlagen, dass alle Berechnungen in einer logarithmisch transformierten Skala durchgeführt werden sollten. Meiner Meinung nach gibt es dafür absolut keinen Grund. IEEE binary64 speichert eine Zahl als Binärzahl 1.b_ {51} b {50} ... b_0 mal 2 ^ {e-1023}. Die relative Genauigkeit hängt nicht vom Exponentenwert ab, da eine Zahl im Bereich [~ 10 ^ -308; 10 ^ 308].

Bei RPKM können wir den Bereich nur verlassen, wenn die Gesamtzahl der Lesevorgänge bei 10 ^ 300 liegt, was überhaupt nicht realistisch ist.

Ein An der hellen Stelle schadet es auch nicht sehr, Berechnungen im Log-Maßstab durchzuführen. Im schlimmsten Fall verlieren Sie etwas an Präzision.

und dann kann die Auswahl der richtigen Genlänge selbst je nach Organismus ein Albtraum sein ... Vielleicht sollten Sie auch Strategien entwickeln, um die richtige auszuwählen.
@Mitra Ich denke, das ist etwas außerhalb des Rahmens. Aber wenn Sie Erfahrung damit haben, können Sie hier eine Antwort hinzufügen? Wäre echt toll!
@Mitra dachte noch einmal darüber nach, dies verdient eine separate Frage. Wenn es eine geben wird, werde ich sie mit dieser Antwort verknüpfen.
Ich denke nicht, dass es außerhalb des Rahmens liegt: Sie benötigen die Genlänge, um das RPKM zu berechnen. Andere Antworten haben dieses Problem umgangen, indem sie Transkript- oder Exonlängen verwendet haben, die leicht aus der Anmerkung abgerufen werden können. Dies ist bei der Genlänge, die Sie in Ihrer Formel verwenden, nicht der Fall.
#3
+1
arupgsh
2017-05-17 20:33:34 UTC
view on stackexchange narkive permalink

Wenn Sie eine differenzielle Expressionsanalyse planen, benötigen Sie wahrscheinlich keine RPKM-Berechnung.

RPK = Anzahl der zugeordneten Lesevorgänge / Transkriptlänge in kb (Transkriptlänge / 1000) )

RPKM = RPK / Gesamtzahl der Lesevorgänge in Millionen (Gesamtzahl der Lesevorgänge / 1000000)

Die gesamte Formel zusammen:

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

C = Anzahl der einem Gen zugeordneten Lesevorgänge

N = Gesamtzahl der zugeordneten Lesevorgänge in der Experiment

L = Exonlänge in Basenpaaren für ein Gen

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

Wenn Sie nach einer visuelleren Lösung suchen (zusätzlich zu den anderen Antworten), bietet NCI Genomic Data Commons (TCGA-Repository) eine schöne Formel:

enter image description here

Wobei:

RCg: Anzahl der dem Gen zugeordneten Lesevorgänge

RCpc: Anzahl der zugeordneten Lesevorgänge alle Protein-kodierenden Gene

RCg75: Der 75. Perzentil-Lesezählwert für Gene in der Probe

L: Länge des Gens in Basenpaaren

Wie andere bereits betont haben, haben FPKMs einige Probleme. GDC berechnet auch FPKM-UQ-Werte, die im oberen Quartil normalisiert sind. Diese werden für den Stichprobenvergleich und die Analyse der differentiellen Expression empfohlen.

Der Nenner für FPKM-UQ ist immer noch probenspezifisch, daher ist dies keine geeignete Normalisierung für die Analyse der differentiellen Expression (er liegt innerhalb der Probennormalisierung, nicht dazwischen).
GDC scheint etwas anderes zu implizieren: https://gdc.cancer.gov/about-data/data-harmonization-and-generation/genomic-data-harmonization/high-level-data-generation/rna-seq-quantification


Diese Fragen und Antworten wurden automatisch aus der englischen Sprache übersetzt.Der ursprüngliche Inhalt ist auf stackexchange verfügbar. Wir danken ihm für die cc by-sa 3.0-Lizenz, unter der er vertrieben wird.
Loading...