利用R手动计算基于BH算法的FDR值

R代码：

p = S2_cluster0_da_peaks$p_val[1:8]
p
FDR1 = p*length(p)/rank(p)
FDR2 = p.adjust(p,method = "BH",n=length(p))
FDR1
FDR2

结果：

> p = S2_cluster0_da_peaks$p_val[1:8]
> p
[1] 7.085652e-13 3.675757e-12 1.973203e-14 5.998462e-16 2.103066e-11 4.881295e-14 1.768549e-18 1.376856e-13
> FDR1 = p*length(p)/rank(p) # fdr = p*m/k
> FDR2 = p.adjust(p,method = "BH",n=length(p))
> FDR1
[1] 9.447536e-13 4.200865e-12 5.261875e-14 2.399385e-15 2.103066e-11 9.762590e-14 1.414839e-17 2.202969e-13
> FDR2
[1] 9.447536e-13 4.200865e-12 5.261875e-14 2.399385e-15 2.103066e-11 9.762590e-14 1.414839e-17 2.202969e-13
> 

可以看到，和内置函数计算的结果是一样的。

踩坑：
最开始搞成order函数一直达不到预期应改用rank（秩）
order：是计算排序后的位置索引，方便取值。
rank：是看原向量该值排第几。

sort(p)等于p[order(p)]
rank(p)等于order(order(p))