wenjias / rnaseq Goto Github PK
View Code? Open in Web Editor NEWRNAseq
rnaseq's People
Watchers
rnaseq's Issues
union.DE_gene_subcluster.R
heatmap.2 <- function (x,
## dendrogram control
Rowv = TRUE,
Colv=if(symm)"Rowv" else TRUE,
distfun = dist,
hclustfun = hclust,
dendrogram = c("both","row","column","none"),
reorderfun = function(d, w) reorder(d, w),
symm = FALSE,
## data scaling
scale = c("none","row", "column"),
na.rm=TRUE,
## image plot
revC = identical(Colv, "Rowv"),
add.expr,
## mapping data to colors
breaks,
symbreaks=any(x < 0, na.rm=TRUE) || scale!="none",
## colors
col="heat.colors",
## block sepration
colsep,
rowsep,
sepcolor="white",
sepwidth=c(0.05,0.05),
## cell labeling
cellnote,
notecex=1.0,
notecol="cyan",
na.color=par("bg"),
## level trace
trace=c("column","row","both","none"),
tracecol="cyan",
hline=median(breaks),
vline=median(breaks),
linecol=tracecol,
## Row/Column Labeling
margins = c(5, 5),
ColSideColors,
RowSideColors,
cexRow = 0.2 + 1/log10(nr),
cexCol = 0.2 + 1/log10(nc),
labRow = NULL,
labCol = NULL,
srtRow = NULL,
srtCol = NULL,
adjRow = c(0,NA),
adjCol = c(NA,0),
offsetRow = 0.5,
offsetCol = 0.5,
colRow = NULL,
colCol = NULL,
fontRow = 1,
fontCol = 1,
## color key + density info
key = TRUE,
keysize = 1.5,
density.info=c("histogram","density","none"),
denscol=tracecol,
symkey = any(x < 0, na.rm=TRUE) || symbreaks,
densadj = 0.25,
key.title = NULL,
key.xlab = NULL,
key.ylab = NULL,
key.xtickfun = NULL,
key.ytickfun = NULL,
key.par=list(),
## plot labels
main = NULL,
xlab = NULL,
ylab = NULL,
## plot layout
lmat = NULL,
lhei = NULL,
lwid = NULL,
## extras
extrafun=NULL,
...
)
{
scale01 <- function(x, low=min(x), high=max(x) )
{
x <- (x-low)/(high - low)
x
}
plot.dendrogram <- stats:::plot.dendrogram
retval <- list()
scale <- if(symm && missing(scale)) "none" else match.arg(scale)
dendrogram <- match.arg(dendrogram)
trace <- match.arg(trace)
density.info <- match.arg(density.info)
if(length(col)==1 && is.character(col) )
col <- get(col, mode="function")
if(!missing(breaks) && any(duplicated(breaks)) )
stop("breaks may not contian duplicate values")
if(!missing(breaks) && (scale!="none"))
warning("Using scale=\"row\" or scale=\"column\" when breaks are",
"specified can produce unpredictable results.",
"Please consider using only one or the other.")
if ( is.null(Rowv) || is.na(Rowv) )
Rowv <- FALSE
if ( is.null(Colv) || is.na(Colv) )
Colv <- FALSE
else if( all(Colv=="Rowv") )
Colv <- Rowv
if(length(di <- dim(x)) != 2 || !is.numeric(x))
stop("`x' must be a numeric matrix")
nr <- di[1]
nc <- di[2]
if(nr <= 1 || nc <= 1)
stop("`x' must have at least 2 rows and 2 columns")
if(!is.numeric(margins) || length(margins) != 2)
stop("`margins' must be a numeric vector of length 2")
if(missing(cellnote))
cellnote <- matrix("", ncol=ncol(x), nrow=nrow(x))
if(!inherits(Rowv, "dendrogram")) {
## Check if Rowv and dendrogram arguments are consistent
if (
(
( is.logical(Rowv) && !isTRUE(Rowv) )
||
( is.null(Rowv) )
)
&&
( dendrogram %in% c("both","row") )
)
{
warning("Discrepancy: Rowv is FALSE, while dendrogram is `",
dendrogram, "'. Omitting row dendogram.")
if (dendrogram=="both")
dendrogram <- "column"
else
dendrogram <- "none"
}
}
if(!inherits(Colv, "dendrogram")) {
## Check if Colv and dendrogram arguments are consistent
if (
(
(is.logical(Colv) && !isTRUE(Colv) )
||
(is.null(Colv))
)
&&
( dendrogram %in% c("both","column")) )
{
warning("Discrepancy: Colv is FALSE, while dendrogram is `",
dendrogram, "'. Omitting column dendogram.")
if (dendrogram=="both")
dendrogram <- "row"
else
dendrogram <- "none"
}
}
## by default order by row/col mean
## if(is.null(Rowv)) Rowv <- rowMeans(x, na.rm = na.rm)
## if(is.null(Colv)) Colv <- colMeans(x, na.rm = na.rm)
## get the dendrograms and reordering indices
## if( dendrogram %in% c("both","row") )
## { ## dendrogram option is used *only* for display purposes
if(inherits(Rowv, "dendrogram"))
{
ddr <- Rowv ## use Rowv 'as-is', when it is dendrogram
rowInd <- order.dendrogram(ddr)
if(length(rowInd)>nr || any(rowInd<1 | rowInd > nr ))
stop("Rowv dendrogram doesn't match size of x")
if (length(rowInd) < nr)
nr <- length(rowInd)
}
else if (is.integer(Rowv))
{
## Compute dendrogram and do reordering based on given vector
distr <- distfun(x)
hcr <- hclustfun(distr)
ddr <- as.dendrogram(hcr)
ddr <- reorderfun(ddr, Rowv)
rowInd <- order.dendrogram(ddr)
if(nr != length(rowInd))
stop("row dendrogram ordering gave index of wrong length")
}
else if (isTRUE(Rowv))
{ ## If TRUE, compute dendrogram and do reordering based on rowMeans
Rowv <- rowMeans(x, na.rm = na.rm)
distr <- distfun(x)
hcr <- hclustfun(distr)
ddr <- as.dendrogram(hcr)
ddr <- reorderfun(ddr, Rowv)
rowInd <- order.dendrogram(ddr)
if(nr != length(rowInd))
stop("row dendrogram ordering gave index of wrong length")
}
else if(!isTRUE(Rowv))
{
rowInd <- nr:1
ddr <- as.dendrogram(hclust(dist(diag(nr))))
}
else
{
rowInd <- nr:1
ddr <- as.dendrogram(Rowv)
}
if(inherits(Colv, "dendrogram"))
{
ddc <- Colv ## use Colv 'as-is', when it is dendrogram
colInd <- order.dendrogram(ddc)
if(length(colInd)>nc || any(colInd<1 | colInd > nc ))
stop("Colv dendrogram doesn't match size of x")
if (length(colInd) < nc)
nc <- length(colInd)
}
else if(identical(Colv, "Rowv")) {
if(nr != nc)
stop('Colv = "Rowv" but nrow(x) != ncol(x)')
if(exists("ddr"))
{
ddc <- ddr
colInd <- order.dendrogram(ddc)
} else
colInd <- rowInd
} else if(is.integer(Colv))
{## Compute dendrogram and do reordering based on given vector
distc <- distfun(if(symm)x else t(x))
hcc <- hclustfun(distc)
ddc <- as.dendrogram(hcc)
ddc <- reorderfun(ddc, Colv)
colInd <- order.dendrogram(ddc)
if(nc != length(colInd))
stop("column dendrogram ordering gave index of wrong length")
}
else if (isTRUE(Colv))
{## If TRUE, compute dendrogram and do reordering based on rowMeans
Colv <- colMeans(x, na.rm = na.rm)
distc <- distfun(if(symm)x else t(x))
hcc <- hclustfun(distc)
ddc <- as.dendrogram(hcc)
ddc <- reorderfun(ddc, Colv)
colInd <- order.dendrogram(ddc)
if(nc != length(colInd))
stop("column dendrogram ordering gave index of wrong length")
}
else if(!isTRUE(Colv))
{
colInd <- 1:nc
ddc <- as.dendrogram(hclust(dist(diag(nc))))
}
else
{
colInd <- 1:nc
ddc <- as.dendrogram(Colv)
}
retval$rowInd <- rowInd
retval$colInd <- colInd
retval$call <- match.call()
## reorder x & cellnote
x <- x[rowInd, colInd]
x.unscaled <- x
cellnote <- cellnote[rowInd, colInd]
if(is.null(labRow))
labRow <- if(is.null(rownames(x))) (1:nr)[rowInd] else rownames(x)
else
labRow <- labRow[rowInd]
if(is.null(labCol))
labCol <- if(is.null(colnames(x))) (1:nc)[colInd] else colnames(x)
else
labCol <- labCol[colInd]
if(!is.null(colRow))
colRow <- colRow[rowInd]
if(!is.null(colCol))
colCol <- colCol[colInd]
if(scale == "row") {
retval$rowMeans <- rm <- rowMeans(x, na.rm = na.rm)
x <- sweep(x, 1, rm)
retval$rowSDs <- sx <- apply(x, 1, sd, na.rm = na.rm)
x <- sweep(x, 1, sx, "/")
}
else if(scale == "column") {
retval$colMeans <- rm <- colMeans(x, na.rm = na.rm)
x <- sweep(x, 2, rm)
retval$colSDs <- sx <- apply(x, 2, sd, na.rm = na.rm)
x <- sweep(x, 2, sx, "/")
}
## Set up breaks and force values outside the range into the endmost bins
if(missing(breaks) || is.null(breaks) || length(breaks)<1 )
{
if( missing(col) || is.function(col) )
breaks <- 16
else
breaks <- length(col)+1
}
if(length(breaks)==1)
{
if(!symbreaks)
breaks <- seq( min(x, na.rm=na.rm), max(x,na.rm=na.rm), length=breaks)
else
{
extreme <- max(abs(x), na.rm=TRUE)
breaks <- seq( -extreme, extreme, length=breaks )
}
}
nbr <- length(breaks)
ncol <- length(breaks)-1
if(class(col)=="function")
col <- col(ncol)
min.breaks <- min(breaks)
max.breaks <- max(breaks)
x[x<min.breaks] <- min.breaks
x[x>max.breaks] <- max.breaks
## Calculate the plot layout
if( missing(lhei) || is.null(lhei) )
lhei <- c(keysize, 4)
if( missing(lwid) || is.null(lwid) )
lwid <- c(keysize, 4)
if( missing(lmat) || is.null(lmat) )
{
lmat <- rbind(4:3, 2:1)
if(!missing(ColSideColors)) { ## add middle row to layout
if(!is.character(ColSideColors) || length(ColSideColors) != nc)
stop("'ColSideColors' must be a character vector of length ncol(x)")
lmat <- rbind(lmat[1,]+1, c(NA,1), lmat[2,]+1)
lhei <- c(lhei[1], 0.2, lhei[2])
}
if(!missing(RowSideColors)) { ## add middle column to layout
if(!is.character(RowSideColors) || length(RowSideColors) != nr)
stop("'RowSideColors' must be a character vector of length nrow(x)")
lmat <- cbind(lmat[,1]+1, c(rep(NA, nrow(lmat)-1), 1), lmat[,2]+1)
lwid <- c(lwid[1], 0.2, lwid[2])
}
lmat[is.na(lmat)] <- 0
}
if(length(lhei) != nrow(lmat))
stop("lhei must have length = nrow(lmat) = ", nrow(lmat))
if(length(lwid) != ncol(lmat))
stop("lwid must have length = ncol(lmat) =", ncol(lmat))
## Graphics `output' -----------------------
op <- par(no.readonly = TRUE)
on.exit(par(op))
layout(lmat, widths = lwid, heights = lhei, respect = FALSE)
plot.index <- 1
## draw the side bars
if(!missing(RowSideColors)) {
par(mar = c(margins[1],0, 0,0.5))
image(rbind(1:nr), col = RowSideColors[rowInd], axes = FALSE)
plot.index <- plot.index + 1
}
if(!missing(ColSideColors)) {
par(mar = c(0.5,0, 0,margins[2]))
image(cbind(1:nc), col = ColSideColors[colInd], axes = FALSE)
plot.index <- plot.index + 1
}
## draw the main carpet
par(mar = c(margins[1], 0, 0, margins[2]))
#if(scale != "none" || !symm)
# {
x <- t(x)
cellnote <- t(cellnote)
# }
if(revC)
{ ## x columns reversed
iy <- nr:1
if(exists("ddr"))
ddr <- rev(ddr)
x <- x[,iy]
cellnote <- cellnote[,iy]
}
else iy <- 1:nr
## display the main carpet
image(1:nc, 1:nr, x, xlim = 0.5+ c(0, nc), ylim = 0.5+ c(0, nr),
axes = FALSE, xlab = "", ylab = "", col=col, breaks=breaks,
...)
retval$carpet <- x
if(exists("ddr"))
retval$rowDendrogram <- ddr
if(exists("ddc"))
retval$colDendrogram <- ddc
retval$breaks <- breaks
retval$col <- col
## fill 'na' positions with na.color
# if(!invalid(na.color) & any(is.na(x)))
if(any(is.na(x)))
{
mmat <- ifelse(is.na(x), 1, NA)
image(1:nc, 1:nr, mmat, axes = FALSE, xlab = "", ylab = "",
col=na.color, add=TRUE)
}
## add column labels
if(is.null(srtCol) && is.null(colCol))
axis(1,
1:nc,
labels= labCol,
las= 2,
line= -0.5 + offsetCol,
tick= 0,
cex.axis= cexCol,
hadj=adjCol[1],
padj=adjCol[2],
font=fontCol
)
else
{
if(is.null(srtCol) || is.numeric(srtCol))
{
if(missing(adjCol) || is.null(adjCol))
adjCol=c(1,NA)
if(is.null(srtCol))
srtCol <- 90
xpd.orig <- par("xpd")
par(xpd=NA)
xpos <- axis(1, 1:nc, labels=rep("", nc), las=2, tick=0)
text(x=xpos,
y=par("usr")[3] - (1.0 + offsetCol) * strheight("M"),
labels=labCol,
##pos=1,
adj=adjCol,
cex=cexCol,
srt=srtCol,
col=colCol,
font=fontCol
)
par(xpd=xpd.orig)
}
else
warning("Invalid value for srtCol ignored.")
}
## add row labels
if(is.null(srtRow) && is.null(colRow))
{
axis(4,
iy,
labels=labRow,
las=2,
line=-0.5+offsetRow,
tick=0,
cex.axis=cexRow,
hadj=adjRow[1],
padj=adjRow[2],
font=fontRow
)
}
else
{
if(is.null(srtRow) || is.numeric(srtRow))
{
xpd.orig <- par("xpd")
par(xpd=NA)
ypos <- axis(4, iy, labels=rep("", nr), las=2, line= -0.5, tick=0)
text(x=par("usr")[2] + (1.0 + offsetRow) * strwidth("M"),
y=ypos,
labels=labRow,
adj=adjRow,
cex=cexRow,
srt=srtRow,
col=colRow,
font=fontRow
)
par(xpd=xpd.orig)
}
else
warning("Invalid value for srtRow ignored.")
}
## add row and column headings (xlab, ylab)
if(!is.null(xlab)) mtext(xlab, side = 1, line = margins[1] - 1.25)
if(!is.null(ylab)) mtext(ylab, side = 4, line = margins[2] - 1.25)
## perform user-specified function
if (!missing(add.expr))
eval(substitute(add.expr))
## add 'background' colored spaces to visually separate sections
if(!missing(colsep))
for(csep in colsep)
rect(xleft =csep+0.5, ybottom=0,
xright=csep+0.5+sepwidth[1], ytop=ncol(x)+1,
lty=1, lwd=1, col=sepcolor, border=sepcolor)
if(!missing(rowsep))
for(rsep in rowsep)
rect(xleft =0, ybottom= (ncol(x)+1-rsep)-0.5,
xright=nrow(x)+1, ytop = (ncol(x)+1-rsep)-0.5 - sepwidth[2],
lty=1, lwd=1, col=sepcolor, border=sepcolor)
## show traces
min.scale <- min(breaks)
max.scale <- max(breaks)
x.scaled <- scale01(t(x), min.scale, max.scale)
if(trace %in% c("both","column") )
{
retval$vline <- vline
vline.vals <- scale01(vline, min.scale, max.scale)
for( i in 1:length(colInd) )
{
if(!is.null(vline))
{
abline(v=i-0.5 + vline.vals, col=linecol, lty=2)
}
xv <- rep(i, nrow(x.scaled)) + x.scaled[,i] - 0.5
xv <- c(xv[1], xv)
yv <- 1:length(xv)-0.5
lines(x=xv, y=yv, lwd=1, col=tracecol, type="s")
}
}
if(trace %in% c("both","row") )
{
retval$hline <- hline
hline.vals <- scale01(hline, min.scale, max.scale)
for( i in 1:length(rowInd) )
{
if(!is.null(hline))
{
abline(h=i - 0.5 + hline.vals, col=linecol, lty=2)
}
yv <- rep(i, ncol(x.scaled)) + x.scaled[i,] - 0.5
yv <- rev(c(yv[1], yv))
xv <- length(yv):1-0.5
lines(x=xv, y=yv, lwd=1, col=tracecol, type="s")
}
}
if(!missing(cellnote))
text(x=c(row(cellnote)),
y=c(col(cellnote)),
labels=c(cellnote),
col=notecol,
cex=notecex)
plot.index <- plot.index + 1
## increment plot.index and then do
## latout_set( lmat, plot.index )
## to set to the correct plot region, instead of
## relying on plot.new().
## the two dendrograms :
par(mar = c(margins[1], 0, 0, 0))
if( dendrogram %in% c("both","row") )
{
flag <- try(
plot.dendrogram(ddr, horiz = TRUE, axes = FALSE, yaxs = "i", leaflab = "none")
)
if("try-error" %in% class(flag))
{
cond <- attr(flag, "condition")
if(!is.null(cond) && conditionMessage(cond)=="evaluation nested too deeply: infinite recursion / options(expressions=)?")
stop('Row dendrogram too deeply nested, recursion limit exceeded. Try increasing option("expressions"=...).')
}
}
else
plot.new()
par(mar = c(0, 0, if(!is.null(main)) 5 else 0, margins[2]))
if( dendrogram %in% c("both","column") )
{
flag <- try(
plot.dendrogram(ddc, axes = FALSE, xaxs = "i", leaflab = "none")
)
if("try-error" %in% class(flag))
{
cond <- attr(flag, "condition")
if(!is.null(cond) && conditionMessage(cond)=="evaluation nested too deeply: infinite recursion / options(expressions=)?")
stop('Column dendrogram too deeply nested, recursion limit exceeded. Try increasing option("expressions"=...).')
}
}
else
plot.new()
## title
if(!is.null(main)) title(main, cex.main = 1.5*op[["cex.main"]])
## Add the color-key
if(key)
{
# mar <- c(5, 4, 2, 1)
# if (!is.null(key.xlab) && is.na(key.xlab))
# mar[1] <- 2
# if (!is.null(key.ylab) && is.na(key.ylab))
# mar[2] <- 2
# if (!is.null(key.title) && is.na(key.title))
# mar[3] <- 1
mar = c(4, 0, 0, margins[2])
par(mar = mar, mgp=c(1.7, 0.7, 0)) #cex=0.75,
if (length(key.par) > 0)
do.call(par, key.par)
tmpbreaks <- breaks
if(symkey)
{
max.raw <- max(abs(c(x,breaks)),na.rm=TRUE)
min.raw <- -max.raw
tmpbreaks[1] <- -max(abs(x), na.rm=TRUE)
tmpbreaks[length(tmpbreaks)] <- max(abs(x), na.rm=TRUE)
}
else
{
min.raw <- min.breaks
max.raw <- max.breaks
}
# z <- seq(min.raw, max.raw, by=min(diff(breaks)/100))
z <- breaks
image(z=matrix(z, ncol=1),
col=col, breaks=tmpbreaks,
xaxt="n", yaxt="n")
box()
# z <- barplot(z, col=col, horiz=T, border=F)
par(usr=c(0,1,0,1))
if (is.null(key.xtickfun)) {
lv <- pretty(breaks)
xv <- scale01(as.numeric(lv), min.raw, max.raw)
xargs <- list(at=xv, labels=lv)
} else {
xargs <- key.xtickfun()
}
xargs$side <- 1
xargs$cex.axis <- 0.7
do.call(axis, xargs)
if (is.null(key.xlab)) {
if(scale=="row")
key.xlab <- "Row Z-Score"
else if(scale=="column")
key.xlab <- "Column Z-Score"
else
key.xlab <- "Value"
}
if (!is.na(key.xlab)) {
mtext(side=1, key.xlab, line=par("mgp")[1], padj=0.5, cex=par("cex") * par("cex.lab"))
}
if(density.info=="density")
{
dens <- density(x, adjust=densadj, na.rm=TRUE,
from=min.scale, to=max.scale)
omit <- dens$x < min(breaks) | dens$x > max(breaks)
dens$x <- dens$x[!omit]
dens$y <- dens$y[!omit]
dens$x <- scale01(dens$x, min.raw, max.raw)
lines(dens$x, dens$y / max(dens$y) * 0.95, col=denscol, lwd=1)
if (is.null(key.ytickfun)) {
yargs <- list(at=pretty(dens$y)/max(dens$y) * 0.95, labels=pretty(dens$y))
} else {
yargs <- key.ytickfun()
}
yargs$side <- 2
do.call(axis, yargs)
if (is.null(key.title))
key.title <- "Color Key\nand Density Plot"
if (!is.na(key.title))
title(key.title)
par(cex=0.5)
if (is.null(key.ylab))
key.ylab <- "Density"
if (!is.na(key.ylab))
mtext(side=2,key.ylab, line=par("mgp")[1], padj=0.5, cex=par("cex") * par("cex.lab"))
}
else if(density.info=="histogram")
{
h <- hist(x, plot=FALSE, breaks=breaks)
hx <- scale01(breaks, min.raw, max.raw)
hy <- c(h$counts, h$counts[length(h$counts)])
lines(hx, hy/max(hy)*0.95, lwd=1, type="s", col=denscol)
if (is.null(key.ytickfun)) {
yargs <- list(at=pretty(hy)/max(hy) * 0.95, labels=pretty(hy))
} else {
yargs <- key.ytickfun()
}
yargs$side <- 2
do.call(axis, yargs)
if (is.null(key.title))
key.title <- "Color Key\nand Histogram"
if (!is.na(key.title))
title(key.title)
par(cex=0.5)
if (is.null(key.ylab))
key.ylab <- "Count"
if (!is.na(key.ylab))
mtext(side=2,key.ylab, line=par("mgp")[1], padj=0.5, cex=par("cex") * par("cex.lab"))
}
else
if (is.null(key.title))
title("Color Key")
if(trace %in% c("both","column") )
{
vline.vals <- scale01(vline, min.raw, max.raw)
if(!is.null(vline))
{
abline(v=vline.vals, col=linecol, lty=2)
}
}
if(trace %in% c("both","row") )
{
hline.vals <- scale01(hline, min.raw, max.raw)
if(!is.null(hline))
{
abline(v=hline.vals, col=linecol, lty=2)
}
}
}
else
{
par(mar=c(0, 0, 0, 0))
plot.new()
}
## Create a table showing how colors match to (transformed) data ranges
retval$colorTable <- data.frame(
low=retval$breaks[-length(retval$breaks)],
high=retval$breaks[-1],
color=retval$col
)
# Store layout information, suggested by Jenny Drnevich
retval$layout <- list(lmat = lmat,
lhei = lhei,
lwid = lwid
)
## If user has provided an extra function, call it.
if(!is.null(extrafun))
extrafun()
invisible( retval )
}
plotheatmap <- function(data, name, ColSideColors, RowSideColors, legends, centerw, centerh, inset=-0.05, distfun=dist, labRow=NULL, fontRow=1, xlab='', xtick=TRUE, breaktag=c('num','value'), dendrogram=c('both','row','column','none'), palette_col=c('green', 'black', 'red')) {
breaktag <- match.arg(breaktag)
dendrogram <- match.arg(dendrogram)
nlabel <- ifelse(ncol(data) < 6, 3, 5)
if (missing(centerw)) {
centerw <- ifelse(ncol(data) < 15, ifelse(ncol(data) < 5, 2, ncol(data)/2), ifelse(ncol(data) > 50, ncol(data)/10, ncol(data)/5))
# centerw <- ifelse(ncol(data)/nrow(data) > 2, centerw/1.5, centerw)
}
if (missing(centerh)) {
centerh <- ifelse(nrow(data) < 15, ifelse(nrow(data) < 5, 2, nrow(data)/2), ifelse(nrow(data) > 50, nrow(data)/10, nrow(data)/5))
centerh <- ifelse(nrow(data)/ncol(data) > 2, centerh/1.5, centerh)
}
cexRow <- ifelse(nrow(data) > 15, 0.6, 0.8)
cexCol <- ifelse(ncol(data) > 15, 0.6, 0.8)
denw <- ifelse(dendrogram %in% c('both', 'row'), 1, 0.5)
denh <- ifelse(dendrogram %in% c('both', 'col'), 1, 0.5)
marginr <- ifelse(is.null(labRow), ifelse(max(nchar(rownames(data))) > 10, 5+max(nchar(rownames(data)))/6, 5), 1)
marginc <- ifelse(max(nchar(colnames(data))) > 10, 5+max(nchar(colnames(data)))/6, 5)
margins <- c(marginc, marginr)
lmat <- rbind(c(NA,3), c(2,1), c(NA,4))
lwid <- c(denw, centerw)
lhei <- c(denh, centerh, 0.75)
if(!missing(ColSideColors)) {
lmat <- rbind(lmat[1,]+1, c(NA, 1), lmat[-1,]+1)
lhei <- c(lhei[1], 0.2, lhei[2], lhei[3])
}
if(!missing(RowSideColors)) {
lmat <- cbind(lmat[,1]+1, c(rep(NA, nrow(lmat)-2), 1, NA), lmat[,-1]+1)
lwid <- c(lwid[1], 0.2, lwid[2])
}
lmat[is.na(lmat)] <- 0
pdf(paste0(name,'_heatmap.pdf'), width=sum(lwid), height=sum(lhei))
if (xtick) {
sort_data <- sort(unlist(data))
breaks_num <- seq(1,nlabel-2)/(nlabel-1)
if (breaktag == "num") {
breaks_label <- sort_data[as.integer(breaks_num*length(sort_data))]
} else {
breaks_label <- (max(sort_data)-min(sort_data))*breaks_num + min(sort_data)
}
breaks_label <- round(c(min(sort_data), breaks_label, max(sort_data)),2)
break_cols <- c()
for (i in 1:(length(breaks_label)-1)) {
break_cols <- c(break_cols, seq(breaks_label[i], breaks_label[i+1], , 100))
}
breaks <- unique(break_cols)
col <- colorRampPalette(palette_col)(n=length(breaks)-1)
xtickfun <- function() {
list(at=seq(0,1,,length(breaks_label)), labels=breaks_label)
}
heatmap.2(as.matrix(data), distfun=distfun, ColSideColors=ColSideColors, RowSideColors=RowSideColors, labRow=labRow, col=col, breaks=breaks, dendrogram=dendrogram, offsetRow=0, fontRow=fontRow, offsetCol=0.5, srtCol=45, scale='none', density.info='none', trace='none', key.title=NA, key.xlab=xlab, key.xtickfun=xtickfun, lmat=lmat, lwid=lwid, lhei=lhei, margins=margins, cexRow=cexRow, cexCol=cexCol)
} else {
col <- colorRampPalette(palette_col)(n=100)
heatmap.2(as.matrix(data), distfun=distfun, ColSideColors=ColSideColors, RowSideColors=RowSideColors, labRow=labRow, col=col, dendrogram=dendrogram, offsetRow=0, fontRow=fontRow, offsetCol=0.5, srtCol=45, scale='none', density.info='none', trace='none', key.title=NA, key.xlab=xlab, lmat=lmat, lwid=lwid, lhei=lhei, margins=margins, cexRow=cexRow, cexCol=cexCol)
}
if (!missing(legends)){
legend('topleft', legend=legends, fill=unique(ColSideColors), cex=0.8, bty='n', ncol=ceiling(length(legends)/5), inset=inset, xpd=T)
}
dev.off()
}
qnetwork <- function(M, v, tag, minimum=0.25, cut=0.5) {
library(qgraph)
pdf(paste0(tag, '_network.pdf'), width=8, height=6)
Q <- qgraph(M, groups=v, minimum=minimum, cut=cut, color=cols, posCol='red', negCol='green', vsize=1.5, borders=FALSE, labels=FALSE, legend=TRUE)
dev.off()
pdf(paste0(tag, '_network_venn.pdf'), width=8, height=6)
Q <- qgraph(Q, layout='spring', overlay=TRUE)
dev.off()
# pdf(paste0(tag, '_network_pca.pdf'), width=8, height=6)
# Q <- qgraph.pca(M, length(v), groups=v, minimum=minimum, cut=cut, color=cols, posCol='red', negCol='green', vsize=c(1,7), borders=FALSE, labels=FALSE, vTrans=200, rotation='promax')
# dev.off()
}
data <- read.table('union.DE_gene_tpm.xls', header=T, row.names=1, sep='\t', check.names=F)
sample_list <- read.table('../Sample_list2', header=F, sep='\t', colClasses='character')
tag <- 'union.DE_gene'
outdir <- '../7_DEA'
hcut <- 0.3
ttag <- 'gene'
sample_list <- sample_list[match(colnames(data), sample_list[, 2]), ]
sample_list <- sample_list[order(sample_list[,1], sample_list[,2]), ]
groupnames <- unique(sample_list[,1])
hmdata <- log2(data+1)
hmdata <- t(scale(t(hmdata), scale=F))
hmdata <- hmdata[, match(sample_list[, 2], colnames(hmdata)), drop=F]
if (ncol(hmdata) < 2) {q()}
dir.create(paste0(outdir, '/subcluster/matrix'), showWarnings=FALSE, recursive=TRUE)
# dir.create(paste0(outdir, '/heatmap'), showWarnings=FALSE, recursive=TRUE)
colors <- c('#CA0000','#0087FF','#00BA1D','#CF00FF','#00DBE2','#FFAF00','#0017F4','#006012','#6F1A1A','#B5CF00','#009E73','#F0E442','#D55E00','#CC79A7','#FF8A8A','#AA6400','#50008A','#00FF58','#E175FF','#FFCC99','#009999','#CC0066','#C0C0C0','#666600','#505050','#56B4E9','#0072B2','#33FF33','#99004c','#CCFF99','#660066','#9370DB','#D8BFD8','#BC8F8F','#2F4F4F','#FF6347','#877878','#CD5C5C','#FF0000','#00FF00','#000080','#0D7339','#3D90D9','#D98236','#DF514F','#8CCDA3','#F5E2D7','#FACBBE','#EA4885','#E69F00')
dt <- dist(hmdata, method='euclidean')
hc <- hclust(dt, method='complete')
memb <- cutree(hc, h=hcut*max(hc$height))
if (max(memb) > 30) {memb <- cutree(hc, k=30)}
rb <- rainbow(max(memb), start=0.4, end=0.95)
# if ((ncol(hmdata) >= 15) && (nrow(hmdata) < 1000)) {
# dir.create(paste0(outdir, '/network'), showWarnings=FALSE, recursive=TRUE)
# # high <- apply(data, 1, function(x) sum(x>30)) > ncol(data)*0.9
# # cordata <- hmdata[high, ]
# # v <- memb[high]
# cordata <- hmdata
# v <- memb
# t <- data.frame(num=1:length(v), v=v)
# g <- aggregate(num ~ v, t, c)
# cols <- rb[g$v]
# v <- g$num
# names(v) <- paste0('subcluster', g$v)
# M <- cor(t(cordata))
# s <- unlist(M)
# minimum <- 0.3
# i <- 1
# while (1) {
# if (length(s[s>minimum]) < 5*nrow(cordata)) {break}
# if (minimum >= 0.89) {
# i <- i+1
# minimum <- minimum + 0.1**i*9
# }else {
# minimum <- minimum + 0.1
# }
# }
# cut <- (minimum + 1)/2
# qnetwork(M, v, paste0(outdir, '/network/', tag), minimum, cut)
# }
pdf(paste0(outdir, '/subcluster/', tag, '_subcluster.pdf'))
par(mfrow=c(2, 2))
par(cex=0.6)
par(mar=c(7,5,4,1))
for (k in 1:max(memb)) {
subdata <- hmdata[memb==k, , drop=F]
write.table(subdata, paste0(outdir, '/subcluster/matrix/', tag, '.subcluster', k, '.xls'), sep='\t', quote=F, col.names=NA)
plot(NULL, xlim=c(1,ncol(subdata)), ylim=c(min(subdata),max(subdata)), main=paste0('subcluster', k, ', ', nrow(subdata), ttag, 's'), xaxt='n', xlab='', ylab=bquote(centered~log[2]*(TPM+1)))
axis(side=1, at=1:ncol(subdata), labels=colnames(subdata), las=2)
for(r in 1:nrow(subdata)) {
points(as.numeric(subdata[r,]), type='l', col='lightgray')
}
points(as.numeric(colMeans(subdata)), type='o', col=rb[k])
}
dev.off()
RowSideColors <- rb[memb]
if (length(groupnames) > 1) {
ColSideColors <- colors[as.factor(sample_list[,1])]
plotheatmap(hmdata, name=paste0(outdir, '/subcluster/', tag), ColSideColors=ColSideColors, RowSideColors=RowSideColors, legends=unique(sample_list[,1]), xlab=bquote(centered~log[2]*(TPM+1)), xtick=FALSE, labRow=NA, centerw=5, centerh=7)
} else {
plotheatmap(hmdata, name=paste0(outdir, '/subcluster/', tag), RowSideColors=RowSideColors, xlab=bquote(centered~log[2]*(TPM+1)), xtick=FALSE, labRow=NA, centerw=5, centerh=7)
}
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.