R：从每个基因的 t.test 中获取 p 值答案

【问题标题】：R: obtaining p-values from t.test for each geneR：从每个基因的 t.test 中获取 p 值
【发布时间】：2018-08-28 21:14:02
【问题描述】：

我正在使用 Bioconductor 套件（所有数据集）并尝试对每个基因进行 t.test。目标是查看两性之间的基因表达差异。我可以通过以下方式获得基本的 t.test：

> males <- exprs[, pData(ALL)$sex == "M"]
> females<-exprs[, pData(ALL)$sex == "F"]
> t.test(males, females)

但是当我尝试使用 apply 函数为每个基因提取 p 值时，命令永远不会结束，只会继续无限循环（我认为）。

pvals=apply(exprs,1,function(x) {t.test(x[males],x[females])$p.value})

这里是男性样本，有 12625 行（即探测 ID）。

> males
                                01005     01010     04006     04007     04008
1000_at                      7.597323  7.479445  7.384684  7.905312  7.065914
1001_at                      5.046194  4.932537  4.922627  4.844565  5.147762
1002_f_at                    3.900466  4.208155  4.206798  3.416923  3.945869
1003_s_at                    5.903856  6.169024  6.116890  5.687997  6.208061

【问题讨论】：

这可能有助于stackoverflow.com/a/52011158
请显示一些示例数据。基因是单独的行还是一列中的指标？ 它分崩离析 也没有帮助。发布实际错误或不良结果。
谢谢。我有表达数据（按行），这被标记为exprs，然后是男性和女性。我相信我的男性和女性已经在单独的 data.frames 中。我可以使用以下内容：'sapply(exprs[-1], function(x) {t.test(x[males[,1] == 1], x[females[,1] == 2])$p .value})' 还是我应该做一些data.frame？原始的 data.frame 已经有两性了（我把它们分开了）。它在pData(ALL)$sex。
@Oars 您对我之前就您之前提出的任何建议以及在进行差异基因表达分析时应考虑的相关问题提出的任何建议都非常抗拒。简而言之： 1. 不要使用 t 检验（或上一个问题中的 ANOVA）来寻找差异表达的基因。 2. 您正在处理微阵列数据，其中值是每个探针而不是每个基因的表达值。您首先需要总结每个基因的探针值（通常使用 Tukey 的稳健中值抛光来完成），然后使用例如limma 寻找差异表达的基因。
Maurits - 一如既往的感谢。我的导师将 exprs(ALL) 调用称为获取“基因表达”。我敢肯定你说得更准确（说真的）。无论出于何种原因，我们都希望我们运行一个循环（使用 apply），为每个循环提供来自 t.test 的 p 值......我不想说基因，但他称它们为基因，即 445_at、40419_at、等等。我真的希望你在教这门课！我希望他能向我们展示最合适的包，为什么使用它，如何应用它以及如何使用它。我正在经历折磨。

标签： r

【解决方案1】：

这里有一些东西可以帮助您入门。（冒着重复自己的风险；-) 请注意，这更像是一项统计/计算练习，而不是您真正应该做的事情；正如我在评论中所解释的，存在描述差异基因表达的复杂方法。相比之下，t 检验（或 ANOVA）是一种非常粗略的方法。

我们加载所有库和数据。
```
library(ALL)
data(ALL)
```
为了表征男性和女性个体之间平均探测强度的差异，我们执行两个样本的双边 t 检验并将结果存储在 list。
```
lst <- apply(exprs(ALL), 1, function(x)
    t.test(x[which(pData(ALL)$sex == "M")], x[which(pData(ALL)$sex == "F")]))
```

我们提取每个探针的 t 统计量、平均探针强度和 p 值的差异，并将结果存储在 data.frame 中。

df <- do.call(rbind, lapply(lst, function(x) c(
    statistic = unname(x$statistic),
    diff = unname(diff(x$estimate)),
    pval = unname(x$p.value))))

我们使用Benjamini and Hochberg 的 FDR 方法更正多个假设检验的 p 值。
```
df <- transform(df, padj = p.adjust(pval, method = "BH"))
```

我们检查 df 的前 10 行，从最小到最大调整后的 p 值。

head(df[order(df$padj), ], n = 10)
#        statistic       diff         pval         padj
#37583_at  18.935092 -1.7717178 1.710570e-36 2.159594e-32
#38355_at  20.542586 -4.9979077 6.129942e-32 3.869526e-28
#41214_at  21.494496 -4.3233221 3.937217e-31 1.656912e-27
#34477_at  14.469711 -1.1639971 2.606867e-28 8.227924e-25
#35885_at  14.417265 -1.4006757 5.806146e-28 1.466052e-24
#38446_at -14.357159  2.3848176 1.956173e-21 4.116115e-18
#38182_at  11.052181 -0.7151076 1.140089e-19 2.056232e-16
#40097_at   9.401626 -0.5798433 8.801566e-16 1.388997e-12
#36321_at   9.208492 -0.6499951 1.823511e-15 2.557981e-12
#31534_at   8.939350 -0.5113203 1.077008e-14 1.359723e-11

我们在火山图中显示结果

ggplot(df, aes(diff, -log10(padj))) +
    geom_point() +
    labs(x = "Difference in mean probe intensity", y = "Adjusted p-value")

【讨论】：

毛里求斯 - 我感激不尽！玩过代码后，我会写更多。

【解决方案2】：

感谢 Maurits，我能够使用他的代码来回答我的问题。我还开发了以下完成任务的片段（我实际上更喜欢 Maurits 的解决方案，但这是完成任务的另一种方法：

> exprs<-exprs(ALL)
> pval<-numeric()
> p.dat<-pData(ALL)$sex
> r.sims<-nrow(exprs)
> for(gene in 1:r.sims) { 
+ gexprs<-exprs[gene,]
+ g.data<-data.frame(gexprs,p.dat)
+ ttest<-t.test(gexprs[p.dat=="M"],gexprs[p.dat=="F"])
+ pval[gene]<-ttest$p.value
+ }

【讨论】：

【解决方案3】：

如果你被允许使用外部包，那么：

library(matrixTests)
row_t_welch(exprs[, pData(ALL)$sex == "M"], exprs[, pData(ALL)$sex == "F"])

这是假设基因是按行写的。

【讨论】：