Thursday, December 24, 2009

R pairs plot

RR <- matrix(rnorm(50), 10, 5)
panel.hist <- function(x, ...)
{ usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col="cyan", ...)
}
panel.blank <- function(x, y)
{ }
panel.cor <- function(x, y, digits=2, prefix="", cex.cor)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits=digits)[1]
txt <- paste(prefix, txt, sep="")
if(missing(cex.cor)) cex <- 0.6/strwidth(txt)
#text(0.5, 0.5, txt, cex = cex * r)
text(0.5, 0.5, txt, cex = cex)
}
pairs(RR, lower.panel=panel.smooth, upper.panel=panel.cor, diag.panel=panel.hist)
pairs(RR, lower.panel=panel.blank, upper.panel=panel.cor, diag.panel=panel.hist)


R animation package

It is really interesting and cool: R animation

=========================
"Description This package consists of various functions for animations in statistics, covering many areas such as probability theory, mathematical statistics, multivariate statistics, nonparametric statistics, sampling survey, linear models, time series, computational statistics, data mining and machine learning. These functions might be of help in teaching statistics and data analysis."

Example:
# Animations inside R windows graphics devices
# Bootstrapping
oopt = ani.options(interval = 0.3, nmax = 50)
boot.iid()
ani.options(oopt)

http://cran.r-project.org/web/packages/animation/animation.pdf

Tuesday, April 28, 2009

Add a table to the Figure

library(gplots);
png(filename = "exprs.png", width = 1024, height = 600)
colnames(tmp3) <- c("ORF", "Gene", "sch9/wt.d3", "ras2/wt.d3", "tor1/wt.d3",
"sch9sir2/wt.d3", "sch9sir2/sch9.d3", "CR 48/24", "CR/wt.24h",
"CR/wt.48h", "sir2/wt.logphase", "sir2/wt.d3", "sir2hmra/wt.d3")
textplot(tmp3, halign="left", valign="top", cex=0.7, show.rownames = FALSE,
show.colnames=TRUE, col.data = col.data3, cmar = 0.9, mar = c(0, 0, 0, 0) + 0.1);

Friday, April 24, 2009

Rwui

Rwui is a nice interface to create a user friendly web interface for an R script.
link: http://rwui.cryst.bbk.ac.uk/


Instructions:

(1) you need a Tomcat installation up and running (http://tomcat.apache.org/). Test Tomcat by typing: http://localhost:8080
(2) copying the file "xxx.war" to the Tomcat webapps directory, " "..\Tomcat\webapps\" where Tomcat will automatically unpack and incorporate it.
(3) Test Tomcat by typing: http://localhost:8080/xxx, then Tomcat will automatically unpack and incorporate it
(4) copy and paste the two Rdata files under the folder "..\Tomcat\webapps\xxx\WEB-INF\"
(5) enjoy to run the data interface by http://localhost:8080/xxx

note: In the 'System Variables' box scroll down and select variable 'Path' and press 'Edit'. Add the path to R's bin directory eg by adding something like C:\Program Files\R\R-2.1.0\bin to the ';' separated list.
any problems, you can check the instruction in http://rwui.cryst.bbk.ac.uk/tutorial/Instructions.html#SECTION000161000000000000000

Tuesday, February 17, 2009

combine redundant rows into one

> exprs_133a[1:4, ]
Representative.Public.ID Gene.Symbol Chromosomal.Location CHP_SKN1.CEL MDA.CEL ratios HGU133.IDs..selected.4..
1007_s_at U48705 DDR1 chr6p21.3 10.265215 10.554566 0.2893516 discoidin domain receptor tyrosine kinase 1
1053_at M87338 RFC2 chr7q11.23 9.305431 9.463867 0.1584354 replication factor C (activator 1) 2, 40kDa
117_at X51757 HSPA6 chr1q23 9.255379 9.053673 -0.2017056 heat shock 70kDa protein 6 (HSP70B)
121_at X69699 PAX8 chr2q12-q14 10.405100 10.522243 0.1171425 paired box 8
> geneID.133a <- as.character(exprs_133a[ ,1]);
> length(geneID.133a)
[1] 22283
> sum(geneID.plus2 %in% geneID.133a);
[1] 22442
> tmp1 <- apply(exprs_133a[ ,4:6], 2, function(v) tapply(v, factor(geneID.133a), mean));
> tmp2 <- apply(exprs_133a[ ,c(1, 2, 3, 7)], 2, function(v) tapply(v, factor(geneID.133a), function(v1) v1[1]));

> tmp1[1:3, ]
CHP_SKN1.CEL MDA.CEL ratios
AA001552 9.414568 9.594727 0.18015902
AA004579 8.302277 8.746708 0.44443149
AA004757 9.328749 9.383631 0.05488192
> tmp2[1:3, ]
Representative.Public.ID Gene.Symbol Chromosomal.Location HGU133.IDs..selected.4..
AA001552 "AA001552" "C19orf54" "chr19q13.2" "chromosome 19 open reading frame 54"
AA004579 "AA004579" "TAF1B" "chr2p25" "TATA box binding protein (TBP)-associated factor, RNA polymerase I, B, 63kDa"
AA004757 "AA004757" "ZNF236" "chr18q22-q23" "zinc finger protein 236"
>
> exprs_133a.unique <- data.frame(tmp1, tmp2);
> exprs_133a.unique[1:3, ];
CHP_SKN1.CEL MDA.CEL ratios Representative.Public.ID Gene.Symbol Chromosomal.Location
AA001552 9.414568 9.594727 0.18015902 AA001552 C19orf54 chr19q13.2
AA004579 8.302277 8.746708 0.44443149 AA004579 TAF1B chr2p25
AA004757 9.328749 9.383631 0.05488192 AA004757 ZNF236 chr18q22-q23
HGU133.IDs..selected.4..
AA001552 chromosome 19 open reading frame 54
AA004579 TATA box binding protein (TBP)-associated factor, RNA polymerase I, B, 63kDa
AA004757 zinc finger protein 236

Monday, February 16, 2009

use mtext() to label sub graph

windows(10, 4)
ll = matrix(c(1:3), nrow=1, byrow=TRUE)
width = c(0.4, 0.4, 0.2);
height = c(1);
layout(ll, width, height);
par(mar=c(4, 4, 2, 1));
cols <- c("green", "blue", "red")

plot(0, type="n", ylim=c(0, 1), xlim=c(0, 10), main=expression("type A"), xlab="Time (hr)",
ylab="log2(expr)", axes = FALSE);
axis(1, 1:10, (1:10)*12, cex.axis=0.7)
axis(2)
for(i in 1:3){
lines(1:10, sort(runif(10)), type = "b", col = cols[i], lty = 1, lwd=2, pch=i);
}
box();
mtext("A", at= -1, line = 0.8);

plot(0, type="n", ylim=c(0, 1), xlim=c(0, 10), main=expression("type B"), xlab="Time (hr)",
ylab="log2(expr)", axes = FALSE);
axis(1, 1:10, (1:10)*12, cex.axis=0.7)
axis(2)
for(i in 1:3){
lines(1:10, sort(runif(10)), type = "b", col = cols[i], lty = 1, lwd=2, pch=i);
}
box();
mtext("B", at= -1, line = 0.8);

genelist <- c("AA", "BB", "CC")
par(mar=c(4, 0, 2, 0));
plot(0, type="n", ylim=c(9, 14), xlim=c(0, 9), xlab="", ylab="", axes = FALSE);
legend("topleft", genelist, col=cols, lty=1, cex=0.9, lwd=2, pch=1:3, box.lwd = 0 , box.lty = 0)

---------------------------------------------

Monday, February 9, 2009

colors()



==============
use colors() to specify the color

boxplot

x1 <- rnorm(20, 1, 2);
x2 <- rnorm(20, 2, 3);
y <- list(x1, x2)
boxplot(y)

boxplot(y, horizontal = TRUE)

boxplot(y, horizontal = TRUE, col=c("red", "green"))

boxplot(y, horizontal = TRUE, col=c("red", "green"), notch=TRUE)

boxplot(y, horizontal = TRUE, col=c("pink", "lightgreen"), notch=TRUE, border = c("darkred", "darkgreen"))