library(tidyverse)
#> ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
#> ✔ dplyr 1.2.1 ✔ readr 2.2.0
#> ✔ forcats 1.0.1 ✔ stringr 1.6.0
#> ✔ ggplot2 4.0.3 ✔ tibble 3.3.1
#> ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
#> ✔ purrr 1.2.2
#> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ dplyr::lag() masks stats::lag()
#> ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsCase Study: Glycomics
An end-to-end glycomics analysis with Glycoverse.
This vignette walks you through a complete glycomics analysis using glycoverse. We’ll explore the full spectrum of glycomics data analysis, from data loading and preprocessing to statistical analysis and visualization. We’ll also dive into advanced glycan structure analysis, including motif quantification and derived trait analysis. Ready to dive in? Let’s go!
Heads up: glycoverse is built on tidy principles throughout. If you’re new to tidyverse data analysis, we highly recommend checking out Hadley Wickham’s excellent R for Data Science. Trust us, it’s worth the investment!
Quick readiness check:
- What’s a
tibble? - How do you filter rows in a
tibble? - What’s the modern alternative to
forloops? - What’s the
|>operator? - What makes data “tidy”?
TL;DR
In case you’re in a hurry…
# Load the packages
library(tidyverse)
library(glycoverse)
library(SummarizedExperiment)
# Preprocess the data
clean_exp <- auto_clean(real_experiment2)
# Perform PCA
pca_res <- gly_pca(clean_exp)
autoplot(pca_res)
# Perform differential expression analysis
anova_res <- gly_anova(clean_exp)
get_tidy_result(anova_res, "main_test")
# Perform motif analysis
motifs <- c(
motif1 = "Neu5Ac(??-?)Gal(??-?)GlcNAc(??-",
motif2 = "Gal(??-?)GlcNAc(??-",
motif3 = "GlcNAc(??-"
)
motif_exp <- quantify_motifs(clean_exp, motifs)
motif_anova_res <- gly_anova(motif_exp)
get_tidy_result(motif_anova_res, "main_test")
# Perform derived trait analysis
trait_exp <- derive_traits(clean_exp)
trait_anova_res <- gly_anova(trait_exp)
get_tidy_result(trait_anova_res, "main_test")Loading the Packages
We first load the tidyverse package, as usual.
Just like tidyverse, glycoverse is a meta-package that loads a collection of specialized packages all at once.
library(glycoverse)
#> ── Attaching core glycoverse packages ────────────────────── glycoverse 0.3.1 ──
#> ✔ glyclean 0.15.1 ✔ glyparse 0.7.1
#> ✔ glydet 0.12.1 ✔ glyread 0.12.0
#> ✔ glydraw 0.6.3 ✔ glyrepr 0.13.0
#> ✔ glyexp 0.16.0 ✔ glystats 0.11.1
#> ✔ glymotif 0.17.1 ✔ glyvis 0.7.1
#> ── Conflicts ───────────────────────────────────────── glycoverse_conflicts() ──
#> ✖ glyclean::aggregate() masks stats::aggregate()
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ lubridate::intersect() masks dplyr::intersect(), base::intersect()
#> ✖ dplyr::lag() masks stats::lag()
#> ✖ lubridate::setdiff() masks dplyr::setdiff(), base::setdiff()
#> ✖ dplyr::setequal() masks base::setequal()
#> ✖ lubridate::union() masks dplyr::union(), base::union()
#> ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsThe Glycoverse data containers extend SummarizedExperiment, so we load it to use its standard accessors.
library(SummarizedExperiment)
#> Loading required package: MatrixGenerics
#> Loading required package: matrixStats
#>
#> Attaching package: 'matrixStats'
#> The following object is masked from 'package:dplyr':
#>
#> count
#>
#> 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 object is masked from 'package:lubridate':
#>
#> as.difftime
#> The following object is masked from 'package:dplyr':
#>
#> explain
#> 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 object is masked from 'package:dplyr':
#>
#> combine
#> The following objects are masked from 'package:stats':
#>
#> IQR, mad, sd, var, xtabs
#> 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, table, tapply, unique,
#> unsplit, which.max, which.min
#> Loading required package: S4Vectors
#>
#> Attaching package: 'S4Vectors'
#> The following objects are masked from 'package:lubridate':
#>
#> second, second<-
#> The following objects are masked from 'package:dplyr':
#>
#> first, rename
#> The following object is masked from 'package:tidyr':
#>
#> expand
#> 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
#>
#> Attaching package: 'IRanges'
#> The following object is masked from 'package:lubridate':
#>
#> %within%
#> The following objects are masked from 'package:dplyr':
#>
#> collapse, desc, slice
#> The following object is masked from 'package:purrr':
#>
#> reduce
#> 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, rowMedians
#> The following object is masked from 'package:glyexp':
#>
#> samplesReading the Data
Data import is typically your first step in any analysis. For this tutorial, we’ll use the real_experiment2 dataset that comes with glyexp. This is a real-world N-glycomics dataset from 144 patients across four liver conditions: healthy controls (H), hepatitis (M), cirrhosis (Y), and hepatocellular carcinoma (C).
real_experiment2
#>
#> ── GlycomicSE ──────────────────────────────────────────────────────────────────
#> ℹ Abundance assay: 144 samples, 67 variables
#> ℹ Glycan type: N
#> ℹ Row data fields: glycan_composition <comp>, glycan_structure <struct>
#> ℹ Column data fields: group <fct>
#> ℹ Metadata fields: exp_type <chr>, glycan_type <chr>For your own glycomics projects, create a GlycomicSE object. See Get Started with glyexp for details.
The real_experiment2 object is a GlycomicSE object, a specialized SummarizedExperiment for glycomics data. It manages three key components:
- Abundance assay: quantitative data with samples as columns and variables as rows
- Column data: sample metadata (group, batch, demographics, etc.)
- Row data: feature metadata (glycan compositions, glycan structures, etc.)
You can access these components with assay(), colData(), and rowData().
assay(real_experiment2)[1:5, 1:5]
#> S1 S2 S3 S4 S5
#> Man(3)GlcNAc(3) 1354.352 1884.387 1389.444 2034.693 1472.504
#> Man(3)GlcNAc(7) 3315.779 2500.308 1247.036 3102.668 2903.602
#> Man(5)GlcNAc(2) 6940.940 5911.016 1686.319 4071.061 4349.991
#> Man(4)Gal(2)GlcNAc(4)Neu5Ac(2) 4437.816 7535.886 2053.077 4773.249 3142.817
#> Man(3)Gal(1)GlcNAc(3) 1346.274 1663.375 1043.464 1765.550 1184.083as_tibble(colData(real_experiment2))
#> # A tibble: 144 × 1
#> group
#> <fct>
#> 1 H
#> 2 H
#> 3 Y
#> 4 C
#> 5 H
#> 6 C
#> 7 M
#> 8 C
#> 9 M
#> 10 M
#> # ℹ 134 more rowsas_tibble(rowData(real_experiment2))
#> # A tibble: 67 × 2
#> glycan_composition glycan_structure
#> <comp> <struct>
#> 1 Man(3)GlcNAc(3) GlcNAc(?1-?)Man(?1-?)[Man(?1-?)]Man(?1-…
#> 2 Man(3)GlcNAc(7) GlcNAc(?1-?)[GlcNAc(?1-?)]Man(?1-?)[Glc…
#> 3 Man(5)GlcNAc(2) Man(?1-?)[Man(?1-?)]Man(?1-?)[Man(?1-?)…
#> 4 Man(4)Gal(2)GlcNAc(4)Neu5Ac(2) Neu5Ac(?2-?)Gal(?1-?)GlcNAc(?1-?)[Neu5A…
#> 5 Man(3)Gal(1)GlcNAc(3) Gal(?1-?)GlcNAc(?1-?)Man(?1-?)[Man(?1-?…
#> 6 Man(3)Gal(2)GlcNAc(4)Fuc(2) Gal(?1-?)GlcNAc(?1-?)Man(?1-?)[Gal(?1-?…
#> 7 Man(3)GlcNAc(3)Fuc(1) GlcNAc(?1-?)Man(?1-?)[Man(?1-?)]Man(?1-…
#> 8 Man(3)GlcNAc(4) GlcNAc(?1-?)Man(?1-?)[GlcNAc(?1-?)Man(?…
#> 9 Man(3)Gal(2)GlcNAc(5)Neu5Ac(1) Neu5Ac(?2-?)Gal(?1-?)GlcNAc(?1-?)Man(?1…
#> 10 Man(3)Gal(1)GlcNAc(5)Fuc(1)Neu5Ac(1) Neu5Ac(?2-?)Gal(?1-?)GlcNAc(?1-?)Man(?1…
#> # ℹ 57 more rowsFor a deeper dive into Glycoverse data containers, check out Get Started with glyexp.
Data Preprocessing
Raw quantification data needs preprocessing before analysis—that’s just a fact of life in omics. Typical steps include normalization, missing value imputation, and batch effect correction. Rather than making you implement these tedious steps manually, glyclean provides a comprehensive preprocessing pipeline. Just call auto_clean() on your GlycomicSE object and you’re good to go.
clean_exp <- auto_clean(real_experiment2)
#>
#> ── Removing variables with too many missing values ──
#>
#> ℹ Applying preset "discovery"...
#> ℹ Total removed: 10 (14.93%) variables.
#> ✔ Variable removal completed.
#>
#> ── Imputing missing values ──
#>
#> ℹ Imputation method: `impute_miss_forest()`
#> ℹ Reason: default for "glycomics" with n_samples > 100.
#> ✔ Imputation completed.
#>
#> ── Normalizing data ──
#>
#> ℹ Normalization method: `normalize_total_area()`
#> ℹ Reason: default for "glycomics".
#> ✔ Normalization completed.
#>
#> ── Correcting batch effects ──
#>
#> ℹ Batch column batch not found in sample_info. Skipping batch correction.
#> ✔ Batch correction completed.Your data is now analysis-ready!
Want to customize the preprocessing steps? See Get Started with glyclean for the full toolkit.
Statistical Analysis and Visualization
Time for the fun part—statistical analysis and visualization! We’ll use glystats for the number crunching and glyvis to make sense of the results visually.
Let’s kick off with PCA to get a bird’s-eye view of our data structure.
plot_pca(clean_exp) # from `glyvis`
glyvis isn’t designed for publication-ready figures, but it’s perfect for quick exploratory visualization. Behind the scenes, plot_pca() calls gly_pca() from glystats and renders the results.
You can also break this down into separate steps:
pca_res <- gly_pca(clean_exp) # from `glystats`
autoplot(pca_res) # from `glyvis`
# you can also use `plot_pca(pca_res)`We actually recommend the two-step approach, since it gives you more flexibility with the results. You can create custom ggplot2 visualizations for publications or extract the underlying data when reviewers ask for it.
glystats covers virtually all standard omics analyses. All functions follow the same naming pattern: gly_xxx()—think gly_anova(), gly_ttest(), gly_roc(), gly_cox(), gly_wgcna(), and so on. They take a GlycomicSE, GlycoproteomicSE, or compatible SummarizedExperiment as their first argument.
The return format is consistent across all functions—a list with two components:
tidy_result: cleaned-up tibbles in tidy format. We’ve done the heavy lifting of organizing messy statistical output for you.raw_result: the original statistical objects. These are available when you need to dig deeper or perform advanced analyses.
glystats provides two helper functions to get the tidy result tibble and the raw result list from a glystats result object: get_tidy_result() and get_raw_result(). Let’s now see what the samples tibble looks like:
get_tidy_result(pca_res, "samples") # many tibbles, so we specify one of them
#> # A tibble: 8,208 × 4
#> sample group PC value
#> <chr> <fct> <dbl> <dbl>
#> 1 S1 H 1 -1.12
#> 2 S1 H 2 -1.55
#> 3 S1 H 3 0.708
#> 4 S1 H 4 1.83
#> 5 S1 H 5 -1.00
#> 6 S1 H 6 2.55
#> 7 S1 H 7 1.57
#> 8 S1 H 8 1.32
#> 9 S1 H 9 1.42
#> 10 S1 H 10 0.424
#> # ℹ 8,198 more rowsNotice the “group” column? That’s glystats being helpful— it automatically pulls relevant metadata from your GlycomicSE object and includes it in the results wherever it makes sense.
Back to that autoplot() magic we saw earlier. It automatically recognizes different glystats result types and plots accordingly— no manual specification needed. The plots won’t win any beauty contests, but they’ll get your data insights across fast.
Now let’s dive into differential expression analysis with ANOVA.
anova_res <- gly_anova(clean_exp, contrasts = "H_vs_C") # from `glystats`
#> ℹ Number of groups: 4
#> ℹ Groups: "H", "M", "Y", and "C"
#> ℹ Pairwise comparisons will be performed, with levels coming first as reference groups.
#> Warning: There were 57 warnings in `dplyr::mutate()`.
#> The first warning was:
#> ℹ In argument: `test_result = list(safe_f(formula, data = .data$data, contrasts
#> = "H_vs_C"))`.
#> ℹ In row 1.
#> Caused by warning in `model.matrix.default()`:
#> ! non-list contrasts argument ignored
#> ℹ Run `dplyr::last_dplyr_warnings()` to see the 56 remaining warnings.
get_tidy_result(anova_res, "main_test") # only one tibble here
#> # A tibble: 57 × 12
#> variable glycan_composition glycan_structure term df sumsq meansq
#> <chr> <comp> <struct> <chr> <dbl> <dbl> <dbl>
#> 1 Man(3)GlcNAc(3) Man(3)GlcNAc(3) GlcNAc(?1-?)Man… group 3 0.579 0.193
#> 2 Man(3)GlcNAc(7) Man(3)GlcNAc(7) GlcNAc(?1-?)[Gl… group 3 1.07 0.355
#> 3 Man(5)GlcNAc(2) Man(5)GlcNAc(2) Man(?1-?)[Man(?… group 3 0.887 0.296
#> 4 Man(4)Gal(2)Glc… Man(4)Gal(2)GlcNA… Neu5Ac(?2-?)Gal… group 3 0.898 0.299
#> 5 Man(3)Gal(1)Glc… Man(3)Gal(1)GlcNA… Gal(?1-?)GlcNAc… group 3 0.316 0.105
#> 6 Man(3)Gal(2)Glc… Man(3)Gal(2)GlcNA… Gal(?1-?)GlcNAc… group 3 1.02 0.341
#> 7 Man(3)GlcNAc(3)… Man(3)GlcNAc(3)Fu… GlcNAc(?1-?)Man… group 3 0.561 0.187
#> 8 Man(3)GlcNAc(4) Man(3)GlcNAc(4) GlcNAc(?1-?)Man… group 3 1.24 0.415
#> 9 Man(3)Gal(2)Glc… Man(3)Gal(2)GlcNA… Neu5Ac(?2-?)Gal… group 3 0.195 0.0649
#> 10 Man(3)Gal(1)Glc… Man(3)Gal(1)GlcNA… Neu5Ac(?2-?)Gal… group 3 3.29 1.10
#> # ℹ 47 more rows
#> # ℹ 5 more variables: statistic <dbl>, p_val <dbl>, p_adj <dbl>,
#> # effect_size <dbl>, post_hoc <chr>Excellent! Now let’s identify significantly differentially expressed glycans between HCC and healthy samples.
anova_res |>
get_tidy_result("main_test") |>
filter(p_adj < 0.05) |>
select(glycan_composition, p_adj, effect_size, post_hoc)
#> # A tibble: 18 × 4
#> glycan_composition p_adj effect_size post_hoc
#> <comp> <dbl> <dbl> <chr>
#> 1 Man(3)Gal(1)GlcNAc(5)Fuc(1)Neu5Ac(1) 0.00625 0.104 M_vs_C
#> 2 Man(3)GlcNAc(4)Fuc(1) 0.000000227 0.247 H_vs_Y;H_vs_C;…
#> 3 Man(3)GlcNAc(5) 0.00179 0.126 H_vs_C;M_vs_C
#> 4 Man(3)Gal(1)GlcNAc(5)Neu5Ac(1) 0.0283 0.0809 M_vs_C
#> 5 Man(4)Gal(1)GlcNAc(3)Neu5Ac(1) 0.0380 0.0750 H_vs_Y;M_vs_Y
#> 6 Man(3)GlcNAc(5)Fuc(1) 0.000607 0.145 H_vs_Y;H_vs_C;…
#> 7 Man(3)Gal(1)GlcNAc(5)Fuc(1) 0.00419 0.112 H_vs_Y;H_vs_C
#> 8 Man(3)Gal(1)GlcNAc(5) 0.00578 0.106 H_vs_C;M_vs_C
#> 9 Man(3)Gal(1)GlcNAc(4)Fuc(1) 0.0380 0.0754 H_vs_Y;H_vs_C
#> 10 Man(3)Gal(2)GlcNAc(4)Neu5Ac(1) 0.00542 0.108 H_vs_C;M_vs_C
#> 11 Man(3)Gal(2)GlcNAc(4)Fuc(1) 0.00179 0.127 M_vs_C;Y_vs_C
#> 12 Man(3)Gal(3)GlcNAc(5)Neu5Ac(2) 0.00338 0.116 H_vs_M;H_vs_Y;…
#> 13 Man(3)Gal(3)GlcNAc(5)Fuc(1)Neu5Ac(2) 0.0000000575 0.266 H_vs_C;M_vs_C;…
#> 14 Man(3)Gal(3)GlcNAc(5)Neu5Ac(3) 0.0283 0.0815 H_vs_M;H_vs_C
#> 15 Man(3)Gal(4)GlcNAc(6)Neu5Ac(2) 0.00179 0.129 H_vs_C;M_vs_C;…
#> 16 Man(3)Gal(3)GlcNAc(5)Fuc(1)Neu5Ac(3) 0.0000000169 0.286 H_vs_C;M_vs_C;…
#> 17 Man(3)Gal(3)GlcNAc(5)Fuc(2)Neu5Ac(3) 0.00000675 0.205 H_vs_C;M_vs_C;…
#> 18 Man(3)Gal(4)GlcNAc(6)Fuc(1)Neu5Ac(2) 0.000607 0.146 M_vs_CFor the full statistical arsenal, check out Get Started with glystats and Get Started with glyvis.
Advanced Motif Analysis
Up to now, we’ve covered standard glycomics workflows. While glycoverse certainly streamlines these analyses, it truly shines when it comes to advanced glycan structure analysis.
Before diving into motifs, let’s get acquainted with glyrepr::glycan_structure() vectors.
clean_exp |>
rowData() |>
as_tibble() |>
pull(glycan_structure)
#> <glycan_structure[57]>
#> [1] GlcNAc(?1-?)Man(?1-?)[Man(?1-?)]Man(?1-?)GlcNAc(?1-?)GlcNAc(?1-
#> [2] GlcNAc(?1-?)[GlcNAc(?1-?)]Man(?1-?)[GlcNAc(?1-?)[GlcNAc(?1-?)]Man(?1-?)][GlcNAc(?1-?)]Man(?1-?)GlcNAc(?1-?)GlcNAc(?1-
#> [3] Man(?1-?)[Man(?1-?)]Man(?1-?)[Man(?1-?)]Man(?1-?)GlcNAc(?1-?)GlcNAc(?1-
#> [4] Neu5Ac(?2-?)Gal(?1-?)GlcNAc(?1-?)[Neu5Ac(?2-?)Gal(?1-?)GlcNAc(?1-?)]Man(?1-?)[Man(?1-?)Man(?1-?)]Man(?1-?)GlcNAc(?1-?)GlcNAc(?1-
#> [5] Gal(?1-?)GlcNAc(?1-?)Man(?1-?)[Man(?1-?)]Man(?1-?)GlcNAc(?1-?)GlcNAc(?1-
#> [6] Gal(?1-?)GlcNAc(?1-?)Man(?1-?)[Gal(?1-?)[Fuc(?1-?)]GlcNAc(?1-?)Man(?1-?)]Man(?1-?)GlcNAc(?1-?)[Fuc(?1-?)]GlcNAc(?1-
#> [7] GlcNAc(?1-?)Man(?1-?)[Man(?1-?)]Man(?1-?)GlcNAc(?1-?)[Fuc(?1-?)]GlcNAc(?1-
#> [8] GlcNAc(?1-?)Man(?1-?)[GlcNAc(?1-?)Man(?1-?)]Man(?1-?)GlcNAc(?1-?)GlcNAc(?1-
#> [9] Neu5Ac(?2-?)Gal(?1-?)GlcNAc(?1-?)Man(?1-?)[Gal(?1-?)GlcNAc(?1-?)Man(?1-?)][GlcNAc(?1-?)]Man(?1-?)GlcNAc(?1-?)GlcNAc(?1-
#> [10] Neu5Ac(?2-?)Gal(?1-?)GlcNAc(?1-?)Man(?1-?)[GlcNAc(?1-?)Man(?1-?)][GlcNAc(?1-?)]Man(?1-?)GlcNAc(?1-?)[Fuc(?1-?)]GlcNAc(?1-
#> ... (47 more not shown)
#> # Unique structures: 57Just like integer() and character(), glycan_structure() is a specialized vector type. Some software (like pGlyco3 and StrucGP) outputs structural information as text strings. When you import this data using glyread, the glyparse package automatically converts these strings into proper glycan_structure() vectors and stores them in the variable information tibble. Note that not all software provides structural data—some only give compositions.
For glycomics data, this information is hard to come by automatically. You can do it manually by parsing the glycan structure strings using glyparse and using left_join_var() to join the parsed structures to the variable information tibble.
Fortunately, our example dataset includes structural information, opening up a world of advanced analytical possibilities. Let’s explore motif analysis.
Quick note: The printed structures use IUPAC-condensed notation, which we’ll also use for defining motifs below. Don’t worry if it looks intimidating—we’ll include visual diagrams to help. That said, if you’re planning to do serious structural analysis, learning IUPAC-condensed notation is worth the investment. Check out this guide to get started—it’s easier than it looks!
Human serum N-glycans can have three types of branch terminals (ignoring a1-3 Fuc):
- A Sialyl-LacNAc motif
- A LacNAc motif without sialic acids
- Only a GlcNAc without further elongation
Here’s how we express these motifs in IUPAC-condensed notation:
motifs <- c(
motif1 = "Neu5Ac(??-?)Gal(??-?)GlcNAc(??-",
motif2 = "Gal(??-?)GlcNAc(??-",
motif3 = "GlcNAc(??-"
)The “??-?” represents unknown linkages—a common limitation in mass spectrometry data.
Here’s our research question: Which branching motif show differential expression across conditions? Without glycoverse, this would be a nightmare to tackle manually. Take a moment to imagine the pain of doing this by hand!
Now, the glycoverse solution:
motif_anova_res <- clean_exp |>
quantify_motifs(motifs, alignments = "terminal") |> # quantify these motifs
gly_anova() # and perform ANOVA
#> ℹ Number of groups: 4
#> ℹ Groups: "H", "M", "Y", and "C"
#> ℹ Pairwise comparisons will be performed, with levels coming first as reference groups.
get_tidy_result(motif_anova_res, "main_test")
#> # A tibble: 3 × 12
#> variable trait motif_structure term df sumsq meansq statistic p_val
#> <chr> <chr> <struct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 motif1 motif1 Neu5Ac(??-?)Gal(?… group 3 0.0548 0.0183 1.14 3.35e-1
#> 2 motif2 motif2 Gal(??-?)GlcNAc(?… group 3 0.419 0.140 2.22 8.86e-2
#> 3 motif3 motif3 GlcNAc(??- group 3 3.11 1.04 5.86 8.48e-4
#> # ℹ 3 more variables: p_adj <dbl>, effect_size <dbl>, post_hoc <chr>quantify_motifs() transforms your data into a new SummarizedExperiment. Instead of quantification of glycans, you now have motif abundances across samples. All glystats functions work seamlessly with this container—including gly_anova().
Now we can answer our question using standard tidyverse operations, since motif_anova_res$tidy_result$main_test is just a regular tibble:
motif_anova_res |>
get_tidy_result("main_test") |>
filter(p_adj < 0.05)
#> # A tibble: 1 × 12
#> variable trait motif_structure term df sumsq meansq statistic p_val
#> <chr> <chr> <struct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 motif3 motif3 GlcNAc(??- group 3 3.11 1.04 5.86 0.000848
#> # ℹ 3 more variables: p_adj <dbl>, effect_size <dbl>, post_hoc <chr>glymotif has much more to offer beyond these examples. Dive deeper with Get Started with glymotif.
Derived Trait Analysis
Let’s wrap up with derived traits—a clever analytical approach developed by the N-glycomics community for glycome characterization. Classic examples include:
- High-mannose glycan proportion
- Core-fucosylation rate within complex glycans
- Average sialylation per galactose residue
glydet calculates derived traits in a flash. Using it couldn’t be simpler:
trait_exp <- derive_traits(clean_exp) # from `glydet`
trait_exp
#> class: SummarizedExperiment
#> dim: 14 144
#> metadata(2): exp_type glycan_type
#> assays(1): abundance
#> rownames(14): TM TH ... AG TS
#> rowData names(2): trait explanation
#> colnames(144): S1 S2 ... S143 S144
#> colData names(1): groupThat’s it! Just like quantify_motifs(), derive_traits() creates a new SummarizedExperiment, but now with trait values per sample.
The variable information shows what we’re working with:
as_tibble(rowData(trait_exp))
#> # 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 “trait” column lists all the derived traits we can analyze.
glydet comes with a comprehensive set of built-in traits:
TM: Proportion of high-mannose glycansTH: Proportion of hybrid glycansTC: Proportion of complex glycansMM: Average number of mannoses within high-mannose glycansCA2: Proportion of bi-antennary glycans within complex glycansCA3: Proportion of tri-antennary glycans within complex glycansCA4: Proportion of tetra-antennary glycans within complex glycansTF: Proportion of fucosylated glycansTFc: Proportion of core-fucosylated glycansTFa: Proportion of arm-fucosylated glycansTB: Proportion of glycans with bisecting GlcNAcSG: Average degree of sialylation per galactoseGA: Average degree of galactosylation per antennaTS: Proportion of sialylated glycans
These represent the most widely used traits in glycomics literature.
Let’s identify traits with significantly different values across conditions:
trait_exp |>
gly_anova() |>
get_tidy_result("main_test") |>
filter(p_adj < 0.05)
#> ℹ Number of groups: 4
#> ℹ Groups: "H", "M", "Y", and "C"
#> ℹ Pairwise comparisons will be performed, with levels coming first as reference groups.
#> # A tibble: 8 × 12
#> variable trait explanation term df sumsq meansq statistic p_val
#> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 CA2 CA2 Proportion of bi-… group 3 0.0406 0.0135 5.73 1.00e-3
#> 2 CA3 CA3 Proportion of tri… group 3 1.76 0.585 5.40 1.52e-3
#> 3 CA4 CA4 Proportion of tet… group 3 1.18 0.393 4.87 2.99e-3
#> 4 TF TF Proportion of fuc… group 3 2.15 0.716 8.44 3.41e-5
#> 5 TFc TFc Proportion of cor… group 3 2.15 0.716 8.44 3.41e-5
#> 6 TFa TFa Proportion of arm… group 3 2.78 0.928 7.37 1.28e-4
#> 7 TB TB Proportion of gly… group 3 1.63 0.543 3.67 1.40e-2
#> 8 AG AG Abundance-weighte… group 3 0.0221 0.00738 3.59 1.53e-2
#> # ℹ 3 more variables: p_adj <dbl>, effect_size <dbl>, post_hoc <chr>Once again, it’s just that straightforward.
This just scratches the surface of glydet’s capabilities. The real power lies in defining custom traits tailored to your research questions. Explore the possibilities in Get Started with glydet.
What’s Next?
This vignette has given you a taste of glycoverse in action through a real-world glycomics workflow. But we’ve barely scratched the surface! Now that you’ve got the basics down, you’re ready to unlock the full potential of each package.
Here’s your roadmap to mastering each component:
- glyexp — Master experiment objects and data manipulation
- glyread — Import and organize glycomics data
- glyclean — Build custom preprocessing pipelines
- glystats — Explore the full statistical toolkit
- glyvis — Create stunning visualizations
- glymotif — Define and analyze custom motifs
- glydet — Create powerful derived traits
- glyenzy — Explore enzyme-substrate relationships (we didn’t cover this one, but it’s fascinating!)
- glyrepr — Master glycan structure representation
- glyparse — Parse and convert structural formats
- glydraw — Draw glycan structures
- glydb — Access glycan databases
- glyanno — Annotate glycan structures
- glyfun — Perform functional enrichment analysis
- glysmith — Master the full analytical pipeline
Happy glycan hunting! 🧬