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

【问题标题】：How to use limma ebayes for small samples? (R programming)如何将 limma ebayes 用于小样本？（R 编程）
【发布时间】：2013-04-16 11:34:39
【问题描述】：

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

> 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))

解决方法是什么？从帮助文件中看不清楚。

【问题讨论】：

使用经验贝叶斯对一个基因没有意义——你必须有多个才能发挥作用。为你的例子模拟矩阵 test_gene 和 control_gene 有意义吗？
@DavidRobinson：谢谢戴夫。我用多个基因更新了帖子。
只是好奇-我的回答中是否有我没有涵盖的内容，您希望我添加？

标签： r bioconductor

【解决方案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 更好的原因有很多（它计算其他参数、提供更多选项并执行有用的更正），但如果您确实想要执行 t.test 来获得 p - 矩阵中每个基因的值，你可以这样做：

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

【讨论】：