Thursday, July 17, 2008

One of my R plot

# a plot I spend one hour to adjust the layout
# Jul 17, 2008;5:00:57 PM
#
# Author: Gary HY Ge, hge AT usc DOT edu; hhggee AT gmail DOT com
# Copyright reserved
###############################################################################

mat <- matrix(c(1, 0, 0, 2:4, 5, 0, 0, 6:8, 0, 9, 10), 5, 3, byrow = TRUE)
nf <- layout(mat, widths = c(3, 2.5, 3), height = c(0.37, 3, 0.37, 3, 1.7), TRUE)
#layout.show(nf)

titlenames = c("Sample1", "Sample2"); #Just a name
data1 <- matrix(rnorm(66, mean = 0, sd = 0.5), 11, 6);
data2 <- matrix(rnorm(66, mean = 0.2, sd = 1), 11, 6);
dataset <- list(data1, data2);

for(j in 1:2){
#j = 1;
titlename = titlenames[j];
data1 <- dataset[[j]];

# the plot section name
par(mar=c(0.5, 0.5, 0.5, 0)); #define the margin of the first plot
plot(0, type="n", ylim=c(6, 15), xlim=c(0, 9), main="", xlab="", ylab="", axes = FALSE);
texts <- paste("(", titlename, ")", sep="")
text(2, 10, texts, font = 2, cex = 1.3); #font 1: plain; 2: bold; 3: italic

## histogram of the data:
par(mar=c(4, 4, 3, 1));
resids <- as.vector(data1);
ylabb = "Frequency";
hist(resids, prob=F, breaks=seq(min(resids), max(resids)+0.0599, by=0.06), ylim = c(0, 10), main = "Residual histogram", xlab = "Residuals", ylab= ylabb, xlim = c(-1.5, 1.5));#
#lines(density(resids, kernel = c("gaussian"),adjust = 0.1), col = "red", lty=1);
box();

## heatmap of the reiduals:
#par(mar=c(3, 4, 0.2, 1));
heatmapcols <- rainbow(31);
namess1 <- paste("P_", 1:nrow(data1), sep="");
namess1 <- c(" ", namess1, " ")

image.data <- data1;
cutoff1 <- 0.1;
cutoff2 <- 0.5;
index1 <- abs(image.data) < cutoff1 #small numbers
index4 <- image.data < -cutoff2 #large outliers
index5 <- image.data > cutoff2

image.data[index1] <- 0;
image.data[index4] <- -cutoff2;
image.data[index5] <- cutoff2;
image.data <- t(image.data);

len <- nrow(image.data);
bounders <- rep(cutoff2, 6);
image.data <- cbind(bounders, image.data, -bounders)
dim(image.data);

image(x=1:nrow(image.data), y=1:ncol(image.data), image.data, col=heatmapcols, axes = FALSE, xlab="", ylab="", main="Residual heatmap");
grid(6, 13)
abline(v = 3.5, col="black", lwd=2);
par(las=2)
axis(2, 1:ncol(image.data), namess1, cex.axis=0.9, tick = FALSE);
par(las=1)
print(rownames(image.data))
axis(1, 1:nrow(image.data), c("", "A", " ", " ", "B", " "), cex.axis=0.9, tick = FALSE);#Exp4_R1~3
#box();
rect(-0.60, -0.6, 3.475, 1.6, col = "white");
rect(3.525, -0.6, 6.60, 1.6, col = "white");
rect(-0.5, 12.4, 6.6, 13.7, col = "white")
rect(0.5, -0.6, 0.55, 1.6, col = "black");
rect(6.45, -0.6, 6.5, 1.6, col = "black");

par(las=1)
###########################################################
## curves
#par(mar=c(4, 4, 0.2, 1));
# good to have a set of plot type controller
setcols <- c(
" ", "blue", "blue", "1", "20",
" ", "darkblue", "blue", "2", "18",

" ", "cadetblue", "red", "1", "20",
" ", "cadetblue4", "red", "2", "20",

" ", "darkcyan", "green", "1", "20",
" ", "red", "green", "2", "20",

" ", "brown", "black", "1", "18",
" ", "darkmagenta", "black", "2", "20",

" ", "darkred", "cyan", "1", "20",
" ", "purple", "cyan", "2", "20",

" ", "green", "orange", "1", "20",
" ", "darkolivegreen", "orange", "2", "18",

" ", "darkgoldenrod1", "grey","1", "20",
" ", "tomato", "grey","2", "18"
)
setcols <- matrix(setcols, ncol=5, byrow=TRUE);
setcols <- rbind(setcols, setcols)
cols <- setcols[ ,3];
lty3 <- as.numeric(setcols[ ,4])
pchs <- as.numeric(setcols[ ,5])
###########################################################

#################################
intensities <- data1
ylabb = "Values";
plot(0, type="n", ylim=c(min(intensities), max(intensities)), xlim=c(1, ncol(intensities)), main = "curves", ylab = ylabb, xlab = "index")
for(i in 1:nrow(intensities)){
lines(1:ncol(intensities), intensities[i, ], col=cols[i], lty=lty3[i], lwd=1)
}
abline(v=3, col="grey");
}

par(mar=c(5.7, 4, 3, 1));
lens <- length(heatmapcols)
ColorLevels <- seq(-5, 5, length=lens)
image(ColorLevels, 1, matrix(data=ColorLevels, ncol=1,nrow=length(ColorLevels)), main="Color scale",
col=heatmapcols, xlab="resid values",ylab="", xaxt="n", axes = FALSE);
qvals <- c(0.001, 0.01, 0.05)
par(las=1)
axis(1, c(-5, -3, -1, 0, 1, 3, 5), c(-0.5, -0.3, -0.1, 0, 0.1, 0.3, 0.5), cex.axis=0.9, tick = TRUE)

par(mar=c(0, 4, 3, 1));
plot(0, type="n", ylim=c(6, 15), xlim=c(0, 9), main="Curve legend",
xlab="", ylab="", axes = FALSE);
box();
legend.list <- paste("P", 1:11, sep="_")
len <- length(legend.list)
selected <- 1:6;
legend("topleft", legend.list[selected], col=cols[selected], lty=lty3[selected], cex=1, lwd=1, box.lwd = 0, box.lty = 0)
selected <- 7:11;
legend("topright", legend.list[selected], col=cols[selected], lty=lty3[selected], cex=1, lwd=1, box.lwd = 0 , box.lty = 0)

###########################################################

No comments: