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

```###################################################
#		Get transformed data
###################################################
Trans.Data <- function(n, div=10, pow.file="")
{
res <- Get.Order.Div.Overdisp(n, div=div)
ord <- res\$ord
x <- res\$mea.log

if (pow.file != "")
{
mat <- cbind(x, ord)
colnames(mat) <- c('mean.log.mu', 'one.over.theta')
write.table(mat, file=pow.file, row.names=F, col.names=T, quote=F)
}

### use linear regression
#	ord.pred <- predict(lm(ord ~ x, data=data.frame(x=x, ord=ord)), data.frame(x=log(rowMeans(n))))

### use natrual cubic spline
ord.pred <- predict(lm(ord ~ ns(x, df=3), data=data.frame(x=x, ord=ord)), data.frame(x=log(rowMeans(n))))

n.trans <- n ^ (1 / ord.pred)

ord <- Search.Ord(n.trans)
cat("The power of the transformed data is:", ord, fill=T)

return(n.trans)
}

###################################################
#		Get the overdispersion coefficient
###################################################
Get.Overdisp <- function(n, ord)
{
SMALL.VAL <- 1e-8

### take the transformation
n1 <- n ^ ord

### estimate the cmeans
res <- Est.Depth(n1)
cmeans <- res\$cmeans
keep <- res\$keep

### calculate the overdispersion
n0 <- rowSums(n1) %*% t(cmeans)
overdisp <- sum(rowSums((n1 - n0) ^ 2 / (n0 + SMALL.VAL))[keep]) -
(nrow(n) - 1) * (ncol(n) - 1) / 2

return(overdisp)
}

###################################################
#		Search for the power
###################################################
Search.Ord <- function(n, st=0.5, ed=5.5, nint=5, sround=4)
{
best.ord <- 1

for (i in 1 : sround)
{
sseq <- seq(from=st, to=ed, by=(ed-st)/nint)
overdisp <- rep(0, nint+1)
for (j in 1 : (nint+1))
{
overdisp[j] <- Get.Overdisp(n, 1/sseq[j])
}

flag <- F
for (j in 1 : nint)
{
if (overdisp[j] >= 0 & overdisp[j+1] <= 0)
{
st <- sseq[j]
ed <- sseq[j + 1]

best.ord <- switch((overdisp[j] + overdisp[j+1] >= 0) + 1, st, ed)
flag <- T
break
}
}

if (!flag)
{
cat("ERROR in searching the power: the initial power range is too small!", fill=T)
}
}

return(best.ord)
}

###################################################
#		Compare them: four divisions
###################################################
Get.Order.Div.Overdisp <- function(n, div)
{
ord <- rep(0, div)
mea.log <- rep(0, div)

n <- n[order(rowSums(n), decreasing=T), ]

nr <- floor(nrow(n) / div)
for (i in 1 : div)
{
n.div <- n[((i - 1) * nr + 1) : (i * nr), ]

mea.log[i] <- mean(log(rowMeans(n.div)))

ord[i] <- Search.Ord(n.div)
}

return(list(ord=ord, mea.log=mea.log))
}
```

## 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.