Nothing
pipeline.summarySheetsSamples <- function()
{
if(ncol(indata) >= 1000) return()
dir.create(output.paths["Summary Sheets Samples"], showWarnings=FALSE)
#### Summary Sheets ####
n.genes.in.genesets <- 0
if (preferences$activated.modules$geneset.analysis)
{
n.genes.in.genesets <- length(intersect(unique(unlist(gs.def.list)), gene.info$ids))
}
ylim.max <- 0
for (m in 1:ncol(indata))
{
h <- hist(p.g.m[,m], bre=20, plot=FALSE)
y.max <- max(h$density)
if (y.max > ylim.max)
{
ylim.max <- y.max
}
}
util.info("Writing:", file.path(output.paths["Summary Sheets Samples"], "*.pdf"))
for (m in 1:ncol(indata))
{
basename <- paste(make.names(make.unique(colnames(indata))[m]), ".pdf", sep="")
pdf(file.path(output.paths["Summary Sheets Samples"], basename), 29.7/2.54, 21/2.54, useDingbats=FALSE)
## Global Sheet
layout(matrix(c(1,2,4,1,3,0,5,5,6,7,7,8), 3, 4), widths=c(1,1,2,2), heights=c(2,1,1))
par(mar=c(0,0,0,0))
plot(0, type="n", axes=FALSE, xlab="", ylab="", xlim=c(0,1), ylim=c(0,1))
text(0.1, 0.94, colnames(indata)[m] , cex=3, adj=0)
text(0.1, 0.8, "Global Summary" , cex=1.8, adj=0)
text(0.1, 0.7, paste("%DE =", round(perc.DE.m[colnames(indata)[m]], 2)), adj=0)
all.fdr.genes <- which(fdr.g.m[,m] < 0.2)
plus.fdr.genes <- which(indata[all.fdr.genes, m] > 0)
minus.fdr.genes <- which(indata[all.fdr.genes, m] <= 0)
text(0.1, 0.65, paste("# genes with fdr < 0.2 =",
length(all.fdr.genes), " (",
length(plus.fdr.genes), "+ /",
length(minus.fdr.genes), " -)"), adj=0)
all.fdr.genes <- which(fdr.g.m[,m] < 0.1)
plus.fdr.genes <- which(indata[all.fdr.genes, m] > 0)
minus.fdr.genes <- which(indata[all.fdr.genes, m] <= 0)
text(0.1, 0.6, paste("# genes with fdr < 0.1 =",
length(all.fdr.genes)," (",
length(plus.fdr.genes), "+ /",
length(minus.fdr.genes), " -)"), adj=0)
all.fdr.genes <- which(fdr.g.m[,m] < 0.05)
plus.fdr.genes <- which(indata[all.fdr.genes, m] > 0)
minus.fdr.genes <- which(indata[all.fdr.genes, m] <= 0)
text(0.1, 0.55, paste("# genes with fdr < 0.05 =",
length(all.fdr.genes)," (",
length(plus.fdr.genes), "+ /",
length(minus.fdr.genes), " -)"), adj=0)
all.fdr.genes <- which(fdr.g.m[,m] < 0.01)
plus.fdr.genes <- which(indata[all.fdr.genes, m] > 0)
minus.fdr.genes <- which(indata[all.fdr.genes, m] <= 0)
text(0.1, 0.5, paste("# genes with fdr < 0.01 =",
length(all.fdr.genes)," (",
length(plus.fdr.genes), "+ /",
length(minus.fdr.genes), " -)"), adj=0)
text(0.1, 0.425, paste("# genes in genesets =", n.genes.in.genesets), adj=0)
text(0.1, 0.35, paste("<FC> =", round(mean(indata[,m]), 2)), adj=0)
text(0.1, 0.3, paste("<t-score> =", round(mean(t.g.m[,m]), 2)), adj=0)
text(0.1, 0.25, paste("<p-value> =", round(10 ^ mean(log10(p.g.m[,m])), 2)), adj=0)
text(0.1, 0.2, paste("<fdr> =", round(mean(fdr.g.m[,m]), 2)), adj=0)
par(mar=c(2,3,3,1))
image(matrix(metadata[,m], preferences$dim.1stLvlSom, preferences$dim.1stLvlSom),
axes=FALSE, col = color.palette.portraits(1000), main="Portrait", cex.main=1.5)
axis(1,
seq(0, 1, length.out=preferences$dim.1stLvlSom/10+1),
c(1, seq(10, preferences$dim.1stLvlSom, length.out=preferences$dim.1stLvlSom/10)),
cex.axis=1.0)
axis(2,
seq(0, 1, length.out=preferences$dim.1stLvlSom/10+1),
c(1, seq(10, preferences$dim.1stLvlSom, length.out=preferences$dim.1stLvlSom/10)),
cex.axis=1.0, las=1)
box()
image(matrix(metadata[,m], preferences$dim.1stLvlSom, preferences$dim.1stLvlSom),
axes=FALSE, col=color.palette.portraits(1000), main="Top 100 DE genes", cex.main=1.5)
par(new=TRUE)
n.genes <- 100
o <- order(p.g.m[,m])[1:n.genes]
mask <- matrix(1, preferences$dim.1stLvlSom, preferences$dim.1stLvlSom)
mask[ som.result$feature.BMU[ rownames(indata)[o] ] ] <- NA
image(mask,axes=FALSE, col = "white")
axis(1,
seq(0, 1, length.out=preferences$dim.1stLvlSom/10+1),
c(1, seq(10, preferences$dim.1stLvlSom, length.out=preferences$dim.1stLvlSom/10)),
cex.axis=1.0)
axis(2,
seq(0, 1, length.out=preferences$dim.1stLvlSom/10+1),
c(1, seq(10, preferences$dim.1stLvlSom, length.out=preferences$dim.1stLvlSom/10)),
cex.axis=1.0, las=1)
box()
n.genes <- 20
o <- o[1:n.genes]
par(mar=c(2,8,3,6))
plot(0, type="n", axes=FALSE, xlab="", ylab="", xlim=c(0,1), ylim=c(0,1))
par(mar=c(0,0,0,0))
x.coords <- c(0, 0.06, 0.2, 0.28, 0.36, 0.44, 0.52)
y.coords <- seq(0.75, 0.02, length.out=n.genes)
plot(0, type="n", axes=FALSE, xlab="", ylab="", xlim=c(0,1), ylim=c(0,1))
text(0, 0.88, "Global Genelist", cex=1.8, adj=0)
text(x.coords, rep(c(0.82, 0.80), 4)[1:7],
c("Rank", "ID", "log(FC)", "p-value", "fdr", "Metagene", "Description"),
cex=1, adj=0)
text(x.coords[1], y.coords, c(1:n.genes), adj=0)
text(x.coords[2], y.coords, rownames(indata)[o], cex=0.6, adj=0)
rect(x.coords[3]-0.02, y.coords[1]+0.01, 1, 0, border="white", col="white")
text(x.coords[3], y.coords, round(indata[o, m], 2), cex=0.6, adj=0)
text(x.coords[4], y.coords, format(p.g.m[o, m], digits=1), cex=0.6, adj=0)
text(x.coords[5], y.coords, format(fdr.g.m[o, m], digits=1), cex=0.6, adj=0)
text(x.coords[6], y.coords, gene.info$coordinates[o], cex=0.6, adj=0)
text(x.coords[7], y.coords, gene.info$descriptions[o], cex=0.6, adj=0)
par(mar=c(3,6,2,6))
hist(p.g.m[,m], bre=20, freq=FALSE, xlab="p-value", ylab="", main="p-values",
ylim=c(0,ylim.max), las=1, cex.main=1.5, cex.lab=1, cex.axis=1)
box()
mtext("Density", side=2, line=3, cex=1)
mtext("FDR", side=4, line=3, cex=1)
mtext(paste("%DE =", round(perc.DE.m[colnames(indata)[m]] ,2)), line=-1.2, cex=0.5)
abline(h=n.0.m[colnames(indata)[m]], col="gray", lwd=2)
par(new=TRUE)
plot(0, type="n", xlim=c(0,1), ylim=c(0,1), xlab="", ylab="", axes=FALSE)
axis(4, seq(0, 1, 0.2), seq(0, 1, 0.2), las=1, cex.axis=1)
o <- order(p.g.m[,m])
lines(p.g.m[o,m], Fdr.g.m[o,m], lty=2, lwd=2)
lines(p.g.m[o,m], fdr.g.m[o,m], lty=3, lwd=3)
legend("topright", c("p", expression(eta[0]), "Fdr", "fdr"),
col=c("black","gray","black","black"), lty=c(1,1,2,3),
lwd=c(1,1,1,2), cex=0.7)
if (preferences$activated.modules$geneset.analysis)
{
n.sets <- 20
top.gs.score <- sort(spot.list.samples[[m]]$GSZ.score, decreasing=TRUE)[1:n.sets]
top.gs.p <- spot.list.samples[[m]]$GSZ.p.value[names(top.gs.score)]
par(mar=c(0,0,0,0))
x.coords <- c(0, 0.1, 0.18, 0.30, 0.39, 0.47)
y.coords <- seq(0.75, 0.4, length.out=n.sets)
plot(0, type="n", axes=FALSE, xlab="", ylab="", xlim=c(0,1), ylim=c(0,1))
text(0, 0.88, "Global Geneset Analysis", cex=1.8, adj=0)
text(x.coords, 0.82, c("Rank", "GSZ", "p-value", "#all", "Geneset", ""), cex=1, adj=0)
text(x.coords[1], 0.77, "Overexpressed", cex=0.8, adj=0, font=3)
text(x.coords[1], y.coords, c(1:n.genes), adj=0)
text(x.coords[2], y.coords, round(top.gs.score, 2), cex=0.6, adj=0)
text(x.coords[3], y.coords, format(top.gs.p, digits=1), cex=0.6, adj=0)
text(x.coords[4], y.coords, sapply(gs.def.list[names(top.gs.score)],
function(x) { length(x$Genes) }), cex=0.6, adj=0)
text(x.coords[5], y.coords, sapply(gs.def.list, function(x) { x$Type })[names(top.gs.score)], cex=0.6, adj=0)
text(x.coords[6], y.coords, names(top.gs.score), cex=0.6, adj=0)
top.gs.score <- sort(spot.list.samples[[m]]$GSZ.score, decreasing=FALSE)[1:n.sets]
top.gs.p <- spot.list.samples[[m]]$GSZ.p.value[names(top.gs.score)]
y.coords <- seq(0.35, 0.02, length.out=n.sets)
text(x.coords[1], 0.37, "Underexpressed", cex=0.8, adj=0, font=3)
text(x.coords[1], y.coords, c(1:n.genes), adj=0)
text(x.coords[2], y.coords, round(top.gs.score, 2), cex=0.6, adj=0)
text(x.coords[3], y.coords, format(top.gs.p, digits=1), cex=0.6, adj=0)
text(x.coords[4], y.coords, sapply(gs.def.list[names(top.gs.score)],
function(x) { length(x$Genes) }), cex=0.6, adj=0)
text(x.coords[5], y.coords, sapply(gs.def.list, function(x) { x$Type })[names(top.gs.score)], cex=0.6, adj=0)
text(x.coords[6], y.coords, names(top.gs.score), cex=0.6, adj=0)
if (preferences$activated.modules$geneset.analysis.exact)
{
p <- spot.list.samples[[m]]$GSZ.p.value
fdrtool.result <- suppressWarnings(fdrtool(p, statistic="pvalue", verbose=FALSE, plot=FALSE))
fdr.spot.list.samples <- fdrtool.result$lfdr
Fdr.spot.list.samples <- fdrtool.result$qval
n.0.spot.list.samples <- fdrtool.result$param[1,"eta0"]
perc.DE.spot.list.samples <- 1 - n.0.spot.list.samples
par(mar=c(3,6,2,6))
hist(p, bre=20, freq=FALSE, xlab="p-value", ylab="", main="p-values",
las=1, cex.main=1.5, cex.lab=1, cex.axis=1)
box()
mtext("Density", side=2, line=3, cex=1)
mtext("FDR", side=4, line=3, cex=1)
mtext(paste("%DE =", round(perc.DE.spot.list.samples ,2)), line=-1.2, cex=0.5)
abline(h=n.0.spot.list.samples , col="gray", lwd=2)
par(new=TRUE)
plot(0, type="n", xlim=c(0,1), ylim=c(0,1), xlab="", ylab="", axes=FALSE)
axis(4, seq(0, 1, 0.2), seq(0, 1, 0.2), las=1, cex.axis=1)
o = order(p)
lines(p[o], Fdr.spot.list.samples[o], lty=2, lwd=2)
lines(p[o], fdr.spot.list.samples[o], lty=3, lwd=3)
legend("topright",
c("p", expression(eta[0]), "Fdr", "fdr"),
col=c("black","gray","black","black"),
lty=c(1,1,2,3), lwd=c(1,1,1,2), cex=0.7)
}
}
dev.off()
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.