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)

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

Wednesday, July 16, 2008

color scale in R

require(graphics)
# A Color Wheel
pie(rep(1,12), col=rainbow(12))



##------ Some palettes ------------##
n = 32;
main.name = paste("color palettes; n=",n)
ch.col = c("rainbow(n, start=.7, end=.1)", "heat.colors(n)", "terrain.colors(n)", "topo.colors(n)", "cm.colors(n)");

nt <- length(ch.col)
i <- 1:n;
j <- n/nt;
d <- j/6;
dy <- 2*d;

plot(i,i+d, type="n", yaxt="n", xaxt="n", ylab="", , xlab ="", main=main.name) #yaxt="n" set no y axie label and tick.
for (k in 1:nt) {
rect(i-.5, (k-1)*j+ dy, i+.4, k*j, col = eval(parse(text=ch.col[k])), border = "grey");
text(2.5*j, k * j + dy/2, ch.col[k])
}

layout and arrange the figures on plot

## divide the device into two rows and two columns
## allocate figure 1 all of row 1
## allocate figure 2 the intersection of column 2 and row 2
# example 1
nf <- layout(mat = matrix(c(2,0,1,3),2,2,byrow=TRUE), widths = c(3,1), height = c(1,3), TRUE)
layout.show(nf)

## divide device into two rows and two columns
## allocate figure 1 and figure 2 as above
##-- Create a scatterplot with marginal histograms -----

x <- rnorm(50);
y <- runif(50);
xhist <- hist(x, plot=FALSE)
yhist <- hist(y, plot=FALSE)
top <- max(c(xhist$counts, yhist$counts))
xrange <- c(-1,1)
yrange <- c(0,1)

par(mar=c(3,3,1,1))
plot(x, y, xlim=xrange, ylim=yrange, xlab="", ylab="")
par(mar=c(0,3,1,1))
barplot(xhist$counts, axes=FALSE, ylim=c(0, top), space=0)
par(mar=c(3,0,1,1))
barplot(yhist$counts, axes=FALSE, xlim=c(0, top), space=0, horiz=TRUE)

###############################################
## example 2
## define the figure size on plot
mat <- matrix(c(1:6, 0, 7, 8), 3, 3, byrow = TRUE)
nf <- layout(mat, widths = c(3, 3, 3), height = c(3,3, 1), TRUE)
layout.show(nf)

###############################################
# example 3
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)

Tuesday, July 8, 2008

submit R jobs to Unix cluster nodes

Sample 1:
--------------------
#!/bin/sh
R --vanilla << EOF
cat("Hello", file="test.txt");
EOF
--------------------

Sample 2:
--------------------
rm(list = ls());
setwd("/auto/cmb-01/hge/");
for(i in 1:3){
filenamess = paste("JOB.", i, ".R", sep="");
cat("
#!/bin/sh
R --vanilla << EOF
#source(\"/auto/cmb-01/hge/somecode.R\");
print(", i, ");
}
EOF"
, file=filenamess);
commandlines = paste("qsub -q cmb -j oe -l walltime=120:00:00,nodes=1:myri:ppn=1 ", filenamess, sep = "");
system(commandlines);
}
-----------------------------------------------------

R 2 Latex

rm(list = ls());
#install.packages("xtable")
library(xtable)
object <- matrix(rnorm(50), 5, 10)
newobject<-xtable(object)
setwd("C:/tmp");
print(newobject, type="latex", file="filename.tex")

-----------------------
OUTPUT:

% latex table generated in R 2.5.1 by xtable 1.5-1 package
% Tue Jul 08 13:46:27 2008
\begin{table}[ht]
\begin{center}
\begin{tabular}{rrrrr}
\hline
& 1 & 2 & 3 & 4 \\
\hline
1 & 1.52 & 1.74 & $-$0.70 & $-$2.77 \\
2 & $-$1.37 & 1.70 & 1.27 & 1.22 \\
3 & 0.83 & 0.43 & $-$0.64 & $-$1.85 \\
\hline
\end{tabular}
\end{center}
\end{table}

R 2 html

Generate HTML files with R package R2html

# Sample Session
# install.packages("R2html")
rm(list = ls());
setwd("C:/tmp");
x = 1:10;
y = 1:10;

mydata = list(x=x, y=y);

library(R2HTML)
HTMLStart(outdir="C:/tmp", file="myreport", extension="html", echo=FALSE, HTMLframe=TRUE)
HTML.title("My Report", HR=1)
HTML.title("Description of my data", HR=3)
summary(mydata)
HTMLhr()
HTML.title("X Y Scatter Plot", HR=2)
plot(mydata$y~mydata$x)
HTMLplot()
HTMLStop()