Leider tritt dabei (wieder einmal) das Problem auf, dass ich damit Werkzeuge verwenden muss, deren Komplexität mein Statistikwissen um ein Vielfaches übersteigt. Ich fürchte, ich brauche Hilfe bei der Spezifikation der Modellformel und auch der Interpretation der Ergebnisse.
Mein Datensatz ist im Prinzip so aufgebaut:
Code: Alles auswählen
mydata <- data.frame(
transcript_ID=as.factor(paste("transcript",rep(c(1:50),each=10))),
peak_status=as.factor(rep(c("none","Regular H3K4me3 peak","H3K4me3 buffer domain"),times=c(160,270,70))), # echtes Verhältnis 106:315:1
sample_ID=as.factor(rep(paste("sample",LETTERS[1:10]),times=50)),
genotype=as.factor(rep(rep(c("wt","cchip"),each=5),times=50)),
sample_expression=unlist(lapply(mapply(c,as.list(c(rnorm(43,15,3),rnorm(7,100,20))),as.list(rnorm(50,6,2)),SIMPLIFY=FALSE),function(x){rnorm(10,x[1],x[2])}))
)
Wie muss denn meine Modellformel und der PostHoc-Test ausehen, wenn ich auf Assozation mit genotype und peak_status testen will?
Code: Alles auswählen
library("nlme")
library("multcomp")
fit <- lme(sample_expression ~ genotype + peak_status, data=mydata,random=list(transcript_ID= ~1,sample_ID= ~ 1))
anova(fit)
summary(glht(fit,linfct=mcp(peak_status="Tukey")))
Code: Alles auswählen
# einer von vielen Google-Funden, wo im Random-Teil kuriose Dinge gemacht wurden...
fit <- lme(sample_expression ~ genotype+peak_status, data=mydata,random = ~1|sample_ID/transcript_ID)
Code: Alles auswählen
aggregate(sample_expression ~ genotype + peak_status, data=plotdata_buffered_genes, FUN=median)
genotype peak_status sample_expression
1 wt none 6.388254
2 chip none 5.344367
3 wt Regular H3K4me3 peak 12.299607
4 chip Regular H3K4me3 peak 12.488074
5 wt H3K4me3 buffer domain 147.485491
6 chip H3K4me3 buffer domain 162.046081
anova(fit)
numDF denDF F-value p-value
(Intercept) 1 85542 408.9787 <.0001
genotype 1 85542 12.5218 4e-04
peak_status 2 85542 15.9353 <.0001
summary(glht(fit,linfct=mcp(peak_status="Tukey",genotype = "Tukey")))
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
peak_status: Regular H3K4me3 peak - none == 0 -1.9040 1.9911 -0.956 0.716415
peak_status: H3K4me3 buffer domain - none == 0 75.0342 13.6525 5.496 < 1e-04 ***
peak_status: H3K4me3 buffer domain - Regular H3K4me3 peak == 0 76.9382 13.6775 5.625 < 1e-04 ***
genotype: Dnmt1 -/chip - Dnmt1 +/+ == 0 2.1140 0.5666 3.731 0.000582 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
Matthias