# R/ps_fdr.R In PoissonSeq: Significance analysis of sequencing data based on a Poisson log linear model

```############################################################
#		Get false discovery genes
############################################################
PS.FDR <- function(ps.obj)
{
SMALL.VAL <- 1e-8

cat("Getting the FDR table...\n")

### order the statistics
tt <- sort(abs(ps.obj\$tt))	# from most insig to most sig
ng <- length(tt)
nvals <- ng

### get the dels
dels <- Get.Dels(tt, nvals)
nv.act <- length(dels)

### get the nc of the original data
nc <- ng - (rank(c(dels, tt), ties.method="min")[1 : nv.act] - 1 : nv.act)
nc[nc > ng] <- ng
nc[nc < 0] <- 0

### get the breaks
npermu.act <- ncol(ps.obj\$ttstar0)
ord <- order(abs(ps.obj\$tt), decreasing=T)	# from most sig to most insig
ttstar <- abs(ps.obj\$ttstar0)[ord, ]
fdr <- pval <- rep(0, nv.act)
fdr.mat <- matrix(0, nv.act, npermu.act)

############################################################
### estimate FDR
cat("Our modified permuatation plug-in method...\n")
### estimate pi0
ind.use <- floor(ng * 0.25 + 1) : floor(ng * 0.75)
thrd <- quantile(ttstar[ind.use, ], 0.5)
cat("\t\tthrd =", thrd, fill=T)
pi0 <- sum(tt <= thrd) / (0.5 * ng)
pi0 <- min(pi0, 1)

### get the fdr
disd <- floor(ng * (1 - pi0))
for (i in 1 : npermu.act)
{
fdr.mat[, i] <- (ng - disd) -
(rank(c(dels, ttstar[(disd+1) : ng, i]), ties.method="min")[1 : nv.act] -
1 : nv.act)
}
num.extr <- rowMeans(fdr.mat)
fdr <- num.extr / (nc + SMALL.VAL)
pval <- num.extr / (ng - disd + SMALL.VAL)

fdr[fdr > 1] <- 1
fdr[fdr < 0] <- 0

pval[pval > 1] <- 1
pval[pval < 0] <- 0

### the results
res <- list(
dels = dels,			# deltas
tt = ps.obj\$tt,	# the statistics
sig.ord = ord,		# the order of the stats, from most sig to most insig
nc = nc,					# number of sig-genes called
fdr = fdr,				# fdr
pi0 = pi0,				# pi0
pval = pval				# permuted p-values; not here are sorted, not the original
)

return(res)
}

############################################################
#		Get the true False discovery rate
#	sig.ord: 	order,
#						from most significant to most insignificant.
############################################################
Get.True.FDR <- function(delta, sig.ord, nc=NULL)
{
cat("Getting the true FDR table...\n")

if (is.null(nc))
{
nc <- 1 : length(sig.ord)
}
fdr <- rep(0, length(nc))

for (i in seq_along(nc))
{
if (nc[i] == 0)
{
fdr[i] <- 0
} else
{
fdr[i] <- sum(delta[sig.ord[1 : nc[i]]] == F) / nc[i]
}
}

return(list(nc=nc, fdr=fdr))
}

############################################################
#		Get the deltas
#	Note: tt must be unsigned, and sorted from
#	small to large beforehand.
#	Note: This function has been updated on Jan 1, 2011
#				to make dels always smaller than the true values.
############################################################
Get.Dels <- function(tt, nvals)
{
SMALL.VAL <- 1e-8
ng <- length(tt)

### get the breaks
if (ng <= nvals)
{
dels <- tt - SMALL.VAL
} else
{
dels <- tt[floor(ng / nvals * ((1 : nvals) - 1)) + 1] - SMALL.VAL
}
nv.act <- length(dels)

### make the dels increasing
dels <- sort(dels)
for (i in nv.act : 2)
{
if (dels[i - 1] > dels[i] - SMALL.VAL)
{
dels[i - 1] <- dels[i] - SMALL.VAL
}
}

return(dels)
}
```

## Try the PoissonSeq package in your browser

Any scripts or data that you put into this service are public.

PoissonSeq documentation built on May 1, 2019, 7:33 p.m.