R tutorial part 2

R_tut_2.r

### R code from vignette source 'Presentation_2.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: init
###################################################
options(width=60)


###################################################
### code chunk number 2: baseScatterPlot
###################################################
fib <- function(n=20){
  if(n<3){
  return(c(1,1))
  } else{
    fib <- rep(0, n)
    for(i in 1:n){
      if(i <= 2){
        fib[i] <- 1
      } else{
        fib[i] <- fib[i-1] + fib[i-2]
      }
    }
  }
  return(fib)
}
fib(20)


###################################################
### code chunk number 3: baseScatterPlot
###################################################
head(cars, 3)
plot(cars$speed, cars$dist)


###################################################
### code chunk number 4: baseScatterPlot (eval = FALSE)
###################################################
 pdf(file='scatterplot.pdf', height=7, width=7)
 plot(cars$speed, cars$dist)
 dev.off()


###################################################
### code chunk number 5: baseScatterPlot
###################################################
plot(cars$speed, cars$dist)
fit <- lm(cars$dist ~ cars$speed)
abline(fit, col='red')
summary(fit)


###################################################
### code chunk number 6: baseBoxPlot (eval = FALSE)
###################################################
 pdf(file='boxplot.pdf', height=7, width=7)
 boxplot(len ~ dose * supp,  data=ToothGrowth)
 dev.off()


###################################################
### code chunk number 7: baseHistogram (eval = FALSE)
###################################################
 normal <- rnorm(1000, mean=1, sd=1)
 png(file='histo.png', height=480, width=480)
 par(mfrow=c(2,1))
 hist(normal, main="HISTOGRAM")
 plot(density(normal), col='red', lwd=2, 
     main="DENSITY")
 dev.off()


###################################################
### code chunk number 8: ggplot1 (eval = FALSE)
###################################################
 library(ggplot2)
 qplot(cyl, cty, data=mpg, geom="jitter")


###################################################
### code chunk number 9: ggplot2 (eval = FALSE)
###################################################
qplot(cyl, cty, data=mpg, geom="jitter") + geom_smooth(method="lm")


###################################################
### code chunk number 10: ggplot2 (eval = FALSE)
###################################################
qplot(dose, len, data=ToothGrowth, geom="jitter", colour=factor(supp))


###################################################
### code chunk number 11: ggplot2_box (eval = FALSE)
###################################################
 p <- ggplot(ToothGrowth, aes(as.factor(dose), len))
 theme_set(theme_bw()) # theme_set(theme_grey())
 p + geom_boxplot() + facet_grid(.~supp) + xlab("DOSE") + ylab('Tooth length')


###################################################
### code chunk number 12: welch (eval = FALSE)
###################################################
 t.test(MYC ~ condition, data = experiment, alternative = "two.sided")


###################################################
### code chunk number 13: vartest (eval = FALSE)
###################################################
 var.test(MYC ~ condition, data = experiment)


###################################################
### code chunk number 14: paired1
###################################################
  library(ISwR)
  attach(intake)
  head(intake)


###################################################
### code chunk number 15: paired2
###################################################
 t.test(pre, post, paired = T)


###################################################
### code chunk number 16: fisher1
###################################################
cohort <- matrix(c(12, 4, 3, 9), ncol=2)
colnames(cohort) <- c("AA", "aa")
rownames(cohort) <- c("t2d", "healthy")


###################################################
### code chunk number 17: fisher2
###################################################
cohort


###################################################
### code chunk number 18: fisher3
###################################################
fisher.test(cohort)


###################################################
### code chunk number 19: chi1
###################################################
cohort
chisq.test(cohort)


###################################################
### code chunk number 20: lm1
###################################################
 fit <- lm(post ~ pre)
 summary(fit)


###################################################
### code chunk number 21: detach
###################################################
detach(intake)


###################################################
### code chunk number 22: lm2
###################################################
  attach(intake)
  par(mar=c(5,4,3,0))
  plot(pre, post, cex=2)
  abline(fit, col='red', lwd=2)
  # visualize residuals
  segments(pre, fitted(fit), pre, post)
  # residuals
  # resid(fit)


###################################################
### code chunk number 23: lm2
###################################################
  attach(intake)
  par(mar=c(5,4,3,0))
  plot(pre, post, cex=2)
  abline(fit, col='red', lwd=2)
  # visualize residuals
  segments(pre, fitted(fit), pre, post)
  # residuals
  # resid(fit)


###################################################
### code chunk number 24: detach
###################################################
  detach(intake)
  detach(intake)


###################################################
### code chunk number 25: pairs (eval = FALSE)
###################################################
 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.cor <- 0.8/strwidth(txt)
   text(0.5, 0.5, txt, cex = cex.cor * r)
 }
 pairs(cystfibr[,3:5], lower.panel=panel.smooth, upper.panel=panel.cor)


###################################################
### code chunk number 26: hclust (eval = FALSE)
###################################################
 library(bioDist)
 d <- cor.dist(t(e)) # transpose matrix
 ## dendrogram
 hc = hclust(d, method = "average")
 plot(hc, labels = colnames(da.e), 
 main = "Hier. clust. Pearson", hang=-1)


###################################################
### code chunk number 27: hierarchicalClustering1 (eval = FALSE)
###################################################
library(RColorBrewer)
hmcol = colorRampPalette(brewer.pal(10, "RdBu"))(256)
library(bioDist)
d <- cor.dist(t(da.e))
library(gplots)
heatmap.2(as.matrix(d), 
  distfun=function(x){as.dist(x)}, 
  hclustfun=function(m){hclust(m, method="average")},
  symm=F, col=hmcol, trace='none', notecol='black', 
  denscol='black', notecex=0.8, dendrogram="column")
## For color annotation
heatmap.2(as.matrix(d), 
  distfun=function(x){as.dist(x)}, 
  hclustfun=function(m){hclust(m, method="average")},
  symm=F, col=hmcol, trace='none', notecol='black', 
  denscol='black', notecex=0.8, dendrogram="column", 
  ColSideColors=c(rep("red", 3), rep("blue", 3), rep("green", 4)))


###################################################
### code chunk number 28: hclustplot
###################################################
library(RColorBrewer)
hmcol = colorRampPalette(brewer.pal(10, "RdBu"))(256)
library(bioDist)
d <- cor.dist(t(da.e))
library(gplots)
heatmap.2(as.matrix(d), 
  distfun=function(x){as.dist(x)}, 
  hclustfun=function(m){hclust(m, method="average")},
  symm=F, col=hmcol, trace='none', notecol='black', 
  denscol='black', notecex=0.8, dendrogram="column")
## For color annotation
heatmap.2(as.matrix(d), 
  distfun=function(x){as.dist(x)}, 
  hclustfun=function(m){hclust(m, method="average")},
  symm=F, col=hmcol, trace='none', notecol='black', 
  denscol='black', notecex=0.8, dendrogram="column", 
  ColSideColors=c(rep("red", 3), rep("blue", 3), rep("green", 4)))


###################################################
### code chunk number 29: VolcanoPlot (eval = FALSE)
###################################################
# We need to compute the t-values for plotting the volcano plot
# Rank Product does not produce the t-value
# so we use the simpleaffy package
library(simpleaffy)
da.rma <- rma(gse12499[[1]])
results <- pairwise.comparison(da.rma , "cell_type", c("NSC", "NSC_1F"), spots = da)
plot(-{da.rp$AveFC}, -log(slot(results, "tt")))
## combine the list of up- and down-regulated genes
g.list <- c(rownames(r.nsc.1f_nsc$Table1), rownames(r.nsc.1f_nsc$Table2))
x <- -{da.rp$AveFC[g.list, ]}
y <- -{log(slot(results, "tt")[g.list])}
z <- getSYMBOL(g.list, "mouse4302")
## display the gene significant;y regulated
points(x, y, col = "green", pch = 19)


###################################################
### code chunk number 30: ScatterPlot (eval = FALSE)
###################################################
plot(slot(results, "means"),
  xlim = c(0, 15), ylim=c(0, 15), 
  xlab = "NSC", ylab = "NSC_1F", 
  main = "Gene differentially expressed between NSC and NSC_1F")
x <- slot(results, "means")[g.list, "NSC"]
y <- slot(results, "means")[g.list, "NSC_1F"]
z <- getSYMBOL(g.list, "mouse4302")
points(x, y, col = "green", pch = 19)
abline(2, 1, col='red')
abline(-2, 1, col='red')


###################################################
### code chunk number 31: VennDiagram (eval = FALSE)
###################################################
# TO INSTALL Vennerable
# install.packages("Vennerable", repos="http://R-Forge.R-project.org")
library(Vennerable)
vv <- list(upIniPS_SOMA=geneList_A, upIniPS_newSOMA=geneList_B)
vv <- Venn(vv)
plot(vv)
intersect(geneList_B, geneList_A)
# alternatively use venn() in library(gplots)


###################################################
### code chunk number 32: PAM (eval = FALSE)
###################################################
library(bioDist)
# dummy data
strength <- data.frame(strenght=c(12,15,13,15,14, 18,17,18,20,19), 
    height=c(120,110,130,125,120, 140,135,130,140,140))
# Euclidean distance
d <- euc(as.matrix(strength))
### K-means ###
# Cluster the genes for k = 2
cl <- kmeans(d, centers = 2, iter.max=1000)
plot(strength, col=cl$cluster)


###################################################
### code chunk number 33: GO graph with label (eval = FALSE)
###################################################
library(Rgraphviz)
g1 <- GOGraph(head(summary(mfhyper))$GOBPID, GOBPPARENTS)
# grabbing the label of the nodes
my.labels <- vector()
for(i in 1:length(slot(g1, "nodes"))){
  my.labels[i] <- Term(get(slot(g1, "nodes")[[i]], GOTERM))
}
my.labels
# determining the nodes attributes and associating the node labels
nodattr <- makeNodeAttrs(g1, label=my.labels, 
  shape = "ellipse", fillcolor = "#f2f2f2", fixedsize = FALSE)
# plotting the graph
plot(g1, nodeAttrs = nodattr)


###################################################
### code chunk number 34: GO analysis (eval = FALSE)
###################################################
library(GOstats)
library(mouse4302.db)
uniqueId <- mouse4302ENTREZID
entrezUniverse <- unique(as.character(uniqueId))
ourlist <- mouse4302ENTREZID[g.list]
ourlist <- unique(as.character(ourlist))
# creating the GOHyperGParams object
params = new("GOHyperGParams", geneIds=ourlist, 
universeGeneIds=entrezUniverse, annotation='mouse4302.db', 
ontology="BP", pvalueCutoff=0.001, conditional=FALSE, testDirection="over")
# running the test
mfhyper = hyperGTest(params)


###################################################
### code chunk number 35: GO analysis 2
###################################################
mfhyper
head(summary(mfhyper), 2)
# how many gene were mapped in the end?
geneMappedCount(mfhyper)


###################################################
### code chunk number 36: oligo (eval = FALSE)
###################################################
 vignette('V3AffySnpGenotype')


###################################################
### code chunk number 37: snpStats (eval = FALSE)
###################################################
 vignette('data-input-vignette')


###################################################
### code chunk number 38: GGtools (eval = FALSE)
###################################################
 ?vcf2sm


###################################################
### code chunk number 39: genetics (eval = FALSE)
###################################################
 library(genetics)
 fmsURL <- 'http://people.umass.edu/foulkes/asg/data/FMS_data.txt'
 fms <- read.delim(file=fmsURL, header=T, sep='\t')
 attach(fms); summary(genotype(actn3_rs540874, sep='')); detach(fms)


###################################################
### code chunk number 40: SNPlogistik (eval = FALSE)
###################################################
for(i in 1:nrow(SNPs)){
    fit <- lm(disease_Y-N ~ SNPs[i,] + sex + age)
    raw.pval[i] <- pf(fit$fstatistic[1], fit$fstatistic[2], 
        fit$fstatistic[3], lower.tail=F)
}


###################################################
### code chunk number 41: multiple_hypthesis (eval = FALSE)
###################################################
library(multtest)
mt.rawp2adjp(raw.pval)


###################################################
### code chunk number 42: qvalue (eval = FALSE)
###################################################
library(qvalue)
qvalue(pval)$qvalues


###################################################
### code chunk number 43: snpBoxPlot (eval = FALSE)
###################################################
boxplot(TRAIT ~ sex * rs229922, data=SNPdata, xlab="Genotype by sex", 
    ylab="TRAIT", main='SNP vs. trait', col=c('white', 'grey'))


###################################################
### code chunk number 44: googleVis_earthquakes (eval = FALSE)
###################################################
## Load the library
library(googleVis)
library(XML)
## Set an option to have the data NOT transform to factor 
## when they are imported into a data.frame
options(stringsAsFactors = FALSE)
## Get earthquake data of the last 30 days from the IRIS website
eq=readHTMLTable("http://www.iris.edu/seismon/last30.html")
## eq is a list() object with 2 elements
## Extract the earthquake data.frame
eq <- eq[[2]] 
## We look for Japan earthquakes using the grep function.
## Look at ?grep for more info
jp <- grep("*japan*", eq$REGION, ignore.case=T)
## Filter the original data.frame for the Japanese earthquakes
eq <- eq[jp,]
## Build the location information in the right format
## googleVis expect LAT:LONG as a string
## look at ?paste for the detail of this step
eq$loc <- paste(eq$LAT, eq$LON, sep=":")
## Plot the googleVis graph
M <- gvisGeoMap(eq, locationvar="loc", numvar="MAG", hovervar="DATE",  options=list(region="JP", dataMode="markers"))
plot(M)