在R中进行任何类型的统计分析时，您通常不需要重新发明轮子。
p.adjust(p, method = 'hochberg', n = length(p))

我同意不需要重新执行BH程序，但不同意以前的答案，即“正确的方法”是用p值调整。人们可能合理地假设“执行Benjamini-Hochberg”意味着找到临界值以保持预期的错误发现率，如Benjamini和Hochberg 1995所述; stats::p.adjust实际上是修改p值本身，以保持预期的错误发现率：https://stats.stackexchange.com/a/402217/171222
sgof::BH看起来是Benjamini and Hochberg 1995 procedure的一个很好的实现
它的worker函数bh返回测试拒绝的临界值（FDR）、调整后的p.值（无论它们值多少）以及p.值集合中的拒绝数。它的参数是p值向量u和期望的错误发现率alpha：

function(u, alpha = 0.05) {
  n = length(u)
  r = rank(u, ties.method = "max")
  bh = max(c(r[u <= (r/n) * alpha], 0), na.rm = T)
  su <- sort(u)
  jj <- which(u == 1)
  if (length(jj) != 0) 
    pi0 <- 1
  else pi0 <- min((-1/n) * sum(log(1 - u)), 1)
  if (bh == 0) {
    FDR_BH <- 0
  }
  else {
    FDR_BH <- round((pi0 * su[bh])/(ecdf(u)(su[bh])), 
                    4)
  }
  ad.p = numeric(n)
  ad.p[n] <- sort(u)[n]
  for (i in (n - 1):1) {
    ad.p[i] <- min(sort(u)[i] * (n/i), ad.p[i + 1])
  }
  return(c(list(Rejections = bh, FDR = min(FDR_BH, 1), 
                Adjusted.pvalues = sort(ad.p))))
}

字符串
有趣的是，bh返回排序后的调整后的p.值（ Package 器sgof::BH也对输入的p.值进行排序）;我认为p.adjust按照接收的顺序返回调整后的值。

set.seed(666)
my_p <- runif(10)
new_p <- p.adjust(my_p, method = "fdr")
new_p

BH_crit <- bh(my_p, alpha = 0.5)
BH_crit$FDR # critical p-value for 0.5 level false 
#discovery rate here would be ~ 0.13 

all.equal(new_p, BH_crit$Adjusted.pvalues) # FALSE
all.equal(sort(new_p), sort(BH_crit$Adjusted.pvalues))

型