###################################################
### chunk number 1: Loading rMAT
###################################################
library(rMAT)


###################################################
### chunk number 2: Reading the BPMAP File Header
###################################################


pwd<-"" #INPUT FILES- BPMAP, ARRAYS, etc.
path<- system.file("doc/Sc03b_MR_v04_10000.bpmap",package="rMAT")

bpmapFile = paste(pwd, path, sep="")
seqHeader <-ReadBPMAPAllSeqHeader(bpmapFile)  #save the list of output contents to seqHeader
#print(seqHeader)


###################################################
### chunk number 3: arrayFile
###################################################

pathCEL<- system.file("doc/Swr1WTIP_Short.CEL",package="rMAT")
arrayFile<-paste(pwd,c(pathCEL),sep="")


###################################################
### chunk number 4: BPMAPCelParser
###################################################
ScSet<-BPMAPCelParser(bpmapFile, arrayFile, verbose=FALSE,groupName="Sc")     


###################################################
### chunk number 5: Summary of Scset
###################################################
summary(ScSet)


###################################################
### chunk number 6: Normalize array by array
###################################################
ScSetNorm<-NormalizeProbes(ScSet,method="MAT",robust=FALSE,all=FALSE,standard=TRUE,verbose=FALSE)  


###################################################
### chunk number 7: Summary of ScSetNorm
###################################################
summary(ScSetNorm)


###################################################
### chunk number 8: COMPUTING rMAT SCORES
###################################################
ScScore<- MATScore(ScSetNorm, cName=NULL, dMax=600,nProbesMin=8, dMerge=300,method="score",threshold=5,verbos=TRUE,bedName="MyBedFile")


###################################################
### chunk number 9: Reading library
###################################################
library(GenomeGraphs)


###################################################
### chunk number 10: Ensembl BioMart
###################################################

mart<-useMart("ensembl",dataset="scerevisiae_gene_ensembl")


###################################################
### chunk number 11: Genome Axis
###################################################
genomeAxis<-makeGenomeAxis(add53 = TRUE, add35=TRUE)
minbase<-1
maxbase<-100000


###################################################
### chunk number 12: Plotting Gene eval=FALSE
###################################################
## genesplus<-makeGeneRegion(start = minbase, end = maxbase, strand = "+", chromosome = "I", biomart = mart) 
## genesmin<-makeGeneRegion(start = minbase, end = maxbase, strand = "-", chromosome = "I", biomart = mart)


###################################################
### chunk number 13: MatScore for chromosome I eval=FALSE
###################################################
## MatScore<-makeGenericArray(intensity=as.matrix(ScScore@score[ScScore@featureChromosome=="chr1"]),  probeStart=ScScore@featurePosition[ScScore@featureChromosome=="chr1"], dp=DisplayPars(size=1, color="black", type="l"))


###################################################
### chunk number 14: Overlays eval=FALSE
###################################################
## featurePositionForRegion<-ScScore@featurePosition[ScScore@featureChromosome=="chr1" & ScScore@featurePosition< maxbase & ScScore@featurePosition> minbase]
## regIndexForRegion<-ScScore@regIndex[ScScore@featureChromosome=="chr1" & ScScore@featurePosition< maxbase & ScScore@featurePosition> minbase]
## RegionUnique<-unique(regIndexForRegion[regIndexForRegion>0])
## rectList<-vector("list",length(RegionUnique))
## for(i in 1:length(RegionUnique))
## {
##   ## Minimum
##   m<-min(featurePositionForRegion[regIndexForRegion==RegionUnique[i]])
##   ## Maximum
##   M<-max(featurePositionForRegion[regIndexForRegion==RegionUnique[i]])
##   rectList[i] <- makeRectangleOverlay(start = m, end = M, region = c(1, 4), dp = DisplayPars(color = "green", alpha = 0.1))
## }


###################################################
### chunk number 15:  eval=FALSE
###################################################
## gdPlot(list("score" = MatScore, "Gene +" = genesplus, Position = genomeAxis, "Gene -" = genesmin), minBase = minbase, maxBase = maxbase, labelCex = 1, overlays=rectList)