# Olink® Analyze Vignette

#### Olink DS team

#### 2025-09-24

Olink® Analyze is an R package that provides a versatile toolbox to enable fast and easy handling of Olink® NPX data for your proteomics research. Olink® Analyze provides functions for using Olink data, including functions for importing Olink® NPX datasets, as well as quality control (QC) plot functions and functions for various statistical tests. This package is meant to provide a convenient pipeline for your Olink NPX data analysis.

# Installation

You can install Olink® Analyze from CRAN.

    install.packages("OlinkAnalyze")

# List of functions

Preprocessing

- read_NPX Function to read NPX data into long format
- olink_plate_randomizer Randomize samples on plate
- olink_bridgeselector Select bridge samples
- olink_normalization Normalization of all proteins (by OlinkID)
- olink_lod Calculation of Limit of Detection for Explore data

Statistical analysis

- olink_ttest Function which performs a t-test per protein
- olink_wilcox Function which performs a Mann-Whitney U Test per protein
- olink_anova Function which performs an ANOVA per protein
- olink_anova_posthoc Function which performs an ANOVA post-hoc test per protein
- olink_one_non_parametric Function which performs a Kruskal-Wallis Test or Friedman Test per protein
- olink_one_non_parametric_posthoc Function which performs post-hoc test for one way non-parametric test
- olink_ordinalRegression Function which performs an ordinal regression per protein
- olink_ordinalRegression_posthoc Function which performs an ordinal regression post-hoc test per protein
- olink_lmer Function which performs a linear mixed model per protein
- olink_lmer_posthoc Function which performs a linear mixed model post-hoc per protein
- olink_pathway_enrichment Function which performs GSEA or ORA pathway enrichment using outcome from other statistical tests

Visualization

- olink_boxplot Function which plots boxplots of a selected variable
- olink_dist_plot Function to plot the NPX distribution by panel
- olink_lmer_plot Function which performs a point-range plot per protein on a linear mixed model
- olink_pathway_visualization Function which plots a bar graph for pathways of interest
- olink_pathway_heatmap Function which plots estimates of proteins associated with pathways of interest
- olink_pca_plot Function to plot a PCA of the data
- olink_qc_plot Function to plot an overview of a sample cohort per Panel
- olink_umap_plot Function to plot a UMAP of the data
- olink_volcano_plot Easy volcano plot with Olink theme
- olink_heatmap_plot Function which generates a heatmap over all proteins
- set_plot_theme Function to set plot theme
- olink_bridgeability_plot Function which generates plots that illustrate the three criteria for determining whether an assay is bridgeable in cross-product bridge normalization.

Sample datasets

- npx_data1 NPX Data in Long format
- npx_data2 NPX Data in Long format, Follow-up
- manifest A sample manifest including Sample ID, Subject ID and clinical variables

# Usage

## Load the library

    # Load OlinkAnalyze
    library(OlinkAnalyze)

    # Load other libraries used in Vignette
    library(dplyr)
    library(ggplot2)
    library(stringr)

## Introduction to Olink NPX data format

The package contains two test data files named npx_data1 and npx_data2. These are synthetic datasets that resemble Olink® data accompanied by clinical variables. Olink® data that is delivered in long format or imported with the function read_NPX (that converts the data into a long format) contain the following columns:

- SampleID <chr>: Sample names or IDs.
- OlinkID <chr>: Unique ID for each assay assigned by Olink. In case the assay is included in more than one panels it will have a different OlinkID in each one.
- UniProt <chr>: UniProt ID.
- Assay <chr>: Common gene name for the assay.
- MissingFreq <dbl>: Missing frequency for the OlinkID, i.e. frequency of samples with NPX value below limit of detection (LOD).
- Panel <chr>: Olink Panel that samples ran on. Read more about Olink Panels here: https://olink.com/products/compare.
- Panel_Version <chr>: Version of the panel. A new panel version might include some different or improved assays.
- PlateID <chr>: Name of the plate.
- QC_Warning or SampleQC <chr>: Indication whether the sample passed Olink QC. More information about Olink quality control metrics can be found in our FAQ by searching “Quality control”.
- LOD <dbl>: Limit of detection (LOD) is the minimum level of an individual protein that can be measured. LOD is defined as 3 times the standard deviation over background.
- NPX <dbl>: Normalized Protein eXpression, is Olink®’s unit of protein expression level in a log2 scale. The majority of the functions of this package use NPX values for calculations. Read more about NPX in the Olink FAQ (Search “What is NPX?”) or in Olink’s Data normalization and standardization white paper.

Note: There are 5 additional variables in the sample datasets npx_data1 and npx_data2 that include clinical or other information, namely: Subject <chr>, Treatment <chr>, Site <chr>, Time <chr>, Project <chr>.

The columns found in an Olink data set may vary based on the version and product.

# Preprocessing

## Read NPX data (read_NPX)

The read_NPX function imports an NPX file into a tidy format to work with in R. This function supports Olink® NPX files generated by Olink® data software in CSV, Excel, and Parquet formats. No prior alterations to the NPX output file should be made for this function to work as expected.

### Function arguments

- filename: Path to the NPX output file.

    data <- read_NPX("~/NPX_file_location.xlsx")

### Function output

A tibble in long format containing:

- SampleID: Sample names or IDs.
- Index: Unique number for each SampleID. It is used to make up for non unique sample IDs.
- OlinkID: Unique ID for each assay assigned by Olink. In case the assay is included in more than one panels it will have a different OlinkID in each one.
- UniProt: UniProt ID.
- Assay: Common gene name for the assay.
- MissingFreq: Missing frequency for the OlinkID, i.e. frequency of samples with NPX value below limit of detection (LOD).
- Panel: Olink Panel that samples ran on. Read more about Olink Panels here: https://olink.com/products/compare
- Panel_Version: Version of the panel. A new panel version might include some different or improved assays.
- PlateID: Name of the plate.
- QC_Warning: Indication whether the sample passed Olink QC. More information about Olink quality control metrics can be found in our FAQ (Search “Quality control”).
- LOD: Limit of detection (LOD) is the minimum level of an individual protein that can be measured. LOD is defined as 3 times the standard deviation over background.
- NPX: Normalized Protein eXpression, is Olink’s unit of protein expression level in a log2 scale. The majority of the functions of this package use NPX values for calculations. Read more about NPX in the Olink FAQ (Search “What is NPX?”) or in Olink’s Data normalization and standardization white paper.

## Read multiple NPX data files (read_NPX)

In order to import multiple NPX data files at once, the read_NPX function can be used in combination with the list.files, lapply and dplyr::bind_rows functions, as seen below. The pattern argument of the list.files function specifies the NPX file format (.csv, .xlsx, .parquet, or any combination of these). This method requires that all NPX files are stored in the same folder and have identical column names. No prior alterations to the NPX output file should be made for this method to work as expected.

    # Read in multiple NPX files in .csv format
    data <-  list.files(path = "path/to/dir/with/NPX/files",
                        pattern = "csv$",
                        full.names = TRUE) |>
             lapply(FUN = function(x){
        OlinkAnalyze::read_NPX(x) |>
          # Optionally add additional columns to add file identifiers
          dplyr::mutate(File = x)
          }  |>
             # optional to return a single data frame of all files instead of a list of data frames
             dplyr::bind_rows()

    # Read in multiple NPX files in .parquet format
    data <-  list.files(path = "path/to/dir/with/NPX/files",
                        pattern = "parquet$",
                        full.names = TRUE) |>
             lapply(OlinkAnalyze::read_NPX)  |>
             dplyr::bind_rows()

    # Read in multiple NPX files in either format
    data <-  list.files(path = "path/to/dir/with/NPX/files",
                        pattern = "parquet$|csv$",
                        full.names = TRUE) |>
             lapply(OlinkAnalyze::read_NPX)  |>
             dplyr::bind_rows()

## Randomize samples on plate (olink_plate_randomizer)

The olink_plate_randomizer function randomly assigns samples to a plate well with the option to keep the same individuals on the same plate. Olink® does not recommend to force balance based on other clinical variables.

For more information on plate randomization, consult the Plate Randomization Vignette.

## Select bridge samples (olink_bridgeselector)

The bridge selection function selects a number of bridge samples based on the input data. Bridge samples are used to normalize two dataframes/projects that have been ran at different time points, hence, a batch effect is expected. It selects samples that have good detectability (if applicable), pass quality control, and cover a wide range of data points.

For more information on bridge sample selection, consult the Introduction to bridging Olink® NPX datasets tutorial.

## Normalizing NPX data (olink_normalization)

The Olink® normalization function normalizes NPX values between two different datasets or one Olink® dataset to a set of reference medians.

The function handles four different types of normalization:

- Bridging normalization: One of the dataframes is adjusted to another using overlapping samples (bridge samples). The overlapping samples should have the same IDs between dataframes. Adjustment is made using the median of the paired differences between the bridge samples. For more information on bridging, consult the Introduction to bridging Olink® NPX datasets tutorial.
- Subset normalization: A subset of samples is used to normalize two dataframes, one of which is used as a reference_project. Adjustment is made using the differences of medians between the sample subsets from the two dataframes. Subset normalization is useful if no bridge samples were included and one can assume that the distribution of the two datasets is very similar.
- Reference median normalization: Works only on one dataframe. This is effectively subset normalization, but using difference of medians to pre-recorded median values. df1, overlapping_samples_df1 and reference_medians need to be specified.
- Cross-product bridging: Similar to bridging normalization but bridging across products, for example bridging Explore 3072 data to Explore HT data. Overlapping samples are run on both products and used to determine which assays are bridgeable and what method should be used to bridge each assay. For more information on the between-product bridging methodology, consult the Bridging across NGS-based Olink® products Tutorial.

## Integrating Explore NPX LOD (olink_lod)

The olink_lod function adds LOD information to an Explore HT or Explore 3072 NPX dataframe. This function can incorporate LOD based on either an Explore dataset’s negative controls or using predetermined fixed LOD values, which can be downloaded from the Document Download Center at olink.com, or using both methods. The default LOD calculation method is based off of the negative controls. If an NPX file is intensity normalized, both intensity normalized and PC normalized LODs are provided.

For more information on calculating LOD, consult the Calculating LOD from Olink® Explore data tutorial.

# Statistical analysis

## T-test analysis (olink_ttest)

The olink_ttest function performs a Welch 2-sample t-test or paired t-test at confidence level 0.95 for every protein (by OlinkID) for a given grouping variable. It corrects for multiple testing using the Benjamini-Hochberg method (“fdr”). Adjusted p-values are logically evaluated towards adjusted p-value < 0.05. The resulting t-test table is arranged by ascending p-values.

### Function arguments

- df: NPX data frame in long format should minimally contain protein name (Assay), OlinkID, UniProt, Panel and an outcome factor with 2 levels.
- variable: Character value that should represent a column in the df to be used as a grouping variable. Needs to have exactly 2 levels.
- pair_id: Character value indicating which column contains the paired sample identifier. Only used for paired t-tests.

    olink_ttest(df = npx_data1,
                variable = 'Treatment')

### Function output

A tibble with the following columns:

- Assay <chr>: Assay name.
- OlinkID <chr>: Unique Olink® ID.
- UniProt <chr>: UniProt ID.
- Panel <chr>: Olink® Panel.
- estimate <dbl>: Difference in mean NPX between groups.
- statistic <dbl>: Value of the t-statistic.
- p.value <dbl>: P-value for the test.
- parameter <dbl>: Degrees of freedom for the t-statistic.
- conf.low <dbl>: Low bound of the confidence interval for the mean.
- conf.high <dbl>: High bound of the confidence interval for the mean.
- method <chr>: Method that was used.
- alternative <chr>: : Description of the alternative hypothesis.
- Adjusted_pval <dbl>: Adjusted p-value for the test (Benjamini & Hochberg).
- Threshold <chr>: Text indication if assay is significant (adjusted p-value < 0.05).

## Mann-Whitney U Test analysis (olink_wilcox)

The olink_wilcox function performs a 2-sample Mann-Whitney U test or paired Mann-Whitney U test at confidence level 0.95 for every protein (by OlinkID) for a given grouping variable. It corrects for multiple testing using the Benjamini-Hochberg method (“fdr”). Adjusted p-values are logically evaluated towards adjusted p-value<0.05. The resulting Mann-Whitney U table is arranged by ascending p-values.

### Function arguments

- df: NPX data frame in long format should minimally contain protein name (Assay), OlinkID, UniProt, Panel and an outcome factor with 2 levels.
- variable: Character value that should represent a column in the df to be used as a grouping variable. Needs to have exactly 2 levels.
- pair_id: Character value indicating which column contains the paired sample identifier. Only used for paired Mann-Whitney U tests.

    olink_wilcox(df = npx_data1,
                 variable = 'Treatment')

### Function output

A tibble with the following columns:

- Assay <chr>: Assay name.
- OlinkID <chr>: Unique Olink® ID.
- UniProt <chr>: UniProt ID.
- Panel <chr>: Olink® Panel.
- statistic <dbl>: Value of the Mann-Whitney U statistic.
- p.value <dbl>: P-value for the test.
- method <chr>: Method that was used.
- alternative <chr>: : Description of the alternative hypothesis.
- Adjusted_pval <dbl>: Adjusted p-value for the test (Benjamini & Hochberg).
- Threshold <chr>: Text indication if assay is significant (adjusted p-value < 0.05).

## Analysis for variance (ANOVA) (olink_anova)

The olink_anova function performs an ANOVA F-test for each assay (by OlinkID) using Type III sum of squares. The function handles both factor and numerical variables, and/or confounding factors.

Samples with missing variable information or factor levels are excluded from the analysis. Character columns in the input data frame are converted to factors.

Control samples and control assays should be removed before using this function.

Crossed/interaction analysis, i.e. A*B formula notation, is inferred from the variable argument in the following cases:

- c(‘A’,‘B’)
- c(‘A:B’)
- c(‘A:B’, ‘B’) or c(‘A:B’, ‘A’)

For covariates, crossed analyses need to be specified explicitly, i.e. two main effects will not be expanded with a c(‘A’,‘B’) notation. Main effects present in the variable take precedence.

Adjusted p-values are calculated using the Benjamini & Hochberg (1995) method (“fdr”). The threshold is determined by logic evaluation of Adjusted_pval < 0.05. Covariates are not included in the p-value adjustment.

### Function arguments

- df: NPX data frame in long format should minimally contain protein name (Assay), OlinkID, UniProt, Panel and an outcome factor with at least 3 levels.
- variable: Single character value or character array. In case of single character then that should represent a column in the df. Otherwise, if length > 1, the included variable names will be used in crossed analyses. It can also accept the notations ‘:’ or ‘*’.
- outcome: Name of the column from df that contains the dependent variable. Default: NPX.
- covariates: Single character value or character array. Default: NULL. Confounding factors to include in the analysis. In case of single character then that should represent a column in the df. It can also accept the notations ‘:’ or ‘*’, while crossed analysis will not be inferred from main effects.
- return.covariates: Logical. Default: False. Returns F-test results for the covariates. Note: Adjusted p-values will be NA for covariates.
- verbose: Logical. Default: True. If information about removed samples, factor conversion and final model formula is to be printed to the console.

    # Remove control samples and assays
    npx_data1_no_controls <- npx_data1 |>
      dplyr::filter(!stringr::str_detect(SampleID,
                                         regex("control|ctrl",
                                               ignore_case = TRUE))) |>
      dplyr::filter(!stringr::str_detect(Assay,
                                         regex("control|ctrl",
                                               ignore_case = TRUE)))

    # One-way ANOVA, no covariates
    anova_results_oneway <- olink_anova(df = npx_data1_no_controls,
                                        variable = 'Site')
    # Two-way ANOVA, no covariates
    anova_results_twoway <- olink_anova(df = npx_data1_no_controls,
                                        variable = c('Site', 'Time'))
    # One-way ANOVA, Treatment as covariates
    anova_results_oneway <- olink_anova(df = npx_data1_no_controls,
                                        variable = 'Site',
                                        covariates = 'Treatment')

### Function output

A tibble with the following columns:

- Assay <chr>: Assay name.
- OlinkID <chr>: Unique Olink ID.
- UniProt <chr>: UniProt ID.
- Panel <chr>: Olink Panel.
- term <chr>: Name of the variable that was used for the p-value calculation. The “:” between variables indicates interaction between variables.
- df <dbl>: Numerator of degrees of freedom.
- sumsq <dbl>: Sum of squares.
- meansq <dbl>: Mean of squares.
- statistic <dbl>: Value of F-statistic.
- p.value <dbl>: P-value for the test.
- Adjusted_pval <dbl>: Adjusted p-value for the test (Benjamini & Hochberg).
- Threshold <chr>: Text indication if assay is significant (adjusted p-value < 0.05).

## Post-hoc ANOVA analysis (olink_anova_posthoc)

olink_anova_posthoc performs a post-hoc ANOVA test with Tukey p-value adjustment per assay (by OlinkID) at confidence level 0.95.

The function handles both factor and numerical variables and/or covariates. The post-hoc test for a numerical variable compares the difference in means of the outcome variable (default: NPX) for 1 standard deviation (SD) difference in the numerical variable, e.g. mean NPX at mean (numerical variable) versus mean NPX at mean (numerical variable) + 1*SD (numerical variable).

Control samples and control assays (AssayType is not “assay”, or Assay contains “control” or “ctrl”) should be removed before using this function.

### Function arguments

- df: NPX data frame in long format should minimally contain protein name (Assay), OlinkID, UniProt, Panel and an outcome factor with at least 3 levels.
- olinkid_list: Character vector of OlinkID’s on which to perform the post-hoc analysis. If not specified, all assays in df are used.
- variable: Single character value or character array. In case of single character then that should represent a column in the df. Otherwise, if length > 1, the included variable names will be used in crossed analyses. It can also accept the notations ‘:’ or ‘*’.
- covariates: Single character value or character array. Default: NULL. Confounding factors to include in the analysis. In case of single character then that should represent a column in the df. It can also accept the notations ‘:’ or ‘*’, while crossed analysis will not be inferred from main effects.
- outcome: Name of the column from df that contains the dependent variable. Default: NPX.
- effect: Term on which to perform the post-hoc analysis. Character vector. Must be subset of or identical to the variable and no adjustment is performed.
- mean_return: Logical. If true, returns the mean of each factor level rather than the difference in means (default). Note that no p-value is returned for mean_return = TRUE.
- verbose: Logical. Default: True. If information about removed samples, factor conversion and final model formula is to be printed to the console.

    # calculate the p-value for the ANOVA
    anova_results_oneway <- olink_anova(df = npx_data1_no_controls,
                                        variable = 'Site')

    # extracting the significant proteins
    anova_results_oneway_significant <- anova_results_oneway %>%
      dplyr::filter(Threshold == 'Significant') %>%
      dplyr::pull(OlinkID)

    anova_posthoc_oneway_results <- olink_anova_posthoc(df = npx_data1_no_controls,
                                                        olinkid_list = anova_results_oneway_significant,
                                                        variable = 'Site',
                                                        effect = 'Site')

### Function output

A tibble with the following columns:

- Assay <chr>: Assay name.
- OlinkID <chr>: Unique Olink ID.
- UniProt <chr>: UniProt ID.
- Panel <chr>: Olink Panel.
- term <chr>: Name of the variable that was used for the p-value calculation. The “:” between variables indicates interaction between variables.
- contrast <chr>: Variables (in term) that are compared.
- estimate <dbl>: Difference in mean NPX between variables (from contrast).
- conf.low <dbl>: Low bound of the confidence interval for the mean.
- conf.high <dbl>: High bound of the confidence interval for the mean.
- Adjusted_pval <dbl>: Adjusted p-value for the test (Benjamini & Hochberg).
- Threshold <chr>: Text indication if assay is significant (adjusted p-value < 0.05).

## Linear mixed effects model analysis (olink_lmer)

The olink_lmer fits a linear mixed effects model for every protein (by OlinkID) in every panel. The function handles both factor and numerical variables and/or covariates.

Samples with missing variable information or factor levels are excluded from the analysis. Character columns in the input data frame are converted to factors.

Crossed/interaction analysis, i.e. A*B formula notation, is inferred from the variable argument in the following cases:

- c(‘A’,‘B’)
- c(‘A:B’)
- c(‘A:B’, ‘B’) or c(‘A:B’, ‘A’)

For covariates, crossed analyses need to be specified explicitly, i.e. two main effects will not be expanded with a c(‘A’,‘B’) notation. Main effects present in the variable take precedence.

Adjusted p-values are calculated using the Benjamini & Hochberg (1995) method (“fdr”). The threshold is determined by logic evaluation of Adjusted_pval < 0.05. Covariates are not included in the p-value adjustment.

### Function arguments

- df: NPX data frame in long format should minimally contain protein name (Assay), OlinkID, UniProt, Panel and 1-2 variables with at least 2 levels and subject ID.
- variable: Single character value or character array. In case of single character then that should represent a column in the df. Otherwise, if length > 1, the included variable names will be used in crossed analyses. It can also accept the notations ‘:’ or ‘*’.
- outcome: Name of the column from df that contains the dependent variable. Default: NPX.
- random: Single character value or character array with random effects.
- covariates: Single character value or character array. Default: NULL. Confounding factors to include in the analysis. In case of single character then that should represent a column in the df. It can also accept the notations ‘:’ or ‘*’, while crossed analysis will not be inferred from main effects.
- return.covariates: Logical. Default: False. Returns F-test results for the covariates. Note: Adjusted p-values will be NA for covariates.
- verbose: Logical. Default: True. If information about removed samples, factor conversion and final model formula is to be printed to the console.

    if (requireNamespace("lme4", quietly = TRUE) & requireNamespace("lmerTest", quietly = TRUE)){
    # Linear mixed model with one variable.
      lmer_results_oneway <- olink_lmer(df = npx_data1,
                                        variable = 'Site',
                                        random = 'Subject')
    # Linear mixed model with two variables.
      lmer_results_twoway <- olink_lmer(df = npx_data1,
                                        variable = c('Site', 'Treatment'),
                                        random = 'Subject')
    }

### Function outcome

A tibble with the following columns:

- Assay <chr>: Assay name.
- OlinkID <chr>: Unique Olink ID.
- UniProt <chr>: UniProt ID.
- Panel <chr>: Olink Panel.
- term <chr>: Name of the variable that was used for the p-value calculation. The “:” between variables indicates interaction between variables.
- sumsq <dbl>: Sum of squares.
- meansq <dbl>: Mean of squares.
- NumDF <dbl>: Numerator of degrees of freedom.
- DenDF <dbl>: Denominator of degrees of freedom.
- statistic <dbl>: Value of F-statistic.
- p.value <dbl>: P-value for the test.
- Adjusted_pval <dbl>: Adjusted p-value for the test (Benjamini & Hochberg).
- Threshold <chr>: Text indication if assay is significant (adjusted p-value < 0.05).

## Post-hoc linear mixed effects model analysis (olink_lmer_posthoc)

The olink_lmer_posthoc function is similar to olink_lmer but performs a post-hoc analysis based on a linear mixed model effects model. The function handles both factor and numerical variables and/or covariates. Differences in estimated marginal means are calculated for all pairwise levels of a given output variable. Degrees of freedom are estimated using Satterthwaite’s approximation. The post-hoc test for a numerical variable compares the difference in means of the outcome variable (default: NPX) for 1 standard deviation difference in the numerical variable, e.g. mean NPX at mean(numerical variable) versus mean NPX at mean(numerical variable) + 1*SD(numerical variable). The output tibble is arranged by ascending adjusted p-values.

### Function arguments

- df: NPX data frame in long format should minimally contain protein name (Assay), OlinkID, UniProt, Panel and 1-2 variables with at least 2 levels and subject ID.
- variable: Single character value or character array. In case of single character then that should represent a column in the df. Otherwise, if length > 1, the included variable names will be used in crossed analyses. It can also accept the notations ‘:’ or ‘*’.
- olinkid_list: Character vector of OlinkID’s on which to perform the post-hoc analysis. If not specified, all assays in df are used.
- effect: Term on which to perform the post-hoc analysis. Character vector. Must be subset of or identical to the variable.
- outcome: Name of the column from df that contains the dependent variable. Default: NPX.
- random: Single character value or character array with random effects.
- covariates: Single character value or character array. Default: NULL. Confounding factors to include in the analysis. In case of single character then that should represent a column in the df. It can also accept the notations ‘:’ or ‘*’, while crossed analysis will not be inferred from main effects.
- mean_return: Logical. If true, returns the mean of each factor level rather than the difference in means (default). Note that no p-value is returned for mean_return = TRUE and no adjustment is performed.
- verbose: Logical. Default: True. If information about removed samples, factor conversion and final model formula is to be printed to the console.

    if (requireNamespace("lme4", quietly = TRUE) & requireNamespace("lmerTest", quietly = TRUE)){
      # Linear mixed model with two variables.
      lmer_results_twoway <- olink_lmer(df = npx_data1,
                                        variable = c('Site', 'Treatment'),
                                        random = 'Subject')
      # extracting the significant proteins
      lmer_results_twoway_significant <- lmer_results_twoway %>%
        dplyr::filter(Threshold == 'Significant', term == 'Treatment') %>%
        dplyr::pull(OlinkID)
      # performing post-hoc analysis
      lmer_posthoc_twoway_results <- olink_lmer_posthoc(df = npx_data1,
                                                        olinkid_list = lmer_results_twoway_significant,
                                                        variable = c('Site', 'Treatment'),
                                                        random = 'Subject',
                                                        effect = 'Treatment')
    }

### Function output

A tibble with the following columns:

- Assay <chr>: Assay name.
- OlinkID <chr>: Unique Olink ID.
- UniProt <chr>: UniProt ID.
- Panel <chr>: Olink Panel.
- term <chr>: Name of the variable that was used for the p-value calculation. The “:” between variables indicates interaction between variables.
- contrast <chr>: Variables (in term) that are compared.
- estimate <dbl>: Difference in mean NPX between variables (from contrast).
- conf.low <dbl>: Low bound of the confidence interval for the mean.
- conf.high <dbl>: High bound of the confidence interval for the mean.
- Adjusted_pval <dbl>: Adjusted p-value for the test (Benjamini & Hochberg).
- Threshold <chr>: Text indication if assay is significant (adjusted p-value < 0.05).

## Additional Statistical Tests

Many other statistical functions can be found within Olink Analyze, including:

- olink_one_non_parametric Function which performs a Kruskal-Wallis Test or Friedman Test per protein
- olink_one_non_parametric_posthoc Function which performs post-hoc test for one way non-parametric test
- olink_ordinalRegression Function which performs an ordinal regression per protein
- olink_ordinalRegression_posthoc Function which performs an ordinal regression post-hoc test per protein

To learn more about these function, consult their help documentation using the `help()` function.

## Pathway Enrichment (olink_pathway_enrichment)

The olink_pathway_enrichment function can be used to perform Gene Set Enrichment Analysis (GSEA) or Over-representation Analysis (ORA) using MSigDB, Reactome, KEGG, or GO. MSigDB includes curated gene sets (C2) and ontology gene sets (C5) which encompasses Reactome, KEGG, and GO. This function performs enrichment using the gsea or enrich functions from clusterProfiler from BioConductor. The function uses the estimate from a previous statistical analysis for one contrast for all proteins. MSigDB is subset if ontology is KEGG, GO, or Reactome. test_results must contain estimates for all assays. Posthoc results can be used but should be filtered for one contrast to improve interpretability.

Alternative statistical results can be used as input as long as they include the columns “OlinkID”, “Assay”, and “estimate”. A column named “Adjusted_pal” is also needed for ORA. Any statistical results that contains one estimate per protein will work as long as the estimates are comparable to each other.

### Function Arguments

- data: NPX data frame in long format with columns Assay, OlinkID, UniProt, SampleID, QC_Warning, NPX, and LOD
- test_results: a data frame of statistical test results including Adjusted_pval and estimate columns
- method: String of method name. Must be either “GSEA” (default) or “ORA”
- ontology: String of database to query. Must be either “MSigDb”, “KEGG”, “GO”, and “Reactome”
- organism: String of name of organism. Must be either “human” or “mouse”

    npx_df <- npx_data1 |>
      dplyr::filter(!grepl("control", SampleID, ignore.case = TRUE))

    ttest_results <- olink_ttest(
      df = npx_df,
      variable = "Treatment",
      alternative = "two.sided")

    try({ # This expression might fail if dependencies are not installed
    gsea_results <- olink_pathway_enrichment(data = npx_data1, test_results = ttest_results)
    ora_results <- olink_pathway_enrichment(
      data = npx_data1,
      test_results = ttest_results, method = "ORA")
    }, silent = TRUE)

### Function Output

A data frame of enrichment results. Columns for ORA include:

- ID <chr>: Pathway ID from MSigDB
- Description <chr>: Description of Pathway from MSigDB
- GeneRatio <chr>: ratio of input proteins that are annotated in a term
- BgRatio <chr>: ratio of all genes that are annotated in this term
- pvalue <dbl>: p-value of enrichment
- p.adjust <dbl>: Adjusted p-value (Benjamini-Hochberg)
- qvalue <dbl>: false discovery rate, the estimated probability that the normalized enrichment score represents a false positive finding
- geneID: <chr> list of input proteins (Gene Symbols) annotated in a term delimited by “/”
- Count <dbl>: Number of input proteins that are annotated in a term

Columns for GSEA:

- ID <chr>: Pathway ID from MSigDB
- Description <chr>: Description of Pathway from MSigDB
- setSize <dbl>: ratio of input proteins that are annotated in a term
- enrichmentScore <dbl>: Enrichment score, degree to which a gene set is over-represented at the top or bottom of the ranked list of genes
- NES <dbl>: Normalized Enrichment Score, normalized to account for differences in gene set size and in correlations between gene sets and expression data sets. NES can be used to compare analysis results across gene sets.
- pvalue <dbl>: p-value of enrichment
- p.adjust <dbl>: Adjusted p-value (Benjamini-Hochberg)
- qvalue <dbl>: false discovery rate, the estimated probability that the normalized enrichment score represents a false positive finding
- rank <dbl>: the position in the ranked list where the maximum enrichment score occurred
- leading_edge <chr>: contains tags, list, and signal. Tags gives an indication of the percentage of genes contributing to the enrichment score. List gives an indication of where in the list the enrichment score is obtained. Signal represents the enrichment signal strength and combines the tag and list.
- core_enrichment <chr>: list of input proteins (Gene Symbols) annotated in a term delimited by “/”

# Exploratory analysis

## Principal components analysis (PCA) plot (olink_pca_plot)

Generates PCA projection of all samples from NPX data along two principal components (default PC2 vs PC1) colored by the variable specified by color_g (default QC_Warning) and including the percentage of explained variance. By default, the values are scaled and centered in the PCA and proteins with missing NPX values removed from the corresponding assay(s). Unique sample names are required. Imputation by median value is done for assays with missingness <10% and for multi-plate projects, and for missingness <5% for single plate projects.

More information about `olink_pca()` can be found in the Outlier Exclusion Vignette.

## Uniform Manifold Approximation and Projection for Dimension Reduction (UMAP) (olink_umap_plot)

Computes a manifold approximation and projection and plots the two specified components. Unique sample names are required and imputation by the median is done for assays with missingness <10% for multi-plate projects and <5% for single plate projects.

The arguments outlierDefX and outlierDefY can be used to identify outliers in the UMAP results. Sample outliers will be labelled.

NOTE: UMAP is a non-linear data transformation that might not accurately preserve the properties of the data. Distances in the UMAP plane should therefore be interpreted with caution

### Function arguments (selection)

- df: NPX data frame in long format should minimally contain SampleID, NPX and column that will be used for grouping/coloring.
- color_g: Character value indicating the column name that should be used as fill color. Default QC_Warning.
- x_val: Integer indicating which principal component to plot along the x-axis. Default 1.
- y_val: Integer indicating which principal component to plot along the y-axis. Default 2.
- config: Object of class umap.config, specifying the parameters for the UMAP algorithm.
- label_samples: Logical. If TRUE, points are replaced with SampleID. Default FALSE.
- drop_assays: Logical. All assays with any missing values will be dropped. Takes precedence over sample drop.
- drop_samples: Logical. All samples with any missing values will be dropped.
- byPanel: Logical. Perform the UMAP per panel (default FALSE)
- outlierDefX: (Optional) The number standard deviations along the UMAP dimension plotted on the x-axis that defines an outlier.
- outlierDefY: (Optional) The number standard deviations along the UMAP dimension plotted on the y-axis that defines an outlier.
- OutlierLines: Logical. Draw dashed lines at +/-outlierDef[X,Y] standard deviations from the mean of the plotted UMAP dimensions (default FALSE)
- verbose: Logical. Default: True. If information about removed samples, factor conversion and final model formula is to be printed to the console.
- quiet: Logical. Default: False. If TRUE, the resulting plot is not printed.

    if (requireNamespace("umap", quietly = TRUE) ){
    npx_data1 %>%
      dplyr::filter(!stringr::str_detect(SampleID, 'CONTROL_SAMPLE')) %>%
      olink_umap_plot(df = .,
                     color_g = "QC_Warning", byPanel = TRUE)
      }

    npx_data1 %>%
      dplyr::filter(!stringr::str_detect(SampleID, 'CONTROL_SAMPLE')) %>%
      olink_umap_plot(df = .,
                     color_g = "QC_Warning", byPanel = TRUE)

### Function output

A list of objects of class ggplot (silently returned). Plots are also printed unless option `quiet = TRUE` is set.

# Visualization

## Boxplots for outcomes (olink_boxplot)

The olink_boxplot function is used to generate boxplots of NPX values stratified on a variable for a given list of proteins. In order to annotate the plot with ANOVA posthoc analysis results (i.e. include statistical asterisks in the plot), control samples and control assays should be removed from the data.

### Function arguments

- df: NPX data frame in long format should minimally contain protein name (Assay), OlinkID, UniProt and a grouping variable.
- variable: Single character value indicating the column name to use as a grouping variable in the x axis.
- olinkid_list: Character vector of OlinkID’s that should be used for the boxplot. If not specified, all assays in df are used.
- posthoc_results: Data frame from ANOVA posthoc analysis. This data frame need to be generated using the olink_anova_posthoc() function.
- ttest_results: Data frame from ttest analysis. This data frame need to be generated using the olink_ttest() function.
- verbose: Logical. Default: False. Flag indicating if plots shall be printed additionally to assigned to a list variable.
- number_of_proteins_per_plot: Number of boxplots to include in the facets plot. Default 6.

    # Remove control samples and assays
    npx_data1_no_controls <- npx_data1 |>
      dplyr::filter(!stringr::str_detect(SampleID,
                                         regex("control|ctrl",
                                               ignore_case = TRUE))) |>
      dplyr::filter(!stringr::str_detect(Assay,
                                         regex("control|ctrl",
                                               ignore_case = TRUE)))

    plot <- npx_data1_no_controls |>
      dplyr::filter(!is.na(Site)) |> # removing missing values which exist for Site
      olink_boxplot(variable = "Site",
                    olinkid_list = c("OID00488", "OID01276"),
                    number_of_proteins_per_plot  = 2)

    plot[[1]]

    anova_posthoc_results<-npx_data1_no_controls %>%
      olink_anova_posthoc(olinkid_list = c("OID00488", "OID01276"),
                          variable = 'Site',
                          effect = 'Site')

    plot2 <- npx_data1_no_controls %>%
      na.omit() %>% # removing missing values which exists for Site
      olink_boxplot(variable = "Site",
                    olinkid_list = c("OID00488", "OID01276"),
                    number_of_proteins_per_plot  = 2,
                    posthoc_results = anova_posthoc_results)

    plot2[[1]]

### Function output

A list of objects of class ggplot.

Note: Please note that plots will not appear in the plots panel of Rstudio if not assigned to a variable and printing it (see sample code above).

## Boxplots for QC (olink_dist_plot)

The olink_dist_plot function generates boxplots of NPX values for each sample, faceted by Olink panel. This is used as an initial QC step to identify potential outliers.

More information about `olink_dist_plot()` can be found in the Outlier Exclusion Vignette.

## Point-range plot for LMER (olink_lmer_plot)

The function olink_lmer_plot generates a point-range plot for a given list of proteins based on linear mixed effect model. The points illustrate the mean NPX level for each group and the error bars illustrate 95% confidence intervals. Facets are labeled by the protein name and corresponding OlinkID for the protein.

### Function arguments

- df: NPX data frame in long format should minimally contain protein name (Assay), OlinkID, UniProt, Panel and 1-2 variables with at least 2 levels and subject ID.
- variable: Single character value or character array. In case of single character then that should represent a column in the df. Otherwise, if length > 1, the included variable names will be used in crossed analyses. It can also accept the notations ‘:’ or ‘*’.
- outcome: Name of the column from df that contains the dependent variable. Default: NPX.
- random: Single character value or character array with random effects.
- covariates: Single character value or character array. Default: NULL. Confounding factors to include in the analysis. In case of single character then that should represent a column in the df. It can also accept the notations ‘:’ or ‘*’, while crossed analysis will not be inferred from main effects.
- x_axis_variable: Character. Which main effect to use as x-axis in the plot.
- col_variable: Character. If provided, the interaction effect col_variable:x_axis_variable will be plotted with x_axis_variable on the x-axis and col_variable as color.
- number_of_proteins_per_plot: Number plots to include in the list of point-range plots. Defaults to 6 plots per figure.
- verbose: Logical. Default: True. If information about removed samples, factor conversion and final model formula is to be printed to the console.

    if (requireNamespace("lme4", quietly = TRUE) & requireNamespace("lmerTest", quietly = TRUE)){
      plot <- olink_lmer_plot(df = npx_data1,
                              olinkid_list = c("OID01216", "OID01217"),
                              variable = c('Site', 'Treatment'),
                              x_axis_variable =  'Site',
                              col_variable = 'Treatment',
                              random = 'Subject')
      plot[[1]]
    }

### Function output

A list of objects of class ggplot.

Note: Please note that plots will not appear in the plots panel of Rstudio if not assigned to a variable and printing it (see sample code above).

## Heatmap for visualizing pathway enrichment (olink_pathway_heatmap)

The olink_pathway_heatmap function generates a heatmap of proteins related to pathways using the enrichment results from the olink_pathway_enrichment function. Either the top terms can be visualized or terms containing a certain keyword. For each term, the proteins in the test_result data frame that are related to that term will be visualized by their estimate. This visualization can be used to determining how many proteins of interest are involved in a particular pathway and in which direction their estimates are.

### Function arguments

- enrich_results: data frame of enrichment results from olink_pathway_enrichment()
- test_results: filtered results from statistical test with Assay, OlinkID, and estimate columns
- method: method used in olink_pathway_enrichment (“GSEA” (default) or “ORA”)
- keyword: (optional) keyword to filter enrichment results on, if not specified, displays top terms
- number_of_terms: number of terms to display, default is 20

    # GSEA Heatmap from t-test results
    try({ # This expression might fail if dependencies are not installed
    olink_pathway_heatmap(enrich_results = gsea_results,
                          test_results = ttest_results)
    })

    # ORA Heatmap from t-test results with cell keyword
    try({ # This expression might fail if dependencies are not installed
    olink_pathway_heatmap(enrich_results = ora_results,
                          test_results = ttest_results,
                          method = "ORA", keyword = "immune")
    })

### Function output

A heatmap as a ggplot object

## Bargraph for visualizing pathway enrichment (olink_pathway_visualization)

The olink_pathway_visualization function generates a bar graph of the top terms or terms related to a certain keyword for results from the olink_pathway_enrichment function. The bar represents either the normalized enrichment score (NES) for GSEA results or counts (number of proteins) for ORA results colored by adjusted p-value. Pathways are ordered by unadjusted p-value. The ORA visualization also contains the number of proteins out of the total proteins in that pathway as a ratio after the bar.

### Function arguments

- enrich_results: data frame of enrichment results from olink_pathway_enrichment()
- method: method used in olink_pathway_enrichment (“GSEA” (default) or “ORA”)
- keyword: (optional) keyword to filter enrichment results on, if not specified, displays top terms
- number_of_terms: number of terms to display, default is 20

### Function output

A bar graph as a ggplot object.

## Scatterplot for QC (olink_qc_plot)

The olink_qc_plot function generates a plot faceted by Panel, plotting IQR vs. median for all samples. This is a good first check to find out if any samples have a tendency to be classified as outliers. Horizontal dashed lines indicate +/-3 standard deviations from the mean IQR. Vertical dashed lines indicate +/-3 standard deviations from the mean sample median.

More information about `olink_qc_plot()` can be found in the Outlier Exclusion Vignette.

## Heatmap (olink_heatmap_plot)

The olink_heatmap_plot function generates a heatmap for all samples and proteins. By default, the heatmap centers and scales NPX across all proteins and clusters samples and proteins using a dendrogram. Unique sample names are required.

The grouping variable(s) are annotated and colored in the left side of the heatmap.

### Function arguments

- df: NPX data frame in long format which should minimally contain SampleID, NPX, OlinkID, Assay. Optionally columns of choice for annotations.
- variable_row_list: Columns in df to be annotated for rows in the heatmap.
- variable_col_list: Columns in df to be annotated for columns in the heatmap.
- center_scale: Logical. Default: True. If data should be centered and scaled across assays.
- cluster_rows: Logical. Default: True. Determining if rows should be clustered.
- cluster_cols: Logical. Default: True. Determining if columns should be clustered.
- show_rownames: Logical. Default: True. Determining if row names are shown.
- show_colnames: Logical. Default: True. Determining if column names are shown.
- annotation_legend: Logical. Default: True. Determining if legend for annotations should be shown.
- fontsize. Default: 10. Fontsize for all text.
- na_col. Default: Black. Color of the cells with NA.

    first10 <- npx_data1 %>%
      dplyr::pull(OlinkID) %>%
      unique() %>%
      head(10)

    first15samples <- npx_data1$SampleID %>%
      unique() %>%
      head(15)

    npx_data_small <- npx_data1 %>%
      dplyr::filter(!str_detect(SampleID, 'CONT')) %>%
      dplyr::filter(OlinkID %in% first10) %>%
      dplyr::filter(SampleID %in% first15samples)

    try({ #This statement might fail if dependencies are not installed
      olink_heatmap_plot(npx_data_small, variable_row_list =  'Treatment')
      })

### Function output

An object of class ggplot.

## Plot results of t-test (olink_volcano_plot)

The olink_volcano_plot function generates a volcano plot using results from the olink_ttest function. The estimated difference is shown in the x-axis and -log10(p-value) in the y-axis. The horizontal dotted line indicates p-value = 0.05. Dots are colored based on significance following Benjamini-Hochberg adjustment with a p-value cutoff of 0.05. Significant assays after adjustment can optionally be annotated by OlinkID.

### Function arguments

- p.val_tbl: a data frame of results generated by olink_ttest.
- x_lab: Optional. Character value to use as the x-axis label.
- olinkid_list: Optional. Character vector of proteins (OlinkID) to label in the plot. If not provided, by default the function will label all significant proteins.

    # perform t-test
    ttest_results <- olink_ttest(df = npx_data1,
                                 variable = 'Treatment')
    # select names of proteins to show
    top_10_name <- ttest_results %>%
      dplyr::slice_head(n = 10) %>%
      dplyr::pull(OlinkID)
    # volcano plot
    olink_volcano_plot(p.val_tbl = ttest_results,
                       x_lab = 'Treatment',
                       olinkid_list = top_10_name)

### Function output

An object of class ggplot.

## Theming function (set_plot_theme)

This function sets a coherent plot theme for plots by adding it to a ggplot object. It is mainly used for aesthetic reasons.

    npx_data1 %>%
      dplyr::filter(!is.na(Treatment)) |>
      dplyr::filter(OlinkID == 'OID01216') %>%
      ggplot(aes(x = Treatment, y = NPX, fill = Treatment)) +
      geom_boxplot() +
      set_plot_theme()

## Color theming (olink_color_discrete, olink_color_gradient, olink_fill_discrete, olink_fill_gradient)

These functions sets a coherent coloring theme for the plots by adding it to a ggplot object. It is mainly used for aesthetic reasons.

    npx_data1 %>%
      dplyr::filter(!is.na(Treatment)) |>
      dplyr::filter(OlinkID == 'OID01216') %>%
      ggplot(aes(x = Treatment, y = NPX, fill = Treatment)) +
      geom_boxplot() +
      set_plot_theme() +
      olink_fill_discrete()

## Visualizing bridgeability criteria for between-product normalization (olink_bridgeability_plot)

The olink_bridgeability_plot function generates a series of plots on a per-assay basis for a data frame generated from between-product bridging. The coloration of the figure headers indicate whether that assay has been defined as bridgeable or not bridgeable. The correlation plot, violin plot, and bar chart figures illustrate the three criteria for determining whether an assay is bridgeable. For assays determined to be bridgeable, the ECDF curve and corresponding KS statistic are used to determine which normalization approach (median centering or quantile smoothing) is most suitable for between-product normalization. For more information on the between-product bridging methodology and bridgeability criteria, consult the Bridging across NGS-based Olink® products Tutorial.

### Function arguments

- df: NPX data frame generated from between-product bridging in long format. Should minimally contain Assay, OlinkID, OlinkID_E3072, and all data points corresponding to bridging samples from the reference project and the new project.
- median_counts_threshold: Integer indicating minimum median counts allowed for each platform. If either platform has median counts below 150 for an assay, the assay fails the counts criteria when evaluating bridgeability. Default: 150.
- min_counts: Integer indicating minimum counts allowed for a data point. If any data point in the bridge normalized dataframe contains fewer than the defined minimum count cutoff, it is excluded from the bridgeability assessment and corresponding figures. Default: 10.
- bridge_sampleid: Character vector containing overlapping SampleIDs between the two bridging projects. If this argument is not provided, the function will look for overlapping SampleID values between the two projects in the bridged dataframe. Default: NULL.

    npx_ht <- data_exploreht |>
      dplyr::filter(SampleType == "SAMPLE") |>
      dplyr::mutate(Project = "data1")

    npx_3072 <- data_explore3072 |>
      dplyr::filter(SampleType == "SAMPLE") |>
      dplyr::mutate(Project = "data2")

    overlapping_samples <- unique(intersect(npx_ht |>
                                              dplyr::distinct(SampleID) |>
                                              dplyr::pull(),
                                            npx_3072 |>
                                              dplyr::distinct(SampleID) |>
                                              dplyr::pull()))

    npx_br_data <- OlinkAnalyze::olink_normalization(
      df1 = npx_ht,
      df2 = npx_3072,
      overlapping_samples_df1 = overlapping_samples,
      df1_project_nr = "Explore HT",
      df2_project_nr = "Explore 3072",
      reference_project = "Explore HT"
      )

    npx_br_data_bridgeable_plt <- olink_bridgeability_plot(
      data = npx_br_data,
      median_counts_threshold = 150,
      min_count = 10
      )

    npx_br_data_bridgeable_plt[[1]]

### Function output

A list of objects of class ggplot.

## Contact Us

We are always happy to help. Email us with any questions:

- biostat@olink.com for statistical services and general stats questions

- support@olink.com for Olink lab product and technical support

- info@olink.com for more information

## Legal Disclaimer

© 2025 Olink Proteomics AB, part of Thermo Fisher Scientific.

Olink products and services are For Research Use Only. Not for use in diagnostic procedures.

All information in this document is subject to change without notice. This document is not intended to convey any warranties, representations and/or recommendations of any kind, unless such warranties, representations and/or recommendations are explicitly stated.

Olink assumes no liability arising from a prospective reader’s actions based on this document.

OLINK, NPX, PEA, PROXIMITY EXTENSION, INSIGHT and the Olink logotype are trademarks registered, or pending registration, by Olink Proteomics AB. All third-party trademarks are the property of their respective owners.

Olink products and assay methods are covered by several patents and patent applications https://www.olink.com/patents/.