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)

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

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)

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()

R Create Venn diagram

Using the package limma to draw the venn diagram.

Code:
-----------------------------------------------
A1 <- c("A", "B");
A2 <- c("A", "B", "C");
A3 <- c("D", "B");

library(limma);
SET.list <- list(A1=A1, A2=A2, A3=A3);
len <- length(SET.list);
all.unique <- unique(unlist(SET.list));
tmp <- matrix(0, ncol=len, nrow=length(all.unique));
colnames(tmp) <- names(SET.list);
rownames(tmp) <- all.unique;
for(ith in 1:len){
tmp[SET.list[[ith]], ith] <- 1;
}
count.sum <- apply(tmp, 1, sum);
count.sum
if(len < 4){
windows();
vennDiagram(vennCounts(tmp), main = "gene.intersection", cex = 1)
}
-----------------------------------------------

R colors

R colors in a way that is intended to aid finding colors by name, or by index in the It contains 657 kinds of colors:
For example:
--------------------------------------------------------------------------
> colors()[c(552,254,26)]
[1] "red" "green" "blue"
> colors()[grep("red",colors())]
[1] "darkred" "indianred" "indianred1" "indianred2"
[5] "indianred3" "indianred4" "mediumvioletred" "orangered"
[9] "orangered1" "orangered2" "orangered3" "orangered4"
[13] "palevioletred" "palevioletred1" "palevioletred2" "palevioletred3"
[17] "palevioletred4" "red" "red1" "red2"
[21] "red3" "red4" "violetred" "violetred1"
[25] "violetred2" "violetred3" "violetred4"
> colors()[grep("sky",colors())]
[1] "deepskyblue" "deepskyblue1" "deepskyblue2" "deepskyblue3"
[5] "deepskyblue4" "lightskyblue" "lightskyblue1" "lightskyblue2"
[9] "lightskyblue3" "lightskyblue4" "skyblue" "skyblue1"
[13] "skyblue2" "skyblue3" "skyblue4"
> col2rgb("yellow")
[,1]
red 255
green 255
blue 0
--------------------------------------------------------------------
A full set of color image is here, download from the reference site.



The chart can be generate by the code:

setwd("C:/");
source("http://research.stowers-institute.org/efg/R/Color/Chart/ColorChart.R")

and check the pdf file under your your C: directory.


Reference:
http://research.stowers-institute.org/efg/R/Color/Chart/index.htm
http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf

Letter counter using R

We can use R to do something that they are not designed to do so, like the letter counting using R.

In R, the package "seqinr" is design to do the data analysis for DNA sequence and protein sequence.

For example, in a txt file "tmp.txt", we have a

R string manipulation

> x <- c("a", "b", "c"); > paste(x, 1:3, sep="..");
[1] "a..1" "b..2" "c..3"
> paste(x, collapse="..");
[1] "a..b..c"
> strsplit("a.b.c", ".", fixed = TRUE)
[[1]]
[1] "a" "b" "c"

> unlist(strsplit("a.b.c", ".", fixed = TRUE));
[1] "a" "b" "c"
> strtrim(c("abcdef", "abcdef", "abcdef"), c(1,5,10))
[1] "a" "abcde" "abcdef"
> substr("abcdef",2,4)
[1] "bcd"
> substring("abcdef",c(1, 3), c(2, 6))
[1] "ab" "cdef"
> strtrim(c("abcdef", "abcdef", "abcdef"), c(1,5,10))
[1] "a" "abcde" "abcdef"

#############################################
x <- c("a", "b", "c");
paste(x, 1:3, sep="..");
paste(x, collapse="..");
strsplit("a.b.c", ".", fixed = TRUE)
unlist(strsplit("a.b.c", ".", fixed = TRUE));
strtrim(c("abcdef", "abcdef", "abcdef"), c(1,5,10))
substr("abcdef",2,4)
substring("abcdef",c(1, 3), c(2, 6))
strtrim(c("abcdef", "abcdef", "abcdef"), c(1,5,10))

R points

R use pch to control the point symbol. This can either be a single character or an integer code for one of a set of graphics symbols. The full set of S symbols is available with pch=0:18.

special symbols on R plot

Sometimes we need put special symbols, like Greek letters and italic type, on the plots. Sometimes, we also want to put mathematical annotation on the plot. In R we can use the function expression() do this job:

Sample codes

xlab.name = expression(paste(italic(vti), Delta, sep=""))
ylab.name = expression(mu * "ml")
main.name = expression(paste(plain(sin) * phi))
plot(0, 0, xlab=xlab.name, ylab=ylab.name, main=main.name, xlim=c(-pi, pi), ylim=c(-1.5, 1.5), axes=FALSE)
axis(1, at = c(-pi, -pi/2, 0, pi/2, pi), labels = expression(-pi, -pi/2, 0, pi/2, pi))
axis(2)
box()
text(-pi/2, 0, expression(hat(alpha) == (X^t * X)^{-1} * X^t * y))
text(pi/2, 0, expression(paste(frac(1, sigma*sqrt(2*pi)), plain(e)^{frac(-(x-mu)^2, 2*sigma^2)}, sep="")), cex = 1.2)

### #######################
PS: for the subscript expression,
We should use
> ylab.name = expression(sigma[21])

############################
There are webpages talking about this issue in details. Mathematical Annotation in R: http://rweb.stat.umn.edu/R/library/grDevices/html/plotmath.html
##############################
others:
main=expression(paste(italic(Sch9), Delta, " (18s and 5.8s)")),

However, if we want to use a variable in the expression, we should use substitute:
e.g.

n <- 20
plot(0, 0, main = substitute(paste(n[i], " = ", k), list(k = n)));
i=2;
range.name = substitute(paste(italic(Sch9), Delta, " ", p, "/", k, " hr"), list(k = i*12, p = (i+1)*12));
text(0, 0, range.name)
# with a variable, we can use for iteration to produce multiple tags
##############################

Plot in details

When we plot, we can assign the plot title, label.
We can also choose to show axis or not. On axis, we can set notations' letter size by cex.axis, direction by las.


Sample codes:

image.data = matrix(rnorm(60), nrow=10, ncol=6)
windows()
par(mar=c(5, 5, 4, 1))
image(x=1:nrow(image.data), y=1:ncol(image.data), image.data, axes = FALSE, xlab="X", ylab = "Y", main = "SAMPLE")
par(las=2) #control word's direction on plot, las = 1 or 2
namess1 = paste("row", 1:ncol(image.data), sep="-")
axis(2, 1:ncol(image.data), namess1, cex.axis=0.9, tick = FALSE)
par(las=2)
namess2 = paste("col", 1:nrow(image.data), sep="-")
axis(1, 1:nrow(image.data), namess2, cex.axis=0.9, tick = TRUE)
box()



windows()
par(mar=c(5, 5, 4, 1))
boxplot(data.frame(image.data), axes = FALSE, xlab="cols", ylab = "Distribution", main = "SAMPLE")
axis(2, tick = TRUE)
par(las=1)
namess2 = paste("col", 1:nrow(image.data), sep="-")
axis(1, 1:nrow(image.data), namess2, cex.axis=0.9, tick = TRUE)
box();

First Post for R

I have been using R since 2005. I love it a lot.
There are several tricks we can manipulate the R syntax to do the plot and calculation.
This blog is my note book. I wanna share my experience with other people.

post some of useful links:

download R:
http://www.r-project.org/
Bioconductor (Bioinf):
http://www.bioconductor.org/download/
R graphics:
http://addictedtor.free.fr/graphiques/thumbs.php
Eclipse (editor platform):
http://www.eclipse.org/downloads/
statET (eclipse plug-in):
http://www.walware.de/goto/statet
Rpy (use R in python):
http://rpy.sourceforge.net/
R note:
1: http://www.math.ncu.edu.tw/~chenwc/R_note/index.php?item=about