Skip to content.

bioconductor.org

Bioconductor is an open source and open development software project
for the analysis and comprehension of genomic data.

Sections

lab2.R

################################################### ### chunk number 1: ################################################### set.seed(1234) x <- rnorm(5) print(x)

################################################### ### chunk number 2: ################################################### print(3:10)

################################################### ### chunk number 3: ################################################### print(x[2:4])

################################################### ### chunk number 4: ################################################### inds <- c(2,3,4) print(x[inds])

################################################### ### chunk number 5: ################################################### inds <- c(low=2,mid=3,hi=4) print(inds) print(x[ inds[c("low","hi")] ])

################################################### ### chunk number 6: ################################################### myList <- list(a=x, b=inds) print(myList[[1]]) print(myList[["a"]]) print(myList$b)

################################################### ### chunk number 7: ################################################### library(cluster) library(Biobase) library(annotate) library(genefilter) library(AggPred) library(sma)

################################################### ### chunk number 8: ################################################### library(golubEsets) data(golubTrain) data(golubTest) print(golubTrain) print(table(golubTrain$ALL.AML))

################################################### ### chunk number 9: ################################################### LS <- exprs(golubTrain) cl <- golubTrain$ALL.AML TS <- exprs(golubTest) clts <- golubTest$ALL.AML

################################################### ### chunk number 10: ################################################### LS[LS<100]<-100 TS[TS<100]<-100 LS[LS>16000]<-16000 TS[TS>16000]<-16000

################################################### ### chunk number 11: ################################################### mmfilt <- function(r=5, d=500, na.rm=TRUE) { function(x) { minval <- min(x, na.rm=na.rm) maxval <- max(x, na.rm=na.rm) (maxval/minval > 5) && (maxval-minval > 500) } }

mmfun <- mmfilt()

ffun <- filterfun(mmfun)

good <- genefilter(cbind(LS, TS), ffun )

sum(good) ##I get 3571 not 3051!

LSsub <- log10(LS[good,]) TSsub <- log10(TS[good,]) print(dim(LSsub))

################################################### ### chunk number 12: ################################################### cor.all<-cor(LSsub) plot.cor(cor.all[order(cl), order(cl)], labels=sort(cl), labcols=unclass(sort(cl)), title="Correlation matrix, all genes",zlim=c(-1,1))

################################################### ### chunk number 13: ################################################### library(ellipse) plotcorr(cor.all[order(cl), order(cl)], labels=sort(cl), labcols=unclass(sort(cl)), title="Correlation matrix, all genes")

################################################### ### chunk number 14: ################################################### gf <- gapFilter(250, 500, .1) ff <- filterfun(gf) good <- genefilter(10^LSsub, gf) sum(good) # should be 1091

cor.gf <- cor(LS[good,])

################################################### ### chunk number 15: ################################################### plot.cor(cor.gf[order(cl), order(cl)], labels=sort(cl), labcols=unclass(sort(cl)), title="Correlation matrix, filtered genes",zlim=c(-1,1))

################################################### ### chunk number 16: ################################################### plotcorr(cor.gf[order(cl), order(cl)], labels=sort(cl), labcols=unclass(sort(cl)), title="Correlation matrix, filtered genes")

################################################### ### chunk number 17: ################################################### library(edd) ALLtrDists <- edd.unsupervised( golubTrain[, golubTrain$ALL=="ALL"] ) print(table(ALLtrDists))

################################################### ### chunk number 18: ################################################### golURL <- url("http://www-genome.wi.mit.edu/mpr/data_set_ALL_AML_train.txt","r")

################################################### ### chunk number 19: ################################################### golVec <- scan(golURL, sep="\t", what="")

################################################### ### chunk number 20: ################################################### print(length(golVec)) print(length(golVec)/79)

################################################### ### chunk number 21: ################################################### golVec <- c(golVec[1:78],"",golVec[-(1:78)]) print(length(golVec)/79) golMat <- t(matrix(golVec,nr=79)) print(dim(golMat)) print(golMat[1:5,1:5])

################################################### ### chunk number 22: ################################################### cExprs <- golMat[-1,seq(3,78,2)] numExprs <- t(apply(cExprs,1,as.numeric)) print(numExprs[1:4,1:4])

################################################### ### chunk number 23: ################################################### accnos <- golMat[-1,2] dimnames(numExprs) <- list(accnos,NULL)

################################################### ### chunk number 24: ################################################### exprSamps <- golMat[1,seq(3,78,2)]

################################################### ### chunk number 25: ################################################### pheu <- url("http://www-genome.wi.mit.edu/mpr/table_ALL_AML_samples.txt","r") PH <- readLines(pheu) print(length(PH)) print(PH[1:5])

################################################### ### chunk number 26: ################################################### print(PH[10]) print(PH[47])

################################################### ### chunk number 27: ################################################### print(strsplit(PH[10],"\t"))

################################################### ### chunk number 28: ################################################### samp <- rep(NA,38) dx <- rep(NA,38) cellType <- rep(NA,38) for (i in 10:47) { tmp <- strsplit(PH[i],"\t")[[1]] samp[i-9] <- tmp[1] dx[i-9] <- tmp[3] cellType[i-9] <- tmp[6] } dx <- substring(dx,1,3) samp <- c(substring(samp[1:9],1,1),substring(samp[10:38],1,2)) phenoDF <- data.frame(samp=as.numeric(samp), dx=dx, cellType=cellType) print(phenoDF[1:4,])

################################################### ### chunk number 29: ################################################### library(Biobase) myPh <- new("phenoData", pData=phenoDF[as.numeric(exprSamps),], varLabels=list(samp="sample code", dx="diagnosis", cellType="cell type")) myES <- new("exprSet", exprs=numExprs, phenoData=myPh, description=names(phenoDF), annotation="affy", notes="from web")

################################################### ### chunk number 30: ################################################### print(myES) print(myES[1:10,]) print(myES$dx[c(1:5,37:38)])

################################################### ### chunk number 31: ################################################### print(geneNames(myES)[c(1,1000,2000)]) print(exprs(myES)[c(1,1000,2000),1:3])

News
2009-10-26

BioC 2.5, consisting of 352 packages and designed to work with R 2.10.z, was released today.

2009-01-07

R, the open source platform used by Bioconductor, featured in a series of articles in the New York Times.