## ----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)
})

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

## ----load-packages------------------------------------------------------------
suppressPackageStartupMessages({
  library(linkSet)
  # Load additional packages only if available
  if (requireNamespace("TxDb.Mmusculus.UCSC.mm10.knownGene", quietly = TRUE)) {
    library(TxDb.Mmusculus.UCSC.mm10.knownGene)
  }
  if (requireNamespace("org.Mm.eg.db", quietly = TRUE)) {
    library(org.Mm.eg.db)
  }
  if (requireNamespace("Organism.dplyr", quietly = TRUE)) {
    library(Organism.dplyr)
  }
  if (requireNamespace("InteractionSet", quietly = TRUE)) {
    library(InteractionSet)
  }
  if (requireNamespace("rtracklayer", quietly = TRUE)) {
    library(rtracklayer)
  }
})

## ----setup-data---------------------------------------------------------------
hic_file <- system.file("extdata", "toyExample.pair.gz",
  package = "linkSet"
)
gi <- readvalidPairs(hic_file, format = "pair")
promoterGr <- tryCatch(
  {
    withTxDb("mm10", function(src) {
      genes <- Organism.dplyr::genes(src, columns = "symbol")
      IRanges::promoters(genes, upstream = 10000)
    })
  },
  error = function(e) {
    # Fallback approach using TxDb directly
    if (requireNamespace("TxDb.Mmusculus.UCSC.mm10.knownGene", quietly = TRUE) &&
      requireNamespace("org.Mm.eg.db", quietly = TRUE)) {
      txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene::TxDb.Mmusculus.UCSC.mm10.knownGene
      gene_ranges <- GenomicFeatures::genes(txdb)

      # Add gene symbols
      gene_ids <- names(gene_ranges)
      symbols <- AnnotationDbi::select(org.Mm.eg.db::org.Mm.eg.db,
        keys = gene_ids,
        columns = "SYMBOL",
        keytype = "ENTREZID"
      )
      symbol_map <- setNames(symbols$SYMBOL, symbols$ENTREZID)
      mcols(gene_ranges)$symbol <- symbol_map[names(gene_ranges)]
      genes <- gene_ranges[!is.na(mcols(gene_ranges)$symbol)]

      IRanges::promoters(genes, upstream = 10000)
    } else {
      stop("Required packages not available for mm10 annotation")
    }
  }
)
file_path <- system.file("extdata", "Embryo_body.bed.gz", package = "linkSet")
enhancer <- rtracklayer::import(file_path, format = "BED")

## ----create-linkset-----------------------------------------------------------
gi <- resize(gi, 10000, fix = "center")
ls <- baitGInteractions(gi, promoterGr, enhancer, geneSymbol = "symbol")
ls

## ----show-linkset-------------------------------------------------------------
showLinkSet(ls, baitRegion = TRUE)

## ----annotate-and-diagnose----------------------------------------------------
ls <- suppressWarnings(annotatePromoter(ls, genome = "mm10", upstream = 1000))
diagnoseLinkSet(ls)

## ----filter-links-------------------------------------------------------------
ls <- countInteractibility(ls)
ls <- filterLinks(ls, filter_intra = TRUE, filter_unannotate = TRUE, distance = 50000000)

## ----count-and-order----------------------------------------------------------
ls <- countInteractions(ls)
orderLinks(ls, by = "count", decreasing = TRUE)

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

## ----plot-genomic-ranges------------------------------------------------------
plot_genomic_ranges(ls, showOE = oe(ls)[1])

## ----plot-bait-region---------------------------------------------------------
# Check if regionsBait exists and has names
if (!is.null(regionsBait(ls)) && !is.null(names(regionsBait(ls)))) {
  # Get available bait names from regionsBait
  available_baits <- unique(names(regionsBait(ls)))
  if (length(available_baits) > 0) {
    # Use the first available bait for demonstration
    plot_genomic_ranges(ls, showBait = available_baits[1])
  } else {
    cat("No named bait regions available for plotting")
  }
} else {
  # If no regionsBait, just plot without showBait
  plot_genomic_ranges(ls)
}

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