Get Started with glyclean

Step 1: Normalization – Getting Everyone on the Same Page 📏

Imagine you’re comparing heights of people measured in different units – some in feet, some in meters, some in furlongs (okay, maybe not furlongs). Normalization does the same thing for your omics data, bringing all intensities to a comparable scale.

Available normalization methods in glyclean:

Pro Tips:

• Notice the by parameter in many functions? This allows stratified normalization within groups – super useful when you have distinct experimental conditions!
• Maximum flexibility: The by parameter accepts both column names from your sample metadata and direct vectors! This gives you complete control over grouping, even for custom groupings that aren’t stored in your experiment object.

# Example: Using by parameter with a custom vector
# Normalize within custom groups (e.g., based on some external criteria)
custom_groups <- c("A", "A", "B", "B", "C", "C", "A", "A", "B", "B", "C", "C")
normalized_exp <- normalize_median(exp, by = custom_groups)

Compositional Data: CLR, ALR, and Median Quotient Normalization

For glycomics data (where intensities represent relative abundances of glycans), standard normalization methods may not be appropriate. Instead, compositional data analysis techniques are recommended:

• normalize_median_quotient(): Median quotient normalization. A robust method specifically designed for compositional data, dividing each sample by its median quotient relative to a reference sample.

• transform_clr(): Centered Log-Ratio transformation with glycowork-compatible internal logic. The transformation is carried out on the log2 scale and then back-transformed to the original ratio space, so zeros remain zeros in the returned result. The stochastic branch uses per-feature noise, matching glycowork.

• transform_alr(): Additive Log-Ratio transformation with glycowork-compatible internal logic. The reference glycan is chosen based on Procrustes correlation to CLR geometry and minimal between-group variance, and the final result is returned in the original ratio space. If no suitable reference is found, ALR falls back to CLR.

The CLR and ALR methods align with the best practices described in DOI: 10.1038/s41467-025-56249-3 for glycomics data preprocessing.

Here we do the median normalization manually:

normed_exp <- normalize_median(exp)

Step 2: Variable Filtering – Saying Goodbye to the Unreliable 🧹

Some variables are like that friend who never shows up to plans – they’re missing most of the time and aren’t very helpful. In omics data, variables with too many missing values are often more noise than signal.

Available variable filtering functions in glyclean:

Here we remove variables with more than 50% missing values in all samples:

filtered_exp <- remove_rare(normed_exp, prop = 0.5)
#> ℹ Removed 137 of 4262 (3.21%) variables.

Step 3: Imputation – Filling in the Blanks Intelligently 🔮

Here’s where things get scientifically interesting! Missing values in mass spectrometry aren’t randomly distributed – they follow patterns based on the physics and chemistry of the measurement process.

The Science: Missing values in MS-based omics data are typically “Missing Not At Random” (MNAR), meaning they’re related to the intensity of the signal. Low-abundance ions are more likely to be missed, either due to poor ionization efficiency or because they fall below the detection threshold.

Your imputation toolkit:

For demonstration, let’s use the simple zero imputation:

imputed_exp <- impute_zero(filtered_exp)

Step 4: Aggregation – From Peptides to Glycoforms 🔄

Here’s where glycoproteomics gets uniquely challenging! Search engines typically report results at the PSM or peptide level, but what we really care about are glycoforms – the specific glycan structures attached to specific protein sites.

The Problem: One glycoform can appear multiple times in your data due to:

• Different charge states
• Post-translational modifications
• Missed protease cleavages
• Different peptide sequences covering the same site

The Solution: Intelligent aggregation!

Aggregation levels available:

• “gps”: Glycopeptides with structures (most detailed)
• “gp”: Glycopeptides with compositions
• “gfs”: Glycoforms with structures (recommended default)
• “gf”: Glycoforms with compositions (most condensed)

aggregated_exp <- aggregate(imputed_exp, to_level = "gf")

Pro move: Re-normalize after aggregation to account for the new intensity distributions:

aggregated_exp2 <- normalize_median(aggregated_exp)

Step 5: Batch Effect Correction – Harmonizing Your Orchestra 🎼

Batch effects are the “different violinists playing the same piece” problem of omics data. Even when following identical protocols, subtle differences in instruments, reagents, or environmental conditions can introduce systematic bias.

The Good News: If your experimental design is well-controlled (conditions distributed across batches), batch effect correction can work wonders!

# To detect batch effects:
p_values <- detect_batch_effect(aggregated_exp2)
#> ℹ Detecting batch effects using ANOVA for 2738 variables...
#> ✔ Batch effect detection completed. 19 out of 2738 variables show significant batch effects (p < 0.05).
p_values[1:5]
#> P08185-176-Hex(5)HexNAc(4)NeuAc(2) P04196-344-Hex(5)HexNAc(4)NeuAc(1) 
#>                          0.6998998                          0.3951691 
#>         P04196-344-Hex(5)HexNAc(4)         P10909-291-Hex(6)HexNAc(5) 
#>                          0.4852698                          0.3329619 
#> P04196-344-Hex(5)HexNAc(4)NeuAc(2) 
#>                          0.4907154

Here we do not have batch effects, but we will correct it anyway for demonstration.

# To correct batch effects:
corrected_exp <- correct_batch_effect(aggregated_exp2)
#> Found3batches
#> Adjusting for0covariate(s) or covariate level(s)
#> Standardizing Data across genes
#> Fitting L/S model and finding priors
#> Finding parametric adjustments
#> Adjusting the Data

Step 6: Protein Expression Adjustment – Separating Signal from Noise 🎯

Here’s a fascinating biological puzzle: when you measure a glycopeptide’s intensity, you’re actually seeing two stories at once!

The Dual Nature Problem: Every glycopeptide intensity is a tale of two factors:

• The protein story: How much of the substrate protein is present
• The glycosylation story: How actively that protein is being glycosylated

It’s like trying to understand applause volume – is the audience bigger, or are they just more enthusiastic?

The Solution: If you’re specifically interested in glycosylation changes (not protein abundance changes), you need to mathematically “subtract out” the protein expression component. Think of it as noise cancellation for biology!

# We don't run the code here because we don't have the protein expression matrix.
adjusted_exp <- adjust_protein(inferred_exp, pro_expr_mat)

Important Prerequisites:

• Your protein expression matrix must have the same samples as your glycopeptide data
• The protein expression matrix should be properly preprocessed (good news: you can use glyclean functions for this too!)
• Sample matching needs to be perfect – no mixing apples with oranges!

Pro Tip: Most glyclean functions happily accept plain matrices as input, so you can preprocess your protein expression data using the same toolkit. Consistency is key!

🎉 Mission Accomplished!

Congratulations! You now have beautifully preprocessed, analysis-ready glycoproteomics data. Your data has been normalized, filtered, imputed, aggregated, and batch-corrected – it’s ready to reveal its biological secrets!

Working with Plain Matrices

Most glyclean functions also support plain matrices as input! This means you can use the package’s powerful preprocessing capabilities even if you’re working with traditional data formats. The functions intelligently detect your input type and handle it appropriately.

# Example: Using glyclean with a plain matrix
my_matrix <- matrix(rnorm(100), nrow = 10)
normalized_matrix <- normalize_median(my_matrix)
imputed_matrix <- impute_sample_min(normalized_matrix)

Flexible Grouping with Custom Vectors

Many functions in glyclean accept a by parameter for stratified processing. This parameter offers maximum flexibility – it accepts both column names from your sample metadata and direct vectors!

# Using column names (standard approach)
normalized_exp <- normalize_median(exp, by = "group")

# Using custom vectors (advanced approach)
custom_groups <- c("A", "A", "B", "B", "C", "C", "A", "A", "B", "B", "C", "C")
normalized_exp <- normalize_median(exp, by = custom_groups)

# This also works with matrices
normalized_matrix <- normalize_median(my_matrix, by = custom_groups)

This gives you complete control over grouping, even for custom groupings that aren’t stored in your experiment object. Perfect for complex experimental designs or external grouping criteria!