https://github.com/tansey/smoothfdr
Raw File
Tip revision: c5b693d0a66e83c9387433b33c0eab481bd4a763 authored by Wesley Tansey on 08 May 2020, 15:42:20 UTC
Fixed bug in easy that created too large a support for the alternative distribution
Tip revision: c5b693d
run_bh.r
# Benjamini-Hochberg with two-sided p-values
BenjaminiHochberg = function(zscores, fdr_level) {
# zscores is a vector of z scores
# fdr_level is the desired level (e.g. 0.1) of control over FDR
# returns a binary vector where 0=nofinding, 1=finding at given FDR level
    N = length(zscores)
    pval2 = 2*pmin(pnorm(zscores), 1- pnorm(zscores))
    cuts = (1:N)*fdr_level/N
    bhdiff = sort(pval2)-cuts
    bhcutind2 = max(which(bhdiff < 0))
    bhcut2 = sort(pval2)[bhcutind2]
    0+{pval2 <= bhcut2}
}

args <- commandArgs(trailingOnly = TRUE)
z = as.numeric(t(as.matrix(read.csv(args[1], header=FALSE))))#testing
discoveries = BenjaminiHochberg(z, 0.1)
write.table(discoveries, args[2], row.names=FALSE, col.names=FALSE, sep=",")
back to top