Package 'glystats'

Description

Filtering the experiment to keep only significant variables is a common task. This function provides a convenient way to do this. It supports results from all glystats DEA functions including gly_anova(), gly_ancova(), gly_kruskal(), gly_ttest(), gly_wilcox(), gly_limma(), and gly_linear_model().

Usage

filter_sig_vars(
  exp,
  res,
  p_adj_cutoff = 0.05,
  p_val_cutoff = NULL,
  fc_cutoff = NULL,
  ...
)

## S3 method for class 'glystats_anova_res'
filter_sig_vars(
  exp,
  res,
  p_adj_cutoff = 0.05,
  p_val_cutoff = NULL,
  fc_cutoff = NULL,
  on = "main_test",
  comparison = NULL,
  ...
)

## S3 method for class 'glystats_ancova_res'
filter_sig_vars(
  exp,
  res,
  p_adj_cutoff = 0.05,
  p_val_cutoff = NULL,
  fc_cutoff = NULL,
  on = "main_test",
  comparison = NULL,
  ...
)

## S3 method for class 'glystats_kruskal_res'
filter_sig_vars(
  exp,
  res,
  p_adj_cutoff = 0.05,
  p_val_cutoff = NULL,
  fc_cutoff = NULL,
  on = "main_test",
  comparison = NULL,
  ...
)

## S3 method for class 'glystats_ttest_res'
filter_sig_vars(
  exp,
  res,
  p_adj_cutoff = 0.05,
  p_val_cutoff = NULL,
  fc_cutoff = NULL,
  ...
)

## S3 method for class 'glystats_wilcox_res'
filter_sig_vars(
  exp,
  res,
  p_adj_cutoff = 0.05,
  p_val_cutoff = NULL,
  fc_cutoff = NULL,
  ...
)

## S3 method for class 'glystats_limma_res'
filter_sig_vars(
  exp,
  res,
  p_adj_cutoff = 0.05,
  p_val_cutoff = NULL,
  fc_cutoff = NULL,
  comparison = NULL,
  ...
)

## S3 method for class 'glystats_linear_model_res'
filter_sig_vars(
  exp,
  res,
  p_adj_cutoff = 0.05,
  p_val_cutoff = NULL,
  fc_cutoff = NULL,
  term = NULL,
  ...
)

Arguments

Value

A filtered object of the same class as exp.

Examples

library(glyexp)
library(glyclean)

exp <- auto_clean(real_experiment) |>
  glyexp::slice_head_row(n = 10)
res <- gly_anova(exp)
sig_exp <- filter_sig_vars(exp, res)

Description

Syntax sugar for ⁠$tidy_result⁠ and ⁠$raw_result⁠ elements of a glystats result object. It's useful to be used in pipes.

Usage

get_tidy_result(res, which = NULL)

get_raw_result(res, which = NULL)

Arguments

Value

A tibble.

Examples

library(glyexp)
library(glyclean)
library(dplyr)

exp <- auto_clean(real_experiment) |>
  glyexp::slice_head_row(n = 10)

# Using a pipe
sig_res <- exp |>
  gly_anova() |>
  get_tidy_result("main_test") |>
  filter(p_adj < 0.05)

# Equivalent to
anova_res <- gly_anova(exp)
sig_res <- anova_res$tidy_result$main_test |>
  filter(p_adj < 0.05)

Description

Perform ANCOVA for glycomics or glycoproteomics data with sample-level covariates. The function supports parametric comparison of multiple groups while adjusting for covariates. For significant results, emmeans post-hoc comparisons (Tukey adjustment) are automatically performed. P-values are adjusted for multiple testing using the method specified by p_adj_method.

Usage

gly_ancova(
  exp,
  group_col = "group",
  covariate_cols = NULL,
  p_adj_method = "BH",
  add_info = TRUE,
  ...,
  log_transform = TRUE
)

Arguments

Details

Tests and effect sizes use the selected analysis scale. Reported fold changes continue to use ratios of group means from the original supplied values.

By default, expression data are transformed using log2(x + 1e-6) before statistical testing. Set log_transform = FALSE to skip this transformation. At least 2 groups and at least 1 covariate are required.

For any variable failed to fit a stats::aov() model, NAs will be assigned to the results in both main test and post-hoc test.

Post-hoc Test: emmeans pairwise comparisons with Tukey adjustment (emmeans::emmeans()) are performed for variables with significant main effects (p_adj < 0.05).

Value

A list containing two elements:

• tidy_result: A list containing:

  • main_test: A tibble with ANCOVA results containing the following columns:

    • variable: Variable name

    • term: ANCOVA term (usually "group")

    • df: Degrees of freedom

    • sumsq: Sum of squares

    • meansq: Mean squares

    • statistic: F-statistic

    • p_val: Raw p-value from ANCOVA

    • p_adj: Adjusted p-value (if p_adj_method is not NULL)

    • post_hoc: Significant group pairs from post-hoc test, in the format of "ref_vs_test".

  • post_hoc_test: A tibble with pairwise comparison results containing the following columns:

    • variable: Variable name

    • ref_group: Reference group

    • test_group: Test/treatment/case group

    • p_val: Adjusted p-value from emmeans pairwise test

    • p_adj: Adjusted p-value from emmeans pairwise test

• raw_result: A list containing:

  • main_test: A list of raw aov model objects.

  • post_hoc_test: A list of raw emmeans results.

Required packages

This function requires the emmeans package for post-hoc comparisons.

See Also

stats::aov(), emmeans::emmeans(), emmeans::contrast()

Description

Perform one-way or repeated-measures ANOVA for glycomics or glycoproteomics data. For significant results, Tukey's HSD post-hoc test is automatically performed. P-values are adjusted for multiple testing using the method specified by p_adj_method.

Usage

gly_anova(
  exp,
  group_col = "group",
  p_adj_method = "BH",
  add_info = TRUE,
  ...,
  subject_col = NULL,
  log_transform = TRUE
)

Arguments

Details

Tests and effect sizes use the selected analysis scale. Reported fold changes continue to use ratios of group means from the original supplied values.

By default, expression data are transformed using log2(x + 1e-6) before statistical testing. Set log_transform = FALSE to skip this transformation. At least 2 groups are required in the grouping variable. In paired analyses, only subjects observed in every group are used, with at most one sample per subject and group. The reported effect size is partial eta-squared for the group term.

For any variable failed to fit a stats::aov() model, NAs will be assigned to the results in both main test and post-hoc test.

Post-hoc Test: Tukey's HSD test for pairwise comparisons (stats::TukeyHSD()) is performed for variables with significant main effects (p_adj < 0.05). Paired analyses instead use paired t-tests with Holm correction.

Value

A list containing three elements:

• tidy_result: A list containing:

  • main_test: A tibble with ANOVA results containing the following columns:

    • variable: Variable name

    • term: ANOVA term (usually "groups")

    • df: Degrees of freedom

    • sumsq: Sum of squares

    • meansq: Mean squares

    • statistic: F-statistic

    • p_val: Raw p-value from ANOVA

    • effect_size: Eta-squared, or partial eta-squared for paired analyses

    • p_adj: Adjusted p-value (if p_adj_method is not NULL)

    • post_hoc: Significant group pairs from post-hoc test, in the format of "ref_vs_test".

  • post_hoc_test: A tibble with pairwise comparison results containing the following columns:

    • variable: Variable name

    • ref_group: Reference group

    • test_group: Test/treatment/case group

    • p_val: Raw post-hoc p-value

    • p_adj: Adjusted post-hoc p-value

• raw_result: A list containing:

  • main_test: A list of raw aov model objects.

  • post_hoc_test: A list of raw TukeyHSD objects or paired t-test summaries.

• meta_data: A list containing metadata from the input experiment.

See Also

stats::aov(), stats::TukeyHSD()

Description

Perform pairwise correlation analysis on variables or samples in the expression data. The function calculates correlation coefficients and p-values for all pairs, with optional multiple testing correction.

Usage

gly_cor(
  exp,
  on = "variable",
  method = "pearson",
  p_adj_method = "BH",
  ...,
  log_transform = TRUE
)

Arguments

Details

Set log_transform = FALSE to use the expression values as supplied, for example when they are already log-transformed.

By default, the function performs log2 transformation on the expression data (log2(x + 1e-6)) before correlation analysis. When on = "variable" (default), correlations are calculated between variables across samples. When on = "sample", correlations are calculated between samples across variables.

Correlation Calculation: Correlation coefficients and p-values are calculated using Hmisc::rcorr() which is more efficient than pairwise stats::cor.test() calls.

Multiple Testing Correction: P-values are adjusted for multiple testing using the method specified by p_adj_method.

Value

A list with three elements:

• tidy_result: A tibble with correlation results containing the following columns:

  • variable1 or sample1: First element of the pair (depending on on parameter)

  • variable2 or sample2: Second element of the pair (depending on on parameter)

  • cor: Correlation coefficient

  • p_val: Raw p-value from correlation test

  • p_adj: Adjusted p-value (if p_adj_method is not NULL)

• raw_result: The raw rcorr object from Hmisc::rcorr()

• meta_data: A list containing metadata from the input experiment The list has classes glystats_cor_res and glystats_res.

Required packages

This function requires the Hmisc package for efficient correlation calculation.

See Also

Hmisc::rcorr(), stats::cor()

Description

Fit a Cox proportional hazards model for each variable in the expression data, and extract p-values and hazard ratios from it.

Usage

gly_cox(
  exp,
  time_col = "time",
  event_col = "event",
  p_adj_method = "BH",
  add_info = TRUE,
  ...,
  log_transform = FALSE
)

Arguments

Details

Expression values are used as supplied by default. Set log_transform = TRUE to apply log2(x + 1e-6) before analysis.

The function fits an Cox proportional hazards model for each variable by:

coxph(Surv(time, event) ~ expression_value)

P-values are adjusted by Benjamini-Hochberg method by default.

Value

A list with three elements:

• tidy_result: A tibble with Cox model results containing the following columns:

  • variable: Variable name

  • coefficient: Regression coefficient (log hazard ratio)

  • std.error: Standard error of the coefficient

  • statistic: Wald test statistic

  • p_val: Raw p-value from Wald test

  • hr: Hazard ratio (exp(coefficient))

  • p_adj: Adjusted p-value (if p_adj_method is not NULL)

• raw_result: A list of raw coxph model objects.

• meta_data: A list containing metadata from the input experiment

See Also

survival::coxph(), survival::Surv()

Description

Calculate fold change for a given expression matrix and group information. It could only be used for 2-group analysis. When you run this function, you will see message about "Ref Group" and "Test Group". "Ref Group" is the reference group, and "Test Group" is the test/treatment/case group.

Usage

gly_fold_change(
  exp,
  group_col = "group",
  add_info = TRUE,
  log_transform = FALSE
)

Arguments

Value

A list with three elements:

• tidy_result: A tibble with fold change results containing the following columns:

  • variable: Variable name

  • log2fc: Log2 ratio of arithmetic group means by default, or difference of mean log2-transformed values when log_transform = TRUE. If more than two groups, two additional columns ref_group and test_group will be added.

• raw_result: The raw result (same as tidy_result for this function)

• meta_data: A list containing metadata from the input experiment

Description

Perform hierarchical clustering on the expression data. The function uses stats::hclust() to perform clustering and provides tidy results including cluster assignments, dendrogram data for plotting, and merge heights.

Usage

gly_hclust(
  exp,
  on = "variable",
  k_values = c(2, 3, 4, 5),
  scale = TRUE,
  add_info = TRUE,
  ...,
  log_transform = TRUE
)

Arguments

Details

Set log_transform = FALSE to use the expression values as supplied, for example when they are already log-transformed.

By default, the function performs log2 transformation on the expression data (log2(x + 1e-6)) before clustering. When on = "variable" (default), variables are clustered based on their expression patterns across samples. When on = "sample", samples are clustered based on their expression profiles across variables.

Distance Calculation: Distance is calculated using stats::dist() with the specified method.

Clustering Method: Hierarchical clustering is performed using stats::hclust() with the specified method.

Cluster Assignment: The dendrogram is cut at different heights to produce cluster assignments for the specified k values using stats::cutree().

Value

A list containing:

• tidy_result: A list of tibbles with clustering results:

  • clusters: Cluster assignments containing the following columns:

    • variable or sample: Variable or sample name (depending on on parameter)

    • cluster_k2, cluster_k3, etc.: Cluster assignments for different k values

  • dendrogram: Dendrogram segment data for plotting (if ggdendro is available) containing:

    • x, y, xend, yend: Segment coordinates for plotting

  • heights: Merge heights and steps containing the following columns:

    • merge_step: Step number in the clustering process

    • height: Height at which clusters are merged

    • n_clusters: Number of clusters remaining after this merge

  • labels: Labels and their positions (if ggdendro is available) containing:

    • x, y: Position coordinates

    • label: Label text

• raw_result: The raw hclust object from stats::hclust()

• meta_data: A list containing metadata from the input experiment

Required packages

This function only uses base R packages and does not require additional dependencies. For enhanced dendrogram plotting capabilities, the ggdendro package is recommended but not required.

See Also

stats::hclust(), stats::dist(), stats::cutree()

Description

Perform k-means clustering on the expression data. The function uses stats::kmeans() to perform clustering and provides tidy results with cluster assignments.

Usage

gly_kmeans(
  exp,
  on = "variable",
  centers = 3,
  scale = TRUE,
  add_info = TRUE,
  ...,
  log_transform = TRUE
)

Arguments

Details

Set log_transform = FALSE to use the expression values as supplied, for example when they are already log-transformed.

By default, the function performs log2 transformation on the expression data (log2(x + 1e-6)) before clustering. When on = "variable" (default), variables are clustered based on their expression patterns across samples. When on = "sample", samples are clustered based on their expression profiles across variables.

Data Preparation: Data is log2-transformed when log_transform = TRUE and optionally scaled before clustering.

Clustering Method: K-means clustering is performed using stats::kmeans() with the specified parameters.

Value

A list with three elements:

• tidy_result: A tibble with cluster assignments containing the following columns:

  • variable or sample: Variable or sample name (depending on on parameter)

  • cluster: Cluster assignment

• raw_result: The raw kmeans object from stats::kmeans().

• meta_data: A list containing metadata from the input experiment

Required packages

This function only uses base R packages and does not require additional dependencies.

See Also

stats::kmeans()

Description

Perform Kruskal-Wallis or Friedman tests for glycomics or glycoproteomics data. For significant results, Dunn's post-hoc test is automatically performed. P-values are adjusted for multiple testing using the method specified by p_adj_method.

Usage

gly_kruskal(
  exp,
  group_col = "group",
  p_adj_method = "BH",
  add_info = TRUE,
  ...,
  subject_col = NULL,
  log_transform = TRUE
)

Arguments

Details

Tests and effect sizes use the selected analysis scale. Reported fold changes continue to use ratios of group means from the original supplied values.

By default, expression data are transformed using log2(x + 1e-6) before statistical testing. Set log_transform = FALSE to skip this transformation. At least 2 groups are required in the grouping variable. In paired analyses, only subjects observed in every group are used, with at most one sample per subject and group. The reported effect size is Kendall's $W$.

For any variable failed to fit a stats::kruskal.test() or stats::friedman.test() model, NAs will be assigned to the results in both main test and post-hoc test.

Post-hoc Test: Dunn's test with Holm correction for multiple comparisons (FSA::dunnTest()) is performed for variables with significant main effects (p_adj < 0.05). Paired analyses instead use paired Wilcoxon signed-rank tests with Holm correction.

Value

A list containing three elements:

• tidy_result: A list containing:

  • main_test: A tibble with Kruskal-Wallis or Friedman test results containing the following columns:

    • variable: Variable name

    • statistic: Kruskal-Wallis test statistic

    • p_val: Raw p-value from Kruskal-Wallis test

    • parameter: Degrees of freedom

    • method: Statistical method used

    • effect_size: Epsilon-squared, or Kendall's W for paired analyses

    • p_adj: Adjusted p-value (if p_adj_method is not NULL)

    • post_hoc: Significant group pairs from post-hoc test, in the format of "ref_vs_test".

  • post_hoc_test: A tibble with pairwise comparison results containing the following columns:

    • variable: Variable name

    • ref_group: Reference group

    • test_group: Test/treatment/case group

    • p_val: Raw post-hoc p-value

    • p_adj: Holm-adjusted post-hoc p-value

• raw_result: A list containing:

  • main_test: A list of raw kruskal.test or friedman.test objects.

  • post_hoc_test: A list of raw Dunn or paired Wilcoxon summaries.

• meta_data: A list containing metadata from the input experiment.

Required packages

Unpaired analyses require the FSA package for Dunn's post-hoc test.

See Also

stats::kruskal.test(), stats::friedman.test(), FSA::dunnTest()

Description

Perform differential expression analysis using linear models with empirical Bayes moderation from the limma package. Supports both two-group and multi-group comparisons.

Usage

gly_limma(
  exp,
  group_col = "group",
  covariate_cols = NULL,
  subject_col = NULL,
  p_adj_method = "BH",
  ref_group = NULL,
  contrasts = NULL,
  add_info = TRUE,
  ...,
  log_transform = TRUE
)

Arguments

Details

By default, expression data are transformed using log2(x + 1e-6) before statistical testing. Set log_transform = FALSE to skip this transformation. The analysis uses linear models with empirical Bayes moderation to improve statistical power, especially for small sample sizes.

The log2fc column always contains the fitted model contrast. When log_transform = FALSE, its units are those of the supplied expression data; it is a log2 fold change only when those values are already on a log2 scale.

For two or more groups, all pairwise comparisons are automatically generated (one contrast for two groups) and performed using contrast matrices unless custom contrasts are specified.

If covariates are provided, they are added as additional terms in the design matrix and are not part of the group contrasts.

If subjects are provided, they are included as a blocking factor in the design matrix using the formula ~ 0 + group + subject (plus any covariates).

When specifying custom contrasts, use either "group1-group2" or "group1_vs_group2" format. If group names contain hyphens and you use the first format, an error will be raised suggesting to use the second format.

Value

A list with three elements:

• tidy_result: A tibble with limma results containing the following columns:

  • variable: Variable name

  • log2fc: Model contrast on the analysis scale (log2 fold change when expression values are on a log2 scale)

  • AveExpr: Average expression level

  • t: t-statistic

  • p_val: Raw p-value

  • p_adj: Adjusted p-value (if p_adj_method is not NULL)

  • b: B-statistic (log-odds of differential expression)

  • ref_group: Reference group

  • test_group: Test group

• raw_result: The raw limma fit object(s).

• meta_data: A list containing metadata from the input experiment.

Required packages

This function requires the following packages to be installed:

• limma for linear model fitting and empirical Bayes moderation

See Also

limma::lmFit(), limma::eBayes(), limma::makeContrasts()

Description

Fit a linear model to every variable in a glycomics or glycoproteomics experiment using limma's empirical Bayes moderation.

Usage

gly_linear_model(
  exp,
  formula,
  contrasts = NULL,
  p_adj_method = "BH",
  add_info = TRUE,
  ...,
  log_transform = TRUE
)

Arguments

Details

gly_linear_model() and gly_limma() both use limma under the hood, but serve complementary use cases. gly_linear_model() provides a general formula interface for specifying a wide range of analysis designs and custom contrasts. gly_limma() is a dedicated wrapper for common differential expression analysis tasks, with convenient group, covariate, subject, and pairwise-comparison arguments.

By default, expression values are transformed using log2(x + 1e-6) before modeling. Set log_transform = FALSE to use the supplied values unchanged. Coefficient estimates and average expression are reported on the selected analysis scale. When contrasts is NULL, all non-intercept coefficients are tested. When contrasts are supplied, only the named contrasts are returned.

Character predictors are handled as factors by stats::model.matrix(). Factor reference levels therefore determine the meaning of coefficients. Set factor levels in the sample information before calling this function when a specific reference level is required.

P-values are adjusted independently across variables for each coefficient or contrast.

Value

A list with classes glystats_linear_model_res and glystats_res containing:

• tidy_result: A tibble with variable, term, estimate, AveExpr, t, p_val, optional p_adj, and b columns.

• raw_result: A list containing the moderated limma fit, design matrix, contrast matrix, and coefficient-name mapping.

• meta_data: Metadata copied from exp.

See Also

gly_limma(), limma::lmFit(), limma::eBayes()

Examples

exp <- glyexp::real_experiment
complete <- rowSums(is.finite(SummarizedExperiment::assay(exp))) == ncol(exp)
exp <- exp[which(complete)[seq_len(10)], ]
SummarizedExperiment::colData(exp)$time <- factor(
  rep(c("early", "late"), length.out = ncol(exp))
)

model_res <- gly_linear_model(exp, ~ group * time)

contrast_res <- gly_linear_model(
  exp,
  ~ group * time,
  contrasts = c(
    C_late_vs_H_late = "groupC + `groupC:timelate`"
  )
)

Description

Perform orthogonal partial least squares discriminant analysis on the expression data. The function uses ropls::opls() to perform OPLS-DA and returns tidy results.

Usage

gly_oplsda(
  exp,
  group_col = "group",
  pred_i = 1,
  ortho_i = NA,
  scale = TRUE,
  add_info = TRUE,
  ...,
  log_transform = TRUE
)

Arguments

Details

By default, expression values are transformed using log(x + 1) before analysis. Set log_transform = FALSE to use the expression values as supplied, for example when they are already log-transformed.

Value

A list containing:

• tidy_result: A list of tibbles with OPLS-DA results:

  • samples: OPLS-DA scores for each sample containing the following columns:

    • sample: Sample name

    • group: Group assignment

    • p1, p2, etc.: Predictive component scores

    • o1, o2, etc.: Orthogonal component scores

  • variables: OPLS-DA loadings for each variable containing the following columns:

    • variable: Variable name

    • p1, p2, etc.: Predictive component loadings

    • o1, o2, etc.: Orthogonal component loadings

    • pcorr1, pcorr2, etc.: Correlation between each variable and the corresponding predictive component

  • variance: OPLS-DA explained variance containing the following columns:

    • component: Component name (p1, o1, etc.)

    • prop_var_explained: Proportion of variance explained by each component

    • cumulative_prop_var: Cumulative proportion of variance explained

  • vip: Variable Importance in Projection scores containing the following columns:

    • variable: Variable name

    • vip: VIP score

  • perm_test: Permutation test results containing the following columns:

    • model: Model type ("Original" for the original model, "Permutation" for permuted models)

    • perm_id: Permutation ID (0 for original model, 1+ for permutations)

    • Additional columns from the permutation test matrix (e.g., R2X, R2Y, Q2, etc.)

• raw_result: The raw ropls opls object from ropls::opls()

• meta_data: A list containing metadata from the input experiment

Required packages

This function requires the following packages to be installed:

• ropls for OPLS-DA analysis

See Also

ropls::opls()

Description

Perform principal component analysis on the expression data. The function uses prcomp() to perform PCA and broom::tidy() to tidy the results. If scale = TRUE, constant variables (zero variance) will be removed before PCA.

Usage

gly_pca(
  exp,
  center = TRUE,
  scale = TRUE,
  add_info = TRUE,
  ...,
  log_transform = TRUE
)

Arguments

Details

By default, the function performs log transformation on the expression data (log(x + 1)) before PCA analysis. Set log_transform = FALSE to use the expression values as supplied, for example when they are already log-transformed.

Value

A list containing:

• tidy_result: A list of tibbles with PCA results:

  • samples: PCA scores for each sample containing the following columns:

    • sample: Sample name

    • PC: Principal component name (PC1, PC2, etc.)

    • value: Score value for the principal component

  • variables: PCA loadings for each variable containing the following columns:

    • variable: Variable name

    • PC: Principal component name (PC1, PC2, etc.)

    • value: Loading value for the principal component

  • eigenvalues: PCA eigenvalues containing the following columns:

    • PC: Principal component name (PC1, PC2, etc.)

    • std.dev: Standard deviation

    • percent: Percentage of variance explained

    • cumulative: Cumulative percentage of variance explained

• raw_result: The raw prcomp object from stats::prcomp()

• meta_data: A list containing metadata from the input experiment

Required packages

This function only uses base R packages and does not require additional dependencies.

See Also

stats::prcomp()

Description

Perform partial least squares discriminant analysis on the expression data. The function uses ropls::opls() to perform PLS-DA and returns tidy results.

Usage

gly_plsda(
  exp,
  group_col = "group",
  ncomp = 2,
  scale = TRUE,
  add_info = TRUE,
  ...,
  log_transform = TRUE
)

Arguments

Details

By default, expression values are transformed using log(x + 1) before analysis. Set log_transform = FALSE to use the expression values as supplied, for example when they are already log-transformed.

Value

A list containing:

• tidy_result: A list of tibbles with PLS-DA results:

  • samples: PLS-DA scores for each sample containing the following columns:

    • sample: Sample name

    • group: Group assignment

    • p1, p2, etc.: PLS-DA component scores

  • variables: PLS-DA loadings for each variable containing the following columns:

    • variable: Variable name

    • p1, p2, etc.: PLS-DA component loadings

    • pcorr1, pcorr2, etc.: Correlation between each variable and component

  • variance: PLS-DA explained variance containing the following columns:

    • component: Component name (p1, p2, etc.)

    • prop_var_explained: Proportion of variance explained by each component

    • cumulative_prop_var: Cumulative proportion of variance explained

  • vip: Variable Importance in Projection scores containing the following columns:

    • variable: Variable name

    • vip: VIP score

  • perm_test: Permutation test results containing the following columns:

    • model: Model type ("Original" for the original model, "Permutation" for permuted models)

    • perm_id: Permutation ID (0 for original model, 1+ for permutations)

    • Additional columns from the permutation test matrix (e.g., R2X, R2Y, Q2, etc.)

• raw_result: The raw ropls opls object from ropls::opls()

• meta_data: A list containing metadata from the input experiment

Required packages

This function requires the following packages to be installed:

• ropls for PLS-DA analysis

See Also

ropls::opls()

Description

Perform Receiver Operating Characteristic (ROC) analysis for binary classification of glycomics or glycoproteomics data. The function calculates ROC curves and Area Under the Curve (AUC) values for each variable to assess their discriminatory power between two groups.

Usage

gly_roc(
  exp,
  group_col = "group",
  pos_class = NULL,
  add_info = TRUE,
  log_transform = FALSE
)

Arguments

Details

Expression values are used as supplied by default. Set log_transform = TRUE to apply log2(x + 1e-6) before analysis.

For each variable, a ROC curve is computed using the expression values as predictor and the binary group labels as response.

The function requires exactly 2 groups in the specified grouping variable. If more than 2 groups are present, an error will be thrown.

Underlying Function:

• ROC analysis is performed using pROC::roc()

• Coordinates are extracted using pROC::coords()

Value

A list with three elements:

• tidy_result: A list containing two tibbles:

  • auc: A tibble containing AUC values for each variable with the following columns:

    • variable: Variable name

    • auc: Area Under the Curve value

    • auc_ci_low: Lower bound of 95% confidence interval for AUC

    • auc_ci_high: Upper bound of 95% confidence interval for AUC

  • coords: A tibble containing ROC curve coordinates with the following columns:

    • variable: Variable name

    • threshold: Threshold value for classification

    • specificity: Specificity (True Negative Rate)

    • sensitivity: Sensitivity (True Positive Rate)

• raw_result: A list of pROC objects

• meta_data: A list containing metadata from the input experiment The list has classes glystats_roc_res and glystats_res.

Required packages

This function requires the pROC package to be installed for ROC curve computation.

See Also

pROC::roc(), pROC::coords()

Description

Construct and jointly test correlated variable sets with Hotelling's $T^2$ statistic.

Usage

gly_set_test(
  exp,
  sets = NULL,
  group_col = "group",
  subject_col = NULL,
  method = "hotelling",
  p_adj_method = "BH",
  threshold = 0.9,
  correlation = "pearson",
  clustering = "connected",
  within = NULL,
  log_transform = TRUE
)

Arguments

Details

By default, expression values are transformed using log2(x + 1e-6). The first factor level of group_col is the reference group and the second is the test group. The estimate is therefore test - reference.

Every usable variable is assigned to a set, including a one-variable set when it is not correlated above threshold with another variable. Structurally rank-deficient covariance matrices are standardized to unit member variance and tested in their nonredundant subspace using a Moore-Penrose inverse when the original set dimension satisfies the classical sample-size requirement and the observed mean contrast lies in the estimable subspace. Other ineligible sets are retained with status = "failed", NA statistics, and an explicit failure_reason.

Value

A list with classes glystats_set_test_res and glystats_res containing:

• tidy_result$sets: One row per set with estimates for every member, the effective covariance rank in test_dimension (NA when rejected before covariance estimation), Hotelling statistic, degrees of freedom, Mahalanobis effect size, p-value, adjusted p-value, and fit status.

• tidy_result$members: One row per set member with its marginal estimate and mean within-set correlation.

• raw_result: Set definitions, their correlation matrix, individual Hotelling test objects, and automatic set-construction details (when sets is NULL).

• meta_data: Metadata copied from exp.

See Also

stats::cor(), stats::hclust(), stats::p.adjust()

Examples

set.seed(1)
expression <- matrix(
  stats::rexp(36),
  nrow = 3,
  dimnames = list(c("A", "B", "C"), paste0("S", 1:12))
)
expression["B", ] <- expression["A", ] * exp(stats::rnorm(12, sd = 0.01))
exp <- SummarizedExperiment::SummarizedExperiment(
  assays = list(expression = expression),
  colData = S4Vectors::DataFrame(
    group = factor(
      rep(c("control", "case"), each = 6),
      levels = c("control", "case")
    )
  )
)
result <- gly_set_test(exp, threshold = 0.8)
result$tidy_result$sets

Description

Perform t-SNE dimensionality reduction on the expression data. The function uses Rtsne::Rtsne() to perform t-SNE analysis.

Usage

gly_tsne(
  exp,
  dims = 2,
  perplexity = 30,
  add_info = TRUE,
  ...,
  log_transform = TRUE
)

Arguments

Details

By default, expression values are transformed using log(x + 1) before analysis. Set log_transform = FALSE to use the expression values as supplied, for example when they are already log-transformed.

Value

A list with three elements:

• tidy_result: A tibble with t-SNE coordinates containing the following columns:

  • sample: Sample name

  • tsne1: First t-SNE dimension

  • tsne2: Second t-SNE dimension

• raw_result: The raw Rtsne object

• meta_data: A list containing metadata from the input experiment The list has classes glystats_tsne_res and glystats_res.

Required packages

This function requires the Rtsne package to be installed for t-SNE analysis.

See Also

Rtsne::Rtsne()

Description

Perform independent or paired t-tests for glycomics or glycoproteomics data. P-values are adjusted for multiple testing using the method specified by p_adj_method.

Usage

gly_ttest(
  exp,
  group_col = "group",
  p_adj_method = "BH",
  ref_group = NULL,
  add_info = TRUE,
  ...,
  subject_col = NULL,
  log_transform = TRUE
)

Arguments

Details

Tests and effect sizes use the selected analysis scale. Reported fold changes continue to use ratios of group means from the original supplied values.

By default, expression data are transformed using log2(x + 1e-6) before statistical testing. Set log_transform = FALSE to skip this transformation. Exactly 2 groups are required in the grouping variable. In paired analyses, only subjects observed in both groups are used and effect_size is Cohen's $d_z$ based on within-subject differences.

Value

A list with three elements:

• tidy_result: A tibble with t-test results containing the following columns:

  • variable: Variable name

  • estimate: Difference in group means (group2 - group1)

  • estimate1: Mean of group 1

  • estimate2: Mean of group 2

  • statistic: t-statistic

  • p_val: Raw p-value from t-test

  • parameter: Degrees of freedom

  • conf_low: Lower bound of 95% confidence interval

  • conf_high: Upper bound of 95% confidence interval

  • effect_size: Cohen's d, or Cohen's dz for paired analyses

  • method: Statistical method used

  • alternative: Alternative hypothesis

  • p_adj: Adjusted p-value (if p_adj_method is not NULL)

  • log2fc: Log2 fold change (log2(group2_mean / group1_mean))

• raw_result: A list of t.test model objects

• meta_data: A list containing metadata from the input experiment The list has classes glystats_ttest_res and glystats_res.

See Also

stats::t.test()

Description

Perform UMAP dimensionality reduction on the expression data. The function uses uwot::umap() to perform UMAP analysis.

Usage

gly_umap(
  exp,
  n_neighbors = 15,
  n_components = 2,
  add_info = TRUE,
  ...,
  log_transform = TRUE
)

Arguments

Details

By default, expression values are transformed using log(x + 1) before analysis. Set log_transform = FALSE to use the expression values as supplied, for example when they are already log-transformed.

Value

A list with three elements:

• tidy_result: A tibble with UMAP coordinates containing the following columns:

  • sample: Sample name

  • umap1: First UMAP dimension

  • umap2: Second UMAP dimension

  • umap3, umap4, etc.: Additional UMAP dimensions (if n_components > 2)

• raw_result: The raw UMAP result matrix

• meta_data: A list containing metadata from the input experiment The list has classes glystats_umap_res and glystats_res.

Required packages

This function requires the uwot package to be installed for UMAP analysis.

See Also

uwot::umap()

Description

Perform Wilcoxon rank-sum or signed-rank tests for glycomics or glycoproteomics data. P-values are adjusted for multiple testing using the method specified by p_adj_method.

Usage

gly_wilcox(
  exp,
  group_col = "group",
  p_adj_method = "BH",
  ref_group = NULL,
  add_info = TRUE,
  ...,
  subject_col = NULL,
  log_transform = TRUE
)

Arguments

Details

Tests and effect sizes use the selected analysis scale. Reported fold changes continue to use ratios of group means from the original supplied values.

By default, expression data are transformed using log2(x + 1e-6) before statistical testing. Set log_transform = FALSE to skip this transformation. Exactly 2 groups are required in the grouping variable. In paired analyses, only subjects observed in both groups are used and effect_size is the matched-pairs rank-biserial correlation.

Value

A list with three elements:

• tidy_result: A tibble with Wilcoxon test results containing the following columns:

  • variable: Variable name

  • statistic: Wilcoxon test statistic

  • p_val: Raw p-value from Wilcoxon test

  • effect_size: Rank-biserial correlation, using matched pairs in paired analyses

  • method: Statistical method used

  • alternative: Alternative hypothesis

  • p_adj: Adjusted p-value (if p_adj_method is not NULL)

  • log2fc: Log2 fold change (log2(group2_mean / group1_mean)) Additional columns from experiment metadata may be included if add_info = TRUE.

• raw_result: A list of wilcox.test model objects

• meta_data: A list containing metadata from the input experiment The list has classes glystats_wilcox_res and glystats_res.

See Also

stats::wilcox.test()