在 Wilcoxon Rank 之后在 R 中使用 Benjamini-Hochberg 错误发现率时出错

Question

我进行了 Wilcoxon 秩和检验，看 598019 个基因在三个疾病样本和三个对照样本之间的表达是否有显着差异。我在R.

当我看到有多少基因的 p 值 < 0.05 时，我总共得到 41913 个。我将Wilcoxon的参数设置如下；

wilcox.test(currRow[4:6], currRow[1:3], paired=F, alternative="two.sided", exact=F, correct=F)$p.value

（这个是在一个apply函数里面，如果有必要我可以提供我的全部代码，我有点不确定alternative="two.sided"是否正确）。

但是，由于我假设使用 Benjamini Hochberg 错误发现率校正多重比较会降低这个数字，然后我通过以下代码调整了 p 值 pvaluesadjust1 <- p.adjust(pvalues_genes, method="BH")

通过以下代码重新评估哪些 p 值小于 0.05，我得到 0！

p_thresh1 <- 0.05
names(pvaluesadjust1) <- rownames(gene_analysis1)
output <- names(pvaluesadjust1)[pvaluesadjust1 < p_thresh1]
length(output)

如果有人能解释一下，或者将我引导到可以帮助我理解发生了什么的地方，我将不胜感激！

谢谢（作为一个额外的问题，由于数据的大小，t 检验是否可以，Anderson-Darling 检验表明基础数据不正常。使用此统计检验，我的基因远少于 0.05 而不是比 Wilcoxon（2000 年左右）。

Answer 1

Wilcoxon 是一种基于等级的参数检验。如果您只有 6 个样本，您可以获得的最佳结果是疾病排名 2、2、2 与对照排名 5、5、5，反之亦然。

例如，在下面的这些随机值上尝试您在测试中使用的参数，您会得到相同的 p.value 0.02534732。

wilcox.test(c(100,100,100),c(1,1,1),exact=F, correct=F)$p.value
wilcox.test(c(5,5,5),c(15,15,15),exact=F, correct=F)$p.value

所以是的，使用 598019 可以得到 41913 < 0.05，这些 p 值不够低并且通过 FDR 调整，none 将永远通过。

您使用了错误的测试。要回答您的第二个问题，t.test 效果不佳，因为您没有足够的样本来正确估计标准偏差。下面我给大家展示一个使用DESeq2寻找差异基因的例子

library(zebrafishRNASeq)
data(zfGenes)
# remove spikeins
zfGenes = zfGenes[-grep("^ERCC", rownames(zfGenes)),]
head(zfGenes)
                   Ctl1 Ctl3 Ctl5 Trt9 Trt11 Trt13
ENSDARG00000000001  304  129  339  102    16   617
ENSDARG00000000002  605  637  406   82   230  1245

前三个是对照，后三个是处理，就像你的数据集一样。为了验证我之前所说的，你可以看到如果你做一个wilcoxon.test，最小值是0.02534732

all_pvalues = apply(zfGenes,1,function(i){
wilcox.test(i[1:3],i[4:6],exact=F, correct=F)$p.value
})
min(all_pvalues,na.rm=T)
# returns 0.02534732

所以我们继续进行 DESeq2

library(DESeq2)
#create a data.frame to annotate your samples
DF = data.frame(id=colnames(zfGenes),type=rep(c("ctrl","treat"),each=3))
# run DESeq2
dds = DESeqDataSetFromMatrix(zfGenes,DF,~type)
dds = DESeq(dds)
summary(results(dds),alpha=0.05)

out of 25839 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up)       : 69, 0.27%
LFC < 0 (down)     : 47, 0.18%
outliers [1]       : 1270, 4.9%
low counts [2]     : 5930, 23%
(mean count < 7)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

因此您确实获得了通过 FDR 截止值的命中。最后我们可以提取重要基因列表

res = results(dds)
res[which(res$padj < 0.05),]

在 Wilcoxon Rank 之后在 R 中使用 Benjamini-Hochberg 错误发现率时出错

Error when using the Benjamini-Hochberg false discovery rate in R after Wilcoxon Rank

r

bioinformatics

p-value