### R code from vignette source 'vignettes/HiTC/inst/doc/HiTC.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: head
###################################################
library(HiTC)
showClass("HTCexp")


###################################################
### code chunk number 2: importSimple
###################################################
## Two genome intervals objects with primers informations
reverse <- new("Genome_intervals", 
               matrix(c(98831149, 98837507, 98834145, 98840771), ncol = 2), 
               closed=c(TRUE, TRUE), 
               annotation=data.frame(seq_name=c("chr1","chr1"), 
               id=c("REV_2","REV_4"), inter_base=FALSE))

forward <- new("Genome_intervals", 
               matrix(c(98834146, 98840772, 98837506, 98841227), ncol = 2), 
               closed=c(TRUE, TRUE), 
               annotation=data.frame(seq_name=c("chr1","chr1"), 
               id=c("FOR_3","FOR_5"), inter_base=FALSE))

## A matrix of interaction counts
interac <- matrix(c(8463, 7144, 2494, 8310), ncol=2)
colnames(interac) <- c("REV_2","REV_4")
rownames(interac) <- c("FOR_3","FOR_5")

z <- HTCexp(interac, xgi=reverse, ygi=forward)
show(z)


###################################################
### code chunk number 3: importMatrix
###################################################
## Load Dekker et al. Chromosome 14 data (from GEO GSE18199)
exDir <- system.file("extdata", package="HiTC")
hiC <- import.my5C(file.path(exDir,"HIC_gm06690_chr14_chr14_1000000_obs.txt"))
detail(hiC$chr14chr14)


###################################################
### code chunk number 4: importCsv
###################################################
## Import/Export of csv format
exportC(hiC$chr14chr14, "HIC_chr14chr14.csv")
importC("HIC_chr14chr14.csv")


###################################################
### code chunk number 5: importNora
###################################################
## Load Nora et al 5C dataset
data(Nora_5C)
show(E14)
show(MEF)


###################################################
### code chunk number 6: qc5c (eval = FALSE)
###################################################
## CQC(E14)


###################################################
### code chunk number 7: QCeval
###################################################
png(file="HiTC-qc5c.png", res=300, units="in", width=5, height=5)
CQC(E14)
graphics.off()


###################################################
### code chunk number 8: mapC5c (eval = FALSE)
###################################################
## mapC(E14$chrXchrX)


###################################################
### code chunk number 9: norm5CEval
###################################################
png(file="HiTC-mapC.png", res=300, units="in", width=5, height=5)
mapC(E14$chrXchrX)
graphics.off()


###################################################
### code chunk number 10: bin5C
###################################################
## Focus on a subset chrX:100295000:102250000
E14subset<-extractRegion(E14$chrXchrX, 
                         chr="chrX", from=100295000, to=102250000)

## Binning of 5C interaction map
E14.subset.binned <- binningC(E14subset, binsize=100000, step=3)
mapC(E14.subset.binned)


###################################################
### code chunk number 11: norm5Cexp (eval = FALSE)
###################################################
## ## Look at exptected counts
## E14exp <- getExpectedCounts(E14subset, stdev=TRUE, plot=TRUE)


###################################################
### code chunk number 12: norm5CexpEval
###################################################
png(file="HiTC-expint.png", res=300, units="in", width=6, height=6)
E14exp <- getExpectedCounts(E14subset, stdev=TRUE, plot=TRUE)
graphics.off()


###################################################
### code chunk number 13: norm5Cznorm
###################################################
E14subsetz <- normPerExpected(E14subset, stdev=TRUE)
E14subsetz.binned <- binningC(E14subsetz, binsize=50000, step=3)
mapC(E14subsetz.binned)


###################################################
### code chunk number 14: annot5C
###################################################
E14.binned <- binningC(E14$chrXchrX, binsize=100000, step=3)
Refgene <- readBED(file.path(exDir,"refseq_mm9_chrX_98831149_103425150.bed"))
CTCF <- readBED(file.path(exDir,"CTCF_chrX_98892125_102969775.bed"))
mapC(E14.binned, 
     giblocs=list(RefSeqGene=Refgene$Refseq_Gene, CTCF=CTCF$CTCF), 
     maxrange=10, view=2)


###################################################
### code chunk number 15: comp5C
###################################################
MEF.binned <- binningC(MEF$chrXchrX, binsize=100000, step=3)
mapC(E14.binned, MEF.binned,
     giblocs=list(RefSeqGene=Refgene$Refseq_Gene, CTCF=CTCF$CTCF), 
     maxrange=10)


###################################################
### code chunk number 16: mapChic
###################################################
## Extract region of interest and plot the interaction map
hiC <- extractRegion(hiC$chr14chr14, 
                     chr="chr14", from=1.8e+07, to=106368584)
mapC(hiC, maxrange=100)


###################################################
### code chunk number 17: mapNormhic
###################################################
## Data Normalization by Expected number of Counts
hiCnorm <- normPerExpected(hiC)
mapC(hiCnorm, log.data=TRUE)


###################################################
### code chunk number 18: mapCorhic
###################################################
## Correlation Map of Chromosome 14
mapC(cor(intdata(hiCnorm)), maxrange=1, minrange=-1, 
     col.pos=c("black", NA, "red"), col.neg=c("black",NA, "blue"))


###################################################
### code chunk number 19: sessionInfo
###################################################
toLatex(sessionInfo(), locale=FALSE)