Serum Proteomics, Longitudinal

A full DIA-MS pipeline for longitudinal treatment response.

Overview

Problem. Find response markers in small, sparse, longitudinal data.

Learning goals

• Missingness itself can be signal (MNAR)
• Non-parametric for small n; mixed models for longitudinal

Figures

Tutorial

When to Use This Skill

Use this skill when ALL of the following are true: - You have serum or plasma proteomics data (DIA-MS, Olink, SomaScan, or similar) - The study has ≥2 longitudinal timepoints (e.g., baseline T1, mid-treatment T2, end T3) - Samples are from ≥2 response groups (e.g., Early Responder, Poor Responder) - You want MNAR detection, longitudinal modelling, or biomarker tiering — not just standard DE

Inputs

Input format notes

• Intensity values: Log2-transformed preferred. If raw intensities (median > 100), scripts auto-apply log2.
• NA encoding: Use NA or empty cells for non-detected proteins. Do NOT use 0 — zero is a valid intensity.
• Sample IDs: Must match between intensity matrix column names and metadata sample_id column.
• Protein names: Use gene symbols (CHGA, SCG2, ITIH1) or UniProt accessions. Consistent naming required.
• Supported platforms: Spectronaut 16 (DIA-MS), MaxQuant (LFQ), Olink NPX, SomaScan RFU — any platform that produces a protein × sample matrix with NA for non-detection.

Outputs

Installation

options(repos = c(CRAN = "https://cloud.r-project.org"))

# Core (required)
install.packages(c("dplyr", "tidyr", "nlme"))

# Trajectory clustering (strongly recommended)
install.packages(c("dtw", "cluster"))

# Downstream skills
if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("limma")   # for proteomics-diff-exp upstream
install.packages(c("glmnet", "pROC"))  # for lasso-biomarker-panel downstream

Standard Workflow

MANDATORY: SOURCE SCRIPTS — DO NOT WRITE INLINE CODE

Step 0 — Load data

# Load your protein intensity matrix (rows = proteins, cols = samples, NA = not detected)
intensity_matrix <- read.csv("your_protein_matrix.csv", row.names = 1, check.names = FALSE)

# Load sample metadata
metadata <- read.csv("your_metadata.csv", stringsAsFactors = FALSE)
# Required columns: sample_id, patient_id, timepoint, response_group

# Optional: clinical score table for trajectory clustering
score_table <- read.csv("your_pasi_scores.csv", stringsAsFactors = FALSE)
# Required columns: patient_id, timepoint, PASI (or your score column)

# Optional: cross-species DE table
xspecies_df <- read.csv("your_mouse_de.csv", stringsAsFactors = FALSE)
# Required columns: protein, log2FC, p_value, direction

cat(sprintf("Matrix: %d proteins × %d samples\n", nrow(intensity_matrix), ncol(intensity_matrix)))
cat(sprintf("Missing values: %.1f%%\n", 100 * mean(is.na(intensity_matrix))))
cat(sprintf("Metadata: %d rows, groups: %s\n",
    nrow(metadata), paste(unique(metadata$response_group), collapse=", ")))

Step 1 — MNAR Detection

source("scripts/01_mnar_detection.R")

mnar_result <- run_mnar_detection(
  intensity_matrix = intensity_matrix,
  metadata         = metadata,
  params = list(
    baseline_timepoint = "T1",       # timepoint for MNAR test
    response_col       = "response_group",
    sample_col         = "sample_id",
    mnar_fisher_p      = 0.05        # Fisher p threshold
  )
)

# Optional: check granin dissociation pattern (CHGA vs SCG2)
check_granin_dissociation(mnar_result$mnar_table,
  params = list(responder_group = "Early"))

DO NOT impute MNAR proteins before this step. Imputation destroys the MNAR signal.

Step 2 — Wilcoxon Differential Expression

source("scripts/02_wilcoxon_de.R")

de_result <- run_wilcoxon_de(
  intensity_matrix = intensity_matrix,
  metadata         = metadata,
  mnar_table       = mnar_result$mnar_table,   # pass MNAR flags
  params = list(
    de_timepoint      = "T1",          # baseline comparison
    response_col      = "response_group",
    sample_col        = "sample_id",
    group_a           = "Early",       # numerator group
    group_b           = "Poor",        # denominator group
    de_adj_p          = 0.05,
    de_log2fc         = 0.58,
    min_detected_frac = 0.5
  )
)

# Quick volcano summary (text)
summarise_volcano(de_result$de_table)

Note: MNAR proteins are automatically appended to the DE table with is_mnar = TRUE and wilcox_p = NA. They are included in the biomarker tiering step.

Step 3 — Longitudinal Mixed-Effects Modelling

source("scripts/03_longitudinal_lme.R")

lme_result <- run_longitudinal_lme(
  intensity_matrix = intensity_matrix,
  metadata         = metadata,
  params = list(
    timepoints        = c("T1", "T2", "T3"),
    timepoint_col     = "timepoint",
    response_col      = "response_group",
    patient_col       = "patient_id",
    covariates        = c("baseline_PASI", "age", "sex"),  # adjust to your metadata
    ref_timepoint     = "T1",
    ref_group         = "Poor",          # reference group for contrasts
    lme_interaction_p = 0.05,
    t2_surge_group    = "Early"          # group expected to show T2 surge
  )
)

Note on convergence: 5–15% of proteins may fail to converge — this is normal. They are flagged as model_status = "model_failed" and excluded from interaction ranking but retained in the DE table.

Skip this step if: You have only 1 timepoint, or fewer than 3 patients with repeated measures.

Step 4 — PASI Trajectory Clustering (Optional)

source("scripts/04_trajectory_clustering.R")

cluster_result <- run_trajectory_clustering(
  score_table = score_table,
  params = list(
    patient_col   = "patient_id",
    timepoint_col = "timepoint",
    score_col     = "PASI",
    timepoints    = c("T1", "T2", "T3"),   # NULL = auto-detect
    k_range       = 2:7,
    flare_ratio   = 1.2,    # T2/T1 > 1.2 = flare
    clear_ratio   = 0.5     # T3/T1 < 0.5 = PASI50 response
  )
)

CRITICAL CLINICAL NOTE: If Flare→Clear patients are identified, do NOT recommend treatment discontinuation at 12 weeks based on T2 PASI alone. These patients show transient worsening before resolution — a pattern that would be misclassified as non-response under a standard 12-week endpoint.

Skip this step if: You have no clinical score data, or only 1 timepoint.

Step 5 — Biomarker Tiering

source("scripts/05_biomarker_tiering.R")

tier_result <- run_biomarker_tiering(
  de_result   = de_result,
  lme_result  = lme_result,    # NULL if Step 3 was skipped
  xspecies_df = xspecies_df,   # NULL if no cross-species data
  params = list(
    tier1_mnar_adj_p      = 0.001,
    tier1_lme_adj_p       = 0.001,
    tier1_require_t2surge = TRUE,
    tier2_de_adj_p        = 0.01,
    tier2_de_log2fc       = 1.0,
    tier3_de_adj_p        = 0.05,
    tier3_de_log2fc       = 0.58
  )
)

Step 6 — Export All Results

source("scripts/06_export_results.R")

export_all(
  mnar_result    = mnar_result,
  de_result      = de_result,
  lme_result     = lme_result,       # NULL if skipped
  cluster_result = cluster_result,   # NULL if skipped
  tier_result    = tier_result,
  output_dir     = "results",
  study_label    = "MMIP_psoriasis"  # change to your study name
)

Parameters Reference

Worked Example: MMIP Psoriasis Cohort

Study design: n=90 psoriasis patients, DIA-MS serum proteomics, 3 timepoints (T1=baseline, T2=12 weeks, T3=1 year), 3 response groups (Early Responder n=10, Late Responder n=10, Poor Responder n=10 for discovery; n=60 for validation).

Expected results from this pipeline:

Clarification Questions

ALWAYS ask Question 1 FIRST.

1. Data availability

• Do you have a protein intensity matrix with NA values for non-detected proteins?
• If yes: proceed to Question 2
• If no (all proteins quantified, no NAs): MNAR step will find nothing — consider proteomics-diff-exp instead

2. Study design

• How many timepoints? (2 = skip LME; 3+ = run full longitudinal model)
• How many patients per response group? (< 5 per group = results unreliable; 8–15 = typical discovery cohort)
• Do you have clinical score data (PASI or equivalent) for trajectory clustering?

3. Response group definition

• Are response groups pre-defined in your metadata, or do you need to derive them from clinical scores?
• Pre-defined: proceed directly
• Need to derive: use PASI50/75/90 thresholds or trajectory clustering first, then re-run DE

4. Covariates

• Which covariates are available in your metadata for the LME model?
• Recommended: baseline clinical score, age, sex
• If none available: run without covariates (set covariates = c())

Common Issues

Interpretation Guidelines

MNAR results

• MNAR_up_in_X: protein is MORE OFTEN detected in group X → group X has higher circulating levels (above LOD)
• MNAR_down_in_X: protein is LESS OFTEN detected in group X → group X has lower levels (below LOD)
• CHGA MNAR in responders = responders have LOWER CHGA (below LOD) = chromaffin granule depletion

DE results

• log2FC > 0 = higher in group_a; log2FC < 0 = higher in group_b
• MNAR proteins have wilcox_p = NA — their significance comes from fisher_adj_p
• Do NOT interpret MNAR proteins using log2FC alone

LME results

• trajectory_type = "T2_surge" = protein rises at T2 specifically in the surge group
• trajectory_type = "T2_dip" = protein falls at T2 specifically in the surge group
• model_status = "model_failed" = convergence failure — not a biological finding

Biomarker tiers

• Tier 1 (★★★): Recommend for immediate ELISA validation in independent cohort
• Tier 2 (★★): Recommend for multiplexed panel (Luminex/Olink) validation
• Tier 3 (★): Hypothesis-generating only; do not prioritise for validation without additional evidence
• Mechanistic annotation: "unknown" is not a negative finding — it may indicate a novel axis

AUC / ROC

• AUC values reported in analysis_summary.txt are from the discovery cohort only
• These are expected to be optimistic; external validation is required before clinical claims
• For validated AUC, use lasso-biomarker-panel with an independent cohort

Suggested Next Steps

Related Skills

References

• MNAR handling: Lazar C, et al. (2016) J Proteome Res 15(4):1116–1125
• Wilcoxon for small proteomics: Kammers K, et al. (2015) EuPA Open Proteomics 7:11–19
• nlme mixed-effects: Pinheiro J, Bates D (2000) Mixed-Effects Models in S and S-PLUS. Springer
• DTW clustering: Sakoe H, Chiba S (1978) IEEE Trans Acoust 26(1):43–49
• BH correction: Benjamini Y, Hochberg Y (1995) J R Stat Soc B 57(1):289–300
• MMIP study: Wan F, et al. (in preparation) — Sympatho-adrenal hypothesis for psoriasis treatment response
• See references/method-notes.md for full statistical rationale

Code preview

scripts/01_mnar_detection.R

# =============================================================================
# 01_mnar_detection.R
# Missing-Not-At-Random (MNAR) protein detection for serum proteomics
#
# Purpose:
#   Identify proteins whose absence from serum is non-random with respect to
#   treatment response group. In DIA-MS data, a protein below the detection
#   limit (NA) carries biological information distinct from a low-abundance
#   quantified protein. This script separates MNAR proteins from quantified
#   proteins before downstream DE analysis.
#
# Inputs (loaded from environment or passed as arguments):
#   intensity_matrix  - data.frame, rows = proteins, cols = sample_ids
#                       NA values indicate non-detection (below LOD)
#   metadata          - data.frame with columns: sample_id, patient_id,
#                       timepoint, response_group
#   params            - list of analysis parameters (see run_mnar_detection())
#
# Outputs:
#   mnar_results.csv  - per-protein MNAR classification table
#   Returns:          - list with $mnar_table and $mnar_proteins vector
# =============================================================================

suppressPackageStartupMessages({
  library(dplyr)
  library(tidyr)
})

# -----------------------------------------------------------------------------
# Main function
# -----------------------------------------------------------------------------

run_mnar_detection <- function(
    intensity_matrix,
    metadata,
    params = list()
) {
  # --- Default parameters ---
  p <- modifyList(
    list(
      baseline_timepoint  = "T1",          # timepoint to use for MNAR test
      timepoint_col       = "timepoint",
      response_col        = "response_group",
      sample_col          = "sample_id",
      mnar_fisher_p       = 0.05,          # Fisher p threshold
      min_detected_any    = 2,             # min detections in any group to test
      responder_label     = NULL           # e.g. "Early"; NULL = first level
    ),
    params
  )

  cat("=== MNAR Detection ===\n")
  cat(sprintf("Baseline timepoint: %s\n", p$baseline_timepoint))
  cat(sprintf("Fisher p threshold: %.3f\n", p$mnar_fisher_p))

  # --- Subset metadata to baseline timepoint ---
  meta_t1 <- metadata %>%
    filter(.data[[p$timepoint_col]] == p$baseline_timepoint)

  if (nrow(meta_t1) == 0) {
    stop(sprintf(
      "No samples found for timepoint '%s'. Check params$baseline_timepoint.",
      p$baseline_timepoint
    ))
  }

  # --- Align matrix to baseline samples ---
  t1_samples <- meta_t1[[p$sample_col]]
  t1_samples <- intersect(t1_samples, colnames(intensity_matrix))

  if (length(t1_samples) == 0) {
    stop("No overlap between metadata sample_ids and intensity_matrix column names.")
  }

  mat_t1 <- intensity_matrix[, t1_samples, drop = FALSE]
  meta_t1 <- meta_t1[meta_t1[[p$sample_col]] %in% t1_samples, ]

  groups <- unique(meta_t1[[p$response_col]])
  cat(sprintf("Response groups: %s\n", paste(groups, collapse = ", ")))
  cat(sprintf("Proteins to test: %d\n", nrow(mat_t1)))

scripts/02_wilcoxon_de.R

# =============================================================================
# 02_wilcoxon_de.R
# Wilcoxon rank-sum differential expression for small serum proteomics cohorts
#
# Purpose:
#   Non-parametric DE for discovery cohorts (typically n=8–15 per group) where
#   normality cannot be assumed. Handles MNAR proteins separately: they are
#   flagged from 01_mnar_detection.R and excluded from the Wilcoxon test.
#   log2FC is computed as the difference of group medians (on log2 scale).
#
# Inputs:
#   intensity_matrix  - data.frame, rows = proteins, cols = sample_ids (log2)
#   metadata          - data.frame: sample_id, patient_id, timepoint, response_group
#   mnar_table        - output$mnar_table from run_mnar_detection() [optional]
#   params            - list of analysis parameters
#
# Outputs:
#   de_results.csv    - full DE table with log2FC, p-value, adj_p, MNAR flag
#   Returns:          - list with $de_table
# =============================================================================

suppressPackageStartupMessages({
  library(dplyr)
  library(tidyr)
})

# -----------------------------------------------------------------------------
# Main function
# -----------------------------------------------------------------------------

run_wilcoxon_de <- function(
    intensity_matrix,
    metadata,
    mnar_table  = NULL,
    params      = list()
) {
  p <- modifyList(
    list(
      de_timepoint        = "T1",          # timepoint for baseline DE
      timepoint_col       = "timepoint",
      response_col        = "response_group",
      sample_col          = "sample_id",
      group_a             = NULL,          # "numerator" group (e.g. "Early")
      group_b             = NULL,          # "denominator" group (e.g. "Poor")
      de_adj_p            = 0.05,
      de_log2fc           = 0.58,          # ~1.5-fold
      min_detected_frac   = 0.5,           # min fraction detected in ≥1 group
      log2_transform      = TRUE           # apply log2 if not already done
    ),
    params
  )

  cat("=== Wilcoxon Differential Expression ===\n")
  cat(sprintf("Timepoint: %s\n", p$de_timepoint))

  # --- Subset to DE timepoint ---
  meta_de <- metadata %>%
    filter(.data[[p$timepoint_col]] == p$de_timepoint)

  de_samples <- intersect(meta_de[[p$sample_col]], colnames(intensity_matrix))
  mat_de <- intensity_matrix[, de_samples, drop = FALSE]
  meta_de <- meta_de[meta_de[[p$sample_col]] %in% de_samples, ]

  # --- Determine groups ---
  all_groups <- unique(meta_de[[p$response_col]])
  if (is.null(p$group_a) || is.null(p$group_b)) {
    if (length(all_groups) < 2) stop("Need at least 2 response groups.")
    p$group_a <- as.character(all_groups[1])
    p$group_b <- as.character(all_groups[2])
    cat(sprintf("Auto-selected groups: A=%s vs B=%s\n", p$group_a, p$group_b))
  }
  cat(sprintf("Comparison: %s (A) vs %s (B)\n", p$group_a, p$group_b))
  cat(sprintf("log2FC = median(A) - median(B)\n"))

  samps_a <- meta_de[[p$sample_col]][meta_de[[p$response_col]] == p$group_a]
  samps_b <- meta_de[[p$sample_col]][meta_de[[p$response_col]] == p$group_b]
  samps_a <- intersect(samps_a, colnames(mat_de))
  samps_b <- intersect(samps_b, colnames(mat_de))

  cat(sprintf("n(A)=%d, n(B)=%d\n", length(samps_a), length(samps_b)))

scripts/03_longitudinal_lme.R

# =============================================================================
# 03_longitudinal_lme.R
# Longitudinal mixed-effects modelling for serum proteomics
#
# Purpose:
#   Fit per-protein linear mixed-effects models to detect proteins with
#   significant time × response_group interaction — the statistical signature
#   of "psychological reactogenicity" (T2 surge in responders, absent in
#   non-responders). Random intercept per patient accounts for repeated measures.
#
# Model:
#   intensity ~ time * response_group + [covariates] + (1 | patient_id)
#   Fitted with nlme::lme(); time and response_group are treated as factors.
#
# Inputs:
#   intensity_matrix  - data.frame, rows = proteins, cols = sample_ids (log2)
#   metadata          - data.frame: sample_id, patient_id, timepoint,
#                       response_group [+ optional covariate columns]
#   params            - list of analysis parameters
#
# Outputs:
#   lme_results.csv   - per-protein interaction p-value, trajectory type
#   Returns:          - list with $lme_table
# =============================================================================

suppressPackageStartupMessages({
  library(nlme)
  library(dplyr)
  library(tidyr)
})

# -----------------------------------------------------------------------------
# Main function
# -----------------------------------------------------------------------------

run_longitudinal_lme <- function(
    intensity_matrix,
    metadata,
    params = list()
) {
  p <- modifyList(
    list(
      timepoints          = c("T1", "T2", "T3"),
      timepoint_col       = "timepoint",
      response_col        = "response_group",
      sample_col          = "sample_id",
      patient_col         = "patient_id",
      covariates          = c(),            # e.g. c("age", "sex", "baseline_PASI")
      ref_timepoint       = "T1",           # reference level for time factor
      ref_group           = NULL,           # reference level for response_group
      lme_interaction_p   = 0.05,
      min_obs_per_protein = 10,             # min non-NA observations to fit model
      t2_surge_group      = NULL,           # group expected to show T2 surge
      log2_transform      = TRUE
    ),
    params
  )

  cat("=== Longitudinal Mixed-Effects Modelling ===\n")
  cat(sprintf("Timepoints: %s\n", paste(p$timepoints, collapse = " → ")))
  cat(sprintf("Interaction p threshold: %.3f\n", p$lme_interaction_p))

  # --- Subset metadata to relevant timepoints ---
  meta_long <- metadata %>%
    filter(.data[[p$timepoint_col]] %in% p$timepoints)

  long_samples <- intersect(meta_long[[p$sample_col]], colnames(intensity_matrix))
  mat_long <- intensity_matrix[, long_samples, drop = FALSE]
  meta_long <- meta_long[meta_long[[p$sample_col]] %in% long_samples, ]

  # --- Optional log2 transform ---
  if (p$log2_transform) {
    med_val <- median(as.matrix(mat_long), na.rm = TRUE)
    if (med_val > 100) {
      mat_long <- log2(mat_long + 1)
    }
  }

  # --- Set factor levels ---
  meta_long[[p$timepoint_col]] <- factor(

Companion files