HiTC.R
In HiTC: High Throughput Chromosome Conformation Capture analysis

### R code from vignette source 'HiTC.Rnw'

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex(use.unsrturl=FALSE)


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


###################################################
### code chunk number 3: importSimple
###################################################
require(Matrix)
## Two genome intervals objects with primers informations
reverse <- GRanges(seqnames=c("chr1","chr1"), 
                   ranges = IRanges(start=c(98831149, 98837507), 
                   end=c(98834145, 98840771), 
                   names=c("REV_2","REV_4")))

forward <- GRanges(seqnames=c("chr1","chr1"), 
                   ranges = IRanges(start=c(98834146, 98840772), 
                   end=c(98837506, 98841227), 
                   names=c("FOR_3","FOR_5")))

## 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)
detail(z)

## Access to the slots
x_intervals(z)
y_intervals(z)
intdata(z)

## Methods
range(z)
isBinned(z)
isIntraChrom(z)
seqlevels(z)


###################################################
### code chunk number 4: importMatrix
###################################################
## Load Lieberman et al. Chromosome 12 to 14 data (from GEO GSE18199)
exDir <- system.file("extdata", package="HiTC")
l <- sapply(list.files(exDir, pattern=paste("HIC_gm06690_"), full.names=TRUE),
            import.my5C)
hiC <- HTClist(l)
show(hiC)
names(hiC)

## Methods
ranges(hiC)
range(hiC)
isBinned(hiC)
isIntraChrom(hiC)
isComplete(hiC)
seqlevels(hiC)
summary(hiC)


###################################################
### code chunk number 5: qcc
###################################################
par(mfrow=c(2,2))
CQC(hiC, winsize = 1e+06, dev.new=FALSE, hist.dist=FALSE)


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


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


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


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

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


###################################################
### code chunk number 10: norm5Cexp (eval = FALSE)
###################################################
## ## Look at exptected counts
## E14exp <- getExpectedCounts(E14subset.binned, method="loess", stdev=TRUE, plot=TRUE)


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


###################################################
### code chunk number 12: norm5Cznorm
###################################################
E14norm.binned <- normPerExpected(E14subset.binned, method="loess", stdev=TRUE)
mapC(E14norm.binned)


###################################################
### code chunk number 13: annot5C
###################################################
E14.binned <- binningC(E14$chrXchrX, binsize=100000, method="median", step=3)
require(rtracklayer)
gene <- import(file.path(exDir,"refseq_mm9_chrX_98831149_103425150.bed"), 
               format="bed")
ctcf <- import(file.path(exDir,"CTCF_chrX_98892125_102969775.bed"), 
               format="bed")
mapC(E14.binned, 
     tracks=list(RefSeqGene=gene, CTCF=ctcf), 
     maxrange=10)


###################################################
### code chunk number 14: comp5C
###################################################
MEF.binned <- binningC(MEF$chrXchrX, binsize=100000, method="median", step=3)
mapC(E14.binned, MEF.binned,
     tracks=list(RefSeqGene=gene, CTCF=ctcf), 
     maxrange=10)


###################################################
### code chunk number 15: mapClist
###################################################
mapC(hiC, maxrange=100)


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


###################################################
### code chunk number 17: mapNormhic
###################################################
## Data Normalization by Expected number of Counts
hiC14norm <- normPerExpected(hiC14, method="loess")
mapC(HTClist(hiC14norm), log.data=TRUE)


###################################################
### code chunk number 18: mapCorhic
###################################################
## Correlation Map of Chromosome 14
#intdata(hiC14norm) <- as(cor(as.matrix(intdata(hiC14norm))),"Matrix")
intdata(hiC14norm) <- HiTC:::sparseCor(intdata(hiC14norm))
mapC(HTClist(hiC14norm), maxrange=1, minrange=-1, 
     col.pos=c("black", "red"), col.neg=c("blue","black"))


###################################################
### code chunk number 19: mapPCAhic
###################################################
## Principal Component Analysis
pc <- pca.hic(hiC14, normPerExpected=TRUE, method="loess", npc=1)
plot(start(pc$PC1), score(pc$PC1), type="h", 
     xlab="chr14", ylab="PC1vec", frame=FALSE)



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