Test auf Homogenität schief verteilter Werte?

[Curnen]

Hallo liebe R-Forum Leser,

Ich brauche eure Hilfe bei der Wahl des richtigen statistischen Tests und dessen korrekter Interpretation. Ausgangspunkt sind numerische Genexpressionsdaten von denen ich vier biologische Replikate habe:

Code: Alles auswählen

str(genes.fpkm)
'data.frame':	97212 obs. of  5 variables:
 $ WT.CKH_0  : num  1.53e+02 0.00 8.33e-08 1.80e+01 4.25e-07 ...
 $ WT.CKH_1  : num  1.68e+02 0.00 0.00 2.30e+01 2.87e-07 ...
 $ WT.CKH_2  : num  1.35e+02 0.00 0.00 2.37e+01 1.31e-08 ...
 $ WT.CKH_3  : num  1.64e+02 0.00 0.00 1.97e+01 1.05e-12 ...
 $ covered_by_enhancer: Factor w/ 2 levels "no significant interaction",..: 2 1 1 1 1 1 1 2 1 1 ...
 
levels(genes.fpkm$covered_by_enhancer)
[1] "no significant interaction" "interacting with enhancer" 

Für rund 10.000 dieser Werte kann ich aus deinem anderen Datensatz ableiten, dass die betreffenden Gene mit einem "Enhancer" interagieren sollten. Dieser tut, was der Name andeutet - er sollte deren Genexpression verstärken. Die Vorhersage wäre also, dass die Werte der interagierenden Gene größer sein sollten als jene, für die keine Interaktion festgestellt werden konnte und das möchte ich formal testen.

Ich habe dazu schon mal mit t.test, ks.test und wilcox.test herumgespielt, aber frage mich, ob die Tests aus folgenden Gründen dafür überhaupt geeignet sind:

• Die Verteilung der Werte ist schief, d.h. viele Gene sind nicht oder praktisch nicht exprimiert, einige wenige extreme Ausreißer. Ein Test auf abweichende Mittelwerte wie der t.test erscheint mir nicht ausreichend.
  Werteverteilung.png
• Der Stichprobenumfang ist für "no significant interaction" ungefähr um den Faktor 8 größer
• Ich habe die Replikate einfach mit rowMeans() auf einen Mittelwert pro Gen eingedampft und dann diese Werte als Stichprobe benutzt. Ist das legitim, oder sollte ich die Einzelwerte der Replikate nutzen (d.h. entweder die Tests für jedes Replikat einzeln machen oder alle Einzelwerte aufnehmen?)

Könnt ihr mir dazu vielleicht den ein oder anderen Tipp geben bzw. den richtigen Test empfehlen?

Vielen Dank schonmal
Matthias

[bigben]

Hallo Matthias,

was den t-Test angeht: Wenn Du aus einer Egal-Wie-Verteilung eine Stichprobe von 97000 ziehst, dann entstammt der Mittelwert dieser Stichprobe einer Normalverteilung. Das sagt der Zentrale Grenzwertsatz und deshalb ist der t-Test bei so großen Stichproben auch recht robust gegen Normalverteilungsverletzungen.

Bei solch riesigen Datensätzen stellt sich eher nicht die Frage nach der Signifikanz sondern nach dem Ausmaß des Unterschieds. Wenn Du einen t-Test rechnest, kannst Du dieses Maß sehr elegant als Differenz der Mittelwerte angeben. Beim Rangsummentest bietet sich kein ähnlich brauchbares Maß an.

Leider würdest Du mit dem t-Test, wie auch mit dem Rangsummentest, die Information der 4 Replikate tatsächlich "eindampfen" müssen und das wäre doch schade um die darin enthaltene Information. Deshalb könnte ich mir vorstellen, dass ein mixed-effects-model oder eine Messwiederholungs-ANOVA das Problem passender angehen könnten. Ich habe allerdings selbst keine Erfahrung mit so umfangreichen Datensätzen.

LG,
Bernhard

[EDi]

Ich hab die Daten/das Setup noch nicht ganz verstanden:

Was sind denn die 97212 Reihen? So wie dein Bild ausschaut, nimmst du diese als unabhängige Beobachten an - ist das auch so haltbar? Oder sind die irgendwie korreliert?

Innerhalb eine Replikats (=Spalten1-4) hast du beide Faktoren drin?

Was sind das für Werte? Scheinen alle >0 zu sein? (oder ist die 0 auch dabei?).


Für mich hört sich das eher nach einem GLMM an.
Das "G" für Generalized um eine anderer Verteilung als die normale zu modellieren (Gamma? Tweedie?)
Das "LM" ist klar
Das "M" für mixed um die hierarchische Struktur (Gene innerhalb eines replikats?) zui brücksichtigen.

Sind nur einpaar Gedanken, die noch nicht fertig gekocht sind....

[Curnen]

Was sind denn die 97212 Reihen? So wie dein Bild ausschaut, nimmst du diese als unabhängige Beobachten an - ist das auch so haltbar? Oder sind die irgendwie korreliert?

Das sind Genexpressionswerte in der Einheit FPKM. Eine Zelle produziert Arbeitskopien von Genen, die gebraucht werden und deren Häufigkeit kann in einem Verfahren namens RNA-seq gemessen werden. Die Idee ist, dass eine größere Zahl von Arbeitskopien dafür spricht, dass diese Gene für die Zelle wirklich entscheidend sind. (Wobei das nicht immer zutreffend ist, wie in Einzelfällen bereits eindeutig nachgewiesen werden konnte)

Die Replikate sind weitgehend unabhängig (urpsrünglich aus den gleichen Zellen erzeugt, aber als eigene Krebszelllinie etabliert und kultiviert), sodass hier eine Korrelation dafür spräche, dass die Genexpression eine biologische Bedeutung hat.

Innerhalb eines Replikats sind die ~97.000 Werte nicht unabhängig, denn Gene sind in Signalkaskaden organisiert, d.h. die Expression eines Gens kann die Expression weiterer Gene anstoßen oder abschalten. Leider sind nicht alle Signalkaskaden bekannt und mitunter im Krebsfall auch von gesunden Zellen abweichend. Die Werte innerhalb eines Replikats als unabhängig anzusehen ist daher leider trotz des Wissens um ein eigentlich gegenteiliges Funktionsprizip unumgänglich.

Innerhalb eine Replikats (=Spalten1-4) hast du beide Faktoren drin?

Nein, es handelt sich dabei ebenfalls um eine sträflich grobe Vereinfachung der tatsächlichen biologischen Komplexität. Bei den Replikaten handelt es sich nämlich noch nicht einmal um die Werte einer einzelnen Zelle sondern um den Durchschnitt durch je eine Population von rund einer Million heterogener Krebszellen. Letztlich könnte es also sein, dass ein mittelmäßig exprimiertes Gen in 10% der Zellen extrem stark und in den anderen 90% gar nicht exprimiert wird.

Genauso verhält es sich mit der Interaktion mit dem Enhancer. Diese Interaktion kann nicht direkt gemessen werden sondern nur aus Daten gemutmaßt werden, die die dreidimensionale Struktur des Genoms (in Chromosomen) berücksichtigt. Aus der Kenntnis wo sich ein Gen und wo sich ein Enhancer befindet kann man darauf schließen, ob es plausibel ist, dass diese sich berühren könnten (so wie es unplausibel ist, dass du dich mit deinem großen Zeh selbst zwischen den Schulterblättern am Rücken kratzen kannst ). Aber nur weil diese Interaktion theoretisch möglich wäre, heißt das natürlich noch nicht, dass sie tatsächlich stattfindet - als Indiz dafür, dass die überwiegende Zahl der vorhergesagten Interaktionen zumindest phasenweise stattfindet, möchte dich den formalen Test machen und sehen, ob es eine "unerklärliche" Erhöhung der Genexpressionswerte gibt, die möglicherweise auf die Interaktion mit Enhancern zurückzuführen ist.

Was sind das für Werte? Scheinen alle >0 zu sein? (oder ist die 0 auch dabei?).

Es gibt auch einige tausend Werte, die exakt 0 sind sowie sehr viele, die nahe 0 sind. Ich habe, weil ich die Log-Skalierung brauchte diese also ausgeschlossen. Ich habe auch einige wenige Ausreißer ganz oben ausgeschlossen, weil es sich dabei sehr sicher um technische Artefakte/Messfehler handelt.

Da gerade der t.test sehr sensitiv je nach ausgeschlossenen Werten mal die Alternativhypothese angenommen oder abgeleht hat, bin ich misstrauisch geworden, ob der so genutzt werden kann.

Für mich hört sich das eher nach einem GLMM an.
Das "G" für Generalized um eine anderer Verteilung als die normale zu modellieren (Gamma? Tweedie?)
Das "LM" ist klar
Das "M" für mixed um die hierarchische Struktur (Gene innerhalb eines replikats?) zui brücksichtigen.

Danke Bernhard und EDi für die Tipps - werde mich mal mit Messwiederholungs-ANOVA und GLMMs beschäftigen. Egal ob fertig gekocht oder nicht - ein paar Stichworte oder Ideen zum weitersuchen sind immer hilfreich. DANKE!

PS: Das Paper mit den Methylierungswerten, wofür ihr mir vor knapp zwei Jahren die GAM-Formeln für ausgearbeitet habt, ist leider im Papierkorb gelandet. Die Hauptautorin hat im Streit das Labor verlassen und das Manuskript wurde nie eingereicht...wobei auch das ganze Projekt schlecht konzipiert war. Das aktuelle ist nicht besser

[EDi]

Ich bin in dem Genetik-Kram leider gar nicht drin, aber ich denke dass die Nullen auch Informationen enthalten können (würed ich also per se nicht ausschließen). Schau dir mal ein Tweedie-GLM an für kontinuierliche Daten >= 0 zu modellieren. Das dürfte auch andere deiner Kollegen deines Faches genutzt habe (?).

PS: Das Paper mit den Methylierungswerten, wofür ihr mir vor knapp zwei Jahren die GAM-Formeln für ausgearbeitet habt, ist leider im Papierkorb gelandet. Die Hauptautorin hat im Streit das Labor verlassen und das Manuskript wurde nie eingereicht...wobei auch das ganze Projekt schlecht konzipiert war. Das aktuelle ist nicht besser

Ich erinnere mich, war mit dem fancy GAM Schade