#config Sys.setenv(R_XMAP_CONF_DIR="~/.exonmap") library(exonmap) #test library(exonmap) xmapConnect() symbol.to.gene("TP53") #load data #library(exonmap) #raw.data <- read.exon() #raw.data@cdfName <- "exon.pmcdf" #x <- rma(raw.data) load("mcf7.mcf10a.R") pData(x) #preprocess-simple r <- pc(x,"group",c("a","b")) plot(fc(r),log10(tt(r)),pch=20) keep <- (abs(fc(r)) > 4) & (tt(r) < 0.001) ps <- names(fc(r))[keep] length(ps) ps[1:100] #using limma library(limma) pairwise <- function(x,group,gp1,gp2,fc=log2(1),p.value=0.001) { pd <- pData(x) grp <- pd[, colnames(pd) == group] c1 <- grp == gp1 c2 <- grp == gp2 design <- as.matrix(cbind(as.numeric(c1),as.numeric(c2))) colnames(design) <- c(gp1,gp2) fit <- lmFit(x,design) ct <- paste(gp1,"-",gp2,sep="") cont.matrix <- makeContrasts(contrasts=ct,levels=design) fit2 <- contrasts.fit(fit,cont.matrix) fit2 <- eBayes(fit2) result <- decideTests(fit2,lfc=fc,p.value=p.value) keep <- result@".Data" != 0 r <- cbind(p.adjust(fit2$"p.value"[keep]),fit2$"coeff"[keep]) rownames(r) <- names(result@".Data"[keep,]) colnames(r) <- c("adj.p.value","fold.change") r } topN <- function(x,N=20,type=c("both","up","down")) { type <- match.arg(type) i <- 1:dim(x)[1] o <- order(x[,2],decreasing=T) i <- i[o] if(type == "both") N <- floor(N/2) u <- i[1:N] d <- i[length(i)- N + 1:N] keep <- switch(type, both=c(u,d), up=u, down=d) x[keep,,drop=FALSE] } r <- pairwise(x,"group","a","b",fc=4) ps.2 <- rownames(r) #first use of exonmap library(exonmap) xmapConnect() g <- symbol.to.gene("TP53") #plot a gene gene.graph(g,x,gps=list(1,2,3,4,5,6),gp.col=c("red","red","red","orange","orange","orange")) #filtering e <- exonic(ps) i <- intronic(ps) ig <- intergenic(ps) #multitarget multitarget(ps) multitarget(ps,exclude=TRUE) #tests is.exonic(ps) is.intronic(ps) is.intergenic(ps) is.multitarget(ps) #different syntax, same operations select.probewise(ps,"exonic") # or intronic, intergenic, multitarget exclude.probewise(ps,"exonic") # or intronic, intergenic, multitarget probeset.stats(ps[40:50]) #mapppings ex <- probeset.to.exon(e[1:100]) ex tr <- probeset.to.transcript(e[1:100]) g <- probeset.to.gene(e[1:100]) exon.to.transcript(ex) transcript.to.gene(tr) #whole gene gene.to.exon.probeset(g[1:2]) gene.to.exon.probeset.expr(x,g[1:2]) #more details gene.details(g) #range queries genes.in.range(7400000,7600000,1,"17") transcripts.in.range(7400000,7600000,1,"17") probesets.in.range(7400000,7600000,1,"17") # Putting it together library(exonmap) xmapConnect() # load("mcf7.mcf10a.R") r <- pc(x,"group",c("a","b")) keep <- (abs(fc(r)) > 4) & (tt(r) < 0.001) ps <- names(fc(r))[keep] ps <- exonic(ps) # also gets rid of multi-targeted probesets g <- probeset.to.gene(ps) gds <- gene.details(g) #filter on variance gepe <- gene.to.exon.probeset.expr(x,g[1:100]) #just do 100 genes fc.var <- function(g) { fcs <- apply(g,1,function(a) { mean(as.numeric(a[1:3])) - mean(as.numeric(a[4:6])) }) var(fcs) } l <- split(gepe,gepe$"gene") vs <- sapply(l,fc.var) o <- order(vs,decreasing=FALSE) g.ordered <- names(l)[o] #gene.strip plots gene.strip(g.ordered[1:40],x,list(1:3,4:6),type="mean-fc") gene.strip(rev(g.ordered)[1:40],x,list(1:3,4:6),type="mean-fc") #splicing index s <- si(x,g,list(1:3,4:6)) keep <- intersect(names(s),g) s <- s[keep] max.si <- sapply(s, function(a) { max(abs(a$si)) }) gene.av <- sapply(s, function(a) { max(a$gene.av) }) up<-gene.av>0 o<-order(max.si[up],decreasing=TRUE) gene.strip(g[up][o[1:40]],x,list(1:3,4:6),type="mean-fc", main="up",probes.min=0) o<-order(max.si[!up],decreasing=TRUE) gene.strip(g[!up][o[1:40]],x,list(1:3,4:6),type="mean-fc", main="down",probes.min=0) #plotting genes plotGene("ENSG00000151914",x,list(1:3,4:6), type="mean-fc", main="up",probes.min=0) 3