## ----options, echo=FALSE-------------------------------------------------
options("scipen"=10, "digits"=5)

## ----imports, warning=F, results="hide", message=FALSE-------------------
library(GenomicInteractions)
library(BSgenome.Mmusculus.UCSC.mm9)

## ------------------------------------------------------------------------
hic_file <- system.file("extdata", "Seitan2013_WT_100kb_interactions.txt", 
                        package="GenomicInteractions")

hic_data <- GenomicInteractions(hic_file, 
                    type="homer", 
                    experiment_name = "HiC 100kb", 
                    description = "HiC 100kb resolution",
                    gname="BSgenome.Mmusculus.UCSC.mm9")

## ------------------------------------------------------------------------
summary(width(anchorOne(hic_data)))
summary(width(anchorTwo(hic_data)))
hic_data

## ------------------------------------------------------------------------
head(count(hic_data))
mean(count(hic_data))

## ------------------------------------------------------------------------
plot(density(pValue(hic_data)))
plot(density(FDR(hic_data)))

## ------------------------------------------------------------------------
sum(FDR(hic_data) < 0.1)
hic_data_subset <- hic_data[FDR(hic_data) < 0.1]

## ------------------------------------------------------------------------
plotCisTrans(hic_data)
plotCisTrans(hic_data_subset)

plotCounts(hic_data, cut=30)
plotCounts(hic_data_subset, cut=30)

## ----eval=FALSE----------------------------------------------------------
#  ## Not run
#  library(GenomicFeatures)
#  mm9.refseq.db <- makeTranscriptDbFromUCSC(genome="mm9", table="refGene")
#  refseq.genes = genes(mm9.refseq.db)
#  refseq.transcripts = transcriptsBy(mm9.refseq.db, by="gene")
#  refseq.transcripts = refseq.transcripts[ names(refseq.transcripts) %in% unlist(refseq.genes$gene_id) ]
#  mm9_refseq_promoters <- promoters(refseq.transcripts, 2500,2500)
#  mm9_refseq_promoters <- unlist(mm9_refseq_promoters[seqnames(mm9_refseq_promoters) %in% c("chr14", "chr15")])

## ------------------------------------------------------------------------
data("mm9_refseq_promoters")

annotation.features <- list(promoter = mm9_refseq_promoters)

annotateInteractions(hic_data_subset, annotation.features)
hic_data_subset

## ------------------------------------------------------------------------
anchorOne(hic_data_subset)
head(anchorOne(hic_data_subset)$promoter.id)

## ------------------------------------------------------------------------
table(anchorOne(hic_data_subset)$node.class)
table(anchorTwo(hic_data_subset)$node.class)

## ------------------------------------------------------------------------
plotInteractionAnnotations(hic_data_subset)

## ------------------------------------------------------------------------
length(hic_data_subset[isInteractionType(hic_data_subset, "promoter", "distal")])

## ------------------------------------------------------------------------
length(hic_data_subset[is.pd(hic_data_subset)])
sum(is.trans(hic_data_subset))

## ------------------------------------------------------------------------
hic_data_pd <- hic_data_subset[is.pd(hic_data_subset)]
max(count(hic_data_pd))
most_counts <- hic_data_pd[which.max(count(hic_data_pd))]
most_counts

## ------------------------------------------------------------------------
min(pValue(hic_data_pd))
min_pval <- hic_data_pd[which.min(pValue(hic_data_pd))]
min_pval

## ------------------------------------------------------------------------
calculateDistances(most_counts, method="midpoint")
calculateDistances(min_pval,method="midpoint")
summary(calculateDistances(hic_data_subset,method="midpoint"))

## ------------------------------------------------------------------------
viewPoint(mm9_refseq_promoters[1], hic_data_subset, leftflank = 2000000, rightflank = 2000000)

set.seed(42)
promoter_subset <- sample(mm9_refseq_promoters, 20)
viewPointAverage(promoter_subset, hic_data_subset, 
                             leftflank = 1000000, rightflank = 1000000)
viewPointAverage(promoter_subset, hic_data_subset, 
                             leftflank = 1000000, rightflank = 1000000, flip = TRUE)

## ----, eval=FALSE--------------------------------------------------------
#  ## Not run
#  export.bed12(hic_data_subset, fn="hic_data_FDR0.1.bed")