-
-
Save glickmac/aaa5a3c6182868972d21 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ### 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