## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.wide = TRUE,
  fig.align = "center",
  fig.retina = 2,
  warning = FALSE,
  message = FALSE,
  cache = FALSE,
  echo = TRUE
)
suppressPackageStartupMessages({
  library(BiocStyle)
})

## ----overview-diagram, out.width = '70%', echo = FALSE, fig.alt="Overview diagram of the linkSet class structure and its components"----
knitr::include_graphics("img/overview.png")

## ----methods-diagram, out.width = '80%', echo = FALSE, fig.alt="Overview diagram of the linkSet methods"----
knitr::include_graphics("img/linksetmethod.png")

## ----data-structure-diagram, out.width = '70%', echo = FALSE, fig.alt="data structure of the linkSet class"----
knitr::include_graphics("img/9.10_data_structure.png")

## ----import-methods-diagram, out.width = '70%', echo = FALSE, fig.alt="Import methods"----
knitr::include_graphics("img/linkSetImport.png")

## ----installation, eval = FALSE-----------------------------------------------
# if (!require("BiocManager", quietly = TRUE)) {
#   install.packages("BiocManager")
# }

## ----load-linkset-------------------------------------------------------------
suppressPackageStartupMessages(library(linkSet))

## ----construct-from-granges---------------------------------------------------
library(GenomicRanges)
gr1 <- GRanges(
  seqnames = c("chr1", "chr1", "chr2"),
  ranges = IRanges(start = c(1, 100, 200), width = 10),
  strand = "+", symbol = c("Gene1", "Gene2", "Gene3")
)
gr2 <- GRanges(
  seqnames = c("chr1", "chr2", "chr2"),
  ranges = IRanges(start = c(50, 150, 250), width = 10),
  strand = "+"
)
## construct linkSet
linkObj <- linkSet(gr1, gr2, specificCol = "symbol")
linkObj

## ----construct-from-ginteractions---------------------------------------------
library(InteractionSet)
gi <- GInteractions(
  anchor1 = c(1, 3, 5),
  anchor2 = c(2, 4, 6),
  regions = GRanges(
    seqnames = c("chr1", "chr1", "chr2", "chr2", "chr3", "chr3"),
    ranges = IRanges(start = c(1, 50, 100, 150, 200, 250), width = 10)
  )
)
mcols(gi)$score <- c(0.1, 0.2, 0.3)
mcols(gi)$gene <- c("Sox2", "Sox9", "Sox10")
## Convert to linkSet
linkObj_from_gi <- Convert(gi, baitCol = "gene")
linkObj_from_gi

## ----construct-with-gene-enhancer---------------------------------------------
geneGr <- GRanges(
  seqnames = c("chr1", "chr2", "chr3"),
  ranges = IRanges(start = c(5, 105, 205), width = 20),
  symbol = c("Gene1", "Gene2", "Gene3")
)

peakGr <- GRanges(
  seqnames = c("chr1", "chr2", "chr3"),
  ranges = IRanges(start = c(45, 145, 245), width = 10)
)

# Run baitGInteractions
linkObj_from_gi_2 <- baitGInteractions(gi, geneGr, peakGr, geneSymbol = "symbol")

linkObj_from_gi_2

## ----construct-from-dataframe-------------------------------------------------
test_df <- data.frame(
  gene = c("gene1", "gene2", "gene3"),
  peak = c("chr1:1000-2000", "chr2:3000-4000", "chr3:5000-6000"),
  score = c(0.5, 0.7, 0.9),
  stringsAsFactors = FALSE
)

# Test successful conversion
linkObj_from_df <- Convert(test_df)
linkObj_from_df

## ----get-regions--------------------------------------------------------------
linkSet::regions(linkObj)

## ----get-regions-bait---------------------------------------------------------
regionsBait(linkObj)

## ----get-oe-------------------------------------------------------------------
oe(linkObj)

## ----get-bait-----------------------------------------------------------------
bait(linkObj)

## ----get-mcols----------------------------------------------------------------
mcols(linkObj)

## ----set-new-bait-oe----------------------------------------------------------
new_bait <- c("Gene40", "Gene41", "Gene42")
new_oe <- GRanges(
  seqnames = c("chr1", "chr2", "chr3"),
  ranges = IRanges(start = c(5, 105, 205), width = 20),
  symbol = c("Gene1", "Gene2", "Gene3")
)
bait(linkObj) <- new_bait
oe(linkObj) <- new_oe
linkObj

## ----subset-by-index----------------------------------------------------------
linkset_sub <- linkObj[1:2]
linkset_sub

## ----subset-by-bait-----------------------------------------------------------
linkset_sub_2 <- subsetBait(linkObj, "Gene1")
linkset_sub_2
linkset_sub_3 <- subsetBaitRegion(linkObj, "chr1:100-200")
linkset_sub_3

## ----concatenate-linksets-----------------------------------------------------
linkset_concat <- c(linkObj, linkObj)
linkset_concat

## ----resize-regions-----------------------------------------------------------
data(linkExample)
linkExample

## resize the bait region
resize_bait <- resizeRegions(linkExample, width = 75, fix = "start", region = "bait")
resize_bait

## resize the oe region
resize_oe <- resizeRegions(linkExample, width = 75, fix = "start", region = "oe")
resize_oe

## ----reduce-linkset-----------------------------------------------------------
reduce_linkset <- reduce(linkset_concat)
reduce_linkset

## ----reduce-no-count----------------------------------------------------------
reduce_linkset_2 <- linkSet::reduceRegions(linkset_concat, countInteractions = FALSE)
reduce_linkset_2

## ----diagnose-linkset---------------------------------------------------------
diagnoseLinkSet(linkExample)

## ----show-linkset-concat------------------------------------------------------
linkset_concat

## ----filter-intra-------------------------------------------------------------
linkset_concat <- filterLinks(linkset_concat, filter_intra = TRUE)
linkset_concat

## ----annotate-promoter--------------------------------------------------------
gr1 <- GRanges(
  seqnames = c("chr1", "chr3", "chr3"),
  ranges = IRanges(start = c(1000, 2000, 3000), width = 100),
  strand = "+", symbol = c("BRCA1", "TP53", "NONEXISTENT")
)
gr2 <- GRanges(
  seqnames = c("chr1", "chr2", "chr3"),
  ranges = IRanges(start = c(5000, 6000, 7000), width = 100),
  strand = "+"
)
linkObj <- linkSet(gr1, gr2, specificCol = "symbol")

# Test annotatePromoter
annotated_linkset <- suppressWarnings(annotatePromoter(linkObj, genome = "hg38", upstream = 500, overwrite = TRUE))
annotated_linkset

## ----cross-gene-enhancer------------------------------------------------------
annotated_linkset <- crossGeneEnhancer(annotated_linkset)
annotated_linkset <- orderLinks(annotated_linkset, by = "crossFreq", decreasing = TRUE)
annotated_linkset

## ----chicane-analysis---------------------------------------------------------
annotated_linkset <- countInteractibility(annotated_linkset)
annotated_linkset <- linkSet::pairdist(annotated_linkset)
annotated_linkset <- run_chicane(annotated_linkset)
annotated_linkset

## ----plot-genomic-ranges, eval=FALSE------------------------------------------
# # Note: visualization requires annotated bait regions
# # linkExample doesn't have regionsBait annotated, so this will be skipped
# plot_genomic_ranges(linkExample, extend.base = 10)

## ----plot-genomic-ranges-custom, eval=FALSE-----------------------------------
# # Note: visualization requires annotated bait regions
# plot_genomic_ranges(linkExample, extend.base = 10, x.range = c(50, 200))

## ----session-info-------------------------------------------------------------
sessionInfo()