R/Output.R
In bioinfo-lab/TADBD: A sensitive and fast tool for detection of TAD boundaries
Defines functions
Output
Coordinate
colorpanel
Documented in
Output
# @fn Output
# @param df_result : dataframe, the return df_result from function TADBD
# @param species : string, species name
# @param chr : string, chromosome number of input data
# @param resolution : int, resolution of input data
# @param outxtfile : string, directory and name of output file where the TAD boundary coordinates are written
# @param bheatmap : string, if draw heatmap, bheatmap = TRUE
# if not draw heatmap, bheatmap = (default)
# @param heatmapfile : string, directory and name of output file where heatmap and TAD boundary tracks are visualized
# @param hicmat : matrix, matrixFile by hicdata
# @param map_start : int, the starting coordinates of the heatmap
# value is an integer multiple of resolution
# @param map_end : int, the ending coordinates of the heatmap,
# value is an integer multiple of resolution
# @param l_color : string, color of the TAD track
# @param l_width : float, width of the TAD track
### Supplementary notes ###
# The source of function heatmap.2 is Rpackage gplots::heatmap.2.
# The source of function colorpanel is Rpackage gplots::colorpanel.
Output <- function(df_result, species = "hg19", chr = "chr1",
resolution = 50000, outxtfile=NULL, bheatmap = FALSE, heatmapfile=NULL, hicmat = NULL,
map_start = 0, map_end = 10000000,l_color="blue",l_width=2.5)
{
#### set range of heatmap ########
if (bheatmap == TRUE)
{
map_up <- floor(map_start/resolution)+1
map_down <- floor(map_end/resolution)+1
if (map_up > nrow(hicmat))
{
map_up <- 1
}
if (map_down > nrow(hicmat))
{
map_down <- nrow(hicmat)
}
}
#######################
if (length(df_result) == 1)
{
if ( !is.null(outxtfile) )
{
############# output all peaks ####################
pathname_df_result <- paste(outxtfile, ".allpeaks", sep="")
start <- (df_result[["allpeaks"]]$from.id3)*resolution
end <- (df_result[["allpeaks"]]$to.id3)*resolution
chrnum <- rep(chr,length(start))
outfile0 <- data.frame(chrnum, start, end)
write.table(outfile0, file = pathname_df_result, append = FALSE,
sep = " ", row.names = FALSE, col.names = TRUE)
}
############ draw heatmap #####################
if (bheatmap == TRUE && !is.null(heatmapfile))
{
x_peaks_resolution <- outfile0$end
x_peaks <- x_peaks_resolution/resolution + 1
hicmat <- log10(hicmat+1)
hicmat <- hicmat[c(map_up:map_down),c(map_up:map_down)]
x_peaks <<- x_peaks
pathname_map <- paste(heatmapfile, "_AllpeaksMap.TIFF", sep="")
tiff(pathname_map, width = 18, height = 18, units = "cm", res=350)
LXY <- Coordinate(x_peaks, map_up, map_down)
Laxx <<- LXY$Laxx
Layy <<- LXY$Layy
Lbxx <<- LXY$Lbxx
Lbyy <<- LXY$Lbyy
l_color <<- l_color
l_width <<- l_width
####################
rownames(hicmat) <- resolution * c(map_up:map_down)
heatmap.2(hicmat, Rowv=FALSE,col=colorpanel(128,"lightyellow","red"),symm=TRUE,dendrogram="none",
trace = "none", offsetRow=2.2,labCol = FALSE, xlab=paste0(chr,":",map_up*resolution,"-",map_down*resolution),
lmat=rbind(c(3,4),c(2,1)),lhei=c(1,5),lwid=c(1,5),keysize=3.5,key.xlab="",
add.expr=for(j in 1:length(Laxx)){lines(c(Laxx[j],Lbxx[j]),c(Layy[j],Lbyy[j]),
col = l_color, lwd = l_width)})
dev.off()
}
}
if (length(df_result) == 2)
{
if ( !is.null(outxtfile) )
{
############# output all peaks ####################
pathname_df_result1 <- paste(outxtfile, ".allpeaks", sep="")
start <- (df_result[["allpeaks"]]$from.id1)*resolution
end <- (df_result[["allpeaks"]]$to.id1)*resolution
pvalue <- df_result[["allpeaks"]]$pvalue
qvalue <- df_result[["allpeaks"]]$qvalue
chrnum <- rep(chr,length(start))
outfile1 <- data.frame(chrnum, start, end, pvalue, qvalue)
write.table(outfile1, file = pathname_df_result1, append = FALSE,
sep = " ", row.names = FALSE, col.names = TRUE)
############ output domains ####################
pathname_df_result2 <- paste(outxtfile, ".domains", sep="")
start <- df_result[["domains"]]$from.id2*resolution
end <- df_result[["domains"]]$to.id2*resolution
outfile2 <- data.frame(start, end)
write.table(outfile2, file = pathname_df_result2, append = FALSE,
sep = " ", row.names = FALSE, col.names = TRUE)
}
########### draw heatmap #####################
if (bheatmap == TRUE && !is.null(heatmapfile))
{
boundary_resolution <- outfile2$end
boundary <- boundary_resolution/resolution
hicmat <- log10(hicmat+1)
hicmat <- hicmat[c(map_up:map_down),c(map_up:map_down)]
pathname_map <- paste(heatmapfile, "_DomainsMap.TIFF", sep="")
tiff(pathname_map, width = 18, height = 18, units = "cm", res=350)
LXY <- Coordinate(boundary, map_up, map_down)
Laxx <<- LXY$Laxx
Layy <<- LXY$Layy
Lbxx <<- LXY$Lbxx
Lbyy <<- LXY$Lbyy
l_color <<- l_color
l_width <<- l_width
############
rownames(hicmat) <- resolution * c(map_up:map_down)
heatmap.2(hicmat, Rowv=FALSE,col=colorpanel(128,"lightyellow","red"),symm=TRUE,dendrogram="none",
trace = "none", offsetRow=2.2,labCol = FALSE, xlab=paste0(chr,":",map_up*resolution,"-",map_down*resolution),
lmat=rbind(c(3,4),c(2,1)),lhei=c(1,5),lwid=c(1,5),keysize=3.5,key.xlab="",
add.expr=for(j in 1:length(Laxx)){lines(c(Laxx[j],Lbxx[j]),c(Layy[j],Lbyy[j]),
col = l_color, lwd = l_width)})
dev.off()
}
}
}
############# Lines X Y ##############
Coordinate <- function(listC, map_up, map_down)
{
ax <- c()
ay <- c()
bx <- c()
by <- c()
j <- 1
for(i in 1:length(listC))
{
if (listC[i] >= map_up)
{
ax[j] <- listC[i]
ay[j] <- listC[i]
if (i == length(listC) | listC[i+1]>=map_down)
{
break
}
else
{
bx[j] <- listC[i+1]
by[j] <- listC[i]
}
j <- j + 1
}
}
axx <- c()
ayy <- c()
for(i in 1:length(ax))
{
axx <- c(axx,ax[i],ax[i])
ayy <- c(ayy,ay[i],ay[i])
}
bxx <- c()
byy <- c()
for(i in 1:length(bx))
{
bxx <- c(bxx,bx[i],bx[i])
byy <- c(byy,by[i],by[i])
}
Laxx <- c(map_up,axx,map_down) - map_up + 1
Layy <- (map_down) - c(map_up,axx,map_down) + map_up- map_up + 1
Lbxx <- c(axx[1],axx[1],bxx,map_down,map_down)- map_up + 1
Lbyy <- (map_down) - c(map_up,map_up,byy,axx[length(axx)],axx[length(axx)]) + map_up- map_up + 1
LXY <- data.frame(Laxx, Layy, Lbxx, Lbyy)
return(LXY)
}
############### heatmap.2 ####################
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,
## 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
}
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) || any(is.na(Rowv)) )
Rowv <- FALSE
if ( is.null(Colv) || any(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)))
# {
# 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]
)
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
)
par(xpd=xpd.orig)
}
# else
# warning("Invalid value for srtCol ignored.")
}
## add row labels
if(is.null(srtRow) && is.null(colRow))
{
axis(2,
iy,
labels=labRow,
las=2,
line=-0.5+offsetRow,
tick=0,
cex.axis=cexRow,
hadj=adjRow[1],
padj=adjRow[2]
)
}
else
{
if(is.null(srtRow) || is.numeric(srtRow))
{
xpd.orig <- par("xpd")
par(xpd=NA)
ypos <- axis(2, 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
)
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 = 2, 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
par(mar = mar, cex=0.75, mgp=c(2, 1, 0))
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))
image(z=matrix(z, ncol=1),
col=col, breaks=tmpbreaks,
xaxt="n", yaxt="n")
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
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 )
}
colorpanel <- function(n,low,mid,high)
{
if(missing(mid) || missing(high) )
{
## convert to rgb
low <- col2rgb(low)
if(missing(high))
high <- col2rgb(mid)
else
high <- col2rgb(high)
red <- seq(low[1,1], high[1,1], length=n)/255
green <- seq(low[3,1], high[3,1], length=n)/255
blue <- seq(low[2,1], high[2,1], length=n)/255
}
else # use a center color
{
isodd <- odd(n)
if(isodd)
{
n <- n+1
}
## convert to rgb
low <- col2rgb(low)
mid <- col2rgb(mid)
high <- col2rgb(high)
## determine length of each component
lower <- floor(n/2)
upper <- n - lower
red <- c(
seq(low[1,1], mid [1,1], length=lower),
seq(mid[1,1], high[1,1], length=upper)
)/255
green <- c(
seq(low[3,1], mid [3,1], length=lower),
seq(mid[3,1], high[3,1], length=upper)
)/255
blue <- c(
seq(low[2,1], mid [2,1], length=lower),
seq(mid[2,1], high[2,1], length=upper)
)/255
if(isodd)
{
red <- red [-(lower+1)]
green <- green[-(lower+1)]
blue <- blue [-(lower+1)]
}
}
rgb(red,blue,green)
}
bioinfo-lab/TADBD documentation built on March 15, 2020, 8:53 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions Output Coordinate colorpanel
Documented in Output
# @fn Output
# @param df_result : dataframe, the return df_result from function TADBD
# @param species : string, species name
# @param chr : string, chromosome number of input data
# @param resolution : int, resolution of input data
# @param outxtfile : string, directory and name of output file where the TAD boundary coordinates are written
# @param bheatmap : string, if draw heatmap, bheatmap = TRUE
# if not draw heatmap, bheatmap = (default)
# @param heatmapfile : string, directory and name of output file where heatmap and TAD boundary tracks are visualized
# @param hicmat : matrix, matrixFile by hicdata
# @param map_start : int, the starting coordinates of the heatmap
# value is an integer multiple of resolution
# @param map_end : int, the ending coordinates of the heatmap,
# value is an integer multiple of resolution
# @param l_color : string, color of the TAD track
# @param l_width : float, width of the TAD track
### Supplementary notes ###
# The source of function heatmap.2 is Rpackage gplots::heatmap.2.
# The source of function colorpanel is Rpackage gplots::colorpanel.
Output <- function(df_result, species = "hg19", chr = "chr1",
resolution = 50000, outxtfile=NULL, bheatmap = FALSE, heatmapfile=NULL, hicmat = NULL,
map_start = 0, map_end = 10000000,l_color="blue",l_width=2.5)
{
#### set range of heatmap ########
if (bheatmap == TRUE)
{
map_up <- floor(map_start/resolution)+1
map_down <- floor(map_end/resolution)+1
if (map_up > nrow(hicmat))
{
map_up <- 1
}
if (map_down > nrow(hicmat))
{
map_down <- nrow(hicmat)
}
}
#######################
if (length(df_result) == 1)
{
if ( !is.null(outxtfile) )
{
############# output all peaks ####################
pathname_df_result <- paste(outxtfile, ".allpeaks", sep="")
start <- (df_result[["allpeaks"]]$from.id3)*resolution
end <- (df_result[["allpeaks"]]$to.id3)*resolution
chrnum <- rep(chr,length(start))
outfile0 <- data.frame(chrnum, start, end)
write.table(outfile0, file = pathname_df_result, append = FALSE,
sep = " ", row.names = FALSE, col.names = TRUE)
}
############ draw heatmap #####################
if (bheatmap == TRUE && !is.null(heatmapfile))
{
x_peaks_resolution <- outfile0$end
x_peaks <- x_peaks_resolution/resolution + 1
hicmat <- log10(hicmat+1)
hicmat <- hicmat[c(map_up:map_down),c(map_up:map_down)]
x_peaks <<- x_peaks
pathname_map <- paste(heatmapfile, "_AllpeaksMap.TIFF", sep="")
tiff(pathname_map, width = 18, height = 18, units = "cm", res=350)
LXY <- Coordinate(x_peaks, map_up, map_down)
Laxx <<- LXY$Laxx
Layy <<- LXY$Layy
Lbxx <<- LXY$Lbxx
Lbyy <<- LXY$Lbyy
l_color <<- l_color
l_width <<- l_width
####################
rownames(hicmat) <- resolution * c(map_up:map_down)
heatmap.2(hicmat, Rowv=FALSE,col=colorpanel(128,"lightyellow","red"),symm=TRUE,dendrogram="none",
trace = "none", offsetRow=2.2,labCol = FALSE, xlab=paste0(chr,":",map_up*resolution,"-",map_down*resolution),
lmat=rbind(c(3,4),c(2,1)),lhei=c(1,5),lwid=c(1,5),keysize=3.5,key.xlab="",
add.expr=for(j in 1:length(Laxx)){lines(c(Laxx[j],Lbxx[j]),c(Layy[j],Lbyy[j]),
col = l_color, lwd = l_width)})
dev.off()
}
}
if (length(df_result) == 2)
{
if ( !is.null(outxtfile) )
{
############# output all peaks ####################
pathname_df_result1 <- paste(outxtfile, ".allpeaks", sep="")
start <- (df_result[["allpeaks"]]$from.id1)*resolution
end <- (df_result[["allpeaks"]]$to.id1)*resolution
pvalue <- df_result[["allpeaks"]]$pvalue
qvalue <- df_result[["allpeaks"]]$qvalue
chrnum <- rep(chr,length(start))
outfile1 <- data.frame(chrnum, start, end, pvalue, qvalue)
write.table(outfile1, file = pathname_df_result1, append = FALSE,
sep = " ", row.names = FALSE, col.names = TRUE)
############ output domains ####################
pathname_df_result2 <- paste(outxtfile, ".domains", sep="")
start <- df_result[["domains"]]$from.id2*resolution
end <- df_result[["domains"]]$to.id2*resolution
outfile2 <- data.frame(start, end)
write.table(outfile2, file = pathname_df_result2, append = FALSE,
sep = " ", row.names = FALSE, col.names = TRUE)
}
########### draw heatmap #####################
if (bheatmap == TRUE && !is.null(heatmapfile))
{
boundary_resolution <- outfile2$end
boundary <- boundary_resolution/resolution
hicmat <- log10(hicmat+1)
hicmat <- hicmat[c(map_up:map_down),c(map_up:map_down)]
pathname_map <- paste(heatmapfile, "_DomainsMap.TIFF", sep="")
tiff(pathname_map, width = 18, height = 18, units = "cm", res=350)
LXY <- Coordinate(boundary, map_up, map_down)
Laxx <<- LXY$Laxx
Layy <<- LXY$Layy
Lbxx <<- LXY$Lbxx
Lbyy <<- LXY$Lbyy
l_color <<- l_color
l_width <<- l_width
############
rownames(hicmat) <- resolution * c(map_up:map_down)
heatmap.2(hicmat, Rowv=FALSE,col=colorpanel(128,"lightyellow","red"),symm=TRUE,dendrogram="none",
trace = "none", offsetRow=2.2,labCol = FALSE, xlab=paste0(chr,":",map_up*resolution,"-",map_down*resolution),
lmat=rbind(c(3,4),c(2,1)),lhei=c(1,5),lwid=c(1,5),keysize=3.5,key.xlab="",
add.expr=for(j in 1:length(Laxx)){lines(c(Laxx[j],Lbxx[j]),c(Layy[j],Lbyy[j]),
col = l_color, lwd = l_width)})
dev.off()
}
}
}
############# Lines X Y ##############
Coordinate <- function(listC, map_up, map_down)
{
ax <- c()
ay <- c()
bx <- c()
by <- c()
j <- 1
for(i in 1:length(listC))
{
if (listC[i] >= map_up)
{
ax[j] <- listC[i]
ay[j] <- listC[i]
if (i == length(listC) | listC[i+1]>=map_down)
{
break
}
else
{
bx[j] <- listC[i+1]
by[j] <- listC[i]
}
j <- j + 1
}
}
axx <- c()
ayy <- c()
for(i in 1:length(ax))
{
axx <- c(axx,ax[i],ax[i])
ayy <- c(ayy,ay[i],ay[i])
}
bxx <- c()
byy <- c()
for(i in 1:length(bx))
{
bxx <- c(bxx,bx[i],bx[i])
byy <- c(byy,by[i],by[i])
}
Laxx <- c(map_up,axx,map_down) - map_up + 1
Layy <- (map_down) - c(map_up,axx,map_down) + map_up- map_up + 1
Lbxx <- c(axx[1],axx[1],bxx,map_down,map_down)- map_up + 1
Lbyy <- (map_down) - c(map_up,map_up,byy,axx[length(axx)],axx[length(axx)]) + map_up- map_up + 1
LXY <- data.frame(Laxx, Layy, Lbxx, Lbyy)
return(LXY)
}
############### heatmap.2 ####################
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,
## 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
}
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) || any(is.na(Rowv)) )
Rowv <- FALSE
if ( is.null(Colv) || any(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)))
# {
# 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]
)
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
)
par(xpd=xpd.orig)
}
# else
# warning("Invalid value for srtCol ignored.")
}
## add row labels
if(is.null(srtRow) && is.null(colRow))
{
axis(2,
iy,
labels=labRow,
las=2,
line=-0.5+offsetRow,
tick=0,
cex.axis=cexRow,
hadj=adjRow[1],
padj=adjRow[2]
)
}
else
{
if(is.null(srtRow) || is.numeric(srtRow))
{
xpd.orig <- par("xpd")
par(xpd=NA)
ypos <- axis(2, 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
)
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 = 2, 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
par(mar = mar, cex=0.75, mgp=c(2, 1, 0))
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))
image(z=matrix(z, ncol=1),
col=col, breaks=tmpbreaks,
xaxt="n", yaxt="n")
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
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 )
}
colorpanel <- function(n,low,mid,high)
{
if(missing(mid) || missing(high) )
{
## convert to rgb
low <- col2rgb(low)
if(missing(high))
high <- col2rgb(mid)
else
high <- col2rgb(high)
red <- seq(low[1,1], high[1,1], length=n)/255
green <- seq(low[3,1], high[3,1], length=n)/255
blue <- seq(low[2,1], high[2,1], length=n)/255
}
else # use a center color
{
isodd <- odd(n)
if(isodd)
{
n <- n+1
}
## convert to rgb
low <- col2rgb(low)
mid <- col2rgb(mid)
high <- col2rgb(high)
## determine length of each component
lower <- floor(n/2)
upper <- n - lower
red <- c(
seq(low[1,1], mid [1,1], length=lower),
seq(mid[1,1], high[1,1], length=upper)
)/255
green <- c(
seq(low[3,1], mid [3,1], length=lower),
seq(mid[3,1], high[3,1], length=upper)
)/255
blue <- c(
seq(low[2,1], mid [2,1], length=lower),
seq(mid[2,1], high[2,1], length=upper)
)/255
if(isodd)
{
red <- red [-(lower+1)]
green <- green[-(lower+1)]
blue <- blue [-(lower+1)]
}
}
rgb(red,blue,green)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.