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

2024-05-09 • 问答

我已经进行了Wilcoxon秩和检验，以查看三个疾病样品与三个对照样品之间598019基因的表达之间是否存在显着差异。我在河里。

当我看到有多少个基因的p值

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

（这在apply函数中，如果有必要，我可以提供全部代码，但我不确定Alternative =“ two.side”是否正确。）

但是，由于我假设使用Benjamini Hochberg错误发现率校正多个比较会降低该数字，因此我通过以下代码调整了p值 pvaluesadjust1

通过下面的代码重新评估哪个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年）大。

Wilcoxon是基于等级的参数测试。如果您只有6个样本，则可获得的最佳结果是疾病的等级2、2、2与对照组的等级5、5、5，或者相反。

例如，在下面的这些随机值上尝试使用您在测试中使用的参数，您将获得相同的p。值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),correct=F)$p.value

所以是的，使用598019，您可以获得41913

您使用了错误的测试。要回答第二个问题，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,function(i){
wilcox.test(i[1:3],i[4:6],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排名后在R中使用Benjamini-Hochberg错误发现率时出错

routernsn 回答：在Wilcoxon排名后在R中使用Benjamini-Hochberg错误发现率时出错

大家都在问