inst/doc/hic_vignette.R

In GenomicInteractions: Utilities for handling genomic interaction data

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

## ----imports, warning=F, results="hide", message=FALSE------------------------
library(Gviz)
library(GenomicInteractions)
library(GenomicRanges)
library(InteractionSet)


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

hic_data <- makeGenomicInteractionsFromFile(hic_file, 
                    type="homer", 
                    experiment_name = "HiC 100kb", 
                    description = "HiC 100kb resolution")
seqlengths(hic_data) <- c(chr15 = 103494974, chr14 = 125194864)

## -----------------------------------------------------------------------------
hic_data
mcols(hic_data)
head(hic_data$LogP)
hic_data$p.value <- exp(hic_data$LogP)

## -----------------------------------------------------------------------------
regions(hic_data)
anchors(hic_data, type = "first")
head(anchors(hic_data, type = "first", id = TRUE))
anchorOne(hic_data)

## -----------------------------------------------------------------------------
summary(width(regions(hic_data)))

## -----------------------------------------------------------------------------
head(interactionCounts(hic_data))
mean(interactionCounts(hic_data))

## -----------------------------------------------------------------------------
plot(density(hic_data$p.value))

hic_data$fdr <- hic_data$FDR.Benjamini..based.on.3.68e.08.total.tests.
plot(density(hic_data$fdr))

## -----------------------------------------------------------------------------
sum(hic_data$fdr < 0.1)
hic_data_subset <- hic_data[hic_data$fdr < 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 <- makeTxDbFromUCSC(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")])
#  mm9_refseq_promoters <- unique(mm9_refseq_promoters) # some duplicate promoters from different transcript isoforms
#  
#  #get gene symbols
#  mart = useMart("ensembl", dataset="mmusculus_gene_ensembl")
#  genes <- getBM(attributes = c("mgi_symbol", "refseq_mrna"), filter = "refseq_mrna",
#                 values = mm9_refseq_promoters$tx_name, mart = mart)
#  mm9_refseq_promoters$geneSymbol <- genes$mgi_symbol[match(mm9_refseq_promoters$tx_name, genes$refseq_mrna)]
#  
#  names(mm9_refseq_promoters) <- mm9_refseq_promoters$geneSymbol
#  na.symbol <- is.na(names(mm9_refseq_promoters))
#  names(mm9_refseq_promoters)[na.symbol] <- mm9_refseq_promoters$tx_name[na.symbol]

## ----eval=FALSE---------------------------------------------------------------
#  #Not run
#  
#  ## get enhancers from http://chromosome.sdsc.edu/mouse/download.html
#  download.file("http://chromosome.sdsc.edu/mouse/download/thymus.zip", "thymus.zip")
#  unzip("thymus.zip")
#  thymus_enh <- read.table("thymus/thymus.enhancer.txt", sep="\t", stringsAsFactors = FALSE)
#  thymus_enh <- GRanges(seqnames=thymus_enh$V1, ranges=IRanges(thymus_enh$V2, width=1))
#  thymus_enh <- resize(thymus_enh, fix="center", width=500)
#  thymus_enh <- thymus_enh[seqnames(thymus_enh) %in% c("chr14", "chr15")]
#  names(thymus_enh) <- paste("ENH", as.character(thymus_enh), sep = "_")

## -----------------------------------------------------------------------------
data("mm9_refseq_promoters")
data("thymus_enhancers")
annotation.features <- list(promoter = mm9_refseq_promoters, enhancer = thymus_enh)
annotateInteractions(hic_data_subset, annotation.features)

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

## -----------------------------------------------------------------------------
table(regions(hic_data_subset)$node.class)

## -----------------------------------------------------------------------------
plotInteractionAnnotations(hic_data_subset, legend = TRUE)

## -----------------------------------------------------------------------------
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_ep <- hic_data_subset[isInteractionType(hic_data_subset, "promoter", "enhancer")]

max(interactionCounts(hic_data_ep))
most_counts <- hic_data_ep[which.max(interactionCounts(hic_data_ep))]
most_counts

## -----------------------------------------------------------------------------
min(hic_data_ep$p.value)
min_pval <- hic_data_ep[which.min(hic_data_ep$p.value)]
min_pval

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

## -----------------------------------------------------------------------------
anchorOne(most_counts)$promoter.id
anchorTwo(most_counts)$enhancer.id

## -----------------------------------------------------------------------------
Trib1_region <- resize(mm9_refseq_promoters["Trib1"], fix = "center", width = 1000000)
interaction_track <- InteractionTrack(hic_data_subset, name = "HiC", chromosome = "chr15")
plotTracks(interaction_track, chromosome="chr15", 
           from=start(Trib1_region), to=end(Trib1_region))

## -----------------------------------------------------------------------------
promoterTrack <- AnnotationTrack(mm9_refseq_promoters, genome="mm9", name="Promoters",
                             id=names(mm9_refseq_promoters),  featureAnnotation="id")
enhTrack <- AnnotationTrack(thymus_enh, genome="mm9", name="Enhancers", stacking = "dense")

displayPars(promoterTrack) <- list(fill = "deepskyblue", col = NA, 
                                   fontcolor.feature = "black", fontsize=8,
                                   just.group="below")
displayPars(enhTrack) <- list(fill = "black", col = NA)
displayPars(interaction_track) = list(col.interactions="red", 
                                      col.anchors.fill ="blue",
                                      col.anchors.line = "black",
                                      interaction.dimension="height", 
                                      interaction.measure ="counts",
                                      plot.trans=FALSE,
                                      plot.outside = TRUE, 
                                      col.outside="lightblue", 
                                      anchor.height = 0.1)

plotTracks(list(interaction_track, promoterTrack, enhTrack),
           chromosome="chr15", from=start(Trib1_region), to=end(Trib1_region), 
           sizes=c(0.6, 0.2, 0.2))

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