r-如何在Bonferroni校正的情况下计算数据帧中每行的超几何检验

Question

我正在计算R.中数据帧中每行的超几何测试，希望是正确的

其中列1是基因(微小RNA(的名称；Total_mRNAs"；是指基因组中总共存在多少信使核糖核酸，所以这不会改变。列"；Total_targets_targets"是如果所有的信使核糖核酸都存在，每个微小核糖核酸可以靶向多少信使核糖核酸。然而对于这个例子；mRNAs"亚基_；存在(这个数字也总是相同的(，并且在这些中，我知道每个微小RNA可以靶向多少mRNA；subset_targets"；。

为了确定每种微小RNA的靶标与背景相比是否富集(总mRNA和靶向它们的总微小RNA(，我正在每行进行超几何测试，如下所示：

phyper(targets-in-subset, targets-in-bkgd, failure-in-bkgd, sample-size-subset, lower.tail= FALSE)

dput(df1)
structure(list(Genes_names = c("microRNA-1", "microRNA-2", "microRNA-3", 
"microRNA-4", "microRNA-5", "microRNA-6", "microRNA-7", "microRNA-8", 
"microRNA-9", "microRNA-10"), Total_mRNAs = c(61064L, 61064L, 
61064L, 61064L, 61064L, 61064L, 61064L, 61064L, 61064L, 61064L
), Total_targets_targets = c(1918L, 7807L, 3969L, 771L, 2850L, 
1355L, 1560L, 2478L, 1560L, 2478L), subset_mRNAs = c(17571L, 
17571L, 17571L, 17571L, 17571L, 17571L, 17571L, 17571L, 17571L, 
17571L), subset_targets = c(544L, 2109L, 1137L, 213L, 793L, 394L, 
430L, 686L, 430L, 686L)), class = "data.frame", row.names = c(NA, 
-10L))

df1$pvalue <- phyper(df1$subset_targets, df1$Total_targets_targets, df1$Total_mRNAs-df1$Total_targets_targets, df1$subset_mRNAs, lower.tail= FALSE)

现在的问题是，我如何才能纠正这些值？这个计算在理论上正确吗？

Answer 1

如果您没有很多样本，请避免这种痛苦，只需使用fisher测试并使用p.adjust:进行bonferroni

library(broom)
result = lapply(1:nrow(df1),function(i){
not_target_subset = df1$Total_targets_targets[i] - df1$subset_targets[i]
not_subset = df1$Total_mRNAs[i] - df1$subset_mRNAs[i] - not_target_subset


M = cbind(c(df1$subset_targets[i],df1$subset_mRNAs[i]-df1$subset_targets[i]),
c(not_target_subset,not_subset))

res = data.frame(Genes_names=df1$Genes_names[i],
tidy(fisher.test(M,alternative="greater")))
return(res)
})
result= do.call(rbind,result)
result$padj = p.adjust(result$p.value,"bonferroni")

你的超几何代码有点偏离。请注意，你正在进行单侧超几何测试。

你可以查看这篇文章，了解如何将表格放入phyper中，以及为什么你需要-1。所以我们计算超几何p值：

result$hyper_p = with(df1, 
phyper(subset_targets-1,subset_mRNAs,Total_mRNAs-subset_mRNAs, Total_targets_targets, lower.tail= FALSE)
)

你可以看到它的计数：

Genes_names  estimate   p.value  conf.low conf.high
1   microRNA-1 0.9793710 0.6655527 0.8984025       Inf
2   microRNA-2 0.9047305 0.9998968 0.8647701       Inf
3   microRNA-3 0.9933480 0.5791759 0.9350214       Inf
4   microRNA-4 0.9441864 0.7722712 0.8229140       Inf
5   microRNA-5 0.9520878 0.8789562 0.8863785       Inf
6   microRNA-6 1.0151760 0.4119600 0.9168998       Inf
7   microRNA-7 0.9404619 0.8641420 0.8539585       Inf
8   microRNA-8 0.9454359 0.8942082 0.8756678       Inf
9   microRNA-9 0.9404619 0.8641420 0.8539585       Inf
10 microRNA-10 0.9454359 0.8942082 0.8756678       Inf
method alternative   hyper_p
1  Fisher's Exact Test for Count Data     greater 0.6655527
2  Fisher's Exact Test for Count Data     greater 0.9998968
3  Fisher's Exact Test for Count Data     greater 0.5791759
4  Fisher's Exact Test for Count Data     greater 0.7722712
5  Fisher's Exact Test for Count Data     greater 0.8789562
6  Fisher's Exact Test for Count Data     greater 0.4119600
7  Fisher's Exact Test for Count Data     greater 0.8641420
8  Fisher's Exact Test for Count Data     greater 0.8942082
9  Fisher's Exact Test for Count Data     greater 0.8641420
10 Fisher's Exact Test for Count Data     greater 0.8942082

Answer 2

警告：提出该问题的用户指出该答案中的计算可能是错误的。请参阅下面的评论。

根据您的编辑，似乎您正在寻找的是

df1$new.column <- apply(df1,
margin = 1,
function(row),
{
return(phyper(row$targets.1, row$targets, sum(row$targets.1, row$targets), row$subset, lower.tail= FALSE))
}

编辑

正如StupidWolf在评论中指出的，phyper是矢量化的。因此，您可以使用(从评论中复制(

with(df1, phyper(targets.1, targets, sum(targets.1, targets), subset, lower.tail= FALSE)

啊！