glyfun is a package
for functional enrichment analysis of glycoproteomics data. Why does
glycoproteomics need a dedicated tool? For transcriptomics and
proteomics, functional enrichment analysis is usually fairly
straightforward. For glycoproteomics, however, we need to account for
both macro- and micro-heterogeneity. Macro-heterogeneity means that one
protein can carry multiple glycosylation sites, while
micro-heterogeneity means that each site can carry multiple glycan
structures. As a result, a single protein may have several dysregulated
glycosylation events across different sites and glycans. That extra
layer of complexity makes traditional enrichment workflows harder to
apply directly.
In general, glyfun answers two questions:
The following sections walk through both perspectives and show when each one is useful.
library(glyfun)
library(SummarizedExperiment)
#> Loading required package: MatrixGenerics
#> Loading required package: matrixStats
#>
#> Attaching package: 'MatrixGenerics'
#> The following objects are masked from 'package:matrixStats':
#>
#> colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
#> colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
#> colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
#> colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
#> colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
#> colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
#> colWeightedMeans, colWeightedMedians, colWeightedSds,
#> colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
#> rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
#> rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
#> rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
#> rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
#> rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
#> rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
#> rowWeightedSds, rowWeightedVars
#> Loading required package: GenomicRanges
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> Loading required package: generics
#>
#> Attaching package: 'generics'
#> The following objects are masked from 'package:base':
#>
#> as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
#> setequal, union
#>
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#>
#> IQR, mad, sd, var, xtabs
#> The following object is masked from 'package:utils':
#>
#> data
#> The following objects are masked from 'package:base':
#>
#> anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#> colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#> get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget,
#> order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
#> rbind, Reduce, rownames, sapply, saveRDS, scale, sequence, table,
#> tapply, transform, unique, unsplit, which.max, which.min
#> Loading required package: S4Vectors
#>
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:utils':
#>
#> findMatches
#> The following objects are masked from 'package:base':
#>
#> expand.grid, I, unname
#> Loading required package: IRanges
#> Loading required package: Seqinfo
#> Loading required package: Biobase
#> Welcome to Bioconductor
#>
#> Vignettes contain introductory material; view with
#> 'browseVignettes()'. To cite Bioconductor, see
#> 'citation("Biobase")', and for packages 'citation("pkgname")'.
#>
#> Attaching package: 'Biobase'
#> The following object is masked from 'package:MatrixGenerics':
#>
#> rowMedians
#> The following objects are masked from 'package:matrixStats':
#>
#> anyMissing, rowMediansPrerequisites: To follow this vignette, you should
be familiar with the core glycoverse packages,
especially glyexp and
glystats. For the
glycan-centric enrichment analysis in the second half, you will also
need some familiarity with glydet. Experience with
clusterProfiler is helpful too, since glyfun
returns standard enrichment result objects that you can inspect and
visualize with that ecosystem.
Let’s start with the first question. We first run differential analysis on the glycoforms.
library(glyexp)
#>
#> Attaching package: 'glyexp'
#> The following object is masked from 'package:Biobase':
#>
#> samples
library(glyclean)
#>
#> Attaching package: 'glyclean'
#> The following object is masked from 'package:S4Vectors':
#>
#> aggregate
#> The following object is masked from 'package:stats':
#>
#> aggregate
library(glystats)
# `real_experiment` is a bundled GlycoproteomicSE.
gp_se <- real_experiment |>
auto_clean() |>
filter_col(group %in% c("H", "C")) |>
mutate_col(group = droplevels(group))
#>
#> ── Removing variables with too many missing values ──
#>
#> ℹ Applying preset "discovery"...
#> ℹ Total removed: 24 (0.56%) variables.
#> ✔ Variable removal completed.
#>
#> ── Normalizing data ──
#>
#> ℹ Normalization method: `normalize_median()`
#> ℹ Reason: default for "glycoproteomics".
#> ✔ Normalization completed.
#>
#> ── Imputing missing values ──
#>
#> ℹ Imputation method: `impute_min_prob()`
#> ℹ Reason: default for "glycoproteomics" with n_samples < 30.
#> ✔ Imputation completed.
#>
#> ── Aggregating data ──
#>
#> ℹ Aggregating to "gfs" level
#> ✔ Aggregation completed.
#>
#> ── Normalizing data again ──
#>
#> ℹ Normalization method: `normalize_median()`
#> ℹ Reason: default for "glycoproteomics".
#> ✔ Normalization completed.
#>
#> ── Correcting batch effects ──
#>
#> ℹ Batch column batch not found in sample_info. Skipping batch correction.
#> ✔ Batch correction completed.
dea_res <- gly_limma(gp_se)
#> ℹ Ref Group: "H"
#> ℹ Test Group: "C"
get_tidy_result(dea_res)
#> # A tibble: 3,979 × 14
#> variable protein glycan_composition glycan_structure protein_site gene
#> <chr> <chr> <comp> <struct> <int> <chr>
#> 1 P08185-176-He… P08185 Hex(5)HexNAc(4)Ne… NeuAc(??-?)Hex(… 176 SERP…
#> 2 P04196-344-He… P04196 Hex(5)HexNAc(4)Ne… NeuAc(??-?)Hex(… 344 HRG
#> 3 P04196-344-He… P04196 Hex(5)HexNAc(4) Hex(??-?)HexNAc… 344 HRG
#> 4 P04196-344-He… P04196 Hex(5)HexNAc(4)Ne… NeuAc(??-?)Hex(… 344 HRG
#> 5 P10909-291-He… P10909 Hex(6)HexNAc(5) Hex(??-?)HexNAc… 291 CLU
#> 6 P04196-344-He… P04196 Hex(5)HexNAc(4)Ne… NeuAc(??-?)Hex(… 344 HRG
#> 7 P04196-345-He… P04196 Hex(5)HexNAc(4) Hex(??-?)HexNAc… 345 HRG
#> 8 P04196-344-He… P04196 Hex(5)HexNAc(4)dH… dHex(??-?)Hex(?… 344 HRG
#> 9 P04196-344-He… P04196 Hex(4)HexNAc(3) Hex(??-?)HexNAc… 344 HRG
#> 10 P04196-344-He… P04196 Hex(4)HexNAc(4)Ne… NeuAc(??-?)Hex(… 344 HRG
#> # ℹ 3,969 more rows
#> # ℹ 8 more variables: log2fc <dbl>, AveExpr <dbl>, t <dbl>, p_val <dbl>,
#> # p_adj <dbl>, b <dbl>, ref_group <chr>, test_group <chr>Now that we have differential results for glycoforms, we can ask which functions are enriched among proteins with dysregulated glycosylation. Here, a protein is included if it has at least one dysregulated glycoform. The workflow has three steps:
glyfun wraps these steps in a set of convenience
functions:
# Gene Ontology (GO) enrichment analysis
enrich_ora_go(
dea_res,
dea_p_cutoff = 0.05,
dea_log2fc_cutoff = c(-1, 1)
)
# KEGG pathway enrichment analysis
enrich_ora_kegg(dea_res)
# Reactome pathway enrichment analysis
enrich_ora_reactome(dea_res)
# WikiPathways enrichment analysis
enrich_ora_wp(dea_res)
# Disease Ontology (DO) enrichment analysis
enrich_ora_do(dea_res)
# Network of Cancer Genes (NCG) enrichment analysis
enrich_ora_ncg(dea_res)Use dea_p_cutoff and dea_log2fc_cutoff to
define which glycoforms count as dysregulated. The two values in
dea_log2fc_cutoff are the lower and upper boundaries;
c(-1, 1) selects glycoforms with log2 fold changes below -1
or above 1. p_cutoff separately filters the enriched terms
returned by the downstream enrichment method.
Under the hood, these functions call the corresponding
clusterProfiler functions and return an
enrichResult object. That means you can convert the result
to a tibble with as_tibble(), or visualize it with
dotplot().
You can also perform Gene Set Enrichment Analysis (GSEA). Unlike ORA, GSEA does not require a hard cutoff to define dysregulated glycoforms. Instead, it works from a ranked protein list, often based on statistics such as log2FC.
This raises one practical question: one protein can have multiple
glycoforms with different log2FC values, so how should it be ranked?
glyfun handles this by aggregating glycoform-level
statistics into one value per protein. The default is the median, but
you can also choose other methods such as the mean or maximum.
# GSEA for Gene Ontology (GO)
enrich_gsea_go(dea_res, aggr = "median")
# GSEA for KEGG pathway
enrich_gsea_kegg(dea_res, aggr = "median")
# GSEA for Reactome pathway
enrich_gsea_reactome(dea_res, aggr = "median")
# GSEA for WikiPathways
enrich_gsea_wp(dea_res, aggr = "median")
# GSEA for Disease Ontology (DO)
enrich_gsea_do(dea_res, aggr = "median")
# GSEA for Network of Cancer Genes (NCG)
enrich_gsea_ncg(dea_res, aggr = "median")These functions also support several ranking methods. By default,
glyfun uses "signed_log10p", which combines
the sign of log2FC with the p-value. Other options include
"log2fc", "abs_log2fc", "log10p",
and "log2fc_log10p". See the GSEA function documentation
for the full details.
Note: You should be familiar with the glydet package to understand this section.
The protein-centric approach is useful and easy to interpret, but it hides important glycan-level information. Saying that a protein has “dysregulated glycosylation” can be too broad. Often, we want a more specific question: which functions are associated with particular glycosylation changes? For example, we might ask which functions are linked to increased fucosylation or decreased sialylation.
To do that, we first calculate glycosylation traits with
glydet.
library(glydet)
trait_se <- derive_traits(gp_se)
trait_info <- tibble::as_tibble(rowData(trait_se), rownames = "variable")
trait_info |>
dplyr::distinct(trait, explanation)
#> # A tibble: 14 × 2
#> trait explanation
#> <chr> <chr>
#> 1 TM Proportion of high-mannose glycans among all glycans.
#> 2 TH Proportion of hybrid glycans among all glycans.
#> 3 TC Proportion of complex glycans among all glycans.
#> 4 MM Abundance-weighted mean of mannose count within high-mannose glycans.
#> 5 CA2 Proportion of bi-antennary glycans within complex glycans.
#> 6 CA3 Proportion of tri-antennary glycans within complex glycans.
#> 7 CA4 Proportion of tetra-antennary glycans within complex glycans.
#> 8 TF Proportion of fucosylated glycans among all glycans.
#> 9 TFc Proportion of core-fucosylated glycans among all glycans.
#> 10 TFa Proportion of arm-fucosylated glycans among all glycans.
#> 11 TB Proportion of glycans with bisecting GlcNAc among all glycans.
#> 12 GS Abundance-weighted mean of degree of sialylation per galactose among a…
#> 13 AG Abundance-weighted mean of degree of galactosylation per antenna among…
#> 14 TS Proportion of sialylated glycans among all glycans.The result contains several trait types, including TF
for fucosylation and GS for sialylation. These traits are
calculated site by site, so each glycosylation site on each protein has
its own trait values.
Next, we run differential analysis on the trait experiment.
trait_dea_res <- gly_limma(trait_se)
#> ℹ Ref Group: "H"
#> ℹ Test Group: "C"
#> Warning in splines::ns(covariate, df = splinedf, intercept = TRUE): shoving
#> 'interior' knots matching boundary knots to insideAs before, ORA uses hard cutoffs to identify dysregulated features.
The difference is that enrichment is now performed separately for each
trait. For example, for the TF trait, we filter
statistically significant glycosites and then run enrichment analysis on
the proteins with dysregulated TF. The same process is
repeated for GS and the other traits. In
glyfun, this is called glycan-centric enrichment.
glyfun provides the following functions for
glycan-centric ORA:
# Glycan-centric ORA for Gene Ontology (GO)
enrich_gc_ora_go(trait_dea_res)
# Glycan-centric ORA for KEGG pathway
enrich_gc_ora_kegg(trait_dea_res)
# Glycan-centric ORA for Reactome pathway
enrich_gc_ora_reactome(trait_dea_res)
# Glycan-centric ORA for WikiPathways
enrich_gc_ora_wp(trait_dea_res)
# Glycan-centric ORA for Disease Ontology (DO)
enrich_gc_ora_do(trait_dea_res)
# Glycan-centric ORA for Network of Cancer Genes (NCG)
enrich_gc_ora_ncg(trait_dea_res)These functions return compareClusterResult objects,
which can also be visualized with dotplot().
Trait differential results are not the only valid input for glycan-centric enrichment. You can also use raw glycoform experiments or motif quantification experiments. With raw glycoform experiments, enrichment is performed for each glycan composition or structure; with motif quantification experiments, it is performed for each motif. Because this can trigger many enrichment analyses, raw glycoform-level input can be slow and is usually not the first option to try.
GSEA can also be performed in a glycan-centric manner. The idea mirrors glycan-centric ORA: rank proteins separately for each trait, then run GSEA for each ranked list.
# Glycan-centric GSEA for Gene Ontology (GO)
enrich_gc_gsea_go(trait_dea_res, aggr = "median")
# Glycan-centric GSEA for KEGG pathway
enrich_gc_gsea_kegg(trait_dea_res, aggr = "median")
# Glycan-centric GSEA for Reactome pathway
enrich_gc_gsea_reactome(trait_dea_res, aggr = "median")
# Glycan-centric GSEA for WikiPathways
enrich_gc_gsea_wp(trait_dea_res, aggr = "median")
# Glycan-centric GSEA for Disease Ontology (DO)
enrich_gc_gsea_do(trait_dea_res, aggr = "median")
# Glycan-centric GSEA for Network of Cancer Genes (NCG)
enrich_gc_gsea_ncg(trait_dea_res, aggr = "median")Try these workflows on your own data and see which biological patterns emerge.