## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = TRUE,
  warning = FALSE,
  message = FALSE,
  fig.align = "center",
  out.width = "80%"
)

## ----logo, echo=FALSE, eval=TRUE, out.width='10%'-----------------------------
knitr::include_graphics("../man/figures/SPICEY_LOGO.svg", dpi = 800)

## ----install, eval=FALSE, echo=TRUE-------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# 
# BiocManager::install("SPICEY")

## ----load-spicey--------------------------------------------------------------
library(SPICEY)

## ----da-atac, message=FALSE, warning=FALSE------------------------------------
data("atac")

## ----da-rna, message=FALSE, warning=FALSE-------------------------------------
data("rna")

## ----links, message=FALSE, warning=FALSE--------------------------------------
data("cicero_links")
head(cicero_links)

## ----quick-example, fig.width=2, fig.height=3---------------------------------
# Load necessary packages
library(dplyr)
library(GenomicRanges)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(org.Hs.eg.db)
library(SPICEY)

data(atac)
data(rna)
data(cicero_links)

# Annotate peaks to genes with coaccessibility
peaks <- unique(unlist(atac)[, c("region_id")])
annotation_coacc <- annotate_with_coaccessibility(
  peaks = peaks,
  txdb = TxDb.Hsapiens.UCSC.hg38.knownGene,
  links_df = cicero_links,
  annot_dbi = org.Hs.eg.db,
  protein_coding_only = TRUE,
  verbose = TRUE,
  add_tss_annotation = FALSE,
  upstream = 2000,
  downstream = 2000
)

# Calculate SPICEY measures and link them with coaccessibility
spicey_coacc <- SPICEY(
  rna = rna,
  atac = atac,
  annotation = annotation_coacc
)

# Plot results
spicey_heatmap(spicey_coacc$linked,
  spicey_measure = "SPICEY",
  combined_zscore = TRUE
)

## ----retsi, message=FALSE, warning=FALSE--------------------------------------
retsi <- SPICEY(atac = atac)
head(retsi)

## ----getsi, message=FALSE, warning=FALSE--------------------------------------
getsi <- SPICEY(rna = rna)
head(getsi)

## ----spicey, message=FALSE, warning=FALSE-------------------------------------
spicey <- SPICEY(
  atac = atac,
  rna = rna
)
lapply(spicey, head)

## ----re-gene-nearest, message=FALSE, warning=FALSE----------------------------
peaks <- unique(unlist(atac)[, c("region_id")])

annotation_near <- annotate_with_nearest(
  peaks = peaks,
  txdb = TxDb.Hsapiens.UCSC.hg38.knownGene,
  annot_dbi = org.Hs.eg.db,
  protein_coding_only = TRUE,
  verbose = TRUE,
  add_tss_annotation = FALSE,
  upstream = 2000,
  downstream = 2000
)

head(annotation_near)

## ----re-gene-coaccessibility, message=FALSE, warning=FALSE--------------------
peaks <- unique(unlist(atac)[, c("region_id")])

annotation_coacc <- annotate_with_coaccessibility(
  peaks = peaks,
  txdb = TxDb.Hsapiens.UCSC.hg38.knownGene,
  links_df = cicero_links,
  annot_dbi = org.Hs.eg.db,
  protein_coding_only = TRUE,
  verbose = TRUE,
  add_tss_annotation = FALSE,
  upstream = 2000,
  downstream = 2000
)

head(annotation_coacc)

## ----spicey_link_annot, message=FALSE, warning=FALSE--------------------------
# Using nearest gene annotation data
spicey_near <- SPICEY(
  rna = rna,
  atac = atac,
  annotation = annotation_near
)


# Using co-accessibility gene annotation data
spicey_coacc <- SPICEY(
  rna = rna,
  atac = atac,
  annotation = annotation_coacc
)

## ----spicey-heatmap-sep, message=FALSE, warning=FALSE, fig.width=2, fig.height=3----
# RETSI and/or GETSI separately
spicey_heatmap(spicey_coacc$GETSI,
  spicey_measure = "GETSI"
)

# If only RETSI is calculated (only ATAC data available), join with coaccessibility or nearest gene annotation to show gene_id in the plots
retsi <- spicey_coacc$RETSI
retsi <- retsi |> left_join(annotation_near, by = c("region_id"))
spicey_heatmap(retsi,
  spicey_measure = "RETSI"
)

## ----spicey-heatmap-pair, message=FALSE, warning=FALSE, fig.width=3.5, fig.height=4.5----
spicey_heatmap(spicey_coacc$linked,
  spicey_measure = "SPICEY",
  combined_zscore = FALSE
)

## ----spicey-heatmap-comb, message=FALSE, warning=FALSE, fig.width=2, fig.height=3----
# SPICEY: combined average RETSI and GETSI scores
spicey_heatmap(spicey_coacc$linked,
  spicey_measure = "SPICEY",
  combined_zscore = TRUE
)

## ----sessionInfo--------------------------------------------------------------
sessionInfo()