## ----style, echo = FALSE, results = 'asis'------------------------------------
BiocStyle::markdown()

## ----install, eval=FALSE------------------------------------------------------
# if (!require("BiocManager", quietly = TRUE)) {
#   install.packages("BiocManager")
# }
# BiocManager::install("matsui-lab/glycoTraitR")

## ----load---------------------------------------------------------------------
library(glycoTraitR)

## ----load-example-data--------------------------------------------------------
# The following is an example file from Zenodo
# url <- "https://zenodo.org/records/17756830/files/pGlycoDB-GP-FDR-Pro-Quant-Site.txt?download=1"
# download.file(url, "pGlyco3_GPSM.txt", mode = "wb")
# gpsm <- read_pGlyco3_gpsm("pGlyco3_GPSM.txt")

# meta_url <- "https://zenodo.org/records/17759790/files/meta_mcp_2022_100433.rds?download=1"
# download.file(meta_url, "meta_example.rds", mode = "wb")
# meta <- readRDS("meta_example.rds")

# The following is a smaller toy example files from package
path <- system.file("extdata", "pGlyco3_gpsm_toyexample.txt", package = "glycoTraitR")
gpsm_toyexample <- read_pGlyco3_gpsm(path)
data("meta_toyexample")

head(gpsm_toyexample)
head(meta_toyexample)

## ----build-trait--------------------------------------------------------------
trait_se <- build_trait_se(
  gpsm_toyexample,
  from = "pGlyco3",
  motifs = NULL,
  level = "protein",
  meta = meta_toyexample
)
trait_se

## ----DA-trait-----------------------------------------------------------------
changed_traits <- analyze_trait_changes(
  trait_se     = trait_se,
  group_col    = "Diagnosis",
  group_levels = c("Normal", "Symptomatic"),
  min_psm      = 20
)
head(changed_traits)

## ----plot-trait---------------------------------------------------------------
p <- plot_trait_distribution(
  trait_se     = trait_se,
  group_col    = "Diagnosis",
  group_levels = c("Normal", "Symptomatic"),
  trait_name   = "Hexose",
  feature      = "MOG"
)

p$p_hist
p$p_box

## ----load-igraph, message=FALSE-----------------------------------------------
library(igraph)

## ----example-pGlyco3----------------------------------------------------------
# Example: parsed glycan trees from existing GPSM data
path <- system.file("extdata", "pGlyco3_gpsm_toyexample.txt", package = "glycoTraitR")
gpsm_toyexample <- read_pGlyco3_gpsm(path)
data("meta_toyexample")
glycans <- lapply(gpsm_toyexample$GlycanStructure, glycoTraitR:::pGlyco3_to_tree)

## ----visualize-tree-----------------------------------------------------------
# visualize example structures
igraph_obj <- build_glycan_igraph(glycans[[1]])
plot_glycan_tree(igraph_obj)

## ----example1-----------------------------------------------------------------
#### Example 1 ####
# (1) Select a glycan from pGlyco3 data
tree_full <- glycans[[24]]
g_full <- build_glycan_igraph(tree_full)

# (2) Define motif manually
motif <- list(
  node = c("H", "H", "H"),
  edge = c("a-b", "b-c")
)
g_motif <- build_glycan_igraph(motif)

# (3) Perform subgraph matching
n_match <- count_subgraph_isomorphisms(
  g_motif, g_full,
  method = "vf2",
  vertex.color1 = factor(V(g_full)$type),
  vertex.color2 = factor(V(g_motif)$type)
)

# (4) Visualize
par(mfrow = c(1, 2))
plot_glycan_tree(g_motif)
plot_glycan_tree(g_full)
print(paste(n_match, "motif(s) found in the glycan"))

## ----example2-----------------------------------------------------------------
#### Example 2 ####
# (1) Select a glycan
tree_full <- glycans[[1]]
g_full <- glycoTraitR::build_glycan_igraph(tree_full)

# (2) Symbolic motif definition
motif <- list(
  node = c("H", "N", "F"),
  edge = c("a-b", "b-c")
)
g_motif <- glycoTraitR::build_glycan_igraph(motif)

# (3) Subgraph matching
n_match <- igraph::count_subgraph_isomorphisms(
  g_motif, g_full,
  method = "vf2",
  vertex.color1 = factor(V(g_full)$type),
  vertex.color2 = factor(V(g_motif)$type)
)

# (4) Visualize motif vs. full glycan
par(mfrow = c(1, 2))
plot_glycan_tree(g_motif)
plot_glycan_tree(g_full)
print(paste(n_match, "motif(s) found in the glycan"))

## ----udmotifs-----------------------------------------------------------------
user_motifs <- list(
  LinearH3 = list(
    node = c("H", "H", "H"),
    edge = c("a-b", "b-c")
  ),
  FucBranch = list(
    node = c("H", "N", "F"),
    edge = c("a-b", "b-c")
  )
)

## ----se-show------------------------------------------------------------------
trait_se <- build_trait_se(
  gpsm_toyexample,
  from = "pGlyco3",
  motifs = user_motifs,
  level = "protein",
  meta = meta_toyexample
)
SummarizedExperiment::assayNames(trait_se)

## ----session-info, echo=FALSE-------------------------------------------------
sessionInfo()