如何将limma ebayes用于小样本？（R编程）

Question

我有以下来自测试和控制的表达数据，每个数据包含两个样本。

> test <- read.table("http://dpaste.com/1059889/plain/")
> control <-read.table("http://dpaste.com/1059891/plain/")
# in reality there are more  genes in each file.


> test
        V1         V2         V3
1   Gene_1 3694.11963 3591.95732
2   Gene_2  791.57558  753.04814
3   Gene_3 2751.34223 2562.46166
4   Gene_4 3068.07188 2651.62906
5   Gene_5  295.00476  247.78944
6   Gene_6 2737.22068 2824.85853
7   Gene_7 1274.54016 1196.54412
8   Gene_8 7011.31102 6703.59308
9   Gene_9  991.71137 1170.66072
10 Gene_10   67.83878   81.69033
11 Gene_11  139.96068  141.97499
12 Gene_12  337.40039  354.96356
13 Gene_13 2861.67548 3132.97426
14 Gene_14 1264.63942 1547.56872    

> control
        V1        V2        V3
1   Gene_1  98.76904 219.95533
2   Gene_2  64.13716 152.69867
3   Gene_3  84.54906 194.95583
4   Gene_4 106.64893 220.18668
5   Gene_5  50.30000  40.20000
6   Gene_6  24.22860  56.13421
7   Gene_7  43.08251  63.50765
8   Gene_8 408.95196 589.50150
9   Gene_9  37.68644  58.33591
10 Gene_10 100.33000 430.00000
11 Gene_11  23.24041  20.00000
12 Gene_12  17.64007  21.34300
13 Gene_13  65.45922  74.02418
14 Gene_14  43.69905  89.19588

我想使用 计算 P 值，以查看它们是否使用ebayes进行差异表达 。

使用标准 t.test 很简单，但我发现它对小样本没有用。

 t.test(c(test[2,]$V1,teste[2,]$V3),c(control[2,]$V1,control[2,]$V3))

有什么办法呢？从帮助文件中看不清楚。

Answer 1

有两个步骤——首先你创建一个 4x14 的矩阵——14 个（在你的情况下）基因的 4 个重复。还要创建一个实验设计：一个包含一列的数据框c(1, 1, 0, 0)，表示矩阵的前两列是测试，后两列是控制）。

library(limma)
m = as.matrix(cbind(test[, 2:3], control[, 2:3]))
rownames(m) = test[, 1]
d = data.frame(condition=c(1, 1, 0, 0))

然后使用lmFit拟合线性模型，eBayes调整它并计算 p 值：

fit = lmFit(m, d)
e = eBayes(fit)

类的e对象MArrayLM包含 p 值和估计系数，以及一些其他信息：

print(e$coefficients)
         condition
# Gene_1  3643.03848
# Gene_2   772.31186
# ...
print(e$p.value)
#            condition
# Gene_1  6.299507e-07
# Gene_2  1.333970e-04
names(e)
#  [1] "coefficients"     "rank"             "assign"           "qr"              
#  [5] "df.residual"      "sigma"            "cov.coefficients" "stdev.unscaled"  
#  [9] "pivot"            "genes"            "Amean"            "method"          
# [13] "design"           "df.prior"         "s2.prior"         "var.prior"       
# [17] "proportion"       "s2.post"          "t"                "df.total"        
# [21] "p.value"          "lods"             "F"                "F.p.value"       

顺便说一句，虽然limma由于多种原因而更好（它计算其他参数，提供更多选项并执行有用的校正），但如果您确实想执行 at.test以获取矩阵中每个基因的 p 值，您可以这样做：

pvals = apply(m, 1, function(r) t.test(r ~ d$condition)$p.value)