DeTSPM <- function(mtrx, grps, paired=FALSE) {
require(DEGandMore);
prepared <- PrepareDe(mtrx, grps, paired);
mtrx <- prepared[[1]];
grps <- prepared[[2]];
paired <- prepared[[3]];
if (paired) warning("Paired test not supported by TSPM; performing unpaired test instead.\n");
#######################################################################
# Downloaded from http://www.stat.purdue.edu/~doerge/software/TSPM.R
#######################################################################
##-------------------------------------------------------------------
## Name: TSPM.R
## R code for the paper by Paul L. Auer and R.W. Doerge:
## "A Two-Stage Poisson Model for Testing RNA-Seq Data"
## Date: February 2011
## Contact: Paul Auer plivermo@fhcrc.org
## R.W. Doerge doerge@purdue.edu
## Example:
## counts <- matrix(0, nrow=1000, ncol=10)
## for(i in 1:1000){
## lambda <- rpois(n=1, lambda=10)
## counts[i,] <- rpois(n=10, lambda=lambda)
## }
## x1 <- gl(n=2, k=5, labels=c("T", "C"))
## x0 <- rep(1, times=10)
## lib.size <- apply(counts,2,sum)
## result <- TSPM(counts, x1, x0, lib.size)
##---------------------------------------------------------------------
#######################################################################
###### The TSPM function ##############################################
#######################################################################
TSPM <- function(counts, x1, x0, lib.size, alpha.wh=0.05){
## Input:
#counts: a matrix of RNA-Seq gene counts (genes are rows, samples are columns)
#x1: a vector of treatment group factors (under the alternative hypothesis)
#x0: a vector of treatment group factors (under the null hypothesis)
#lib.size: a vector of RNA-Seq library sizes. This could simply be obtained
# by specifying lib.size <- apply(counts,2,sum). It may also be any other
# appropriate scaling factor.
#alpha.wh: the significance threshold to use for deciding whether a gene is overdispersed.
# Defaults to 0.05.
## Output:
#log.fold.change: a vector containing the estimated log fold changes for each gene
#pvalues: a vector containing the raw p-values testing differential expression for each gene.
#index.over.disp: a vector of integer values containing the indices of the over-dispersed genes.
#index.not.over.disp: a vector of integer values containing the indices of the non-over-dispersed genes.
#padj: a vector containing the p-values after adjusting for multiple testing using the
# method of Benjamini-Hochberg
######## The main loop that fits the GLMs to each gene ########################
### Initializing model parameters ####
n <- dim(counts)[1]
per.gene.disp <- NULL
LRT <- NULL
score.test <- NULL
LFC <- NULL
###### Fitting the GLMs for each gene #################
for(i in 1:n){
### Fit full and reduced models ###
model.1 <- glm(as.numeric(counts[i,]) ~ x1, offset=log(lib.size), family=poisson)
model.0 <- glm(as.numeric(counts[i,]) ~ x0, offset=log(lib.size), family=poisson)
### Obtain diagonals of Hat matrix from the full model fit ###
hats <- hatvalues(model.1)
### Obtain Pearson overdispersion estimate ####
per.gene.disp[i] <- sum(residuals(model.1, type="pearson")^2)/model.1$df.residual
### Obtain Likelihood ratio statistic ####
LRT[i] <- deviance(model.0)-deviance(model.1)
### Obtain score test statistic ####
score.test[i] <- 1/(2*length(counts[i,])) * sum(residuals(model.1, type="pearson")^2 - ((counts[i,] - hats*model.1$fitted.values)/model.1$fitted.values))^2
### Obtain the estimated log fold change ###
LFC[i] <- -model.1$coef[2]
}
## Initialize parameters for Working-Hotelling bands around the score TSs ###
qchi <- qchisq(df=1, (1:n-0.5)/n)
MSE <- 2
UL <- NULL
#### Obtain the upper boundary of the WH bands #######################################
xbar <- mean(qchi)
bottom <- sum((qchi-xbar)^2)
top <- (qchi-xbar)^2
s <- sqrt(MSE*(1/n) + (top/bottom))
W <- sqrt(2*qf(df1=1, df2=n-1, p=1-(alpha.wh/n)))
UL <- pmax(qchi + W*s,1)
###### Obtain the indices of the over-dispersed and not-over-dispersed genes, respectively ##########
cutoff <- min(which(sort(score.test)-UL > 0))
temp <- cutoff-1 + seq(cutoff:length(score.test))
over.disp <- which(score.test %in% sort(score.test)[temp])
not.over.disp <- setdiff(1:length(score.test), over.disp)
###### Compute p-values ####################################
p.f <- pf(LRT[over.disp]/per.gene.disp[over.disp], df1=1, df2=model.1$df.residual, lower.tail=FALSE)
p.chi <- pchisq(LRT[not.over.disp], df=1, lower.tail=FALSE)
p <- NULL
p[over.disp] <- p.f
p[not.over.disp] <- p.chi
##### Adjust the p-values using the B-H method ####################
p.bh.f <- p.adjust(p.f, method="BH")
p.bh.chi <- p.adjust(p.chi, method="BH")
final.p.bh.tagwise <- NULL
final.p.bh.tagwise[over.disp] <- p.bh.f
final.p.bh.tagwise[not.over.disp] <- p.bh.chi
### Output ###
list(log.fold.change=LFC, pvalues=p, index.over.disp=over.disp, index.not.over.disp=not.over.disp,
padj=final.p.bh.tagwise)
}
n <- sapply(grps, length);
x1 <- as.factor(rep(names(grps), n));
x0 <- rep(1, sum(n));
##########################################
res <- TSPM(mtrx, x1, x0, colSums(mtrx));
##########################################
l2 <- -res$log.fold.change;
l2 <- log2(exp(l2));
pv <- res$pvalues;
qv <- p.adjust(pv, method='BH');
aj <- colSums(mtrx)/mean(colSums(mtrx));
m1 <- rowMeans(mtrx[, grps[[1]]])/mean(aj[grps[[1]]]);
m2 <- rowMeans(mtrx[, grps[[2]]])/mean(aj[grps[[2]]]);
m1[is.na(m1)] <- 0;
m2[is.na(m2)] <- 0;
l2[is.na(l2)] <- 0;
pv[is.na(pv)] <- 1;
qv[is.na(qv)] <- 1;
s <- cbind(m1, m2, m2-m1, l2, pv, qv);
colnames(s) <- c(paste('Mean', names(grps), sep='_'), 'Mean_Change', 'LogFC', 'Pvalue', 'FDR');
rownames(s) <- rownames(mtrx);
list(stat=s, group=grps, TSPM=res);
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.