#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