2016-12-07 79 views
2

我的數據集由3個時間點記錄的基因表達值構成 我正在嘗試應用tukey校正的anova測試來查找跨時間點基因的差異表達。因此,對於每一個基因,我想像的比較: 基因中的時間點1比2 基因中的時間點2比3 基因中的時間點3比1數據子集內的方差分析

我的數據是在以下格式:

> head(rf) 
       gene  expn timepoint rep 
2 EG620009 // EG620009 8.428851  x0hr 0 
3     LYPLA1 10.386500  x0hr 0 
21 EG620009 // EG620009 8.582346  x0hr 1 
31    LYPLA1 10.379710  x0hr 1 
22 EG620009 // EG620009 8.566248  x0hr 2 
32    LYPLA1 10.399080  x0hr 2 
    > tail(rf) 
       gene  expn timepoint rep 
23 EG620009 // EG620009 8.561409  x24hr 0 
33    LYPLA1 10.233400  x24hr 0 
24 EG620009 // EG620009 8.750639  x24hr 1 
34    LYPLA1 10.023780  x24hr 1 
25 EG620009 // EG620009 8.560267  x24hr 2 
35    LYPLA1 10.025980  x24hr 2 

如果我是這樣做:

TukeyHSD(aov(rf$expn ~ rf$timepoint * rf$gene)) 

這會給我所有的基因 即每一個時間點之間的比較。包括比較如 基因a時間點1對基因b時間點2

我一直在努力研究如何將aov函數應用於基因數據子集化。我已經定義了一個函數,它給出了p值作爲輸出,並嘗試使用下面的by函數分別將其應用於每個基因;

> gene.aov = function(x) {TukeyHSD(aov(expn ~ timepoint, data = x))} 
> aov.pval = function(y) {y$timepoint[,4]} 
> gene.pval = function(z) {aov.pval(gene.aov(z))} 
> pvals = by(rf$expn,list(rf$gene),gene.pval) 
> Error in eval(predvars, data, env) : 
    numeric 'envir' arg not of length one 

任何暗示爲什麼這不起作用?或者我應該以完全不同的方式來解決這個問題? 謝謝!

回答

1

它不工作,因爲by預計它的第一個參數是data.frame或矩陣,你傳遞的是rf$exp這是一個numeric向量。你可以做到這一點,它會正常工作(我放棄了多個功能的可讀性)。

by(rf, rf$gene, function(x) {TukeyHSD(aov(expn ~ timepoint, data = x))}, simplify = F) 

rf$gene: EG620009 
    Tukey multiple comparisons of means 
    95% family-wise confidence level 

Fit: aov(formula = expn ~ timepoint, data = x) 

$timepoint 
       diff  lwr  upr  p adj 
x24hr-x0hr 0.09829 -0.123391 0.319971 0.2857424 

--------------------------------------------------------------------------- 
rf$gene: LYPLA1 
    Tukey multiple comparisons of means 
    95% family-wise confidence level 

Fit: aov(formula = expn ~ timepoint, data = x) 

$timepoint 
       diff  lwr  upr  p adj 
x24hr-x0hr -0.2940433 -0.4876756 -0.100411 0.0135193 
+0

太好了,非常感謝您的幫助:) – Sarah