Skip to content

Instantly share code, notes, and snippets.

@glickmac
Forked from bobthecat/R_tut_2.r
Created July 9, 2014 19:14
Show Gist options
  • Save glickmac/aaa5a3c6182868972d21 to your computer and use it in GitHub Desktop.
Save glickmac/aaa5a3c6182868972d21 to your computer and use it in GitHub Desktop.
### 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment